|
1 | 1 | import macro from 'vtk.js/Sources/macros';
|
2 | 2 | import vtkSpline1D from 'vtk.js/Sources/Common/DataModel/Spline1D';
|
3 | 3 |
|
4 |
| -const { vtkErrorMacro } = macro; |
| 4 | +import { BoundaryCondition } from 'vtk.js/Sources/Common/DataModel/Spline1D/Constants'; |
| 5 | + |
| 6 | +const VTK_EPSILON = 0.0001; |
5 | 7 |
|
6 | 8 | // ----------------------------------------------------------------------------
|
7 | 9 | // vtkCardinalSpline1D methods
|
@@ -46,9 +48,9 @@ function vtkCardinalSpline1D(publicAPI, model) {
|
46 | 48 | const dN = work[N];
|
47 | 49 |
|
48 | 50 | // solve resulting set of equations.
|
49 |
| - model.coefficients[0 * 4 + 2] = 0; |
| 51 | + model.coefficients[4 * 0 + 2] = 0; |
50 | 52 | work[0] = 0;
|
51 |
| - model.coefficients[0 * 4 + 3] = 1; |
| 53 | + model.coefficients[4 * 0 + 3] = 1; |
52 | 54 |
|
53 | 55 | for (let k = 1; k <= N; k++) {
|
54 | 56 | model.coefficients[4 * k + 1] -=
|
@@ -114,8 +116,145 @@ function vtkCardinalSpline1D(publicAPI, model) {
|
114 | 116 |
|
115 | 117 | // --------------------------------------------------------------------------
|
116 | 118 |
|
117 |
| - publicAPI.computeOpenCoefficients = (size, work, x, y) => { |
118 |
| - vtkErrorMacro('Open splines are not implemented yet!'); |
| 119 | + publicAPI.computeOpenCoefficients = (size, work, x, y, options = {}) => { |
| 120 | + if (!model.coefficients || model.coefficients.length !== 4 * size) { |
| 121 | + model.coefficients = new Float32Array(4 * size); |
| 122 | + } |
| 123 | + const N = size - 1; |
| 124 | + // develop constraint at leftmost point. |
| 125 | + switch (options.leftConstraint) { |
| 126 | + case BoundaryCondition.DERIVATIVE: |
| 127 | + // desired slope at leftmost point is leftValue. |
| 128 | + model.coefficients[4 * 0 + 1] = 1.0; |
| 129 | + model.coefficients[4 * 0 + 2] = 0.0; |
| 130 | + work[0] = options.leftValue; |
| 131 | + break; |
| 132 | + case BoundaryCondition.SECOND_DERIVATIVE: |
| 133 | + // desired second derivative at leftmost point is leftValue. |
| 134 | + model.coefficients[4 * 0 + 1] = 2.0; |
| 135 | + model.coefficients[4 * 0 + 2] = 1.0; |
| 136 | + work[0] = |
| 137 | + 3.0 * ((y[1] - y[0]) / (x[1] - x[0])) - |
| 138 | + 0.5 * (x[1] - x[0]) * options.leftValue; |
| 139 | + break; |
| 140 | + case BoundaryCondition.SECOND_DERIVATIVE_INTERIOR_POINT: |
| 141 | + // desired second derivative at leftmost point is |
| 142 | + // leftValue times second derivative at first interior point. |
| 143 | + model.coefficients[4 * 0 + 1] = 2.0; |
| 144 | + |
| 145 | + if (Math.abs(options.leftValue + 2) > VTK_EPSILON) { |
| 146 | + model.coefficients[4 * 0 + 2] = |
| 147 | + 4.0 * ((0.5 + options.leftValue) / (2.0 + options.leftValue)); |
| 148 | + work[0] = |
| 149 | + 6.0 * |
| 150 | + ((1.0 + options.leftValue) / (2.0 + options.leftValue)) * |
| 151 | + ((y[1] - y[0]) / (x[1] - x[0])); |
| 152 | + } else { |
| 153 | + model.coefficients[4 * 0 + 2] = 0; |
| 154 | + work[0] = 0; |
| 155 | + } |
| 156 | + break; |
| 157 | + case BoundaryCondition.DEFAULT: |
| 158 | + default: |
| 159 | + // desired slope at leftmost point is derivative from two points |
| 160 | + model.coefficients[4 * 0 + 1] = 1.0; |
| 161 | + model.coefficients[4 * 0 + 2] = 0.0; |
| 162 | + work[0] = y[2] - y[0]; |
| 163 | + break; |
| 164 | + } |
| 165 | + |
| 166 | + for (let k = 1; k < N; k++) { |
| 167 | + const xlk = x[k] - x[k - 1]; |
| 168 | + const xlkp = x[k + 1] - x[k]; |
| 169 | + |
| 170 | + model.coefficients[4 * k + 0] = xlkp; |
| 171 | + model.coefficients[4 * k + 1] = 2 * (xlkp + xlk); |
| 172 | + model.coefficients[4 * k + 2] = xlk; |
| 173 | + work[k] = |
| 174 | + 3.0 * |
| 175 | + ((xlkp * (y[k] - y[k - 1])) / xlk + (xlk * (y[k + 1] - y[k])) / xlkp); |
| 176 | + } |
| 177 | + |
| 178 | + // develop constraint at rightmost point. |
| 179 | + switch (options.rightConstraint) { |
| 180 | + case BoundaryCondition.DERIVATIVE: |
| 181 | + // desired slope at rightmost point is rightValue |
| 182 | + model.coefficients[4 * N + 0] = 0.0; |
| 183 | + model.coefficients[4 * N + 1] = 1.0; |
| 184 | + work[N] = options.rightValue; |
| 185 | + break; |
| 186 | + case BoundaryCondition.SECOND_DERIVATIVE: |
| 187 | + // desired second derivative at rightmost point is rightValue. |
| 188 | + model.coefficients[4 * N + 0] = 1.0; |
| 189 | + model.coefficients[4 * N + 1] = 2.0; |
| 190 | + work[N] = |
| 191 | + 3.0 * ((y[N] - y[N - 1]) / (x[N] - x[N - 1])) + |
| 192 | + 0.5 * (x[N] - x[N - 1]) * options.rightValue; |
| 193 | + break; |
| 194 | + case BoundaryCondition.SECOND_DERIVATIVE_INTERIOR_POINT: |
| 195 | + // desired second derivative at rightmost point is |
| 196 | + // rightValue times second derivative at last interior point. |
| 197 | + model.coefficients[4 * N + 1] = 2.0; |
| 198 | + if (Math.abs(options.rightValue + 2) > VTK_EPSILON) { |
| 199 | + model.coefficients[4 * N + 0] = |
| 200 | + 4.0 * ((0.5 + options.rightValue) / (2.0 + options.rightValue)); |
| 201 | + work[N] = |
| 202 | + 6.0 * |
| 203 | + ((1.0 + options.rightValue) / (2.0 + options.rightValue)) * |
| 204 | + ((y[N] - y[size - 2]) / (x[N] - x[size - 2])); |
| 205 | + } else { |
| 206 | + model.coefficients[4 * N + 0] = 0; |
| 207 | + work[N] = 0; |
| 208 | + } |
| 209 | + break; |
| 210 | + case BoundaryCondition.DEFAULT: |
| 211 | + default: |
| 212 | + // desired slope at rightmost point is derivative from two points |
| 213 | + model.coefficients[4 * N + 0] = 0.0; |
| 214 | + model.coefficients[4 * N + 1] = 1.0; |
| 215 | + work[N] = y[N] - y[N - 2]; |
| 216 | + break; |
| 217 | + } |
| 218 | + |
| 219 | + // solve resulting set of equations. |
| 220 | + model.coefficients[4 * 0 + 2] /= model.coefficients[4 * 0 + 1]; |
| 221 | + work[0] /= model.coefficients[4 * N + 1]; |
| 222 | + model.coefficients[4 * N + 3] = 1; |
| 223 | + |
| 224 | + for (let k = 1; k <= N; k++) { |
| 225 | + model.coefficients[4 * k + 1] -= |
| 226 | + model.coefficients[4 * k + 0] * model.coefficients[4 * (k - 1) + 2]; |
| 227 | + model.coefficients[4 * k + 2] /= model.coefficients[4 * k + 1]; |
| 228 | + work[k] = |
| 229 | + (work[k] - model.coefficients[4 * k + 0] * work[k - 1]) / |
| 230 | + model.coefficients[4 * k + 1]; |
| 231 | + } |
| 232 | + |
| 233 | + for (let k = N - 1; k >= 0; k--) { |
| 234 | + work[k] -= model.coefficients[4 * k + 2] * work[k + 1]; |
| 235 | + } |
| 236 | + |
| 237 | + // the column vector work now contains the first |
| 238 | + // derivative of the spline function at each joint. |
| 239 | + // compute the coefficients of the cubic between |
| 240 | + // each pair of joints. |
| 241 | + for (let k = 0; k < N; k++) { |
| 242 | + const b = x[k + 1] - x[k]; |
| 243 | + model.coefficients[4 * k + 0] = y[k]; |
| 244 | + model.coefficients[4 * k + 1] = work[k]; |
| 245 | + model.coefficients[4 * k + 2] = |
| 246 | + (3 * (y[k + 1] - y[k])) / (b * b) - (work[k + 1] + 2 * work[k]) / b; |
| 247 | + model.coefficients[4 * k + 3] = |
| 248 | + (2 * (y[k] - y[k + 1])) / (b * b * b) + |
| 249 | + (work[k + 1] + work[k]) / (b * b); |
| 250 | + } |
| 251 | + |
| 252 | + // the coefficients of a fictitious nth cubic |
| 253 | + // are the same as the coefficients in the first interval |
| 254 | + model.coefficients[4 * N + 0] = y[N]; |
| 255 | + model.coefficients[4 * N + 1] = work[N]; |
| 256 | + model.coefficients[4 * N + 2] = model.coefficients[4 * 0 + 2]; |
| 257 | + model.coefficients[4 * N + 3] = model.coefficients[4 * 0 + 3]; |
119 | 258 | };
|
120 | 259 |
|
121 | 260 | // --------------------------------------------------------------------------
|
|
0 commit comments