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 #ifndef MIRTK_BSplineFreeFormTransformationSV_H
21 #define MIRTK_BSplineFreeFormTransformationSV_H
22
23 #include "mirtk/BSplineFreeFormTransformation3D.h"
24
25 #include "mirtk/Math.h"
26 #include "mirtk/ImageFunction.h"
27 #include "mirtk/VoxelFunction.h"
28 #include "mirtk/FFDIntegrationMethod.h"
29
30
31 namespace mirtk {
32
33
34 /**
35 * Free-form transformation parameterized by a stationary velocity field.
36 *
37 * This class implements a free-form transformation which is represented
38 * by a stationary velocity field (2D/3D). The displacement field is obtained
39 * by exponentiating the velocity field, resulting in a diffeomorphic
40 * transformation. A stationary velocity field can be further approximated
41 * from a given displacement field using the group logarithm of diffeomorphisms.
42 */
43 class BSplineFreeFormTransformationSV : public BSplineFreeFormTransformation3D
44 {
45 mirtkTransformationMacro(BSplineFreeFormTransformationSV);
46
47 // ---------------------------------------------------------------------------
48 // Attributes
49
50 /// Cross-sectional upper integration limit, i.e., time difference between
51 /// target and source if both volumes have the same temporal origin.
52 /// Otherwise, the difference between the temporal origins as given by the
53 /// t and t0 parameters of LocalTransform et al. is used instead.
54 ///
55 /// \note Can be set to zero for the generation of a deformation movie using
56 /// the transformation tool based on ImageTransformation by varying
57 /// the -Tt parameter from zero to one and keeping -St fixed at one.
58 /// Use the doftool to set the --cross-sectional-time-interval to zero.
59 mirtkPublicAttributeMacro(double, T);
60
61 /// Temporal unit used to determine actual number of integration steps.
62 /// If zero, the number of steps is independent of the temporal interval.
63 mirtkPublicAttributeMacro(double, TimeUnit);
64
65 /// Number of integration steps per _TimeUnit
66 mirtkPublicAttributeMacro(int, NumberOfSteps);
67
68 /// Maximum absolute value of scaled velocity when using the
69 /// scaling-and-squaring for computing the displacement field.
70 /// This introduces an adapative choice of the number of integration steps
71 /// needed. At a minimum, NumberOfSteps are taken, but if needed, more
72 /// steps added. If this attribute is set to zero, however, the number
73 /// of integration steps always corresponds to the specified NumberOfSteps.
74 mirtkPublicAttributeMacro(double, MaxScaledVelocity);
75
76 /// Integration method used to obtain displacement from velocities
77 mirtkPublicAttributeMacro(FFDIntegrationMethod, IntegrationMethod);
78
79 /// Whether to evaluate BCH formulat for parametric gradient calculation
80 /// on dense image lattice of non-parametric gradient using linear interpolation
81 mirtkPublicAttributeMacro(bool, UseDenseBCHGrid);
82
83 /// Use Lie derivative definition of Lie bracket based on vector field
84 /// Jacobian matrices instead of difference of vector field composition
85 mirtkPublicAttributeMacro(bool, LieDerivative);
86
87 /// Number of BCH terms used by ParametricGradient. When set to an invalid
88 /// value, i.e., no. of terms < 2, the gradient is computed using the specified
89 /// integration method instead.
90 mirtkPublicAttributeMacro(int, NumberOfBCHTerms);
91
92 // ---------------------------------------------------------------------------
93 // Construction/Destruction
94
95 protected:
96
97 /// Initialize extrapolation function
98 void InitializeExtrapolator();
99
100 public:
101
102 /// Default constructor
103 BSplineFreeFormTransformationSV();
104
105 /// Construct free-form transformation for given image domain and lattice spacing
106 explicit BSplineFreeFormTransformationSV(const ImageAttributes &,
107 double = -1, double = -1, double = -1);
108
109 /// Construct free-form transformation for given target image and lattice spacing
110 explicit BSplineFreeFormTransformationSV(const BaseImage &,
111 double, double, double);
112
113 /// Construct free-form transformation from existing 3D+t deformation field
114 explicit BSplineFreeFormTransformationSV(const GenericImage<double> &, bool = false);
115
116 /// Copy constructor
117 BSplineFreeFormTransformationSV(const BSplineFreeFormTransformationSV &);
118
119 /// Destructor
120 virtual ~BSplineFreeFormTransformationSV();
121
122 // Import other Initialize overloads
123 using BSplineFreeFormTransformation3D::Initialize;
124
125 /// Initialize free-form transformation constructed by default constructor
126 virtual void Initialize(const ImageAttributes &);
127
128 /// Initialize free-form transformation which approximates a given transformation
129 virtual void Initialize(const ImageAttributes &, double, double, double,
130 const Transformation *);
131
132 /// Subdivide FFD lattice
133 virtual void Subdivide(bool = true, bool = true, bool = true, bool = true);
134
135 // ---------------------------------------------------------------------------
136 // Approximation/Interpolation
137
138 using BSplineFreeFormTransformation3D::ApproximateAsNew;
139
140 /// Approximate displacements: This function takes a set of points and a set
141 /// of displacements and finds !new! parameters such that the resulting
142 /// transformation approximates the displacements as good as possible.
143 virtual void ApproximateDOFs(const double *, const double *, const double *, const double *,
144 const double *, const double *, const double *, int);
145
146 /// Finds gradient of approximation error: This function takes a set of points
147 /// and a set of errors. It finds a gradient w.r.t. the transformation parameters
148 /// which minimizes the L2 norm of the approximation error and adds it to the
149 /// input gradient with the given weight.
150 virtual void ApproximateDOFsGradient(const double *, const double *, const double *, const double *,
151 const double *, const double *, const double *, int,
152 double *, double = 1.0) const;
153
154 /// Get image domain on which this free-form transformation should be defined
155 /// in order to reduce the error when the given transformation is approximated
156 /// and applied to resample an image on the specified image grid.
157 virtual ImageAttributes ApproximationDomain(const ImageAttributes &,
158 const Transformation *);
159
160 /// Approximate another transformation and return approximation error
161 virtual double ApproximateAsNew(const ImageAttributes &, const Transformation *, int = 1, double = .0);
162
163 /// Approximate displacements: This function takes a 3D displacement field
164 /// and finds a stationary velocity field which, when exponentiated, approximates
165 /// these displacements by computing the group logarithm of diffeomorphisms.
166 /// The displacements are replaced by the residual displacements of the newly
167 /// approximated transformation. Returns the approximation error of the resulting FFD.
168 virtual double ApproximateAsNew(GenericImage<double> &, int = 1, double = .0);
169
170 /// Approximate displacements: This function takes a 3D displacement field
171 /// and finds a stationary velocity field which, when exponentiated, approximates
172 /// these displacements by computing the group logarithm of diffeomorphisms.
173 /// The displacements are replaced by the residual displacements of the newly
174 /// approximated transformation. Returns the approximation error of the resulting FFD.
175 virtual double ApproximateAsNew(GenericImage<double> &, bool, int = 3, int = 8);
176
177 /// Interpolates displacements: This function takes a set of displacements defined at
178 /// the control points and finds a stationary velocity field such that the
179 /// resulting transformation interpolates these displacements.
180 virtual void Interpolate(const double *, const double *, const double *);
181
182 /// Interpolates velocities: This function takes a set of velocites defined at
183 /// the control points and interpolates these velocities.
184 virtual void InterpolateVelocities(const double *, const double *, const double *);
185
186 /// Approximate velocities: This function takes a velocity field and initializes
187 /// the cubic B-spline coefficients of this transformation such that the spline
188 /// function approximates the given velocity field. Returns the approximation error
189 /// of the velocity spline function.
190 virtual double ApproximateVelocitiesAsNew(GenericImage<double> &);
191
192 /// Combine transformations: This function takes a transformation and finds
193 /// a stationary velocity field which, when exponentiated, approximates the
194 /// composition of this and the given transformation.
195 virtual void CombineWith(const Transformation *);
196
197 /// Combine transformations: This function takes a transformation and finds
198 /// a stationary velocity field which, when exponentiated, approximates the
199 /// composition of this and the given transformation.
200 virtual void CombineWith(const BSplineFreeFormTransformationSV *);
201
202 /// Invert this transformation
203 virtual void Invert();
204
205 // ---------------------------------------------------------------------------
206 // Parameters (non-DoFs)
207
208 using BSplineFreeFormTransformation3D::Parameter;
209
210 /// Set named (non-DoF) parameter from value as string
211 virtual bool Set(const char *, const char *);
212
213 /// Get (non-DoF) parameters as key/value as string map
214 virtual ParameterList Parameter() const;
215
216 // ---------------------------------------------------------------------------
217 // Evaluation (for use FreeFormTransformationRungeKutta integration)
218
219 using BSplineFreeFormTransformation3D::Evaluate;
220 using BSplineFreeFormTransformation3D::EvaluateJacobianWorld;
221 using BSplineFreeFormTransformation3D::EvaluateJacobianDOFs;
222
223 /// Evaluates the FFD at a point in lattice coordinates
224 void Evaluate(double &, double &, double &, double) const;
225
226 /// Calculates the Jacobian of the FFD at a point in lattice coordinates
227 /// and converts the resulting Jacobian to derivatives w.r.t world coordinates
228 void EvaluateJacobianWorld(Matrix &, double, double, double, double) const;
229
230 /// Calculates the Jacobian of the 2D transformation w.r.t. the transformation parameters
231 void EvaluateJacobianDOFs(TransformationJacobian &,
232 double, double) const;
233
234 /// Calculates the Jacobian of the 3D transformation w.r.t. the transformation parameters
235 void EvaluateJacobianDOFs(TransformationJacobian &,
236 double, double, double) const;
237
238 /// Calculates the Jacobian of the transformation w.r.t. the transformation parameters
239 void EvaluateJacobianDOFs(TransformationJacobian &,
240 double, double, double, double) const;
241
242 // ---------------------------------------------------------------------------
243 // Point transformation
244
245 public:
246
247 using BSplineFreeFormTransformation3D::Displacement;
248
249 protected:
250
251 /// Get number of integration steps for a given temporal integration interval
252 int NumberOfStepsForIntervalLength(double) const;
253
254 /// Set number of integration steps for a given temporal integration interval
255 void NumberOfStepsForIntervalLength(double, int) const;
256
257 /// Get length of steps for a given temporal integration interval
258 double StepLengthForIntervalLength(double) const;
259
260 public:
261
262 /// Scale velocities
263 ///
264 /// Alternatively, change integration inveral _T.
265 void ScaleVelocities(double);
266
267 /// Upper integration limit used given the temporal origin of both target
268 /// and source images. If both images have the same temporal origin, the value
269 /// of the member variable @c T is returned. Note that @c T may be set to zero
270 /// in order to force no displacement between images located at the same point
271 /// in time to be zero. This is especially useful when animating the
272 /// deformation between two images with successively increasing upper
273 /// integration limit. Otherwise, if the temporal origin of the two images
274 /// differs, the signed difference between these is returned, i.e., t - t0,
275 /// which corresponds to a forward integration of the target point (x, y, z)
276 /// to the time point of the source image.
277 double UpperIntegrationLimit(double t, double t0) const;
278
279 /// Transforms a single point using a Runge-Kutta integration
280 ///
281 /// \param[in] T Upper integration limit, i.e., length of time interval.
282 void IntegrateVelocities(double &, double &, double &, double T = 1.0) const;
283
284 /// Compute displacement field using the scaling and squaring (SS) method
285 ///
286 /// \attention If an input/output displacement field is provided, the resulting
287 /// transformation is the composition exp(v) o d.
288 ///
289 /// \attention The scaling and squaring cannot be used to obtain the displacements
290 /// generated by a 3D velocity field on a 2D slice only. It always
291 /// must be applied to the entire domain of the velocity field.
292 ///
293 /// \param[in,out] d Displacement field generated by stationary velocity field.
294 /// \param[in] T Upper integration limit, i.e., length of time interval.
295 /// \param[in] wc Pre-computed world coordinates of image voxels.
296 template <class VoxelType>
297 void ScalingAndSquaring(GenericImage<VoxelType> *d,
298 double T = 1.0, const WorldCoordsImage *wc = NULL) const;
299
300 /// Compute displacement field and derivatives using the scaling and squaring (SS) method
301 ///
302 /// \attention If an input/output displacement field is provided, the resulting
303 /// transformation is the composition exp(v) o d.
304 ///
305 /// \attention The scaling and squaring cannot be used to obtain the displacements
306 /// generated by a 3D velocity field on a 2D slice only. It always
307 /// must be applied to the entire domain of the velocity field.
308 ///
309 /// \param[in] attr Attributes of output images.
310 /// \param[in,out] d Displacement field generated by stationary velocity field.
311 /// \param[out] dx Partial derivatives of corresponding transformation w.r.t. x.
312 /// \param[out] dj Determinant of Jacobian w.r.t. x, i.e., det(jac(\p dx)).
313 /// \param[out] lj Log of determinant of Jacobian w.r.t. x, i.e., log(det(jac(\p dx))).
314 /// \param[in] T Upper integration limit, i.e., length of time interval.
315 /// \param[in] wc Pre-computed world coordinates of image voxels.
316 template <class VoxelType>
317 void ScalingAndSquaring(const ImageAttributes &attr,
318 GenericImage<VoxelType> *d,
319 GenericImage<VoxelType> *dx,
320 GenericImage<VoxelType> *dj = NULL,
321 GenericImage<VoxelType> *lj = NULL,
322 double T = 1.0, const WorldCoordsImage *wc = NULL) const;
323
324 /// Whether the caching of the transformation displacements is required
325 /// (or preferred) by this transformation. For some transformations such as
326 /// those parameterized by velocities, caching of the displacements for
327 /// each target voxel results in better performance or is needed for example
328 /// for the scaling and squaring method.
329 virtual bool RequiresCachingOfDisplacements() const;
330
331 /// Transforms a single point using the local transformation component only
332 virtual void LocalTransform(double &, double &, double &, double = 0, double = NaN) const;
333
334 /// Transforms a single point using the inverse of the local transformation only
335 virtual bool LocalInverse(double &, double &, double &, double = 0, double = NaN) const;
336
337 /// Get stationary velocity field
338 void Velocity(GenericImage<float> &) const;
339
340 /// Get stationary velocity field
341 void Velocity(GenericImage<double> &) const;
342
343 /// Calculates the displacement vectors for a whole image domain
344 ///
345 /// \attention The displacements are computed at the positions after applying the
346 /// current displacements at each voxel. These displacements are then
347 /// added to the current displacements. Therefore, set the input
348 /// displacements to zero if only interested in the displacements of
349 /// this transformation at the voxel positions.
350 virtual void Displacement(GenericImage<double> &, double, double = NaN,
351 const WorldCoordsImage * = NULL) const;
352
353 /// Calculates the displacement vectors for a whole image domain
354 ///
355 /// \attention The displacements are computed at the positions after applying the
356 /// current displacements at each voxel. These displacements are then
357 /// added to the current displacements. Therefore, set the input
358 /// displacements to zero if only interested in the displacements of
359 /// this transformation at the voxel positions.
360 virtual void Displacement(GenericImage<float> &, double, double = NaN,
361 const WorldCoordsImage * = NULL) const;
362
363 /// Calculates the inverse displacement vectors for a whole image domain
364 ///
365 /// \attention The displacements are computed at the positions after applying the
366 /// current displacements at each voxel. These displacements are then
367 /// added to the current displacements. Therefore, set the input
368 /// displacements to zero if only interested in the displacements of
369 /// this transformation at the voxel positions.
370 ///
371 /// \returns Always zero.
372 virtual int InverseDisplacement(GenericImage<double> &, double, double = NaN,
373 const WorldCoordsImage * = NULL) const;
374
375 /// Calculates the inverse displacement vectors for a whole image domain
376 ///
377 /// \attention The displacements are computed at the positions after applying the
378 /// current displacements at each voxel. These displacements are then
379 /// added to the current displacements. Therefore, set the input
380 /// displacements to zero if only interested in the displacements of
381 /// this transformation at the voxel positions.
382 ///
383 /// \returns Always zero.
384 virtual int InverseDisplacement(GenericImage<float> &, double, double = NaN,
385 const WorldCoordsImage * = NULL) const;
386
387 // ---------------------------------------------------------------------------
388 // Derivatives
389
390 using BSplineFreeFormTransformation3D::JacobianDOFs;
391 using BSplineFreeFormTransformation3D::ParametricGradient;
392
393 /// Calculates the Jacobian of the local transformation w.r.t world coordinates
394 virtual void LocalJacobian(Matrix &, double, double, double, double = 0, double = NaN) const;
395
396 /// Calculates the Hessian for each component of the local transformation w.r.t world coordinates
397 virtual void LocalHessian(Matrix [3], double, double, double, double = 0, double = NaN) const;
398
399 /// Calculates the Jacobian of the transformation w.r.t a control point
400 virtual void JacobianDOFs(Matrix &, int, double, double, double, double = 0, double = NaN) const;
401
402 /// Calculates the Jacobian of the transformation w.r.t a transformation parameters
403 ///
404 /// \attention This overriden base class function returns one column of the 3x3 Jacobian
405 /// matrix, i.e., only the partial derivatives of the transformation
406 /// w.r.t. the specified transformation parameter (not control point).
407 /// This is in accordance to the Transformation::JacobianDOFs
408 /// definition, which is violated, however, by most other FFDs.
409 virtual void JacobianDOFs(double [3], int, double, double, double, double = 0, double = NaN) const;
410
411 /// Calculates the Jacobian of the transformation w.r.t the transformation parameters
412 virtual void JacobianDOFs(TransformationJacobian &, double, double, double, double = 0, double = NaN) const;
413
414 /// Calculates derivatives of the Jacobian determinant at world point w.r.t. DoFs of a control point
415 ///
416 /// \param[out] dJ Partial derivatives of Jacobian determinant at (x, y, z) w.r.t. DoFs of control point.
417 /// \param[in] cp Index of control point w.r.t. whose DoFs the derivatives are computed.
418 /// \param[in] x World coordinate along x axis at which to evaluate derivatives.
419 /// \param[in] y World coordinate along y axis at which to evaluate derivatives.
420 /// \param[in] z World coordinate along z axis at which to evaluate derivatives.
421 /// \param[in] adj Adjugate of Jacobian matrix evaluated at (x, y, z).
422 virtual void JacobianDetDerivative(double dJ[3], const Matrix &adj,
423 int cp, double x, double y, double z, double t = 0, double t0 = NaN,
424 bool wrt_world = true, bool use_spacing = true) const;
425
426 protected:
427
428 /// Evaluates the BCH formula s.t. u = log(exp(tau * v) o exp(eta * w))
429 /// \sa Bossa, M., Hernandez, M., & Olmos, S. Contributions to 3D diffeomorphic
430 /// atlas estimation: application to brain images. MICCAI 2007, 10(Pt 1), 667–74.
431 void EvaluateBCHFormula(int, CPImage &, double, const CPImage &, double, const CPImage &, bool = false) const;
432
433 /// Applies the chain rule to convert spatial non-parametric gradient
434 /// to a gradient w.r.t the parameters of this transformation given the
435 /// temporal interval which the displacement gradient corresponds to
436 void ParametricGradient(const GenericImage<double> *, double *,
437 const WorldCoordsImage *,
438 const WorldCoordsImage *,
439 double, double, double) const;
440
441 public:
442
443 /// Applies the chain rule to convert spatial non-parametric gradient
444 /// to a gradient w.r.t the parameters of this transformation
445 virtual void ParametricGradient(const GenericImage<double> *, double *,
446 const WorldCoordsImage *,
447 const WorldCoordsImage *,
448 double = -1, double = 1) const;
449
450 /// Applies the chain rule to convert point-wise non-parametric gradient
451 /// to a gradient w.r.t the parameters of this transformation.
452 virtual void ParametricGradient(const PointSet &, const Vector3D<double> *,
453 double *, double = 0, double = -1, double = 1) const;
454
455 // ---------------------------------------------------------------------------
456 // I/O
457
458 // Do not hide methods of base class
459 using BSplineFreeFormTransformation3D::Print;
460
461 /// Prints the parameters of the transformation
462 virtual void Print(ostream &, Indent = 0) const;
463
464 /// Whether this transformation can read a file of specified type (i.e. format)
465 virtual bool CanRead(TransformationType) const;
466
467 protected:
468
469 /// Reads transformation parameters from a file stream
470 virtual Cifstream &ReadDOFs(Cifstream &, TransformationType);
471
472 /// Writes transformation parameters to a file stream
473 virtual Cofstream &WriteDOFs(Cofstream &) const;
474
475 // ---------------------------------------------------------------------------
476 // Friends
477 friend class EvaluateBSplineSVFFD2D;
478 friend class EvaluateBSplineSVFFD3D;
479 friend class PartialBSplineFreeFormTransformationSV;
480 friend class MultiLevelStationaryVelocityTransformation;
481 };
482
483 ////////////////////////////////////////////////////////////////////////////////
484 // Auxiliary voxel functions for subclasses
485 ////////////////////////////////////////////////////////////////////////////////
486
487 // -----------------------------------------------------------------------------
488 /// Evaluate global SV FFD at image voxels of vector field
489 class EvaluateGlobalSVFFD : public VoxelFunction
490 {
491 public:
492
EvaluateGlobalSVFFD(const Matrix & logA,BaseImage * output)493 EvaluateGlobalSVFFD(const Matrix &logA, BaseImage *output)
494 :
495 _LogA(logA), _Output(output)
496 {}
497
498 // template <class T>
499 // void Evaluate(Vector2D<T> &v, double x, double y, double z)
500 // {
501 // v._x = static_cast<T>(_LogA(0, 0) * x + _LogA(0, 1) * y + _LogA(0, 2) * z + _LogA(0, 3));
502 // v._y = static_cast<T>(_LogA(1, 0) * x + _LogA(1, 1) * y + _LogA(1, 2) * z + _LogA(1, 3));
503 // }
504
505 template <class T>
Evaluate(Vector3D<T> & v,double x,double y,double z)506 void Evaluate(Vector3D<T> &v, double x, double y, double z)
507 {
508 v._x = static_cast<T>(_LogA(0, 0) * x + _LogA(0, 1) * y + _LogA(0, 2) * z + _LogA(0, 3));
509 v._y = static_cast<T>(_LogA(1, 0) * x + _LogA(1, 1) * y + _LogA(1, 2) * z + _LogA(1, 3));
510 v._z = static_cast<T>(_LogA(2, 0) * x + _LogA(2, 1) * y + _LogA(2, 2) * z + _LogA(2, 3));
511 }
512
513 template <class VectorType>
operator()514 void operator()(int i, int j, int k, int, VectorType *v)
515 {
516 double x = i, y = j, z = k;
517 _Output->ImageToWorld(x, y, z);
518 Evaluate(*v, x, y, z);
519 }
520
521 private:
522 const Matrix _LogA; ///< Input transformation
523 BaseImage *_Output; ///< Output image
524 };
525
526 // -----------------------------------------------------------------------------
527 /// Evaluate global SV FFD at image voxels of 3D vector field
528 class EvaluateGlobalSVFFD3D : public VoxelFunction
529 {
530 public:
531
EvaluateGlobalSVFFD3D(const Matrix & logA,BaseImage * output)532 EvaluateGlobalSVFFD3D(const Matrix &logA, BaseImage *output)
533 :
534 _LogA (logA),
535 _Output(output),
536 _y (output->NumberOfSpatialVoxels()),
537 _z (2 * _y)
538 {}
539
540 template <class T>
operator()541 void operator()(int i, int j, int k, int, T *out)
542 {
543 double x = i, y = j, z = k;
544 _Output->ImageToWorld(x, y, z);
545 out[_x] = static_cast<T>(_LogA(0, 0) * x + _LogA(0, 1) * y + _LogA(0, 2) * z + _LogA(0, 3));
546 out[_y] = static_cast<T>(_LogA(1, 0) * x + _LogA(1, 1) * y + _LogA(1, 2) * z + _LogA(1, 3));
547 out[_z] = static_cast<T>(_LogA(2, 0) * x + _LogA(2, 1) * y + _LogA(2, 2) * z + _LogA(2, 3));
548 }
549
550 private:
551 const Matrix _LogA; ///< Input transformation
552 BaseImage *_Output; ///< Output image
553
554 static const int _x = 0; ///< Offset of x component
555 int _y; ///< Offset of y component
556 int _z; ///< Offset of z component
557 };
558
559 // -----------------------------------------------------------------------------
560 /// Evaluate B-spline SV FFD at image voxels
561 class EvaluateBSplineSVFFD : public VoxelFunction
562 {
563 public:
564
EvaluateBSplineSVFFD(const BSplineFreeFormTransformationSV * input,BaseImage * output)565 EvaluateBSplineSVFFD(const BSplineFreeFormTransformationSV *input,
566 BaseImage *output)
567 :
568 _Input(input), _Output(output)
569 {}
570
571 // template <class T>
572 // void Put(Vector2D<T> *v, double vx, double vy, double)
573 // {
574 // v->_x = vx, v->_y = vy;
575 // }
576
577 template <class T>
Put(Vector3D<T> * v,double vx,double vy,double vz)578 void Put(Vector3D<T> *v, double vx, double vy, double vz)
579 {
580 v->_x = static_cast<T>(vx);
581 v->_y = static_cast<T>(vy);
582 v->_z = static_cast<T>(vz);
583 }
584
585 template <class VectorType>
operator()586 void operator()(int i, int j, int k, int, VectorType *v)
587 {
588 double vx = i, vy = j, vz = k;
589 _Output->ImageToWorld (vx, vy, vz);
590 _Input ->WorldToLattice(vx, vy, vz);
591 _Input ->Evaluate (vx, vy, vz);
592 Put(v, vx, vy, vz);
593 }
594
595 protected:
596 const BSplineFreeFormTransformationSV *_Input; ///< Input transformation
597 BaseImage *_Output; ///< Output image
598 };
599
600 // -----------------------------------------------------------------------------
601 /// Evaluate B-spline SV FFD at image voxels
602 class EvaluateBSplineSVFFD3D : public VoxelFunction
603 {
604 public:
605
EvaluateBSplineSVFFD3D(const BSplineFreeFormTransformationSV * input,BaseImage * output)606 EvaluateBSplineSVFFD3D(const BSplineFreeFormTransformationSV *input, BaseImage *output)
607 :
608 _Input (input),
609 _Output(output),
610 _y (output->NumberOfSpatialVoxels()),
611 _z (2 * _y)
612 {}
613
614 template <class T>
operator()615 void operator()(int i, int j, int k, int, T *out)
616 {
617 double x = i, y = j, z = k;
618 _Output->ImageToWorld (x, y, z);
619 _Input ->WorldToLattice(x, y, z);
620 _Input ->Evaluate (x, y, z);
621 out[_x] = static_cast<T>(x);
622 out[_y] = static_cast<T>(y);
623 out[_z] = static_cast<T>(z);
624 }
625
626 private:
627 const BSplineFreeFormTransformationSV *_Input; ///< Input transformation
628 BaseImage *_Output; ///< Output image
629
630 static const int _x = 0; ///< Offset of x component
631 int _y; ///< Offset of y component
632 int _z; ///< Offset of z component
633 };
634
635 // -----------------------------------------------------------------------------
636 /// Evaluate and add B-spline SV FFD at image voxels
637 class AddBSplineSVFFD : public VoxelFunction
638 {
639 public:
640
AddBSplineSVFFD(const BSplineFreeFormTransformationSV * input,BaseImage * output)641 AddBSplineSVFFD(const BSplineFreeFormTransformationSV *input,
642 BaseImage *output)
643 :
644 _Input(input), _Output(output)
645 {}
646
647 // template <class T>
648 // void Add(Vector2D<T> *v, double vx, double vy, double)
649 // {
650 // v->_x += static_cast<T>(vx);
651 // v->_y += static_cast<T>(vy);
652 // }
653
654 template <class T>
Add(Vector3D<T> * v,double vx,double vy,double vz)655 void Add(Vector3D<T> *v, double vx, double vy, double vz)
656 {
657 v->_x += static_cast<T>(vx);
658 v->_y += static_cast<T>(vy);
659 v->_z += static_cast<T>(vz);
660 }
661
662 template <class VectorType>
operator()663 void operator()(int i, int j, int k, int, VectorType *v)
664 {
665 double vx = i, vy = j, vz = k;
666 _Output->ImageToWorld (vx, vy, vz);
667 _Input ->WorldToLattice(vx, vy, vz);
668 _Input ->Evaluate (vx, vy, vz);
669 Add(v, vx, vy, vz);
670 }
671
672 protected:
673 const BSplineFreeFormTransformationSV *_Input; ///< Input transformation
674 BaseImage *_Output; ///< Output image
675 };
676
677 // -----------------------------------------------------------------------------
678 /// Evaluate and add 3D B-spline SV FFD at image voxels
679 class AddBSplineSVFFD3D : public VoxelFunction
680 {
681 public:
682
AddBSplineSVFFD3D(const BSplineFreeFormTransformationSV * input,BaseImage * output)683 AddBSplineSVFFD3D(const BSplineFreeFormTransformationSV *input, BaseImage *output)
684 :
685 _Input (input),
686 _Output(output),
687 _y (output->NumberOfSpatialVoxels()),
688 _z (2 * _y)
689 {}
690
691 template <class T>
operator()692 void operator()(int i, int j, int k, int, T *out)
693 {
694 double x = i, y = j, z = k;
695 _Output->ImageToWorld (x, y, z);
696 _Input ->WorldToLattice(x, y, z);
697 _Input ->Evaluate (x, y, z);
698 out[_x] += static_cast<T>(x);
699 out[_y] += static_cast<T>(y);
700 out[_z] += static_cast<T>(z);
701 }
702
703 private:
704 const BSplineFreeFormTransformationSV *_Input; ///< Input transformation
705 BaseImage *_Output; ///< Output image
706
707 static const int _x = 0; ///< Offset of x component
708 int _y; ///< Offset of y component
709 int _z; ///< Offset of z component
710 };
711
712 ////////////////////////////////////////////////////////////////////////////////
713 // Inline definitions
714 ////////////////////////////////////////////////////////////////////////////////
715
716 // =============================================================================
717 // Evaluation
718 // =============================================================================
719
720 // -----------------------------------------------------------------------------
Evaluate(double & x,double & y,double & z,double)721 inline void BSplineFreeFormTransformationSV::Evaluate(double &x, double &y, double &z, double) const
722 {
723 BSplineFreeFormTransformation3D::Evaluate(x, y, z);
724 }
725
726 // -----------------------------------------------------------------------------
EvaluateJacobianWorld(Matrix & jac,double x,double y,double z,double)727 inline void BSplineFreeFormTransformationSV::EvaluateJacobianWorld(Matrix &jac, double x, double y, double z, double) const
728 {
729 if (_z == 1) BSplineFreeFormTransformation3D::EvaluateJacobianWorld(jac, x, y);
730 else BSplineFreeFormTransformation3D::EvaluateJacobianWorld(jac, x, y, z);
731 }
732
733 // -----------------------------------------------------------------------------
734 inline void BSplineFreeFormTransformationSV
EvaluateJacobianDOFs(TransformationJacobian & jac,double x,double y,double z,double)735 ::EvaluateJacobianDOFs(TransformationJacobian &jac, double x, double y, double z, double) const
736 {
737 if (_z == 1) EvaluateJacobianDOFs(jac, x, y);
738 else EvaluateJacobianDOFs(jac, x, y, z);
739 }
740
741 // -----------------------------------------------------------------------------
Velocity(GenericImage<float> & v)742 inline void BSplineFreeFormTransformationSV::Velocity(GenericImage<float> &v) const
743 {
744 if (v.T() != 3) {
745 Throw(ERR_InvalidArgument, __FUNCTION__, "Output image must have 3 channels!");
746 }
747 v.PutTSize(0.);
748 ParallelForEachVoxel(EvaluateBSplineSVFFD3D(this, &v), v.Attributes(), v);
749 }
750
751 // -----------------------------------------------------------------------------
Velocity(GenericImage<double> & v)752 inline void BSplineFreeFormTransformationSV::Velocity(GenericImage<double> &v) const
753 {
754 if (v.T() != 3) {
755 Throw(ERR_InvalidArgument, __FUNCTION__, "Output image must have 3 channels!");
756 }
757 v.PutTSize(0.);
758 ParallelForEachVoxel(EvaluateBSplineSVFFD3D(this, &v), v.Attributes(), v);
759 }
760
761 // =============================================================================
762 // Point transformation
763 // =============================================================================
764
765 // -----------------------------------------------------------------------------
RequiresCachingOfDisplacements()766 inline bool BSplineFreeFormTransformationSV::RequiresCachingOfDisplacements() const
767 {
768 return (_IntegrationMethod == FFDIM_SS || _IntegrationMethod == FFDIM_FastSS);
769 }
770
771 // -----------------------------------------------------------------------------
UpperIntegrationLimit(double t,double t0)772 inline double BSplineFreeFormTransformationSV::UpperIntegrationLimit(double t, double t0) const
773 {
774 double T = t - t0;
775 return !IsNaN(T) && !AreEqual(T, 0.) ? T : _T; // if zero/NaN, return cross-sectional interval _T instead
776 }
777
778 // -----------------------------------------------------------------------------
NumberOfStepsForIntervalLength(double T)779 inline int BSplineFreeFormTransformationSV::NumberOfStepsForIntervalLength(double T) const
780 {
781 return (_TimeUnit > .0) ? static_cast<int>((abs(T) / _TimeUnit) * _NumberOfSteps + 0.5) : _NumberOfSteps;
782 }
783
784 // -----------------------------------------------------------------------------
NumberOfStepsForIntervalLength(double T,int nsteps)785 inline void BSplineFreeFormTransformationSV::NumberOfStepsForIntervalLength(double T, int nsteps) const
786 {
787 if (_TimeUnit > .0) nsteps = static_cast<int>((_TimeUnit / abs(T)) * nsteps + 0.5);
788 const_cast<BSplineFreeFormTransformationSV *>(this)->_NumberOfSteps = nsteps;
789 }
790
791 // -----------------------------------------------------------------------------
StepLengthForIntervalLength(double T)792 inline double BSplineFreeFormTransformationSV::StepLengthForIntervalLength(double T) const
793 {
794 if (_NumberOfSteps == 0) return .0;
795 return ((_TimeUnit > .0) ? (_TimeUnit / static_cast<double>(_NumberOfSteps))
796 : ( abs(T) / static_cast<double>(_NumberOfSteps)));
797 }
798
799 // =============================================================================
800 // Derivatives
801 // =============================================================================
802
803 // -----------------------------------------------------------------------------
804 inline void BSplineFreeFormTransformationSV
ParametricGradient(const GenericImage<double> * in,double * out,const WorldCoordsImage * i2w,const WorldCoordsImage * wc,double t0,double w)805 ::ParametricGradient(const GenericImage<double> *in, double *out,
806 const WorldCoordsImage *i2w, const WorldCoordsImage *wc,
807 double t0, double w) const
808 {
809 this->ParametricGradient(in, out, i2w, wc, in->GetTOrigin(), t0, w);
810 }
811
812
813 } // namespace mirtk
814
815 #endif // MIRTK_BSplineFreeFormTransformationSV_H
816