1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2013-2015 Imperial College London
5 * Copyright 2013-2015 Andreas Schuh
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/BSplineFreeFormTransformationSV.h"
21
22 #include "mirtk/Math.h"
23 #include "mirtk/Memory.h"
24 #include "mirtk/Profiling.h"
25 #include "mirtk/ScalingAndSquaring.h"
26 #include "mirtk/DisplacementToVelocityFieldBCH.h"
27 #include "mirtk/VelocityToDisplacementFieldSS.h"
28 #include "mirtk/LieBracketImageFilter.h"
29 #include "mirtk/VoxelFunction.h"
30 #include "mirtk/NaryVoxelFunction.h"
31 #include "mirtk/ImageToInterpolationCoefficients.h"
32 #include "mirtk/HomogeneousTransformation.h"
33 #include "mirtk/MultiLevelTransformation.h"
34 #include "mirtk/InterpolateImageFunction.hxx"
35
36 #include "FreeFormTransformationIntegration.h"
37
38
39 namespace mirtk {
40
41
42 // =============================================================================
43 // Integration methods
44 // =============================================================================
45
46 // Tolerance of embedded Runge-Kutta methods with automatic step length control
47 static const double SVFFD_RKTOL = 1.e-3;
48
49 MIRTK_FFDIM2(RKE1, BSplineFreeFormTransformationSV);
50 MIRTK_FFDIM2(RKEH12, BSplineFreeFormTransformationSV);
51 MIRTK_FFDIM2(RKE2, BSplineFreeFormTransformationSV);
52 MIRTK_FFDIM2(RKH2, BSplineFreeFormTransformationSV);
53 MIRTK_FFDIM2(RKBS23, BSplineFreeFormTransformationSV);
54 MIRTK_FFDIM2(RK4, BSplineFreeFormTransformationSV);
55 MIRTK_FFDIM2(RKF45, BSplineFreeFormTransformationSV);
56 MIRTK_FFDIM2(RKCK45, BSplineFreeFormTransformationSV);
57 MIRTK_FFDIM2(RKDP45, BSplineFreeFormTransformationSV);
58
59 // =============================================================================
60 // Construction/Destruction
61 // =============================================================================
62
63 // -----------------------------------------------------------------------------
DefaultMaximumScaledVelocity(double dx,double dy,double dz)64 inline double DefaultMaximumScaledVelocity(double dx, double dy, double dz)
65 {
66 double ds = .0;
67 if (dx > .0 && (ds == .0 || dx < ds)) ds = dx;
68 if (dy > .0 && (ds == .0 || dy < ds)) ds = dy;
69 if (dz > .0 && (ds == .0 || dz < ds)) ds = dz;
70 return 0.4 * ds;
71 }
72
73 // -----------------------------------------------------------------------------
74 BSplineFreeFormTransformationSV
BSplineFreeFormTransformationSV()75 ::BSplineFreeFormTransformationSV()
76 :
77 _T (1.0),
78 _TimeUnit (1.0),
79 _NumberOfSteps (64),
80 _MaxScaledVelocity(-1.0),
81 _IntegrationMethod(FFDIM_FastSS),
82 _UseDenseBCHGrid (false),
83 _LieDerivative (false),
84 _NumberOfBCHTerms (4)
85 {
86 _ExtrapolationMode = Extrapolation_NN;
87 }
88
89 // -----------------------------------------------------------------------------
90 BSplineFreeFormTransformationSV
BSplineFreeFormTransformationSV(const ImageAttributes & attr,double dx,double dy,double dz)91 ::BSplineFreeFormTransformationSV(const ImageAttributes &attr,
92 double dx, double dy, double dz)
93 :
94 _T (1.0),
95 _TimeUnit (1.0),
96 _NumberOfSteps (64),
97 _MaxScaledVelocity(-1.0),
98 _IntegrationMethod(FFDIM_FastSS),
99 _UseDenseBCHGrid (false),
100 _LieDerivative (false),
101 _NumberOfBCHTerms (4)
102 {
103 _ExtrapolationMode = Extrapolation_NN;
104 Initialize(attr, dx, dy, dz);
105 }
106
107 // -----------------------------------------------------------------------------
108 BSplineFreeFormTransformationSV
BSplineFreeFormTransformationSV(const BaseImage & target,double dx,double dy,double dz)109 ::BSplineFreeFormTransformationSV(const BaseImage &target,
110 double dx, double dy, double dz)
111 :
112 _T (1.0),
113 _TimeUnit (1.0),
114 _NumberOfSteps (64),
115 _MaxScaledVelocity(-1.0),
116 _IntegrationMethod(FFDIM_FastSS),
117 _UseDenseBCHGrid (false),
118 _LieDerivative (false),
119 _NumberOfBCHTerms (4)
120 {
121 _ExtrapolationMode = Extrapolation_NN;
122 Initialize(target.Attributes(), dx, dy, dz);
123 }
124
125 // -----------------------------------------------------------------------------
126 BSplineFreeFormTransformationSV
BSplineFreeFormTransformationSV(const GenericImage<double> & image,bool disp)127 ::BSplineFreeFormTransformationSV(const GenericImage<double> &image, bool disp)
128 :
129 _T (1.0),
130 _TimeUnit (1.0),
131 _NumberOfSteps (64),
132 _MaxScaledVelocity(-1.0),
133 _IntegrationMethod(FFDIM_FastSS),
134 _UseDenseBCHGrid (false),
135 _LieDerivative (false),
136 _NumberOfBCHTerms (4)
137 {
138 Initialize(image, disp);
139 }
140
141 // -----------------------------------------------------------------------------
142 BSplineFreeFormTransformationSV
BSplineFreeFormTransformationSV(const BSplineFreeFormTransformationSV & ffd)143 ::BSplineFreeFormTransformationSV(const BSplineFreeFormTransformationSV &ffd)
144 :
145 BSplineFreeFormTransformation3D(ffd),
146 _T (ffd._T),
147 _TimeUnit (ffd._TimeUnit),
148 _NumberOfSteps (ffd._NumberOfSteps),
149 _MaxScaledVelocity(ffd._MaxScaledVelocity),
150 _IntegrationMethod(ffd._IntegrationMethod),
151 _UseDenseBCHGrid (ffd._UseDenseBCHGrid),
152 _LieDerivative (ffd._LieDerivative),
153 _NumberOfBCHTerms (ffd._NumberOfBCHTerms)
154 {
155 }
156
157 // -----------------------------------------------------------------------------
158 BSplineFreeFormTransformationSV
~BSplineFreeFormTransformationSV()159 ::~BSplineFreeFormTransformationSV()
160 {
161 }
162
163 // -----------------------------------------------------------------------------
Initialize(const ImageAttributes & attr)164 void BSplineFreeFormTransformationSV::Initialize(const ImageAttributes &attr)
165 {
166 ImageAttributes domain = attr;
167 domain._t = 1; // SV FFD can be used to register time series, but not regular FFD
168 BSplineFreeFormTransformation3D::Initialize(domain);
169 if (attr._t > 0 && !AreEqual(attr._dt, 0.)) {
170 _TimeUnit = attr._dt;
171 _T = attr._t * attr._dt;
172 }
173 if (_MaxScaledVelocity < .0) _MaxScaledVelocity = DefaultMaximumScaledVelocity(_dx, _dy, _dz);
174 }
175
176 // -----------------------------------------------------------------------------
Initialize(const ImageAttributes & attr,double dx,double dy,double dz,const Transformation * dof)177 void BSplineFreeFormTransformationSV::Initialize(const ImageAttributes &attr,
178 double dx, double dy, double dz,
179 const Transformation *dof)
180 {
181 // Initialize free-form deformation (for extended image grid)
182 //
183 // Ensure that for all target voxels the displacement can be recovered
184 // without requiring any extrapolation of the velocity field during
185 // computation of the trajectory (integration, i.e., exponentiation)
186 this->Initialize(ApproximationDomain(attr, dof), dx, dy, dz);
187
188 // Approximate given transformation
189 this->ApproximateAsNew(dof);
190 }
191
192 // -----------------------------------------------------------------------------
Subdivide(bool subdivide_x,bool subdivide_y,bool subdivide_z,bool subdivide_t)193 void BSplineFreeFormTransformationSV::Subdivide(bool subdivide_x, bool subdivide_y, bool subdivide_z, bool subdivide_t)
194 {
195 BSplineFreeFormTransformation3D::Subdivide(subdivide_x, subdivide_y, subdivide_z, subdivide_t);
196 if (_MaxScaledVelocity > .0) _MaxScaledVelocity /= 2;
197 }
198
199 // =============================================================================
200 // Auxiliary functions
201 // =============================================================================
202
203 // -----------------------------------------------------------------------------
204 /// Voxel function used to evaluate Lie bracket at each lattice point using
205 /// the Lie derivative definition which is based on the Jacobian matrices
206 /// of the vector fields. Uses nearest neighbor extrapolation of the velocity field.
207 class SVFFDEvaluateLieBracket : public VoxelFunction
208 {
209 private:
210
211 typedef BSplineFreeFormTransformationSV::Vector Vector;
212 typedef BSplineFreeFormTransformationSV::Kernel Kernel;
213
214 const BSplineFreeFormTransformationSV *_FFD;
215
216 double _tau; ///< Scaling of left/first vector field
217 double _eta; ///< Scaling of right/second vector field
218 const GenericImage<Vector> &_v; ///< Left/First vector field
219 const GenericImage<Vector> &_w; ///< Right/Second vector field
220
221 public:
222
223 // ---------------------------------------------------------------------------
224 /// Constructor
SVFFDEvaluateLieBracket(const BSplineFreeFormTransformationSV * ffd,const GenericImage<Vector> & v,const GenericImage<Vector> & w)225 SVFFDEvaluateLieBracket(const BSplineFreeFormTransformationSV *ffd,
226 const GenericImage<Vector> &v,
227 const GenericImage<Vector> &w)
228 :
229 _FFD(ffd), _tau(1.0), _eta(1.0), _v(v), _w(w)
230 {}
231
232 // ---------------------------------------------------------------------------
233 /// Constructor
SVFFDEvaluateLieBracket(const BSplineFreeFormTransformationSV * ffd,double tau,const GenericImage<Vector> & v,double eta,const GenericImage<Vector> & w)234 SVFFDEvaluateLieBracket(const BSplineFreeFormTransformationSV *ffd,
235 double tau,
236 const GenericImage<Vector> &v,
237 double eta,
238 const GenericImage<Vector> &w)
239 :
240 _FFD(ffd), _tau(tau), _eta(eta), _v(v), _w(w)
241 {}
242
243 // ---------------------------------------------------------------------------
244 /// Evaluate velocity at lattice point
Evaluate(Vector & x,double s,const GenericImage<Vector> & v,int i,int j,int k) const245 void Evaluate(Vector &x, double s, const GenericImage<Vector> &v, int i, int j, int k) const
246 {
247 double B_I, B_J, B_K;
248 int II, JJ, KK;
249
250 x = .0;
251 for (int K = k-1; K <= k+1; K++) {
252 B_K = Kernel::LatticeWeights[K - (k-1)];
253 if (K < 0) KK = 0;
254 else if (K >= v.Z()) KK = v.Z()-1;
255 else KK = K;
256 for (int J = j-1; J <= j+1; J++) {
257 B_J = Kernel::LatticeWeights[J - (j-1)];
258 if (J < 0) JJ = 0;
259 else if (J >= v.Y()) JJ = v.Y() - 1;
260 else JJ = J;
261 for (int I = i-1; I <= i+1; I++) {
262 B_I = Kernel::LatticeWeights[I - (i-1)];
263 if (I < 0) II = 0;
264 else if (I >= v.X()) II = v.X() - 1;
265 else II = I;
266 x += B_I * B_J * B_K * s * v(II, JJ, KK);
267 }
268 }
269 }
270 }
271
272 // ---------------------------------------------------------------------------
273 /// Evaluate Jacobian of velocity field at lattice point
Jacobian(Matrix & jac,double s,const GenericImage<Vector> & v,int i,int j,int k) const274 void Jacobian(Matrix &jac, double s, const GenericImage<Vector> &v, int i, int j, int k) const
275 {
276 int I, J, K, II, JJ, KK;
277 double B_I, B_J, B_K, B_I_I, B_J_I, B_K_I;
278 Vector dx, dy, dz;
279
280 for (K = k-1; K <= k+1; K++) {
281 B_K = Kernel::LatticeWeights [K - (k-1)];
282 B_K_I = Kernel::LatticeWeights_I[K - (k-1)];
283 if (K < 0) KK = 0;
284 else if (K >= v.Z()) KK = v.Z()-1;
285 else KK = K;
286 for (J = j-1; J <= j+1; J++) {
287 B_J = Kernel::LatticeWeights [J - (j-1)];
288 B_J_I = Kernel::LatticeWeights_I[J - (j-1)];
289 if (J < 0) JJ = 0;
290 else if (J >= v.Y()) JJ = v.Y()-1;
291 else JJ = J;
292 for (I = i-1; I <= i+1; I++) {
293 B_I = Kernel::LatticeWeights [I - (i-1)];
294 B_I_I = Kernel::LatticeWeights_I[I - (i-1)];
295 if (I < 0) II = 0;
296 else if (I >= v.X()) II = v.X()-1;
297 else II = I;
298 const Vector &coeff = v(II, JJ, KK);
299 dx += B_I_I * B_J * B_K * s * coeff;
300 dy += B_I * B_J_I * B_K * s * coeff;
301 dz += B_I * B_J * B_K_I * s * coeff;
302 }
303 }
304 }
305
306 jac.Initialize(3, 3);
307 jac(0, 0) = dx._x;
308 jac(0, 1) = dy._x;
309 jac(0, 2) = dz._x;
310 jac(1, 0) = dx._y;
311 jac(1, 1) = dy._y;
312 jac(1, 2) = dz._y;
313 jac(2, 0) = dx._z;
314 jac(2, 1) = dy._z;
315 jac(2, 2) = dz._z;
316
317 _FFD->JacobianToWorld(jac);
318 }
319
320 // ---------------------------------------------------------------------------
321 /// Compute product of 3x3 matrix and 3D column vector
MatrixProduct(const Matrix & jac,const Vector & vel)322 static Vector MatrixProduct(const Matrix &jac, const Vector &vel)
323 {
324 Vector v;
325 v._x = jac(0, 0) * vel._x + jac(0, 1) * vel._y + jac(0, 2) * vel._z;
326 v._y = jac(1, 0) * vel._x + jac(1, 1) * vel._y + jac(1, 2) * vel._z;
327 v._z = jac(2, 0) * vel._x + jac(2, 1) * vel._y + jac(2, 2) * vel._z;
328 return v;
329 }
330
331 // ---------------------------------------------------------------------------
332 /// Evaluate Lie bracket at given lattice point, u = [v, w]
operator ()(int i,int j,int k,int,Vector * u) const333 void operator ()(int i, int j, int k, int, Vector *u) const
334 {
335 Matrix jac(3, 3);
336 Vector vel;
337 // u = J_w * v
338 Jacobian(jac, _eta, _w, i, j, k);
339 Evaluate(vel, _tau, _v, i, j, k);
340 *u = MatrixProduct(jac, vel);
341 // u = J_w * v - J_v * w
342 Jacobian(jac, _tau, _v, i, j, k);
343 Evaluate(vel, _eta, _w, i, j, k);
344 *u -= MatrixProduct(jac, vel);
345 }
346 };
347
348 // -----------------------------------------------------------------------------
349 template <class TReal>
EvaluateBCHFormulaDense(int nterms,GenericImage<TReal> & u,const GenericImage<TReal> & v,const GenericImage<TReal> & w,bool minus_v,bool usejac)350 void EvaluateBCHFormulaDense(int nterms, GenericImage<TReal> &u,
351 const GenericImage<TReal> &v,
352 const GenericImage<TReal> &w,
353 bool minus_v, bool usejac)
354 {
355 MIRTK_START_TIMING();
356
357 GenericImage<TReal> l1, l2, l3, l4;
358 if (nterms >= 3) {
359 // - [v, w]
360 liebracket(&l1, &v, &w, usejac);
361 if (nterms >= 4) {
362 // - [v, [v, w]]
363 liebracket(&l2, &v, &l1, usejac);
364 if (nterms >= 5) {
365 // - [w, [w, v]] = [[v, w], w]
366 liebracket(&l3, &l1, &w, usejac);
367 if (nterms >= 6) {
368 // - [[v, [v, w]], w]
369 liebracket(&l4, &l2, &w, usejac);
370 // - [[w, [v, w]], v] == [[v, [v, w]], w]
371 }
372 }
373 }
374 }
375
376 // Evaluate BCH formula given all pre-computed terms and their respective weights
377 const double weight1[] = {0.0, 1.0, 1.0/2.0, 1.0/12.0, 1.0/12.0, 1.0/48.0, 1.0/48.0};
378 const double weight2[] = {1.0, 1.0, 1.0/2.0, 1.0/12.0, 1.0/12.0, 1.0/48.0, 1.0/48.0};
379
380 NaryVoxelFunction::VoxelWiseWeightedSum bch;
381 if (minus_v) bch._Weight = weight1;
382 else bch._Weight = weight2;
383
384 switch (nterms) {
385 case 1: { if (minus_v) { u = .0; } else { u = v; } break; }
386 case 2: { ParallelForEachScalar(v, w, u, bch); break; }
387 case 3: { ParallelForEachScalar(v, w, l1, u, bch); break; }
388 case 4: { ParallelForEachScalar(v, w, l1, l2, u, bch); break; }
389 case 5: { ParallelForEachScalar(v, w, l1, l2, l3, u, bch); break; }
390 case 6: { ParallelForEachScalar(v, w, l1, l2, l3, l4, u, bch); break; }
391 case 7: { ParallelForEachScalar(v, w, l1, l2, l3, l4, l4, u, bch); break; }
392 default:
393 cerr << "BSplineFreeFormTransformationSV::EvaluateBCHFormulaDense: Invalid number of terms " << nterms << endl;
394 exit(1);
395 };
396 MIRTK_DEBUG_TIMING(3, "evaluation of BCH formula (dense)");
397 }
398
399 // -----------------------------------------------------------------------------
400 void BSplineFreeFormTransformationSV
EvaluateBCHFormula(int nterms,CPImage & u,double tau,const CPImage & v,double eta,const CPImage & w,bool minus_v) const401 ::EvaluateBCHFormula(int nterms, CPImage &u,
402 double tau, const CPImage &v,
403 double eta, const CPImage &w,
404 bool minus_v) const
405 {
406 MIRTK_START_TIMING();
407 GenericImage<Vector> l1, l2, l3, l4;
408 const ImageAttributes &lattice = u.Attributes();
409
410 // Calculate required Lie brackets...
411 if (_LieDerivative) {
412 // ... using Lie derivative
413 if (nterms >= 3) {
414 // - [v, w]
415 l1.Initialize(lattice, 3);
416 ParallelForEachVoxel(SVFFDEvaluateLieBracket(this, tau, v, eta, w), lattice, l1);
417 ConvertToCubicBSplineCoefficients(l1);
418 if (nterms >= 4) {
419 // - [v, [v, w]]
420 l2.Initialize(lattice, 3);
421 ParallelForEachVoxel(SVFFDEvaluateLieBracket(this, tau, v, 1., l1), lattice, l2);
422 ConvertToCubicBSplineCoefficients(l2);
423 if (nterms >= 5) {
424 // - [w, [w, v]] = [[v, w], w]
425 l3.Initialize(lattice, 3);
426 ParallelForEachVoxel(SVFFDEvaluateLieBracket(this, 1., l1, eta, w), lattice, l3);
427 ConvertToCubicBSplineCoefficients(l3);
428 if (nterms >= 6) {
429 // - [[v, [v, w]], w]
430 l4.Initialize(lattice, 3);
431 ParallelForEachVoxel(SVFFDEvaluateLieBracket(this, 1., l2, eta, w), lattice, l4);
432 ConvertToCubicBSplineCoefficients(l4);
433 // - [[w, [v, w]], v] == [[v, [v, w]], w]
434 }
435 }
436 }
437 }
438 } else {
439 // ... using composition of vector fields
440 if (nterms >= 3) {
441 DifferenceOfCompositionLieBracketImageFilter3D<Vector> lb;
442 lb.Interpolation(Interpolation_CubicBSpline);
443 lb.Extrapolation(Extrapolation_NN);
444 lb.ComputeInterpolationCoefficients(false);
445 // - [v, w]
446 lb.Input(0, &v), lb.Scaling(0, tau);
447 lb.Input(1, &w), lb.Scaling(1, eta);
448 lb.Output(&l1);
449 lb.Run();
450 ConvertToCubicBSplineCoefficients(l1);
451 if (nterms >= 4) {
452 // - [v, [v, w]]
453 lb.Input(0, &v ), lb.Scaling(0, tau);
454 lb.Input(1, &l1), lb.Scaling(1, 1.0);
455 lb.Output(&l2);
456 lb.Run();
457 ConvertToCubicBSplineCoefficients(l2);
458 if (nterms >= 5) {
459 // - [[v, w], w]
460 lb.Input(0, &l1), lb.Scaling(0, 1.0);
461 lb.Input(1, &w ), lb.Scaling(1, eta);
462 lb.Output(&l3);
463 lb.Run();
464 ConvertToCubicBSplineCoefficients(l3);
465 if (nterms >= 6) {
466 // - [[v, [v, w]], w]
467 lb.Input(0, &l2), lb.Scaling(0, 1.0);
468 lb.Input(1, &w ), lb.Scaling(1, eta);
469 lb.Output(&l4);
470 lb.Run();
471 ConvertToCubicBSplineCoefficients(l4);
472 // - [[w, [v, w]], v] == [[v, [v, w]], w]
473 }
474 }
475 }
476 }
477 }
478
479 // Evaluate BCH formula given all pre-computed terms and their respective weights
480 const double weight1[] = {0.0, eta, 1.0/2.0, 1.0/12.0, 1.0/12.0, 1.0/48.0, 1.0/48.0};
481 const double weight2[] = {tau, eta, 1.0/2.0, 1.0/12.0, 1.0/12.0, 1.0/48.0, 1.0/48.0};
482
483 NaryVoxelFunction::VoxelWiseWeightedSum bch;
484 if (minus_v) bch._Weight = weight1;
485 else bch._Weight = weight2;
486
487 switch (nterms) {
488 case 1: { if (minus_v) { u = .0; } else { u = v; } break; }
489 case 2: { ParallelForEachScalar(v, w, u, bch); break; }
490 case 3: { ParallelForEachScalar(v, w, l1, u, bch); break; }
491 case 4: { ParallelForEachScalar(v, w, l1, l2, u, bch); break; }
492 case 5: { ParallelForEachScalar(v, w, l1, l2, l3, u, bch); break; }
493 case 6: { ParallelForEachScalar(v, w, l1, l2, l3, l4, u, bch); break; }
494 case 7: { ParallelForEachScalar(v, w, l1, l2, l3, l4, l4, u, bch); break; }
495 default:
496 cerr << "BSplineFreeFormTransformationSV::EvaluateBCHFormula: Invalid number of terms " << nterms << endl;
497 exit(1);
498 };
499 MIRTK_DEBUG_TIMING(3, "evaluation of BCH formula");
500 }
501
502 // =============================================================================
503 // Approximation/Interpolation
504 // =============================================================================
505
506 // -----------------------------------------------------------------------------
507 void BSplineFreeFormTransformationSV
ApproximateDOFs(const double * x,const double * y,const double * z,const double * t,const double * dx,const double * dy,const double * dz,int no)508 ::ApproximateDOFs(const double *x, const double *y, const double *z, const double *t,
509 const double *dx, const double *dy, const double *dz, int no)
510 {
511 // FIXME The initial approximation of the displacements on the control point grid
512 // smoothes the displacement field too much and hence introduces quite some
513 // error. Use the overloaded ApproximateAsNew(disp) method when possible.
514 // -as12312
515
516 // Approximate displacements at control points
517 double *rx = Allocate<double>(no);
518 double *ry = Allocate<double>(no);
519 double *rz = Allocate<double>(no);
520
521 memcpy(rx, dx, no * sizeof(double));
522 memcpy(ry, dy, no * sizeof(double));
523 memcpy(rz, dz, no * sizeof(double));
524
525 BSplineFreeFormTransformation3D::ApproximateDOFs(x, y, z, t, rx, ry, rz, no);
526
527 Deallocate(rx);
528 Deallocate(ry);
529 Deallocate(rz);
530
531 // Find stationary velocity field which approximates the displacements
532 GenericImage<double> disp(this->Attributes(), 3);
533 for (int k = 0; k < _z; ++k) {
534 for (int j = 0; j < _y; ++j) {
535 for (int i = 0; i < _x; ++i) {
536 disp(i, j, k, 0) = _CPImage(i, j, k)._x;
537 disp(i, j, k, 1) = _CPImage(i, j, k)._y;
538 disp(i, j, k, 2) = _CPImage(i, j, k)._z;
539 }
540 }
541 }
542 this->ApproximateAsNew(disp);
543 }
544
545 // -----------------------------------------------------------------------------
546 void BSplineFreeFormTransformationSV
ApproximateDOFsGradient(const double *,const double *,const double *,const double *,const double *,const double *,const double *,int,double *,double) const547 ::ApproximateDOFsGradient(const double *, const double *, const double *, const double *,
548 const double *, const double *, const double *, int,
549 double *, double) const
550 {
551 cerr << this->NameOfClass() << "::ApproximateDOFsGradient: Not implemented" << endl;
552 exit(1);
553 }
554
555 // -----------------------------------------------------------------------------
556 ImageAttributes BSplineFreeFormTransformationSV
ApproximationDomain(const ImageAttributes & attr,const Transformation * dof)557 ::ApproximationDomain(const ImageAttributes &attr, const Transformation *dof)
558 {
559 if (!dof) return attr;
560
561 ImageAttributes grid(attr);
562
563 // Ensure that for all target voxels the displacement can be recovered
564 // without requiring any extrapolation of the velocity field during
565 // computation of the trajectory (integration, i.e., exponentiation)
566 const Matrix i2w = grid.GetImageToWorldMatrix();
567 const Matrix w2i = grid.GetWorldToImageMatrix();
568
569 double margin_top = .0;
570 double margin_bottom = .0;
571 double margin_left = .0;
572 double margin_right = .0;
573 double margin_front = .0;
574 double margin_back = .0;
575
576 double x, y, z, wx, wy, wz;
577 for (int k = 0; k < attr._z; ++k)
578 for (int j = 0; j < attr._y; ++j)
579 for (int i = 0; i < attr._x; ++i) {
580 wx = i2w(0, 0) * i + i2w(0, 1) * j + i2w(0, 2) * k + i2w(0, 3);
581 wy = i2w(1, 0) * i + i2w(1, 1) * j + i2w(1, 2) * k + i2w(1, 3);
582 wz = i2w(2, 0) * i + i2w(2, 1) * j + i2w(2, 2) * k + i2w(2, 3);
583 dof->Transform(wx, wy, wz);
584 x = w2i(0, 0) * wx + w2i(0, 1) * wy + w2i(0, 2) * wz + w2i(0, 3);
585 y = w2i(1, 0) * wx + w2i(1, 1) * wy + w2i(1, 2) * wz + w2i(1, 3);
586 z = w2i(2, 0) * wx + w2i(2, 1) * wy + w2i(2, 2) * wz + w2i(2, 3);
587 if (x < 0 && -x > margin_left) margin_left = -x;
588 if (y < 0 && -y > margin_bottom) margin_bottom = -y;
589 if (z < 0 && -z > margin_front) margin_front = -z;
590 if (x >= grid._x - 1 && x - grid._x - 1 > margin_right) margin_right = x - grid._x - 1;
591 if (y >= grid._y - 1 && y - grid._y - 1 > margin_top) margin_top = y - grid._y - 1;
592 if (z >= grid._z - 1 && z - grid._z - 1 > margin_back) margin_back = z - grid._z - 1;
593 }
594
595 // Disregard small precision errors
596 const double eps = 1e-6;
597 margin_left = round(margin_left / eps) * eps;
598 margin_right = round(margin_right / eps) * eps;
599 margin_bottom = round(margin_bottom / eps) * eps;
600 margin_top = round(margin_top / eps) * eps;
601 margin_front = round(margin_front / eps) * eps;
602 margin_back = round(margin_back / eps) * eps;
603
604 // Account for inter-/extrapolation error on boundary of FFD lattice and
605 // therefore make lattice a bit bigger than otherwise needed
606 const double margin_safety = 1.5;
607 margin_left = ceil(margin_left * margin_safety);
608 margin_right = ceil(margin_right * margin_safety);
609 margin_bottom = ceil(margin_bottom * margin_safety);
610 margin_top = ceil(margin_top * margin_safety);
611 margin_front = ceil(margin_front * margin_safety);
612 margin_back = ceil(margin_back * margin_safety);
613
614 // Compute of offsets by which lattice origin must be moved such that
615 // the lattice origin is the center of the extended lattice again
616 const double ox = (margin_right - margin_left ) * grid._dx / 2.0;
617 const double oy = (margin_top - margin_bottom) * grid._dy / 2.0;
618 const double oz = (margin_back - margin_front ) * grid._dz / 2.0;
619
620 // Initialize free-form deformation (for extended image grid)
621 grid._x += static_cast<int>(margin_left + margin_right);
622 grid._y += static_cast<int>(margin_bottom + margin_top );
623 grid._z += static_cast<int>(margin_front + margin_back );
624 grid._xorigin += grid._xaxis[0] * ox + grid._yaxis[0] * oy + grid._zaxis[0] * oz;
625 grid._yorigin += grid._xaxis[1] * ox + grid._yaxis[1] * oy + grid._zaxis[1] * oz;
626 grid._zorigin += grid._xaxis[2] * ox + grid._yaxis[2] * oy + grid._zaxis[2] * oz;
627
628 return grid;
629 }
630
631 // -----------------------------------------------------------------------------
632 double BSplineFreeFormTransformationSV
ApproximateAsNew(const ImageAttributes & domain,const Transformation * dof,int niter,double max_error)633 ::ApproximateAsNew(const ImageAttributes &domain, const Transformation *dof,
634 int niter, double max_error)
635 {
636 const HomogeneousTransformation *lin = nullptr;
637 const BSplineFreeFormTransformationSV *svffd = nullptr;
638 const MultiLevelTransformation *mffd = nullptr;
639
640 (lin = dynamic_cast<const HomogeneousTransformation *>(dof)) ||
641 (svffd = dynamic_cast<const BSplineFreeFormTransformationSV *>(dof)) ||
642 (mffd = dynamic_cast<const MultiLevelTransformation *>(dof));
643
644 if (mffd) {
645 if (mffd->NumberOfLevels() == 0) {
646 lin = mffd->GetGlobalTransformation();
647 } else if (mffd->NumberOfLevels() == 1) {
648 if (mffd->GetGlobalTransformation()->IsIdentity()) {
649 svffd = dynamic_cast<const BSplineFreeFormTransformationSV *>(mffd->GetLocalTransformation(0));
650 }
651 }
652 }
653
654 // Approximate any other transformation using the base class implementation
655 // which simply evaluates the displacement of the transformation at each
656 // control point and then calls Interpolate in order to interpolate these
657 // control point displacements
658 if (!lin && !svffd) {
659 return BSplineFreeFormTransformation3D::ApproximateAsNew(domain, dof, niter, max_error);
660 }
661 // Check if input is the identity transformation
662 if ((lin && lin->IsIdentity()) || (svffd && svffd->IsIdentity())) {
663 _CPImage = 0.;
664 return 0.;
665 }
666
667 double error = inf;
668 if (niter >= 0) {
669 // Compute velocities at control points using log map of affine matrix
670 if (lin) {
671 EvaluateGlobalSVFFD logA(logm(lin->GetMatrix()), &_CPImage);
672 ParallelForEachVoxel(_CPImage.Attributes(), &_CPImage, logA);
673 // Evaluate velocities of other SV FFD at control points of this SV FFD
674 } else {
675 ParallelForEachVoxel(EvaluateBSplineSVFFD(svffd, &_CPImage), _attr, _CPImage);
676 }
677 // Convert velocities to B-spline coefficients
678 ConvertToSplineCoefficients(3, _CPImage);
679 // Evaluate approximation error
680 error = EvaluateRMSError(domain, dof);
681 }
682
683 return error;
684 }
685
686 // -----------------------------------------------------------------------------
687 //double BSplineFreeFormTransformationSV
688 //::ApproximateFFDAsNew(GenericImage<double> &disp, double T, int niter)
689 //{
690 // GenericImage<double> dv;
691 // WorldCoordinatesImage wc;
692 // disp.ImageToWorld(wc);
693 // for (int iter = 0; iter < niter; ++iter) {
694 // // Interpolate velocity field
695 // const ImageAttributes &grid = d.GetImageAttributes();
696 // v->Initialize(grid, d.GetT());
697 // switch (d.GetT()) {
698 // case 2: { ParallelForEachVoxel(EvaluateBSplineSVFFD2D(this, v), grid, v); break; }
699 // case 3: { ParallelForEachVoxel(EvaluateBSplineSVFFD3D(this, v), grid, v); break; }
700 // default: {
701 // cerr << "BSplineFreeFormTransformationSV::ScalingAndSquaring: Vector field must have 2 or 3 components (_t)" << endl;
702 // exit(1);
703 // }
704 // }
705 // // Exponentiate velocity field
706 // VelocityToDisplacementFieldSS<double> exp;
707 // exp.T(T);
708 // exp.NumberOfSteps (NumberOfStepsForIntervalLength(T));
709 // exp.MaxScaledVelocity(_MaxScaledVelocity);
710 // exp.Interpolation (Interpolation_BSpline);
711 // exp.Upsample (false); // better, but too expensive
712 // exp.ExternalCache (cache); // avoid frequent allocation/deallocation
713 // exp.SetInput (0, v); // velocity field to be exponentiated
714 // exp.SetInput (1, &d); // input displacement field (may be zero)
715 // exp.SetOutput (v); // result is exp(v) o d
716 // exp.Run();
717 // d.CopyFrom(*v);
718 // // dv = exp(-v) o disp
719 // dv = disp;
720 // ScalingAndSquaring(dv, T, &wc);
721 // // Approximate dv by B-spline coefficients at control points
722 // BSplineFreeFormTransformation3D ffd(this->GetFFDAttributes());
723 // ffd.Approximate(wc.GetPointerToVoxels(0, 0, 0, 0),
724 // wc.GetPointerToVoxels(0, 0, 0, 1),
725 // wc.GetPointerToVoxels(0, 0, 0, 2),
726 // dv.GetPointerToVoxels(0, 0, 0, 0),
727 // dv.GetPointerToVoxels(0, 0, 0, 1),
728 // dv.GetPointerToVoxels(0, 0, 0, 2));
729 // dv.Initialize(_x, _y, _z, 1, 3);
730 // for (int k = 0; k < _z; ++k) {
731 // for (int j = 0; j < _y; ++j) {
732 // for (int i = 0; i < _x; ++i) {
733 // dv(i, j, k, 0) = ffd._data[k][j][i]._x;
734 // dv(i, j, k, 1) = ffd._data[k][j][i]._x;
735 // dv(i, j, k, 2) = ffd._data[k][j][i]._x;
736 // }
737 // }
738 // }
739 // // Calculate required Lie brackets
740 // if (_NumberOfBCHTerms >= 3) {
741 // DifferenceOfCompositionLieBracketImageFilter3D<double> lb;
742 // lb.Interpolation(Interpolation_CubicBSpline);
743 // lb.Extrapolation(Extrapolation_NN);
744 // lb.ComputeInterpolationCoefficients(false);
745 // // - [v, dv]
746 // lb.SetInput (0, v);
747 // lb.SetInput (1, d);
748 // lb.SetOutput(l1);
749 // lb.Run();
750 // ConvertToCubicBSplineCoefficients(*l1);
751 // if (_NumberOfBCHTerms >= 4) {
752 // // - [v, [v, dv]]
753 // lb.SetInput (0, v);
754 // lb.SetInput (1, l1);
755 // lb.SetOutput(l2);
756 // lb.Run();
757 // ConvertToCubicBSplineCoefficients(*l2);
758 // if (_NumberOfBCHTerms >= 5) {
759 // // - [dv, [v, dv]]
760 // lb.SetInput (0, d);
761 // lb.SetInput (1, l1);
762 // lb.SetOutput(l3);
763 // lb.Run();
764 // ConvertToCubicBSplineCoefficients(*l3);
765 // if (_NumberOfBCHTerms >= 6) {
766 // // - [dv, [v, [v, dv]]]
767 // lb.SetInput (0, d);
768 // lb.SetInput (1, l2);
769 // lb.SetOutput(l4);
770 // lb.Run();
771 // ConvertToCubicBSplineCoefficients(*l4);
772 // }
773 // }
774 // }
775 // }
776 // // Update velocity coefficients using BCH formula
777 // NaryVoxelFunction::EvaluateBCHFormula bch;
778 // switch (_NumberOfBCHTerms) {
779 // case 2: { /* nothing else to do */ break; }
780 // case 3: { ParallelForEachScalar(v, dv, l1, v, bch); break; }
781 // case 4: { ParallelForEachScalar(v, dv, l1, l2, v, bch); break; }
782 // case 5: { ParallelForEachScalar(v, dv, l1, l2, l3, v, bch); break; }
783 // case 6: { ParallelForEachScalar(v, dv, l1, l2, l3, l4, v, bch); break; }
784 // };
785 // }
786 //}
787
788 // -----------------------------------------------------------------------------
789 double BSplineFreeFormTransformationSV
ApproximateAsNew(GenericImage<double> & disp,int niter,double max_error)790 ::ApproximateAsNew(GenericImage<double> &disp, int niter, double max_error)
791 {
792 return this->ApproximateAsNew(disp, false, 3, niter * 8);
793 }
794
795 // -----------------------------------------------------------------------------
796 double BSplineFreeFormTransformationSV
ApproximateAsNew(GenericImage<double> & disp,bool smooth,int nterms,int niter)797 ::ApproximateAsNew(GenericImage<double> &disp, bool smooth, int nterms, int niter)
798 {
799 // TODO: Refactor/review implementation again after update of velocities
800 // from gradient is now implemented and working.
801
802 ImageAttributes grid = this->Attributes();
803 grid._t = 3;
804 grid._dt = disp.GetTSize(); // ignore difference in _dt
805
806 // Sample displacement field at control points using linear interpolation
807 GenericImage<double> *d = NULL;
808
809 if (disp.Attributes() == grid) {
810 d = &disp;
811 } else {
812 double x, y, z, vec[3] = {.0, .0, .0};
813
814 UniquePtr<InterpolateImageFunction> f;
815 f.reset(InterpolateImageFunction::New(Interpolation_Linear, Extrapolation_NN, &disp));
816 f->Input(&disp);
817 f->Initialize();
818
819 d = new GenericImage<double>(grid);
820
821 for (int k = 0; k < grid._z; ++k)
822 for (int j = 0; j < grid._y; ++j)
823 for (int i = 0; i < grid._x; ++i) {
824 x = i, y = j, z = k;
825 d->ImageToWorld(x, y, z);
826 disp.WorldToImage(x, y, z);
827 f->Evaluate(vec, x, y, z);
828 d->Put(i, j, k, 0, vec[0]);
829 d->Put(i, j, k, 1, vec[1]);
830 d->Put(i, j, k, 2, vec[2]);
831 }
832 }
833
834 // Compute stationary velocity field at control points
835 GenericImage<double> v;
836 DisplacementToVelocityFieldBCH<double> dtov;
837
838 dtov.Input (d);
839 dtov.Output(&v);
840
841 dtov.UpperIntegrationLimit(UpperIntegrationLimit(0, 1));
842 dtov.NumberOfIterations(niter);
843 dtov.NumberOfTerms (nterms);
844 dtov.NumberOfSteps (NumberOfStepsForIntervalLength(dtov.UpperIntegrationLimit()));
845 dtov.SmoothVelocities (smooth);
846
847 dtov.Run();
848
849 // Free temporary displacement field
850 if (d != &disp) {
851 delete d;
852 d = NULL;
853 }
854
855 // Interpolate velocities by B-spline function
856 BSplineFreeFormTransformation3D::Interpolate(v.Data(0, 0, 0, 0),
857 v.Data(0, 0, 0, 1),
858 v.Data(0, 0, 0, 2));
859
860 // Evaluate RMS of approximation error
861 double error = .0;
862
863 v = .0;
864 this->Displacement(v, 0, 1);
865
866 for (int k = 0; k < disp.Z(); ++k)
867 for (int j = 0; j < disp.Y(); ++j)
868 for (int i = 0; i < disp.X(); ++i) {
869 disp(i, j, k, 0) -= v(i, j, k, 0);
870 disp(i, j, k, 1) -= v(i, j, k, 1);
871 disp(i, j, k, 2) -= v(i, j, k, 2);
872 error += sqrt(disp(i, j ,k, 0) * disp(i, j, k, 0) +
873 disp(i, j ,k, 1) * disp(i, j, k, 1) +
874 disp(i, j ,k, 2) * disp(i, j, k, 2));
875 }
876 error /= static_cast<double>(disp.X() * disp.Y() * disp.Z());
877
878 return error;
879 }
880
881 // -----------------------------------------------------------------------------
882 void BSplineFreeFormTransformationSV
Interpolate(const double *,const double *,const double *)883 ::Interpolate(const double *, const double *, const double *)
884 {
885 cerr << this->NameOfClass() << "::Interpolate: Not implemented" << endl;
886 exit(1);
887 }
888
889 // -----------------------------------------------------------------------------
890 void BSplineFreeFormTransformationSV
InterpolateVelocities(const double * vx,const double * vy,const double * vz)891 ::InterpolateVelocities(const double *vx, const double *vy, const double *vz)
892 {
893 BSplineFreeFormTransformation3D::Interpolate(vx, vy, vz);
894 }
895
896 // -----------------------------------------------------------------------------
897 double BSplineFreeFormTransformationSV
ApproximateVelocitiesAsNew(GenericImage<double> & v)898 ::ApproximateVelocitiesAsNew(GenericImage<double> &v)
899 {
900 double error = 0.;
901
902 auto * const vx = v.Data(0, 0, 0, 0);
903 auto * const vy = v.Data(0, 0, 0, 1);
904 auto * const vz = v.Data(0, 0, 0, 2);
905
906 if (_attr.EqualInSpace(v.Attributes())) {
907
908 BSplineFreeFormTransformation3D::Interpolate(vx, vy, vz);
909
910 } else {
911
912 // Compute world coordinates of lattice points
913 const int no = v.NumberOfSpatialVoxels();
914 double *x = Allocate<double>(no);
915 double *y = Allocate<double>(no);
916 double *z = Allocate<double>(no);
917 double *t = Allocate<double>(no);
918 v.Attributes().LatticeToWorld(x, y, z, t);
919
920 // Approximate velocities
921 BSplineFreeFormTransformation3D::ApproximateDOFs(x, y, z, t, vx, vy, vz, no);
922
923 // Evaluate error of approximation
924 GenericImage<double> u(v.Attributes(), 3);
925 this->Velocity(u);
926 for (int k = 0; k < v.Z(); ++k)
927 for (int j = 0; j < v.Y(); ++j)
928 for (int i = 0; i < v.X(); ++i) {
929 const int idx = v.VoxelToIndex(i, j, k);
930 vx[idx] -= u(i, j, k, 0);
931 vy[idx] -= u(i, j, k, 1);
932 vz[idx] -= u(i, j, k, 2);
933 error += pow(vx[idx], 2) + pow(vy[idx], 2) + pow(vz[idx], 2);
934 }
935 error = sqrt(error / no);
936
937 // Free memory
938 Deallocate(x);
939 Deallocate(y);
940 Deallocate(z);
941 Deallocate(t);
942 }
943
944 return error;
945 }
946
947 // -----------------------------------------------------------------------------
CombineWith(const Transformation * dof)948 void BSplineFreeFormTransformationSV::CombineWith(const Transformation *dof)
949 {
950 // Convert transformation into SV FFD
951 const BSplineFreeFormTransformationSV *svffd = NULL;
952 svffd = dynamic_cast<const BSplineFreeFormTransformationSV *>(dof);
953 if (!svffd) {
954 BSplineFreeFormTransformationSV *tmp;
955 tmp = new BSplineFreeFormTransformationSV();
956 tmp->Initialize(this->Attributes());
957 tmp->ApproximateAsNew(dof);
958 svffd = tmp;
959 }
960 // Compute coefficients of composite SV FFD using BCH formula
961 this->CombineWith(svffd);
962 // Delete temporary SV FFD
963 if (svffd != dof) delete svffd;
964 }
965
966 // -----------------------------------------------------------------------------
CombineWith(const BSplineFreeFormTransformationSV * svffd)967 void BSplineFreeFormTransformationSV::CombineWith(const BSplineFreeFormTransformationSV *svffd)
968 {
969 EvaluateBCHFormula(_NumberOfBCHTerms, _CPImage, _T, _CPImage, svffd->T(), svffd->_CPImage);
970 if (_T < 0.) this->Invert();
971 }
972
973 // -----------------------------------------------------------------------------
ScaleVelocities(double e)974 void BSplineFreeFormTransformationSV::ScaleVelocities(double e)
975 {
976 _CPImage *= e;
977 }
978
979 // -----------------------------------------------------------------------------
Invert()980 void BSplineFreeFormTransformationSV::Invert()
981 {
982 this->ScaleVelocities(-1.);
983 }
984
985 // =============================================================================
986 // Parameters (non-DoFs)
987 // =============================================================================
988
989 // -----------------------------------------------------------------------------
Set(const char * name,const char * value)990 bool BSplineFreeFormTransformationSV::Set(const char *name, const char *value)
991 {
992 if (strcmp(name, "Cross-sectional time interval") == 0 ||
993 strcmp(name, "Cross sectional time interval") == 0) {
994 return FromString(value, _T);
995 } else if (strcmp(name, "Time unit of integration interval") == 0) {
996 return FromString(value, _TimeUnit);
997 } else if (strcmp(name, "No. of integration steps") == 0 ||
998 strcmp(name, "Number of integration steps") == 0) {
999 return FromString(value, _NumberOfSteps) && _NumberOfSteps > 0;
1000 } else if (strcmp(name, "No. of squaring steps") == 0 ||
1001 strcmp(name, "Number of squaring steps") == 0) {
1002 if (!FromString(value, _NumberOfSteps) || _NumberOfSteps <= 0) return false;
1003 _NumberOfSteps = static_cast<int>(pow(2, _NumberOfSteps));
1004 if (_IntegrationMethod != FFDIM_SS && _IntegrationMethod != FFDIM_FastSS) {
1005 _IntegrationMethod = FFDIM_FastSS;
1006 }
1007 } else if (strcmp(name, "Maximum scaled velocity") == 0) {
1008 return FromString(value, _MaxScaledVelocity);
1009 } else if (strcmp(name, "Use Lie derivative") == 0) {
1010 return FromString(value, _LieDerivative);
1011 } else if (strcmp(name, "Use dense BCH lattice") == 0) {
1012 return FromString(value, _UseDenseBCHGrid);
1013 } else if (strcmp(name, "No. of BCH terms") == 0 ||
1014 strcmp(name, "Number of BCH terms") == 0) {
1015 return FromString(value, _NumberOfBCHTerms) && _NumberOfBCHTerms <= 6;
1016 } else if (strcmp(name, "Integration method") == 0) {
1017 return FromString(value, _IntegrationMethod) && _IntegrationMethod != FFDIM_Unknown;
1018 // deprecated parameters
1019 } else if (strcmp(name, "Use scaling and squaring") == 0) {
1020 bool useSS = false;
1021 if (!FromString(value, useSS)) return false;
1022 if (useSS) {
1023 if (_IntegrationMethod != FFDIM_SS && _IntegrationMethod != FFDIM_FastSS) {
1024 _IntegrationMethod = FFDIM_FastSS;
1025 }
1026 } else {
1027 if (_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) {
1028 _IntegrationMethod = FFDIM_RKE1;
1029 }
1030 }
1031 return true;
1032 } else if (strcmp(name, "Fast scaling and squaring") == 0) {
1033 bool fastSS;
1034 if (!FromString(value, fastSS)) return false;
1035 if (_IntegrationMethod == FFDIM_SS && fastSS) _IntegrationMethod = FFDIM_FastSS;
1036 return true;
1037 }
1038 return BSplineFreeFormTransformation3D::Set(name, value);
1039 }
1040
1041 // -----------------------------------------------------------------------------
Parameter() const1042 ParameterList BSplineFreeFormTransformationSV::Parameter() const
1043 {
1044 ParameterList params = BSplineFreeFormTransformation3D::Parameter();
1045 Insert(params, "Integration method", _IntegrationMethod);
1046 Insert(params, "Cross-sectional time interval", _T);
1047 Insert(params, "Time unit of integration interval", _TimeUnit);
1048 Insert(params, "No. of integration steps", _NumberOfSteps);
1049 Insert(params, "Maximum scaled velocity", _MaxScaledVelocity);
1050 Insert(params, "Use Lie derivative", _LieDerivative);
1051 Insert(params, "Use dense BCH lattice", _UseDenseBCHGrid);
1052 Insert(params, "No. of BCH terms", _NumberOfBCHTerms);
1053 return params;
1054 }
1055
1056 // =============================================================================
1057 // Point transformation
1058 // =============================================================================
1059
1060 // -----------------------------------------------------------------------------
IntegrateVelocities(double & x,double & y,double & z,double T) const1061 void BSplineFreeFormTransformationSV::IntegrateVelocities(double &x, double &y, double &z, double T) const
1062 {
1063 const double dt = StepLengthForIntervalLength(T);
1064 if (dt) {
1065 if (_IntegrationMethod == FFDIM_FastSS ||
1066 _IntegrationMethod == FFDIM_SS ||
1067 _IntegrationMethod == FFDIM_RKE1) RKE1 ::Transform(this, x, y, z, .0, T, dt);
1068 else if (_IntegrationMethod == FFDIM_RKE2) RKE2 ::Transform(this, x, y, z, .0, T, dt);
1069 else if (_IntegrationMethod == FFDIM_RKH2) RKH2 ::Transform(this, x, y, z, .0, T, dt);
1070 else if (_IntegrationMethod == FFDIM_RK4) RK4 ::Transform(this, x, y, z, .0, T, dt);
1071 else if (_IntegrationMethod == FFDIM_RKEH12) RKEH12::Transform(this, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1072 else if (_IntegrationMethod == FFDIM_RKBS23) RKBS23::Transform(this, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1073 else if (_IntegrationMethod == FFDIM_RKF45) RKF45 ::Transform(this, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1074 else if (_IntegrationMethod == FFDIM_RKCK45) RKCK45::Transform(this, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1075 else if (_IntegrationMethod == FFDIM_RKDP45) RKDP45::Transform(this, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1076 else {
1077 cerr << "BSplineFreeFormTransformationSV::IntegrateVelocities: Unknown integration method: " << _IntegrationMethod << endl;
1078 exit(1);
1079 }
1080 }
1081 }
1082
1083 // -----------------------------------------------------------------------------
1084 void BSplineFreeFormTransformationSV
LocalTransform(double & x,double & y,double & z,double t,double t0) const1085 ::LocalTransform(double &x, double &y, double &z, double t, double t0) const
1086 {
1087 IntegrateVelocities(x, y, z, + UpperIntegrationLimit(t, t0));
1088 }
1089
1090 // -----------------------------------------------------------------------------
1091 bool BSplineFreeFormTransformationSV
LocalInverse(double & x,double & y,double & z,double t,double t0) const1092 ::LocalInverse(double &x, double &y, double &z, double t, double t0) const
1093 {
1094 IntegrateVelocities(x, y, z, - UpperIntegrationLimit(t, t0));
1095 return true;
1096 }
1097
1098 // -----------------------------------------------------------------------------
1099 template <class VoxelType>
1100 void BSplineFreeFormTransformationSV
ScalingAndSquaring(GenericImage<VoxelType> * d,double T,const WorldCoordsImage * wc) const1101 ::ScalingAndSquaring(GenericImage<VoxelType> *d,
1102 double T, const WorldCoordsImage *wc) const
1103 {
1104 ScalingAndSquaring<VoxelType>(d->Attributes(), d, nullptr, nullptr, nullptr, T, wc);
1105 }
1106
1107 template void BSplineFreeFormTransformationSV
1108 ::ScalingAndSquaring(GenericImage<float> *d,
1109 double T, const WorldCoordsImage *) const;
1110
1111 template void BSplineFreeFormTransformationSV
1112 ::ScalingAndSquaring(GenericImage<double> *d,
1113 double T, const WorldCoordsImage *) const;
1114
1115 // -----------------------------------------------------------------------------
1116 template <class VoxelType>
1117 void BSplineFreeFormTransformationSV
ScalingAndSquaring(const ImageAttributes & a,GenericImage<VoxelType> * d,GenericImage<VoxelType> * dx,GenericImage<VoxelType> * dj,GenericImage<VoxelType> * lj,double T,const WorldCoordsImage *) const1118 ::ScalingAndSquaring(const ImageAttributes &a,
1119 GenericImage<VoxelType> *d,
1120 GenericImage<VoxelType> *dx,
1121 GenericImage<VoxelType> *dj,
1122 GenericImage<VoxelType> *lj,
1123 double T, const WorldCoordsImage *) const
1124 {
1125 // Attributes of output images
1126 ImageAttributes attr(a);
1127 if (!attr) {
1128 if (d ) attr = d ->Attributes();
1129 else if (dx) attr = dx->Attributes();
1130 else if (dj) attr = dj->Attributes();
1131 else if (lj) attr = lj->Attributes();
1132 }
1133 attr._t = 1, attr._dt = .0;
1134 if (!attr) return;
1135 // Copy input displacement field
1136 UniquePtr<GenericImage<VoxelType> > din(d ? new GenericImage<VoxelType>(*d) : nullptr);
1137 // TODO: The runtime of the ScalingAndSquaring filter has been greatly improved
1138 // to almost match the old VelocityToDisplacementFieldSS implementation.
1139 // However, the latter is still about 15% faster...
1140 if (d && !dx && !dj && !lj) {
1141 GenericImage<VoxelType> v;
1142 if (_IntegrationMethod == FFDIM_FastSS) {
1143 v.Initialize(this->Attributes(), 3);
1144 ParallelForEachVoxel(EvaluateBSplineSVFFD3D(this, &v), this->Attributes(), v);
1145 } else {
1146 v.Initialize(attr, 3);
1147 ParallelForEachVoxel(EvaluateBSplineSVFFD3D(this, &v), attr, v);
1148 }
1149 // Exponentiate velocity field
1150 VelocityToDisplacementFieldSS<VoxelType> exp;
1151 exp.UpperIntegrationLimit(T);
1152 exp.NumberOfSteps(NumberOfStepsForIntervalLength(T));
1153 exp.MaxScaledVelocity(static_cast<VoxelType>(_MaxScaledVelocity));
1154 exp.Interpolation(Interpolation_Linear);
1155 exp.Upsample(false); // better, but too expensive
1156 exp.Input(0, &v); // velocity field to be exponentiated
1157 exp.Input(1, din.get()); // input displacement field (may be nullptr)
1158 exp.Output(d); // result is exp(v) o d
1159 exp.Run();
1160 } else {
1161 // Copy B-spline coefficients of velocity field
1162 GenericImage<VoxelType> v(this->Attributes(), 3);
1163 VoxelType *vx = v.Data(0, 0, 0, 0);
1164 VoxelType *vy = v.Data(0, 0, 0, 1);
1165 VoxelType *vz = v.Data(0, 0, 0, 2);
1166 const Vector *vp = _CPImage.Data();
1167 for (int idx = 0; idx < _CPImage.NumberOfVoxels(); ++idx, ++vx, ++vy, ++vz, ++vp) {
1168 *vx = static_cast<VoxelType>(vp->_x);
1169 *vy = static_cast<VoxelType>(vp->_y);
1170 *vz = static_cast<VoxelType>(vp->_z);
1171 }
1172 // Exponentiate velocity field
1173 mirtk::ScalingAndSquaring<VoxelType> exp;
1174 exp.UpperIntegrationLimit(T);
1175 exp.NumberOfSteps(NumberOfStepsForIntervalLength(T));
1176 exp.MaxScaledVelocity(_MaxScaledVelocity);
1177 exp.InterimAttributes(_IntegrationMethod == FFDIM_FastSS ? this->Attributes() : attr);
1178 exp.OutputAttributes(attr);
1179 exp.Upsample(false); // better, but too expensive
1180 exp.ComputeInterpolationCoefficients(false); // v contains B-spline coefficients
1181 exp.InputVelocity(&v); // velocity field to be exponentiated
1182 exp.InputDisplacement(din.get()); // input displacement field (may be zero)
1183 exp.OutputDisplacement(d); // i.e., d = exp(v) o din
1184 exp.OutputJacobian(dx); // i.e., Jacobian
1185 exp.OutputDetJacobian(dj); // i.e., det(Jacobian)
1186 exp.OutputLogJacobian(lj); // i.e., log(det(Jacobian)
1187 exp.Run();
1188 }
1189 }
1190
1191 template void BSplineFreeFormTransformationSV
1192 ::ScalingAndSquaring(const ImageAttributes &,
1193 GenericImage<float> *,
1194 GenericImage<float> *,
1195 GenericImage<float> *,
1196 GenericImage<float> *,
1197 double, const WorldCoordsImage *) const;
1198
1199 template void BSplineFreeFormTransformationSV
1200 ::ScalingAndSquaring(const ImageAttributes &,
1201 GenericImage<double> *,
1202 GenericImage<double> *,
1203 GenericImage<double> *,
1204 GenericImage<double> *,
1205 double, const WorldCoordsImage *) const;
1206
1207 // -----------------------------------------------------------------------------
1208 void BSplineFreeFormTransformationSV
Displacement(GenericImage<float> & d,double t,double t0,const WorldCoordsImage * wc) const1209 ::Displacement(GenericImage<float> &d, double t, double t0, const WorldCoordsImage *wc) const
1210 {
1211 double T;
1212 if ((T = + UpperIntegrationLimit(t, t0))) {
1213 MIRTK_START_TIMING();
1214 // Use scaling and squaring method to efficiently compute displacements when possible
1215 if ((_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) &&
1216 ((_z <= 1 && d.Z() <= 1) || (_z > 1 && d.Z() > 1))) {
1217 ScalingAndSquaring(&d, T, wc);
1218 // Evaluate transformation at each voxel separately using numerical integration
1219 } else {
1220 Transformation::Displacement(d, t, t0, wc);
1221 }
1222 MIRTK_DEBUG_TIMING(3, "computation of exp(" << T << "*v)");
1223 }
1224 }
1225
1226 // -----------------------------------------------------------------------------
1227 void BSplineFreeFormTransformationSV
Displacement(GenericImage<double> & d,double t,double t0,const WorldCoordsImage * wc) const1228 ::Displacement(GenericImage<double> &d, double t, double t0, const WorldCoordsImage *wc) const
1229 {
1230 double T;
1231 if ((T = + UpperIntegrationLimit(t, t0))) {
1232 MIRTK_START_TIMING();
1233 // Use scaling and squaring method to efficiently compute displacements when possible
1234 if ((_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) &&
1235 ((_z <= 1 && d.Z() <= 1) || (_z > 1 && d.Z() > 1))) {
1236 ScalingAndSquaring(&d, T, wc);
1237 // Evaluate transformation at each voxel separately using numerical integration
1238 } else {
1239 Transformation::Displacement(d, t, t0, wc);
1240 }
1241 MIRTK_DEBUG_TIMING(3, "computation of exp(" << T << "*v)");
1242 }
1243 }
1244
1245 // -----------------------------------------------------------------------------
1246 int BSplineFreeFormTransformationSV
InverseDisplacement(GenericImage<float> & d,double t,double t0,const WorldCoordsImage * wc) const1247 ::InverseDisplacement(GenericImage<float> &d, double t, double t0, const WorldCoordsImage *wc) const
1248 {
1249 double T;
1250 if ((T = - UpperIntegrationLimit(t, t0))) {
1251 MIRTK_START_TIMING();
1252 // Use scaling and squaring method to efficiently compute displacements when possible
1253 if ((_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) &&
1254 ((_z <= 1 && d.Z() <= 1) || (_z > 1 && d.Z() > 1))) {
1255 ScalingAndSquaring(&d, T, wc);
1256 // Evaluate transformation at each voxel separately using numerical integration
1257 } else {
1258 Transformation::InverseDisplacement(d, t, t0, wc);
1259 }
1260 MIRTK_DEBUG_TIMING(3, "computation of exp(" << T << "*v)");
1261 }
1262 return 0;
1263 }
1264
1265 // -----------------------------------------------------------------------------
1266 int BSplineFreeFormTransformationSV
InverseDisplacement(GenericImage<double> & d,double t,double t0,const WorldCoordsImage * wc) const1267 ::InverseDisplacement(GenericImage<double> &d, double t, double t0, const WorldCoordsImage *wc) const
1268 {
1269 double T;
1270 if ((T = - UpperIntegrationLimit(t, t0))) {
1271 MIRTK_START_TIMING();
1272 // Use scaling and squaring method to efficiently compute displacements when possible
1273 if ((_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) &&
1274 ((_z <= 1 && d.Z() <= 1) || (_z > 1 && d.Z() > 1))) {
1275 ScalingAndSquaring(&d, T, wc);
1276 // Evaluate transformation at each voxel separately uing numerical integration
1277 } else {
1278 Transformation::InverseDisplacement(d, t, t0, wc);
1279 }
1280 MIRTK_DEBUG_TIMING(3, "computation of exp(" << T << "*v)");
1281 }
1282 return 0;
1283 }
1284
1285 // =============================================================================
1286 // Derivatives
1287 // =============================================================================
1288
1289 // -----------------------------------------------------------------------------
LocalJacobian(Matrix & jac,double x,double y,double z,double t,double t0) const1290 void BSplineFreeFormTransformationSV::LocalJacobian(Matrix &jac, double x, double y, double z, double t, double t0) const
1291 {
1292 jac.Initialize(3, 3);
1293 jac.Ident();
1294 double dt, T;
1295 if ((dt = StepLengthForIntervalLength(T = UpperIntegrationLimit(t, t0)))) {
1296 if (_IntegrationMethod == FFDIM_SS ||
1297 _IntegrationMethod == FFDIM_FastSS ||
1298 _IntegrationMethod == FFDIM_RKE1) RKE1 ::Jacobian(this, jac, x, y, z, .0, T, dt);
1299 else if (_IntegrationMethod == FFDIM_RKE2) RKE2 ::Jacobian(this, jac, x, y, z, .0, T, dt);
1300 else if (_IntegrationMethod == FFDIM_RKH2) RKH2 ::Jacobian(this, jac, x, y, z, .0, T, dt);
1301 else if (_IntegrationMethod == FFDIM_RK4) RK4 ::Jacobian(this, jac, x, y, z, .0, T, dt);
1302 else if (_IntegrationMethod == FFDIM_RKEH12) RKEH12::Jacobian(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1303 else if (_IntegrationMethod == FFDIM_RKBS23) RKBS23::Jacobian(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1304 else if (_IntegrationMethod == FFDIM_RKF45) RKF45 ::Jacobian(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1305 else if (_IntegrationMethod == FFDIM_RKCK45) RKCK45::Jacobian(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1306 else if (_IntegrationMethod == FFDIM_RKDP45) RKDP45::Jacobian(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1307 else {
1308 cerr << "BSplineFreeFormTransformationSV::Jacobian: Unknown integration method: " << _IntegrationMethod << endl;
1309 exit(1);
1310 }
1311 }
1312 }
1313
1314 // -----------------------------------------------------------------------------
LocalHessian(Matrix[3],double,double,double,double,double) const1315 void BSplineFreeFormTransformationSV::LocalHessian(Matrix [3], double, double, double, double, double) const
1316 {
1317 cerr << this->NameOfClass() << "::LocalHessian: Not implemented" << endl;
1318 exit(1);
1319 }
1320
1321 // -----------------------------------------------------------------------------
JacobianDOFs(Matrix & jac,int cp,double x,double y,double z,double t,double t0) const1322 void BSplineFreeFormTransformationSV::JacobianDOFs(Matrix &jac, int cp, double x, double y, double z, double t, double t0) const
1323 {
1324 jac.Initialize(3, 3);
1325 jac.Ident();
1326 double dt, T;
1327 if ((dt = StepLengthForIntervalLength(T = UpperIntegrationLimit(t, t0)))) {
1328 int ci, cj, ck, cl;
1329 this->IndexToLattice(cp, ci, cj, ck, cl);
1330 if (_IntegrationMethod == FFDIM_SS ||
1331 _IntegrationMethod == FFDIM_FastSS ||
1332 _IntegrationMethod == FFDIM_RKE1) RKE1 ::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, dt);
1333 else if (_IntegrationMethod == FFDIM_RKE2) RKE2 ::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, dt);
1334 else if (_IntegrationMethod == FFDIM_RKH2) RKH2 ::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, dt);
1335 else if (_IntegrationMethod == FFDIM_RK4) RK4 ::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, dt);
1336 else if (_IntegrationMethod == FFDIM_RKEH12) RKEH12::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1337 else if (_IntegrationMethod == FFDIM_RKBS23) RKBS23::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1338 else if (_IntegrationMethod == FFDIM_RKF45) RKF45 ::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1339 else if (_IntegrationMethod == FFDIM_RKCK45) RKCK45::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1340 else if (_IntegrationMethod == FFDIM_RKDP45) RKDP45::JacobianDOFs(this, jac, ci, cj, ck, cl, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1341 else {
1342 cerr << "BSplineFreeFormTransformationSV::JacobianDOFs: Unknown integration method: " << _IntegrationMethod << endl;
1343 exit(1);
1344 }
1345 }
1346 }
1347
1348 // -----------------------------------------------------------------------------
JacobianDOFs(double jac[3],int dof,double x,double y,double z,double t,double t0) const1349 void BSplineFreeFormTransformationSV::JacobianDOFs(double jac[3], int dof, double x, double y, double z, double t, double t0) const
1350 {
1351 Matrix dTdp(3, 3);
1352 this->JacobianDOFs(dTdp, dof / 3, x, y, z, t, t0);
1353 const int c = dof % 3;
1354 jac[0] = dTdp(0, c);
1355 jac[1] = dTdp(1, c);
1356 jac[2] = dTdp(2, c);
1357 }
1358
1359 // -----------------------------------------------------------------------------
1360 void BSplineFreeFormTransformationSV
EvaluateJacobianDOFs(TransformationJacobian & jac,double x,double y) const1361 ::EvaluateJacobianDOFs(TransformationJacobian &jac, double x, double y) const
1362 {
1363 int i = ifloor(x);
1364 int j = ifloor(y);
1365
1366 const int A = Kernel::VariableToIndex(x - i);
1367 const int B = Kernel::VariableToIndex(y - j);
1368
1369 --i, --j;
1370
1371 double wxy, wy;
1372 int ci, cj, xdof, ydof;
1373
1374 for (int b = 0; b <= 3; ++b) {
1375 cj = j + b;
1376 if (cj < 0 || cj >= _y) continue;
1377 wy = Kernel::LookupTable[B][b];
1378 for (int a = 0; a <= 3; ++a) {
1379 ci = i + a;
1380 if (ci < 0 || ci >= _x) continue;
1381 wxy = Kernel::LookupTable[A][a] * wy;
1382 IndexToDOFs(LatticeToIndex(ci, cj), xdof, ydof);
1383 jac(xdof)._x = jac(ydof)._y = wxy;
1384 }
1385 }
1386 }
1387
1388 // -----------------------------------------------------------------------------
1389 void BSplineFreeFormTransformationSV
EvaluateJacobianDOFs(TransformationJacobian & jac,double x,double y,double z) const1390 ::EvaluateJacobianDOFs(TransformationJacobian &jac, double x, double y, double z) const
1391 {
1392 int i = ifloor(x);
1393 int j = ifloor(y);
1394 int k = ifloor(z);
1395
1396 const int A = Kernel::VariableToIndex(x - i);
1397 const int B = Kernel::VariableToIndex(y - j);
1398 const int C = Kernel::VariableToIndex(z - k);
1399
1400 --i, --j, --k;
1401
1402 double wxyz, wyz, wz;
1403 int ci, cj, ck, xdof, ydof, zdof;
1404
1405 for (int c = 0; c <= 3; ++c) {
1406 ck = k + c;
1407 if (ck < 0 || ck >= _z) continue;
1408 wz = Kernel::LookupTable[C][c];
1409 for (int b = 0; b <= 3; ++b) {
1410 cj = j + b;
1411 if (cj < 0 || cj >= _y) continue;
1412 wyz = Kernel::LookupTable[B][b] * wz;
1413 for (int a = 0; a <= 3; ++a) {
1414 ci = i + a;
1415 if (ci < 0 || ci >= _x) continue;
1416 wxyz = Kernel::LookupTable[A][a] * wyz;
1417 IndexToDOFs(LatticeToIndex(ci, cj, ck), xdof, ydof, zdof);
1418 jac(xdof)._x = jac(ydof)._y = jac(zdof)._z = wxyz;
1419 }
1420 }
1421 }
1422 }
1423
1424 // -----------------------------------------------------------------------------
JacobianDOFs(TransformationJacobian & jac,double x,double y,double z,double t,double t0) const1425 void BSplineFreeFormTransformationSV::JacobianDOFs(TransformationJacobian &jac, double x, double y, double z, double t, double t0) const
1426 {
1427 jac.Clear();
1428 double dt, T;
1429 if ((dt = StepLengthForIntervalLength(T = UpperIntegrationLimit(t, t0)))) {
1430 if (_IntegrationMethod == FFDIM_SS ||
1431 _IntegrationMethod == FFDIM_FastSS ||
1432 _IntegrationMethod == FFDIM_RKE1) RKE1 ::JacobianDOFs(this, jac, x, y, z, .0, T, dt);
1433 else if (_IntegrationMethod == FFDIM_RKE2) RKE2 ::JacobianDOFs(this, jac, x, y, z, .0, T, dt);
1434 else if (_IntegrationMethod == FFDIM_RKH2) RKH2 ::JacobianDOFs(this, jac, x, y, z, .0, T, dt);
1435 else if (_IntegrationMethod == FFDIM_RK4) RK4 ::JacobianDOFs(this, jac, x, y, z, .0, T, dt);
1436 else if (_IntegrationMethod == FFDIM_RKEH12) RKEH12::JacobianDOFs(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1437 else if (_IntegrationMethod == FFDIM_RKBS23) RKBS23::JacobianDOFs(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1438 else if (_IntegrationMethod == FFDIM_RKF45) RKF45 ::JacobianDOFs(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1439 else if (_IntegrationMethod == FFDIM_RKCK45) RKCK45::JacobianDOFs(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1440 else if (_IntegrationMethod == FFDIM_RKDP45) RKDP45::JacobianDOFs(this, jac, x, y, z, .0, T, 0.5 * dt, 2.0 * dt, SVFFD_RKTOL);
1441 else {
1442 cerr << "BSplineFreeFormTransformationSV::JacobianDOFs: Unknown integration method: " << _IntegrationMethod << endl;
1443 exit(1);
1444 }
1445 }
1446 }
1447
1448 // -----------------------------------------------------------------------------
1449 inline void BSplineFreeFormTransformationSV
JacobianDetDerivative(double[3],const Matrix &,int,double,double,double,double,double,bool,bool) const1450 ::JacobianDetDerivative(double [3], const Matrix &, int, double, double, double, double, double, bool, bool) const
1451 {
1452 Throw(ERR_NotImplemented, __FUNCTION__, "Not implemented");
1453 }
1454
1455 namespace BSplineFreeFormTransformationSVUtils {
1456
1457 // -----------------------------------------------------------------------------
1458 struct MultiplyDerivatives : public VoxelFunction
1459 {
1460 const int x, y, z, xx, xy, xz, yx, yy, yz, zx, zy, zz; // offsets
1461
MultiplyDerivativesmirtk::BSplineFreeFormTransformationSVUtils::MultiplyDerivatives1462 MultiplyDerivatives(int n)
1463 :
1464 x(0), y(x+n), z(y+n),
1465 xx( 0), xy(xx+n), xz(xy+n),
1466 yx(xz+n), yy(yx+n), yz(yy+n),
1467 zx(yz+n), zy(zx+n), zz(zy+n)
1468 {}
1469
operator ()mirtk::BSplineFreeFormTransformationSVUtils::MultiplyDerivatives1470 void operator ()(const GenericImage<double> &, int, const double *in, const double *d, double *out) const
1471 {
1472 // Attention: out can be equal d, therefore use temporary variables
1473 double gx = in[x] * d[xx] + in[x] * d[xy] + in[x] * d[xz];
1474 double gy = in[y] * d[yx] + in[y] * d[yy] + in[y] * d[yz];
1475 double gz = in[z] * d[zx] + in[z] * d[zy] + in[z] * d[zz];
1476 out[x] = gx, out[y] = gy, out[z] = gz;
1477 }
1478 };
1479
1480
1481 } // namespace BSplineFreeFormTransformationSVUtils
1482 using namespace BSplineFreeFormTransformationSVUtils;
1483
1484 // -----------------------------------------------------------------------------
1485 void BSplineFreeFormTransformationSV
ParametricGradient(const GenericImage<double> * in,double * out,const WorldCoordsImage * i2w,const WorldCoordsImage * wc,double t,double t0,double w) const1486 ::ParametricGradient(const GenericImage<double> *in, double *out,
1487 const WorldCoordsImage *i2w, const WorldCoordsImage *wc,
1488 double t, double t0, double w) const
1489 {
1490 // Upper integration limit for given interval
1491 const double T = UpperIntegrationLimit(t, t0);
1492 if (IsZero(T)) return;
1493
1494 // ---------------------------------------------------------------------------
1495 // BCH based velocity update computation
1496 if (_NumberOfBCHTerms > 1) {
1497
1498 MIRTK_START_TIMING();
1499 if (_UseDenseBCHGrid && _NumberOfBCHTerms > 2) {
1500 GenericImage<double> v(in->Attributes(), 3);
1501 this->Velocity(v);
1502 EvaluateBCHFormulaDense<double>(_NumberOfBCHTerms, v, v, *in, true, _LieDerivative);
1503 BSplineFreeFormTransformation3D::ParametricGradient(&v, out, i2w, wc, t0, w / T);
1504 } else {
1505 // Compute spline coefficients of update velocity field
1506 CPImage u(this->Attributes());
1507 DOFValue * const grd = reinterpret_cast<DOFValue *>(u.Data());
1508 BSplineFreeFormTransformation3D::ParametricGradient(in, grd, i2w, wc, t0, 1.);
1509 // Approximate velocity spline coefficients of composite transformation
1510 // using Baker-Campbell-Hausdorff (BCH) formula and subtract current coefficients
1511 EvaluateBCHFormula(_NumberOfBCHTerms, u, T, _CPImage, 1., u, true);
1512 // Adjust weight as update field is computed for tau * v, i.e.,
1513 // exp(tau * v_{i+1}) = exp(tau v_i) o exp(\delta u)
1514 // ==> v_{i+1} = log(exp(tau * v_{i+1})) / tau
1515 w /= T;
1516 // Add weighted gradient to total energy gradient
1517 for (int dof = 0; dof < this->NumberOfDOFs(); ++dof) {
1518 if (_Status[dof] == Active) out[dof] += w * grd[dof];
1519 }
1520 }
1521 MIRTK_DEBUG_TIMING(2, "parametric gradient computation (BCH, dense=" << (_UseDenseBCHGrid ? "yes" : "no") << ")");
1522
1523 // ---------------------------------------------------------------------------
1524 // Scaling and squaring based gradient computation
1525 } else if (_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS) {
1526
1527 MIRTK_START_TIMING();
1528
1529 // Copy B-spline coefficients of velocity field
1530 GenericImage<double> v(this->Attributes(), 3);
1531 double *vx = v.Data(0, 0, 0, 0);
1532 double *vy = v.Data(0, 0, 0, 1);
1533 double *vz = v.Data(0, 0, 0, 2);
1534 const Vector *vp = _CPImage.Data();
1535 for (int idx = 0; idx < _CPImage.NumberOfVoxels(); ++idx, ++vx, ++vy, ++vz, ++vp) {
1536 *vx = vp->_x;
1537 *vy = vp->_y;
1538 *vz = vp->_z;
1539 }
1540
1541 // Exponentiate non-parametric gradient using inverse mapping
1542 // (cf. Modat et al., "Parametric non-rigid registration using a stationary velocity field", MMBIA, 145–150, 2012)
1543 GenericImage<double> dv;
1544 mirtk::ScalingAndSquaring<double> exp;
1545 exp.UpperIntegrationLimit(T);
1546 exp.ComputeInverse(true);
1547 exp.NumberOfSteps(NumberOfStepsForIntervalLength(T));
1548 exp.MaxScaledVelocity(_MaxScaledVelocity);
1549 exp.InterimAttributes(_IntegrationMethod == FFDIM_FastSS ? this->Attributes() : in->Attributes());
1550 exp.OutputAttributes(exp.InterimAttributes());
1551 exp.InputVelocity(&v);
1552 exp.ComputeInterpolationCoefficients(false);
1553 exp.InputGradient(in);
1554 exp.OutputGradient(&dv);
1555 exp.Run();
1556
1557 // Free copy of velocity coefficients
1558 v.Clear();
1559
1560 // Multiply resulting vectors by derivative of v w.r.t. the DoFs
1561 BSplineFreeFormTransformation3D::ParametricGradient(&dv, out, i2w, wc, t0, w);
1562
1563 MIRTK_DEBUG_TIMING(2, "parametric gradient computation" << (_IntegrationMethod == FFDIM_FastSS ? " (FastSS)" : " (SS)"));
1564
1565 // ---------------------------------------------------------------------------
1566 // Runge-Kutta integration based gradient computation similar to TD FFD
1567 // transformation parameterized by non-stationary velocity field
1568 } else {
1569
1570 // Note: T = in->GetTOrigin() - t0
1571 FreeFormTransformation::ParametricGradient(in, out, i2w, wc, in->GetTOrigin() - T, w);
1572
1573 }
1574 }
1575
1576 // -----------------------------------------------------------------------------
1577 void BSplineFreeFormTransformationSV
ParametricGradient(const PointSet & pos,const Vector3D<double> * in,double * out,double t,double t0,double w) const1578 ::ParametricGradient(const PointSet &pos, const Vector3D<double> *in,
1579 double *out, double t, double t0, double w) const
1580 {
1581 // Runge-Kutta integration based gradient computation similar to TD FFD
1582 // transformation parameterized by non-stationary velocity field
1583 FreeFormTransformation::ParametricGradient(pos, in, out, t, t0, w);
1584 }
1585
1586 // =============================================================================
1587 // I/O
1588 // =============================================================================
1589
1590 // -----------------------------------------------------------------------------
Print(ostream & os,Indent indent) const1591 void BSplineFreeFormTransformationSV::Print(ostream &os, Indent indent) const
1592 {
1593 os << indent << "B-spline SV FFD:" << endl;
1594 indent++;
1595 // Print FFD attributes
1596 FreeFormTransformation3D::Print(os, indent);
1597 // Change output stream settings
1598 const streamsize w = os.width (0);
1599 const streamsize p = os.precision(2);
1600 const ios::fmtflags f = os.flags();
1601 cout.unsetf(ios::floatfield);
1602 // Print SV FFD parameters
1603 os << indent << "Integration method: " << setw(6) << ToString(_IntegrationMethod) << endl;
1604 os << indent << "Cross-sectional time interval: " << setw(6) << _T << endl;
1605 os << indent << "Time unit of integration interval: " << setw(6) << _TimeUnit << endl;
1606 os << indent << "Maximum scaled velocity: " << setw(6) << _MaxScaledVelocity << endl;
1607 os << indent << "No. of integration steps per unit: " << setw(6) << _NumberOfSteps << endl;
1608 os << indent << "No. of cross-sectional steps: " << setw(6) << NumberOfStepsForIntervalLength(_T) << endl;
1609 os << indent << "No. of BCH terms: " << setw(6) << _NumberOfBCHTerms << endl;
1610 os << indent << "Use Lie derivative: " << setw(6) << ToString(_LieDerivative) << endl;
1611 // Restore output stream settings
1612 os.width (w);
1613 os.precision(p);
1614 os.flags (f);
1615 }
1616
1617 // -----------------------------------------------------------------------------
CanRead(TransformationType format) const1618 bool BSplineFreeFormTransformationSV::CanRead(TransformationType format) const
1619 {
1620 switch (format) {
1621 case TRANSFORMATION_BSPLINE_FFD_SV_v1:
1622 case TRANSFORMATION_BSPLINE_FFD_SV_v2:
1623 case TRANSFORMATION_BSPLINE_FFD_SV_v3:
1624 case TRANSFORMATION_BSPLINE_FFD_SV_v4:
1625 case TRANSFORMATION_BSPLINE_FFD_SV_v5:
1626 case TRANSFORMATION_BSPLINE_FFD_SV_v6:
1627 case TRANSFORMATION_BSPLINE_FFD_SV_v7:
1628 case TRANSFORMATION_BSPLINE_FFD_SV_v8:
1629 return true;
1630 default:
1631 return false;
1632 }
1633 }
1634
1635 // -----------------------------------------------------------------------------
ReadDOFs(Cifstream & from,TransformationType format)1636 Cifstream &BSplineFreeFormTransformationSV::ReadDOFs(Cifstream &from, TransformationType format)
1637 {
1638 // Read FFD data
1639 switch (format) {
1640 case TRANSFORMATION_BSPLINE_FFD_SV_v1:
1641 case TRANSFORMATION_BSPLINE_FFD_SV_v2:
1642 case TRANSFORMATION_BSPLINE_FFD_SV_v3:
1643 case TRANSFORMATION_BSPLINE_FFD_SV_v4:
1644 case TRANSFORMATION_BSPLINE_FFD_SV_v5:
1645 case TRANSFORMATION_BSPLINE_FFD_SV_v6:
1646 BSplineFreeFormTransformation3D::ReadDOFs(from, TRANSFORMATION_BSPLINE_FFD_3D_v2);
1647 case TRANSFORMATION_BSPLINE_FFD_SV_v7:
1648 BSplineFreeFormTransformation3D::ReadDOFs(from, TRANSFORMATION_BSPLINE_FFD_3D_v3);
1649 default:
1650 BSplineFreeFormTransformation3D::ReadDOFs(from, TRANSFORMATION_BSPLINE_FFD_3D);
1651 }
1652
1653 // Read number of integration steps
1654 from.ReadAsInt(&_NumberOfSteps, 1);
1655
1656 if (format == TRANSFORMATION_BSPLINE_FFD_SV_v1) return from;
1657
1658 // Read upper integration limit
1659 from.ReadAsDouble(&_T, 1);
1660 // Read number of BCH terms to use for update
1661 from.ReadAsInt(&_NumberOfBCHTerms, 1);
1662
1663 if (format == TRANSFORMATION_BSPLINE_FFD_SV_v2) return from;
1664
1665 // Read time unit of integration interval
1666 from.ReadAsDouble(&_TimeUnit, 1);
1667
1668 if (format == TRANSFORMATION_BSPLINE_FFD_SV_v3) return from;
1669
1670 if (format <= TRANSFORMATION_BSPLINE_FFD_SV_v6) {
1671
1672 // Whether to use scaling and squaring
1673 char useSS;
1674 from.ReadAsChar(&useSS, 1);
1675
1676 // Maximum scaled velocity
1677 from.ReadAsDouble(&_MaxScaledVelocity, 1);
1678
1679 if (format == TRANSFORMATION_BSPLINE_FFD_SV_v4) return from;
1680
1681 // Whether to use fast scaling and squaring
1682 char fastSS;
1683 from.ReadAsChar(&fastSS, 1);
1684
1685 // Set integration method
1686 _IntegrationMethod = (useSS ? (fastSS ? FFDIM_FastSS : FFDIM_SS) : FFDIM_RKE1);
1687
1688 } else {
1689
1690 // Integration method
1691 unsigned int integration_method;
1692 from.ReadAsUInt(&integration_method, 1);
1693 _IntegrationMethod = static_cast<FFDIntegrationMethod>(integration_method);
1694
1695 // Maximum scaled velocity
1696 from.ReadAsDouble(&_MaxScaledVelocity, 1);
1697
1698 }
1699
1700 return from;
1701 }
1702
1703 // -----------------------------------------------------------------------------
WriteDOFs(Cofstream & to) const1704 Cofstream &BSplineFreeFormTransformationSV::WriteDOFs(Cofstream &to) const
1705 {
1706 // Write FFD data
1707 BSplineFreeFormTransformation3D::WriteDOFs(to);
1708
1709 // Write number of integration steps
1710 to.WriteAsInt(&_NumberOfSteps, 1);
1711 // Write upper integration limit
1712 to.WriteAsDouble(&_T, 1);
1713 // Write number of BCH terms to use for update
1714 to.WriteAsInt(&_NumberOfBCHTerms, 1);
1715 // Write time unit of integration interval
1716 to.WriteAsDouble(&_TimeUnit, 1);
1717 // Integration method
1718 const unsigned int integration_method = _IntegrationMethod;
1719 to.WriteAsUInt(&integration_method, 1);
1720 // Maximum scaled velocity
1721 to.WriteAsDouble(&_MaxScaledVelocity, 1);
1722
1723 return to;
1724 }
1725
1726
1727 } // namespace mirtk
1728