1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2008-2017 Imperial College London
5 * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel
6 * Copyright 2013-2017 Andreas Schuh
7 *
8 * Licensed under the Apache License, Version 2.0 (the "License");
9 * you may not use this file except in compliance with the License.
10 * You may obtain a copy of the License at
11 *
12 * http://www.apache.org/licenses/LICENSE-2.0
13 *
14 * Unless required by applicable law or agreed to in writing, software
15 * distributed under the License is distributed on an "AS IS" BASIS,
16 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 * See the License for the specific language governing permissions and
18 * limitations under the License.
19 */
20
21 #ifndef MIRTK_FreeFormTransformation_H
22 #define MIRTK_FreeFormTransformation_H
23
24 #include "mirtk/Transformation.h"
25
26 #include "mirtk/Math.h"
27 #include "mirtk/Memory.h"
28 #include "mirtk/Vector3D.h"
29 #include "mirtk/InterpolateImageFunction.h"
30 #include "mirtk/ExtrapolateImageFunction.h"
31 #include "mirtk/TransformationJacobian.h"
32 #include "mirtk/ImageFunction.h"
33
34
35 namespace mirtk {
36
37
38 /**
39 * Base class for free-form transformations.
40 */
41 class FreeFormTransformation : public Transformation
42 {
43 mirtkAbstractTransformationMacro(FreeFormTransformation);
44
45 public:
46
47 /// Default attributes of free-form transformation lattice
48 ///
49 /// \param[in] attr (Target) image attributes.
50 /// \param[in] dx Control point spacing in x.
51 /// \param[in] dy Control point spacing in y.
52 /// \param[in] dz Control point spacing in z.
53 /// \param[in] dt Control point spacing in t.
54 ///
55 /// \note If a control point number/spacing is negative, it is set to the specified number/size.
56 static ImageAttributes DefaultAttributes(const ImageAttributes &attr,
57 double dx = -1.0,
58 double dy = -1.0,
59 double dz = -1.0,
60 double dt = -1.0);
61
62 /// Default attributes of free-form transformation lattice
63 static ImageAttributes DefaultAttributes(double, double, double, double,
64 double, double, double, double,
65 double, double, double, double,
66 const double *, const double *,
67 const double *);
68
69 // ---------------------------------------------------------------------------
70 // Types
71
72 /// Type of vector storing the coefficients at a control point
73 typedef Vector3D<DOFValue> CPValue;
74
75 /// Type of image representation of free-form transformation
76 typedef GenericImage<CPValue> CPImage;
77
78 /// Type of vector storing status of control point coefficients
79 typedef Vector3D<DOFStatus> CPStatus;
80
81 /// Base type used for interpolation of control point data
82 /// \note Subclasses may use a particular specialized interpolation type.
83 typedef GenericInterpolateImageFunction<CPImage> CPInterpolator;
84
85 /// Base type used for extrapolation of control point data
86 typedef GenericExtrapolateImageFunction<CPImage> CPExtrapolator;
87
88 /// Alias for CPValue type which is preferred when referring to any FFD vector,
89 /// not only those at control point locations
90 typedef CPValue Vector;
91
92 // ---------------------------------------------------------------------------
93 // Attributes
94
95 /// Mode used for extrapolation of control point image
96 mirtkReadOnlyAttributeMacro(enum ExtrapolationMode, ExtrapolationMode);
97
98 /// Speedup factor for gradient computation
99 mirtkPublicAttributeMacro(double, SpeedupFactor);
100
101 public:
102
103 /// Set extrapolation mode
104 virtual void ExtrapolationMode(enum ExtrapolationMode);
105
106 /// Get access to control point coefficients extrapolator
107 /// Intended for use by non-member auxiliary functions of subclass implementations.
108 const CPExtrapolator *Extrapolator() const;
109
110 protected:
111
112 /// Finite discrete representation of free-form transformation coefficients
113 ///
114 /// Note that this image instance stores all the attributes of the control
115 /// point lattice. The actual control point data is stored in the parameters
116 /// of the transformation, the memory of which is managed by the base class.
117 /// The image only wraps this memory and interprets it as an image. This allows
118 /// the use of image functions such as interpolators and efficient image filters.
119 CPImage _CPImage;
120
121 /// Infinite discrete representation of free-form transformation coefficients
122 CPExtrapolator *_CPValue;
123
124 /// Infinite continuous function of free-form transformation
125 ///
126 /// Note that the actual interpolator object is (at the moment) attribute of
127 /// the specific FFD implementation. The base class pointer here is only used
128 /// to update the interpolator when necessary, e.g., when the extrapolation
129 /// mode has changed. It is set by the base class constructor to the
130 /// respective interpolate image function attribute of the subclass.
131 CPInterpolator *_CPFunc;
132 CPInterpolator *_CPFunc2D; // TODO: Remove when 2D FFDs are separate classes
133
134 /// Status of control points
135 ///
136 /// Note that the actual status values are stored by the base class member
137 /// Transformation::_Status in a contiguous memory block.
138 CPStatus ****_CPStatus;
139
140 const ImageAttributes &_attr; ///< Control point lattice attributes
141
142 const int &_x; ///< Read-only reference to _x attribute of _CPImage
143 const int &_y; ///< Read-only reference to _y attribute of _CPImage
144 const int &_z; ///< Read-only reference to _z attribute of _CPImage
145 const int &_t; ///< Read-only reference to _t attribute of _CPImage
146
147 const double &_dx; ///< Read-only reference to _dx attribute of _CPImage
148 const double &_dy; ///< Read-only reference to _dy attribute of _CPImage
149 const double &_dz; ///< Read-only reference to _dz attribute of _CPImage
150 const double &_dt; ///< Read-only reference to _dt attribute of _CPImage
151
152 const Matrix &_matW2L; ///< Read-only reference to _matW2I of _CPImage
153 const Matrix &_matL2W; ///< Read-only reference to _matI2W of _CPImage
154
155 // ---------------------------------------------------------------------------
156 // Construction/Destruction
157
158 /// Initialize interpolator of control points
159 void InitializeInterpolator();
160
161 /// Initialize status of control points
162 void InitializeStatus();
163
164 /// Initialize control points
165 void InitializeCPs(const ImageAttributes &, bool = true);
166
167 /// Copy control points from other transformation
168 void InitializeCPs(const FreeFormTransformation &, bool = true);
169
170 /// Default constructor
171 FreeFormTransformation(CPInterpolator &, CPInterpolator * = NULL);
172
173 /// Copy constructor
174 FreeFormTransformation(const FreeFormTransformation &,
175 CPInterpolator &, CPInterpolator * = NULL);
176
177 public:
178
179 /// Destructor
180 virtual ~FreeFormTransformation();
181
182 /// Initialize free-form transformation
183 virtual void Initialize(const ImageAttributes &) = 0;
184
185 /// Initialize free-form transformation
186 void Initialize(const ImageAttributes &, double, double, double = -1.0, double = -1.0);
187
188 /// Initialize transformation from existing vector field
189 void Initialize(const CPImage &, bool = false);
190
191 /// Initialize transformation from existing 3D+t vector field
192 void Initialize(const GenericImage<double> &, bool = false);
193
194 // ---------------------------------------------------------------------------
195 // Parameters (non-DoFs)
196
197 // Import other Parameter overloads
198 using Transformation::Parameter;
199
200 /// Set named (non-DoF) parameter from value as string
201 virtual bool Set(const char *, const char *);
202
203 /// Get (non-DoF) parameters as key/value as string map
204 virtual ParameterList Parameter() const;
205
206 // ---------------------------------------------------------------------------
207 // Approximation/Interpolation
208
209 using Transformation::EvaluateRMSError;
210 using Transformation::Approximate;
211 using Transformation::ApproximateAsNew;
212 using Transformation::ApproximateGradient;
213
214 /// Get image domain on which this free-form transformation should be defined
215 /// in order to reduce the error when the given transformation is approximated
216 /// and the resulting free-form transformation is applied to resample an image
217 /// that is defined on the specified lattice.
218 virtual ImageAttributes ApproximationDomain(const ImageAttributes &,
219 const Transformation *);
220
221 /// Evaluates RMS error of transformation at control points compared to another.
222 double EvaluateRMSError(const Transformation *) const;
223
224 /// Approximate transformation: This function takes a transformation and finds
225 /// a FFD which approximates this transformation at the control point locations.
226 /// Returns the approximation error of the resulting FFD at the control point locations.
227 virtual double Approximate(const Transformation *, int = 1, double = .0);
228
229 /// Approximate displacements: This function takes a set of points and a set
230 /// of displacements and finds a transformation which approximates these
231 /// displacements. After approximation, the displacements are replaced by the
232 /// residual displacement errors at the points.
233 virtual double Approximate(const ImageAttributes &, double *, double *, double *,
234 int = 1, double = .0);
235
236 /// Approximate displacements: This function takes a set of points and a set
237 /// of displacements and finds a transformation which approximates these
238 /// displacements. After approximation, the displacements are replaced by the
239 /// residual displacement errors at the points.
240 virtual double Approximate(const double *, const double *, const double *,
241 double *, double *, double *, int,
242 int = 1, double = .0);
243
244 /// Approximate displacements: This function takes a set of points and a set
245 /// of displacements and finds a transformation which approximates these
246 /// displacements. After approximation, the displacements are replaced by the
247 /// residual displacement errors at the points.
248 virtual double Approximate(const double *, const double *, const double *, const double *,
249 double *, double *, double *, int,
250 int = 1, double = .0);
251
252 /// Approximate another transformation and return approximation error
253 virtual double ApproximateAsNew(const Transformation *, int = 1, double = .0);
254
255 /// Approximate displacements: This function takes a set of points and a set
256 /// of displacements and finds a !new! transformation which approximates these
257 /// displacements. After approximation, the displacements are replaced by the
258 /// residual displacement errors at the points.
259 virtual double ApproximateAsNew(const ImageAttributes &, double *, double *, double *,
260 int = 1, double = .0);
261
262 /// Interpolates displacements: This function takes a set of displacements defined
263 /// at the control points and finds a FFD which interpolates these displacements.
264 virtual void Interpolate(const double *, const double *, const double *) = 0;
265
266 // ---------------------------------------------------------------------------
267 // Lattice
268
269 /// Returns attributes of control point grid
270 const ImageAttributes &Attributes() const;
271
272 /// Actual number of DoFs, i.e.,
273 /// 2 * NumberOfDOFs for 2D FFD or 3 * NumberOfDOFs for 3D(+t) FFD
274 int ActualNumberOfDOFs() const;
275
276 /// Number of control points
277 int NumberOfCPs() const;
278
279 /// Number of active control points
280 int NumberOfActiveCPs() const;
281
282 /// Number of non-active control points
283 int NumberOfPassiveCPs() const;
284
285 /// Returns the number of control points in x
286 int X() const;
287
288 /// Returns the number of control points in y
289 int Y() const;
290
291 /// Returns the number of control points in z
292 int Z() const;
293
294 /// Returns the number of control points in t
295 int T() const;
296
297 /// Returns the number of control points in x
298 int GetX() const;
299
300 /// Returns the number of control points in y
301 int GetY() const;
302
303 /// Returns the number of control points in z
304 int GetZ() const;
305
306 /// Returns the number of control points in t
307 int GetT() const;
308
309 /// Returns the of control point spacing in x
310 double GetXSpacing() const;
311
312 /// Returns the of control point spacing in y
313 double GetYSpacing() const;
314
315 /// Returns the of control point spacing in z
316 double GetZSpacing() const;
317
318 /// Returns the of control point spacing in t
319 double GetTSpacing() const;
320
321 /// Gets the control point spacing (in mm)
322 void GetSpacing(double &, double &, double &) const;
323
324 /// Gets the control point spacing (in mm)
325 void GetSpacing(double &, double &, double &, double &) const;
326
327 /// Puts the orientation of the free-form deformation lattice
328 void PutOrientation(double *, double *, double *);
329
330 /// Gets the orientation of the free-form deformation lattice
331 void GetOrientation(double *, double *, double *) const;
332
333 /// Get indices of transformation parameters (DoFs)
334 void IndexToDOFs(int, int &, int &) const;
335
336 /// Get indices of transformation parameters (DoFs)
337 void IndexToDOFs(int, int &, int &, int &) const;
338
339 /// Get index of control point corresponding to transformation parameter (DoFs)
340 int DOFToIndex(int) const;
341
342 /// Get index of dimension corresponding to transformation parameter (DoFs)
343 int DOFToDimension(int) const;
344
345 /// Get control point index from lattice coordinates
346 int LatticeToIndex(int, int, int = 0, int = 0) const;
347
348 /// Get control point lattice coordinates from index
349 void IndexToLattice(int, int &, int &) const;
350
351 /// Get control point lattice coordinates from index
352 void IndexToLattice(int, int &, int &, int &) const;
353
354 /// Get control point lattice coordinates from index
355 void IndexToLattice(int, int &, int &, int &, int &) const;
356
357 /// Get world coordinates (in mm) of control point
358 void IndexToWorld(int, double &, double &) const;
359
360 /// Get world coordinates (in mm) of control point
361 void IndexToWorld(int, double &, double &, double &) const;
362
363 /// Get world coordinates (in mm) of control point
364 void IndexToWorld(int, Point &) const;
365
366 /// Get world coordinates (in mm) of control point
367 Point IndexToWorld(int) const;
368
369 /// Transforms world coordinates (in mm) to lattice coordinates
370 virtual void WorldToLattice(double &, double &) const;
371
372 /// Transforms world coordinates (in mm) to lattice coordinates
373 virtual void WorldToLattice(double &, double &, double &) const;
374
375 /// Transforms world coordinates (in mm) to lattice coordinates
376 virtual void WorldToLattice(Point &) const;
377
378 /// Transforms lattice coordinates to world coordinates (in mm)
379 virtual void LatticeToWorld(double &, double &) const;
380
381 /// Transforms lattice coordinates to world coordinates (in mm)
382 virtual void LatticeToWorld(double &, double &, double &) const;
383
384 /// Transforms lattice coordinates to world coordinates (in mm)
385 virtual void LatticeToWorld(Point &) const;
386
387 /// Gets the location of the given control point (in mm)
388 void ControlPointLocation(int, double &, double &) const;
389
390 /// Gets the location of the given control point (in mm)
391 void ControlPointLocation(int, double &, double &, double &) const;
392
393 /// Returns the location of the given control point (in mm)
394 Point ControlPointLocation(int) const;
395
396 /// Transforms time (in ms) to temporal lattice coordinate
397 virtual double TimeToLattice(double) const;
398
399 /// Transforms temporal lattice coordinate to time (in ms)
400 virtual double LatticeToTime(double) const;
401
402 /// Returns the number of control points in x after subdivision
403 virtual int GetXAfterSubdivision() const;
404
405 /// Returns the number of control points in y after subdivision
406 virtual int GetYAfterSubdivision() const;
407
408 /// Returns the number of control points in z after subdivision
409 virtual int GetZAfterSubdivision() const;
410
411 /// Returns the number of control points in t after subdivision
412 virtual int GetTAfterSubdivision() const;
413
414 /// Returns the control point spacing in x after the subdivision
415 virtual double GetXSpacingAfterSubdivision() const;
416
417 /// Returns the control point spacing in y after the subdivision
418 virtual double GetYSpacingAfterSubdivision() const;
419
420 /// Returns the control point spacing in z after the subdivision
421 virtual double GetZSpacingAfterSubdivision() const;
422
423 /// Returns the control point spacing in t after the subdivision
424 virtual double GetTSpacingAfterSubdivision() const;
425
426 /// Subdivide lattice of free-form transformation
427 virtual void Subdivide(bool = true, bool = true, bool = true, bool = true);
428
429 /// Subdivide lattice in first two dimensions
430 void Subdivide2D();
431
432 /// Subdivide lattice in first three dimensions
433 void Subdivide3D();
434
435 /// Subdivide lattice in all four dimensions
436 void Subdivide4D();
437
438 /// Crop/pad lattice to discard passive control points at the boundary,
439 /// keeping only a layer of passive control points of given width.
440 /// The DoF values of passive control points are optionally reset to zero.
441 bool CropPadPassiveCPs(int = 0, bool = false);
442
443 /// Crop/pad lattice to discard passive control points at the boundary,
444 /// keeping only a layer of passive control points of given width.
445 /// The DoF values of passive control points are optionally reset to zero.
446 virtual bool CropPadPassiveCPs(int, int, int = 0, int = 0, bool = false);
447
448 // ---------------------------------------------------------------------------
449 // Bounding box
450
451 /// Size of support region of the used kernel
452 virtual int KernelSize() const = 0;
453
454 /// Radius of support region of the used kernel
455 int KernelRadius() const;
456
457 /// Puts the spatial bounding box for the free-form deformation (in mm)
458 void PutBoundingBox(double, double, double,
459 double, double, double);
460
461 /// Puts the temporal bounding box for the free-form deformation (in mm)
462 void PutBoundingBox(const Point &, const Point &);
463
464 /// Puts the temporal bounding box of the free-form deformation (in ms)
465 void PutBoundingBox(double, double);
466
467 /// Puts the spatio-temporal bounding box for the free-form deformation (in mm)
468 void PutBoundingBox(double, double, double, double,
469 double, double, double, double);
470
471 /// Gets the temporal bounding box of the free-form deformation (in ms)
472 void BoundingBox(double &, double &) const;
473
474 /// Gets the spatial bounding box of the free-form deformation (in mm)
475 void BoundingBox(double &, double &, double &,
476 double &, double &, double &) const;
477
478 /// Gets the spatial bounding box of the free-form deformation (in mm)
479 void BoundingBox(Point &, Point &) const;
480
481 /// Gets the spatio-temporal bounding box of the free-form deformation (in mm and ms)
482 void BoundingBox(double &, double &, double &, double &,
483 double &, double &, double &, double &) const;
484
485 /// Gets the spatio-temporal bounding box of the free-form deformation (in mm and ms)
486 void BoundingBox(Point &, double &, Point &, double &) const;
487
488 /// Gets the temporal bounding box for a control point. The last parameter
489 /// specifies what fraction of the bounding box to return. The default
490 /// is 1 which equals 100% of the bounding box.
491 virtual void BoundingBox(int, double &, double &, double = 1) const;
492
493 /// Gets the spatial bounding box for a control point. The last parameter
494 /// specifies what fraction of the bounding box to return. The default
495 /// is 1 which equals 100% of the bounding box.
496 virtual void BoundingBox(int, double &, double &, double &,
497 double &, double &, double &, double = 1) const = 0;
498
499 /// Gets the spatio-temporal bounding box for a control point. The last parameter
500 /// specifies what fraction of the bounding box to return. The default
501 /// is 1 which equals 100% of the bounding box.
502 virtual void BoundingBox(int, double &, double &, double &, double &,
503 double &, double &, double &, double &, double = 1) const;
504
505 /// Gets the spatial bounding box for a control point. The last parameter
506 /// specifies what fraction of the bounding box to return. The default
507 /// is 1 which equals 100% of the bounding box.
508 void BoundingBox(int, Point &, Point &, double = 1) const;
509
510 /// Gets the spatial bounding box for a control point in lattice coordinates.
511 /// The last parameter specifies what fraction of the bounding box to return.
512 /// The default is 1 which equals 100% of the bounding box.
513 bool BoundingBox(const ImageAttributes &, int, int &, int &, int &,
514 int &, int &, int &, double = 1) const;
515
516 /// Gets the spatio-temporal bounding box for a control point in lattice coordinates.
517 /// The last parameter specifies what fraction of the bounding box to return.
518 /// The default is 1 which equals 100% of the bounding box.
519 bool BoundingBox(const ImageAttributes &, int, int &, int &, int &, int &,
520 int &, int &, int &, int &, double = 1) const;
521
522 /// Gets the spatial bounding box for a transformation parameter in lattice coordinates.
523 /// The last parameter specifies what fraction of the bounding box to return.
524 /// The default is 1 which equals 100% of the bounding box.
525 bool DOFBoundingBox(const ImageAttributes &, int, int &, int &, int &,
526 int &, int &, int &, double = 1) const;
527
528 /// Gets the spatial bounding box for a control point in image coordinates.
529 /// The last parameter specifies what fraction of the bounding box to return.
530 /// The default is 1 which equals 100% of the bounding box.
531 bool BoundingBox(const Image *, int, int &, int &, int &,
532 int &, int &, int &, double = 1) const;
533
534 /// Gets the spatio-temporal bounding box for a control point in image coordinates.
535 /// The last parameter specifies what fraction of the bounding box to return.
536 /// The default is 1 which equals 100% of the bounding box.
537 bool BoundingBox(const Image *, int, int &, int &, int &, int &,
538 int &, int &, int &, int &, double = 1) const;
539
540 /// Gets the spatial bounding box for a transformation parameter in image coordinates.
541 /// The last parameter specifies what fraction of the bounding box to return.
542 /// The default is 1 which equals 100% of the bounding box.
543 bool DOFBoundingBox(const Image *, int, int &, int &, int &,
544 int &, int &, int &, double = 1) const;
545
546 // ---------------------------------------------------------------------------
547 // Transformation parameters (DoFs)
548 using Transformation::Put;
549 using Transformation::Get;
550 using Transformation::PutStatus;
551 using Transformation::GetStatus;
552
553 /// Get norm of the gradient vector
554 virtual double DOFGradientNorm(const double *) const;
555
556 /// Copy active transformation parameters (DoFs) from given
557 /// transformation if possible and return \c false, otherwise
558 virtual bool CopyFrom(const Transformation *);
559
560 /// Puts values of the parameters at a control point
561 void Put(int, const Vector &);
562
563 /// Puts values of the parameters at a control point
564 void Put(int, double, double, double);
565
566 /// Puts values of the parameters at a control point
567 void Put(int, int, int, double, double, double);
568
569 /// Puts values of the parameters at a control point
570 void Put(int, int, int, int, double, double, double);
571
572 /// Gets values of the parameters at a control point
573 void Get(int, Vector &) const;
574
575 /// Gets values of the parameters at a control point
576 void Get(int, double &, double &, double &) const;
577
578 /// Gets values of the parameters at a control point
579 void Get(int, int, int, double &, double &, double &) const;
580
581 /// Gets values of the parameters at a control point
582 void Get(int, int, int, int, double &, double &, double &) const;
583
584 /// Puts status of the parameters at a control point
585 void PutStatus(int, const CPStatus &);
586
587 /// Puts status of the parameters at a control point
588 void PutStatus(int, int, int, DOFStatus, DOFStatus, DOFStatus);
589
590 /// Puts status of the parameters at a control point
591 void PutStatus(int, int, int, int, DOFStatus, DOFStatus, DOFStatus);
592
593 /// Gets status of the parameters at a control point
594 void GetStatus(int, CPStatus &) const;
595
596 /// Gets status of the parameters at a control point
597 void GetStatus(int, int, int, DOFStatus &, DOFStatus &, DOFStatus &) const;
598
599 /// Gets status of the parameters at a control point
600 void GetStatus(int, int, int, int, DOFStatus &, DOFStatus &, DOFStatus &) const;
601
602 /// Whether the control point at given lattice index is active
603 virtual bool IsActive(int) const;
604
605 /// Whether the control point at given lattice coordinates is active
606 virtual bool IsActive(int, int, int = 0, int = 0) const;
607
608 // ---------------------------------------------------------------------------
609 // Point transformation
610 using Transformation::Transform;
611 using Transformation::Inverse;
612
613 /// Transforms a single point using the global transformation component only
614 virtual void GlobalTransform(double &, double &, double &, double = 0, double = NaN) const;
615
616 /// Transforms a single point
617 virtual void Transform(double &, double &, double &, double = 0, double = NaN) const;
618
619 /// Transforms a single point using the inverse of the global transformation only
620 virtual void GlobalInverse(double &, double &, double &, double = 0, double = NaN) const;
621
622 /// Transforms a single point using the inverse of the transformation
623 virtual bool Inverse(double &, double &, double &, double = 0, double = NaN) const;
624
625 // ---------------------------------------------------------------------------
626 // Derivatives
627 using Transformation::Jacobian;
628 using Transformation::GlobalJacobian;
629 using Transformation::JacobianDOFs;
630 using Transformation::ParametricGradient;
631
632 /// Convert 1st order derivatives computed w.r.t 2D lattice coordinates to
633 /// derivatives w.r.t world coordinates
634 void JacobianToWorld(double &, double &) const;
635
636 /// Convert 1st order derivatives computed w.r.t 3D lattice coordinates to
637 /// derivatives w.r.t world coordinates
638 void JacobianToWorld(double &, double &, double &) const;
639
640 /// Convert 1st order derivatives computed w.r.t lattice coordinates to
641 /// derivatives w.r.t world coordinates
642 void JacobianToWorld(Matrix &) const;
643
644 /// Reorient 1st order derivatives computed w.r.t 2D lattice coordinates
645 void JacobianToWorldOrientation(double &, double &) const;
646
647 /// Reorient 1st order derivatives computed w.r.t 2D lattice coordinates
648 void JacobianToWorldOrientation(double &, double &, double &) const;
649
650 /// Convert 2nd order derivatives computed w.r.t 2D lattice coordinates to
651 /// derivatives w.r.t world coordinates
652 void HessianToWorld(double &, double &, double &) const;
653
654 /// Convert 2nd order derivatives computed w.r.t 3D lattice coordinates to
655 /// derivatives w.r.t world coordinates
656 void HessianToWorld(double &, double &, double &, double &, double &, double &) const;
657
658 /// Convert 2nd order derivatives of single transformed coordinate computed
659 /// w.r.t lattice coordinates to derivatives w.r.t world coordinates
660 void HessianToWorld(Matrix &) const;
661
662 /// Convert 2nd order derivatives computed w.r.t lattice coordinates to
663 /// derivatives w.r.t world coordinates
664 void HessianToWorld(Matrix [3]) const;
665
666 /// Calculates the Jacobian of the transformation w.r.t either control point displacements or velocities
667 virtual void FFDJacobianWorld(Matrix &, double, double, double, double = 0, double = NaN) const;
668
669 /// Calculates the Jacobian of the global transformation w.r.t world coordinates
670 virtual void GlobalJacobian(Matrix &, double, double, double, double = 0, double = NaN) const;
671
672 /// Calculates the Jacobian of the transformation w.r.t world coordinates
673 virtual void Jacobian(Matrix &, double, double, double, double = 0, double = NaN) const;
674
675 /// Calculates the Hessian for each component of the global transformation w.r.t world coordinates
676 virtual void GlobalHessian(Matrix [3], double, double, double, double = 0, double = NaN) const;
677
678 /// Calculates the Hessian for each component of the transformation w.r.t world coordinates
679 virtual void Hessian(Matrix [3], double, double, double, double = 0, double = NaN) const;
680
681 /// Calculates the Jacobian of the transformation w.r.t the transformation parameters of a control point
682 virtual void JacobianDOFs(Matrix &, int, double, double, double, double = 0, double = NaN) const;
683
684 /// Calculates the Jacobian of the transformation w.r.t the transformation parameters
685 virtual void JacobianDOFs(TransformationJacobian &, double, double, double, double = 0, double = NaN) const;
686
687 /// Calculates derivatives of the Jacobian determinant of spline function w.r.t. DoFs of a control point
688 ///
689 /// This function is identical to JacobianDetDerivative when the DoFs of the control points are displacements.
690 /// When the DoFs are velocities, however, this function computes the derivatives of the Jacobian determinant
691 /// of the velocity field instead.
692 ///
693 /// \param[out] dJ Partial derivatives of Jacobian determinant at (x, y, z) w.r.t. DoFs of control point.
694 /// \param[in] cp Index of control point w.r.t. whose DoFs the derivatives are computed.
695 /// \param[in] x World coordinate along x axis at which to evaluate derivatives.
696 /// \param[in] y World coordinate along y axis at which to evaluate derivatives.
697 /// \param[in] z World coordinate along z axis at which to evaluate derivatives.
698 /// \param[in] adj Adjugate of Jacobian matrix evaluated at (x, y, z).
699 /// \param[in] wrt_world Whether derivatives are computed w.r.t. world coordinate system.
700 /// \param[in] use_spacing Whether to use grid spacing when \p wrt_world is \c true.
701 virtual void FFDJacobianDetDerivative(double dJ[3], const Matrix &adj,
702 int cp, double x, double y, double z, double = 0, double = NaN,
703 bool wrt_world = true, bool use_spacing = true) const;
704
705 /// Calculates derivatives of the Jacobian determinant at world point w.r.t. DoFs of a control point
706 ///
707 /// \param[out] dJ Partial derivatives of Jacobian determinant w.r.t. DoFs of control point.
708 /// \param[in] adj Pre-computed adjugate of Jacobian matrix at this world point.
709 /// \param[in] cp Index of control point w.r.t. whose DoFs the derivatives are computed.
710 /// \param[in] x World coordinate along x axis at which to evaluate derivatives.
711 /// \param[in] y World coordinate along y axis at which to evaluate derivatives.
712 /// \param[in] z World coordinate along z axis at which to evaluate derivatives.
713 /// \param[in] t Temporal coordinate of point at which to evaluate derivatives.
714 /// \param[in] t0 Temporal coordinate of co-domain (target). Used by velocity-based models.
715 /// \param[in] wrt_world Whether derivatives are computed w.r.t. world coordinate system.
716 /// \param[in] use_spacing Whether to use grid spacing when \p wrt_world is \c true.
717 virtual void JacobianDetDerivative(double dJ[3], const Matrix &adj,
718 int cp, double x, double y, double z, double t = 0, double t0 = NaN,
719 bool wrt_world = true, bool use_spacing = true) const;
720
721 /// Applies the chain rule to convert spatial non-parametric gradient
722 /// to a gradient w.r.t the parameters of this transformation.
723 ///
724 /// If the transformation itself is non-parametric, the gradient will be passed through
725 /// unchanged. The default implementation uses the full Jacobian matrix computed for each
726 /// DoF separately (i.e., calls JacobianDOFs for each DoF).
727 ///
728 /// For 4D transformations, the temporal coordinate t used for the computation of the Jacobian
729 /// of the transformation w.r.t the transformation parameters for all spatial (x, y, z) voxel
730 /// coordinates in world units, is assumed to correspond to the temporal origin of the given
731 /// gradient image. For 4D transformations parameterized by velocities, a second time for the
732 /// upper integration bound can be provided as last argument to this method. This last
733 /// argument is ignored by transformations parameterized by displacements.
734 ///
735 /// \sa ImageSimilarityMetric::EvaluateGradient
736 virtual void ParametricGradient(const GenericImage<double> *, double *,
737 const WorldCoordsImage *,
738 const WorldCoordsImage *,
739 double = NaN, double = 1) const;
740
741 /// Applies the chain rule to convert point-wise non-parametric gradient
742 /// to a gradient w.r.t the parameters of this transformation.
743 virtual void ParametricGradient(const PointSet &, const Vector3D<double> *,
744 double *, double = 0, double = NaN, double = 1) const;
745
746 // ---------------------------------------------------------------------------
747 // Properties
748
749 /// Calculates the bending of the transformation given the 2nd order derivatives
750 static double Bending3D(const Matrix [3]);
751
752 /// Calculates the bending of the transformation
753 virtual double BendingEnergy(double, double, double, double = 0, double = NaN, bool = true) const;
754
755 /// Approximates the bending energy on the control point lattice
756 virtual double BendingEnergy(bool = false, bool = true) const;
757
758 /// Approximates the bending energy on the specified discrete domain
759 virtual double BendingEnergy(const ImageAttributes &attr, double = NaN, bool = true) const;
760
761 /// Approximates the gradient of the bending energy on the control point
762 /// lattice w.r.t the transformation parameters and adds it with the given weight
763 virtual void BendingEnergyGradient(double *, double = 1.0, bool = false, bool = true, bool = true) const;
764
765 // ---------------------------------------------------------------------------
766 // I/O
767
768 // Do not hide methods of base class
769 using Transformation::Print;
770
771 /// Prints the parameters of the transformation
772 virtual void Print(ostream &, Indent = 0) const;
773
774 protected:
775
776 /// Writes the control point and status information to a file stream
777 Cofstream &WriteCPs(Cofstream &) const;
778
779 };
780
781 ////////////////////////////////////////////////////////////////////////////////
782 // Inline definitions
783 ////////////////////////////////////////////////////////////////////////////////
784
785 // -----------------------------------------------------------------------------
Extrapolator()786 inline const FreeFormTransformation::CPExtrapolator *FreeFormTransformation::Extrapolator() const
787 {
788 return _CPValue;
789 }
790
791 // =============================================================================
792 // Lattice
793 // =============================================================================
794
795 // -----------------------------------------------------------------------------
KernelRadius()796 inline int FreeFormTransformation::KernelRadius() const
797 {
798 return this->KernelSize() / 2;
799 }
800
801 // -----------------------------------------------------------------------------
Attributes()802 inline const ImageAttributes &FreeFormTransformation::Attributes() const
803 {
804 return _CPImage.Attributes();
805 }
806
807 // -----------------------------------------------------------------------------
NumberOfCPs()808 inline int FreeFormTransformation::NumberOfCPs() const
809 {
810 return _x * _y * _z * _t;
811 }
812
813 // -----------------------------------------------------------------------------
NumberOfActiveCPs()814 inline int FreeFormTransformation::NumberOfActiveCPs() const
815 {
816 int nactive = 0;
817 for (int cp = 0; cp < NumberOfCPs(); ++cp) {
818 if (this->IsActive(cp)) ++nactive;
819 }
820 return nactive;
821 }
822
823 // -----------------------------------------------------------------------------
NumberOfPassiveCPs()824 inline int FreeFormTransformation::NumberOfPassiveCPs() const
825 {
826 return NumberOfCPs() - NumberOfActiveCPs();
827 }
828
829 // -----------------------------------------------------------------------------
ActualNumberOfDOFs()830 inline int FreeFormTransformation::ActualNumberOfDOFs() const
831 {
832 if (_z == 1) return 2 * NumberOfCPs();
833 else return 3 * NumberOfCPs();
834 }
835
836 // -----------------------------------------------------------------------------
X()837 inline int FreeFormTransformation::X() const
838 {
839 return _x;
840 }
841
842 // -----------------------------------------------------------------------------
Y()843 inline int FreeFormTransformation::Y() const
844 {
845 return _y;
846 }
847
848 // -----------------------------------------------------------------------------
Z()849 inline int FreeFormTransformation::Z() const
850 {
851 return _z;
852 }
853
854 // -----------------------------------------------------------------------------
T()855 inline int FreeFormTransformation::T() const
856 {
857 return _t;
858 }
859
860 // -----------------------------------------------------------------------------
GetX()861 inline int FreeFormTransformation::GetX() const
862 {
863 return X();
864 }
865
866 // -----------------------------------------------------------------------------
GetY()867 inline int FreeFormTransformation::GetY() const
868 {
869 return Y();
870 }
871
872 // -----------------------------------------------------------------------------
GetZ()873 inline int FreeFormTransformation::GetZ() const
874 {
875 return Z();
876 }
877
878 // -----------------------------------------------------------------------------
GetT()879 inline int FreeFormTransformation::GetT() const
880 {
881 return T();
882 }
883
884 // -----------------------------------------------------------------------------
GetXSpacing()885 inline double FreeFormTransformation::GetXSpacing() const
886 {
887 return _dx;
888 }
889
890 // -----------------------------------------------------------------------------
GetYSpacing()891 inline double FreeFormTransformation::GetYSpacing() const
892 {
893 return _dy;
894 }
895
896 // -----------------------------------------------------------------------------
GetZSpacing()897 inline double FreeFormTransformation::GetZSpacing() const
898 {
899 return _dz;
900 }
901
902 // -----------------------------------------------------------------------------
GetTSpacing()903 inline double FreeFormTransformation::GetTSpacing() const
904 {
905 return _dt;
906 }
907
908 // ---------------------------------------------------------------------------
GetSpacing(double & dx,double & dy,double & dz)909 inline void FreeFormTransformation::GetSpacing(double &dx, double &dy, double &dz) const
910 {
911 dx = _dx;
912 dy = _dy;
913 dz = _dz;
914 }
915
916 // ---------------------------------------------------------------------------
GetSpacing(double & dx,double & dy,double & dz,double & dt)917 inline void FreeFormTransformation::GetSpacing(double &dx, double &dy, double &dz, double &dt) const
918 {
919 dx = _dx;
920 dy = _dy;
921 dz = _dz;
922 dt = _dt;
923 }
924
925 // ---------------------------------------------------------------------------
PutOrientation(double * xaxis,double * yaxis,double * zaxis)926 inline void FreeFormTransformation::PutOrientation(double *xaxis, double *yaxis, double *zaxis)
927 {
928 _CPImage.PutOrientation(xaxis, yaxis, zaxis);
929 }
930
931 // ---------------------------------------------------------------------------
GetOrientation(double * xaxis,double * yaxis,double * zaxis)932 inline void FreeFormTransformation::GetOrientation(double *xaxis, double *yaxis, double *zaxis) const
933 {
934 _CPImage.GetOrientation(xaxis, yaxis, zaxis);
935 }
936
937 // -----------------------------------------------------------------------------
IndexToDOFs(int cp,int & x,int & y)938 inline void FreeFormTransformation::IndexToDOFs(int cp, int &x, int &y) const
939 {
940 x = 3 * cp;
941 y = x + 1;
942 }
943
944 // -----------------------------------------------------------------------------
IndexToDOFs(int cp,int & x,int & y,int & z)945 inline void FreeFormTransformation::IndexToDOFs(int cp, int &x, int &y, int &z) const
946 {
947 x = 3 * cp;
948 y = x + 1;
949 z = x + 2;
950 }
951
952 // -----------------------------------------------------------------------------
DOFToIndex(int dof)953 inline int FreeFormTransformation::DOFToIndex(int dof) const
954 {
955 return dof / 3;
956 }
957
958 // -----------------------------------------------------------------------------
DOFToDimension(int dof)959 inline int FreeFormTransformation::DOFToDimension(int dof) const
960 {
961 return dof % 3;
962 }
963
964 // -----------------------------------------------------------------------------
LatticeToIndex(int i,int j,int k,int l)965 inline int FreeFormTransformation::LatticeToIndex(int i, int j, int k, int l) const
966 {
967 return _CPImage.VoxelToIndex(i, j, k, l);
968 }
969
970 // -----------------------------------------------------------------------------
IndexToLattice(int index,int & i,int & j)971 inline void FreeFormTransformation::IndexToLattice(int index, int &i, int &j) const
972 {
973 _CPImage.IndexToVoxel(index, i, j);
974 }
975
976 // -----------------------------------------------------------------------------
IndexToLattice(int index,int & i,int & j,int & k)977 inline void FreeFormTransformation::IndexToLattice(int index, int &i, int &j, int &k) const
978 {
979 _CPImage.IndexToVoxel(index, i, j, k);
980 }
981
982 // -----------------------------------------------------------------------------
IndexToLattice(int index,int & i,int & j,int & k,int & l)983 inline void FreeFormTransformation::IndexToLattice(int index, int &i, int &j, int &k, int &l) const
984 {
985 _CPImage.IndexToVoxel(index, i, j, k, l);
986 }
987
988 // -----------------------------------------------------------------------------
IndexToWorld(int index,double & x,double & y)989 inline void FreeFormTransformation::IndexToWorld(int index, double &x, double &y) const
990 {
991 _CPImage.IndexToWorld(index, x, y);
992 }
993
994 // -----------------------------------------------------------------------------
IndexToWorld(int index,double & x,double & y,double & z)995 inline void FreeFormTransformation::IndexToWorld(int index, double &x, double &y, double &z) const
996 {
997 _CPImage.IndexToWorld(index, x, y, z);
998 }
999
1000 // -----------------------------------------------------------------------------
IndexToWorld(int index,Point & p)1001 inline void FreeFormTransformation::IndexToWorld(int index, Point &p) const
1002 {
1003 _CPImage.IndexToWorld(index, p);
1004 }
1005
1006 // -----------------------------------------------------------------------------
IndexToWorld(int index)1007 inline Point FreeFormTransformation::IndexToWorld(int index) const
1008 {
1009 return _CPImage.IndexToWorld(index);
1010 }
1011
1012 // -----------------------------------------------------------------------------
WorldToLattice(double & x,double & y)1013 inline void FreeFormTransformation::WorldToLattice(double &x, double &y) const
1014 {
1015 _CPImage.WorldToImage(x, y);
1016 }
1017
1018 // -----------------------------------------------------------------------------
WorldToLattice(double & x,double & y,double & z)1019 inline void FreeFormTransformation::WorldToLattice(double &x, double &y, double &z) const
1020 {
1021 _CPImage.WorldToImage(x, y, z);
1022 }
1023
1024 // -----------------------------------------------------------------------------
WorldToLattice(Point & p)1025 inline void FreeFormTransformation::WorldToLattice(Point &p) const
1026 {
1027 this->WorldToLattice(p._x, p._y, p._z);
1028 }
1029
1030 // -----------------------------------------------------------------------------
LatticeToWorld(double & x,double & y)1031 inline void FreeFormTransformation::LatticeToWorld(double &x, double &y) const
1032 {
1033 _CPImage.ImageToWorld(x, y);
1034 }
1035
1036 // -----------------------------------------------------------------------------
LatticeToWorld(double & x,double & y,double & z)1037 inline void FreeFormTransformation::LatticeToWorld(double &x, double &y, double &z) const
1038 {
1039 _CPImage.ImageToWorld(x, y, z);
1040 }
1041
1042 // -----------------------------------------------------------------------------
LatticeToWorld(Point & p)1043 inline void FreeFormTransformation::LatticeToWorld(Point &p) const
1044 {
1045 this->LatticeToWorld(p._x, p._y, p._z);
1046 }
1047
1048 // -----------------------------------------------------------------------------
ControlPointLocation(int cp)1049 inline Point FreeFormTransformation::ControlPointLocation(int cp) const
1050 {
1051 int i, j, k;
1052 IndexToLattice(cp, i, j, k);
1053 Point p(i, j, k);
1054 this->LatticeToWorld(p);
1055 return p;
1056 }
1057
1058 // -----------------------------------------------------------------------------
ControlPointLocation(int cp,double & x,double & y)1059 inline void FreeFormTransformation::ControlPointLocation(int cp, double &x, double &y) const
1060 {
1061 int i, j;
1062 IndexToLattice(cp, i, j);
1063 x = i, y = j;
1064 this->LatticeToWorld(x, y);
1065 }
1066
1067 // -----------------------------------------------------------------------------
ControlPointLocation(int cp,double & x,double & y,double & z)1068 inline void FreeFormTransformation::ControlPointLocation(int cp, double &x, double &y, double &z) const
1069 {
1070 int i, j, k;
1071 IndexToLattice(cp, i, j, k);
1072 x = i, y = j, z = k;
1073 this->LatticeToWorld(x, y, z);
1074 }
1075
1076 // -----------------------------------------------------------------------------
TimeToLattice(double t)1077 inline double FreeFormTransformation::TimeToLattice(double t) const
1078 {
1079 return _CPImage.TimeToImage(t);
1080 }
1081
1082 // -----------------------------------------------------------------------------
LatticeToTime(double t)1083 inline double FreeFormTransformation::LatticeToTime(double t) const
1084 {
1085 return _CPImage.ImageToTime(t);
1086 }
1087
1088 // -----------------------------------------------------------------------------
GetXAfterSubdivision()1089 inline int FreeFormTransformation::GetXAfterSubdivision() const
1090 {
1091 return this->GetX();
1092 }
1093
1094 // -----------------------------------------------------------------------------
GetYAfterSubdivision()1095 inline int FreeFormTransformation::GetYAfterSubdivision() const
1096 {
1097 return this->GetY();
1098 }
1099
1100 // -----------------------------------------------------------------------------
GetZAfterSubdivision()1101 inline int FreeFormTransformation::GetZAfterSubdivision() const
1102 {
1103 return this->GetZ();
1104 }
1105
1106 // -----------------------------------------------------------------------------
GetTAfterSubdivision()1107 inline int FreeFormTransformation::GetTAfterSubdivision() const
1108 {
1109 return this->GetT();
1110 }
1111
1112 // -----------------------------------------------------------------------------
GetXSpacingAfterSubdivision()1113 inline double FreeFormTransformation::GetXSpacingAfterSubdivision() const
1114 {
1115 return this->GetXSpacing();
1116 }
1117
1118 // -----------------------------------------------------------------------------
GetYSpacingAfterSubdivision()1119 inline double FreeFormTransformation::GetYSpacingAfterSubdivision() const
1120 {
1121 return this->GetYSpacing();
1122 }
1123
1124 // -----------------------------------------------------------------------------
GetZSpacingAfterSubdivision()1125 inline double FreeFormTransformation::GetZSpacingAfterSubdivision() const
1126 {
1127 return this->GetZSpacing();
1128 }
1129
1130 // -----------------------------------------------------------------------------
GetTSpacingAfterSubdivision()1131 inline double FreeFormTransformation::GetTSpacingAfterSubdivision() const
1132 {
1133 return this->GetTSpacing();
1134 }
1135
1136 // -----------------------------------------------------------------------------
Subdivide(bool subdivide_x,bool subdivide_y,bool subdivide_z,bool subdivide_t)1137 inline void FreeFormTransformation::Subdivide(bool subdivide_x, bool subdivide_y, bool subdivide_z, bool subdivide_t)
1138 {
1139 if (subdivide_x || subdivide_y || subdivide_z || subdivide_t) {
1140 cerr << this->NameOfClass() << "::Subdivide: Not implemented (or not possible?)" << endl;
1141 exit(1);
1142 }
1143 }
1144
1145 // -----------------------------------------------------------------------------
Subdivide2D()1146 inline void FreeFormTransformation::Subdivide2D()
1147 {
1148 this->Subdivide(true, true, false, false);
1149 }
1150
1151 // -----------------------------------------------------------------------------
Subdivide3D()1152 inline void FreeFormTransformation::Subdivide3D()
1153 {
1154 this->Subdivide(true, true, true, false);
1155 }
1156
1157 // -----------------------------------------------------------------------------
Subdivide4D()1158 inline void FreeFormTransformation::Subdivide4D()
1159 {
1160 this->Subdivide(true, true, true, true);
1161 }
1162
1163 // =============================================================================
1164 // Bounding box
1165 // =============================================================================
1166
1167 // -----------------------------------------------------------------------------
BoundingBox(double & t1,double & t2)1168 inline void FreeFormTransformation::BoundingBox(double &t1, double &t2) const
1169 {
1170 t1 = this->LatticeToTime(.0);
1171 t2 = this->LatticeToTime(_t - 1);
1172 if (t1 > t2) swap(t1, t2);
1173 }
1174
1175 // -----------------------------------------------------------------------------
BoundingBox(double & x1,double & y1,double & z1,double & x2,double & y2,double & z2)1176 inline void FreeFormTransformation::BoundingBox(double &x1, double &y1, double &z1,
1177 double &x2, double &y2, double &z2) const
1178 {
1179 x1 = y1 = z1 = .0;
1180 this->LatticeToWorld(x1, y1, z1);
1181 x2 = _x - 1, y2 = _y - 1, z2 = _z - 1;
1182 this->LatticeToWorld(x2, y2, z2);
1183 if (x1 > x2) swap(x1, x2);
1184 if (y1 > y2) swap(y1, y2);
1185 if (z1 > z2) swap(z1, z2);
1186 }
1187
1188 // -----------------------------------------------------------------------------
BoundingBox(Point & p1,Point & p2)1189 inline void FreeFormTransformation::BoundingBox(Point &p1, Point &p2) const
1190 {
1191 BoundingBox(p1._x, p1._y, p1._z, p2._x, p2._y, p2._z);
1192 }
1193
1194 // -----------------------------------------------------------------------------
BoundingBox(double & x1,double & y1,double & z1,double & t1,double & x2,double & y2,double & z2,double & t2)1195 inline void FreeFormTransformation::BoundingBox(double &x1, double &y1, double &z1, double &t1,
1196 double &x2, double &y2, double &z2, double &t2) const
1197 {
1198 BoundingBox(x1, y1, z1, x2, y2, z2);
1199 BoundingBox(t1, t2);
1200 }
1201
1202 // -----------------------------------------------------------------------------
BoundingBox(Point & p1,double & t1,Point & p2,double & t2)1203 inline void FreeFormTransformation::BoundingBox(Point &p1, double &t1, Point &p2, double &t2) const
1204 {
1205 BoundingBox(p1._x, p1._y, p1._z, t1, p2._x, p2._y, p2._z, t2);
1206 }
1207
1208 // -----------------------------------------------------------------------------
BoundingBox(int,double & t1,double & t2,double)1209 inline void FreeFormTransformation::BoundingBox(int, double &t1, double &t2, double) const
1210 {
1211 t1 = this->LatticeToTime(0);
1212 t2 = this->LatticeToTime(_t - 1);
1213 if (t1 > t2) swap(t1, t2);
1214 }
1215
1216 // -----------------------------------------------------------------------------
BoundingBox(int,double & x1,double & y1,double & z1,double & x2,double & y2,double & z2,double)1217 inline void FreeFormTransformation::BoundingBox(int, double &x1, double &y1, double &z1,
1218 double &x2, double &y2, double &z2, double) const
1219 {
1220 x1 = 0, y1 = 0, z1 = 0;
1221 x2 = _x - 1, y2 = _y - 1, z2 = _z - 1;
1222 this->LatticeToWorld(x1, y1, z1);
1223 this->LatticeToWorld(x2, y2, z2);
1224 if (x1 > x2) swap(x1, x2);
1225 if (y1 > y2) swap(y1, y2);
1226 if (z1 > z2) swap(z1, z2);
1227 }
1228
1229 // -----------------------------------------------------------------------------
BoundingBox(int cp,double & x1,double & y1,double & z1,double & t1,double & x2,double & y2,double & z2,double & t2,double fraction)1230 inline void FreeFormTransformation::BoundingBox(int cp, double &x1, double &y1, double &z1, double &t1,
1231 double &x2, double &y2, double &z2, double &t2,
1232 double fraction) const
1233 {
1234 this->BoundingBox(cp, x1, y1, z1, x2, y2, z2, fraction);
1235 this->BoundingBox(cp, t1, t2, fraction);
1236 }
1237
1238 // -----------------------------------------------------------------------------
BoundingBox(int cp,Point & p1,Point & p2,double fraction)1239 inline void FreeFormTransformation::BoundingBox(int cp, Point &p1, Point &p2, double fraction) const
1240 {
1241 BoundingBox(cp, p1._x, p1._y, p1._z, p2._x, p2._y, p2._z, fraction);
1242 }
1243
1244 // -----------------------------------------------------------------------------
BoundingBox(const ImageAttributes & domain,int cp,int & i1,int & j1,int & k1,int & i2,int & j2,int & k2,double fraction)1245 inline bool FreeFormTransformation::BoundingBox(const ImageAttributes &domain, int cp,
1246 int &i1, int &j1, int &k1,
1247 int &i2, int &j2, int &k2,
1248 double fraction) const
1249 {
1250 // Calculate bounding box in world coordinates parallel to world axes
1251 double x[2], y[2], z[2];
1252 this->BoundingBox(cp, x[0], y[0], z[0], x[1], y[1], z[1], fraction);
1253
1254 // Map bounding box into image space and calculate minimal axes-alinged
1255 // bounding box parallel to image axes which need not coincide with world axes
1256 Point p;
1257 double x1, y1, z1, x2, y2, z2;
1258 x1 = y1 = z1 = + inf;
1259 x2 = y2 = z2 = - inf;
1260
1261 for (int c = 0; c <= 1; ++c)
1262 for (int b = 0; b <= 1; ++b)
1263 for (int a = 0; a <= 1; ++a) {
1264 p = Point(x[a], y[b], z[c]);
1265 domain.WorldToLattice(p);
1266 if (p._x < x1) x1 = p._x;
1267 if (p._x > x2) x2 = p._x;
1268 if (p._y < y1) y1 = p._y;
1269 if (p._y > y2) y2 = p._y;
1270 if (p._z < z1) z1 = p._z;
1271 if (p._z > z2) z2 = p._z;
1272 }
1273
1274 // Round to nearest voxel in image domain
1275 i1 = iround(x1);
1276 i2 = iround(x2);
1277 j1 = iround(y1);
1278 j2 = iround(y2);
1279 k1 = iround(z1);
1280 k2 = iround(z2);
1281
1282 // When both indices are outside in opposite directions,
1283 // use the full range [0, N[. If they are both outside in
1284 // the same direction, the condition i1 <= i2 is false which
1285 // indicates that the bounding box is empty in this case
1286 i1 = (i1 < 0 ? 0 : (i1 >= domain.X() ? domain.X() : i1));
1287 i2 = (i2 < 0 ? -1 : (i2 >= domain.X() ? domain.X() - 1 : i2));
1288 j1 = (j1 < 0 ? 0 : (j1 >= domain.Y() ? domain.Y() : j1));
1289 j2 = (j2 < 0 ? -1 : (j2 >= domain.Y() ? domain.Y() - 1 : j2));
1290 k1 = (k1 < 0 ? 0 : (k1 >= domain.Z() ? domain.Z() : k1));
1291 k2 = (k2 < 0 ? -1 : (k2 >= domain.Z() ? domain.Z() - 1 : k2));
1292 return i1 <= i2 && j1 <= j2 && k1 <= k2;
1293 }
1294
1295 // -----------------------------------------------------------------------------
BoundingBox(const Image * image,int cp,int & i1,int & j1,int & k1,int & i2,int & j2,int & k2,double fraction)1296 inline bool FreeFormTransformation::BoundingBox(const Image *image, int cp,
1297 int &i1, int &j1, int &k1,
1298 int &i2, int &j2, int &k2,
1299 double fraction) const
1300 {
1301 return BoundingBox(image->Attributes(), cp, i1, j1, k1, i2, j2, k2, fraction);
1302 }
1303
1304 // -----------------------------------------------------------------------------
BoundingBox(const ImageAttributes & domain,int cp,int & i1,int & j1,int & k1,int & l1,int & i2,int & j2,int & k2,int & l2,double fraction)1305 inline bool FreeFormTransformation::BoundingBox(const ImageAttributes &domain, int cp,
1306 int &i1, int &j1, int &k1, int &l1,
1307 int &i2, int &j2, int &k2, int &l2,
1308 double fraction) const
1309 {
1310 // Calculate spatial bounding box in image coordinates
1311 bool bbvalid = BoundingBox(domain, cp, i1, j1, k1, i2, j2, k2, fraction);
1312
1313 // Calculate temporal bounding box
1314 double t1, t2;
1315 this->BoundingBox(cp, t1, t2, fraction);
1316
1317 // Convert to image coordinates
1318 t1 = domain.TimeToLattice(t1);
1319 t2 = domain.TimeToLattice(t2);
1320 if (t2 < t1) swap(t1, t2);
1321
1322 // Round to nearest voxel in image domain
1323 l1 = iround(t1);
1324 l2 = iround(t2);
1325
1326 // When both indices are outside in opposite directions,
1327 // use the full range [0, N[. If they are both outside in
1328 // the same direction, the condition l1 <= l2 is false which
1329 // indicates that the bounding box is empty in this case
1330 l1 = (l1 < 0 ? 0 : (l1 >= domain.T() ? domain.T() : l1));
1331 l2 = (l2 < 0 ? -1 : (l2 >= domain.T() ? domain.T() - 1 : l2));
1332 return bbvalid && l1 <= l2;
1333 }
1334
1335 // -----------------------------------------------------------------------------
BoundingBox(const Image * image,int cp,int & i1,int & j1,int & k1,int & l1,int & i2,int & j2,int & k2,int & l2,double fraction)1336 inline bool FreeFormTransformation::BoundingBox(const Image *image, int cp,
1337 int &i1, int &j1, int &k1, int &l1,
1338 int &i2, int &j2, int &k2, int &l2,
1339 double fraction) const
1340 {
1341 return BoundingBox(image->Attributes(), cp, i1, j1, k1, l1, i2, j2, k2, l2, fraction);
1342 }
1343
1344 // -----------------------------------------------------------------------------
DOFBoundingBox(const ImageAttributes & domain,int dof,int & i1,int & j1,int & k1,int & i2,int & j2,int & k2,double fraction)1345 inline bool FreeFormTransformation::DOFBoundingBox(const ImageAttributes &domain, int dof,
1346 int &i1, int &j1, int &k1,
1347 int &i2, int &j2, int &k2,
1348 double fraction) const
1349 {
1350 return BoundingBox(domain, this->DOFToIndex(dof), i1, j1, k1, i2, j2, k2, fraction);
1351 }
1352
1353 // -----------------------------------------------------------------------------
DOFBoundingBox(const Image * image,int dof,int & i1,int & j1,int & k1,int & i2,int & j2,int & k2,double fraction)1354 inline bool FreeFormTransformation::DOFBoundingBox(const Image *image, int dof,
1355 int &i1, int &j1, int &k1,
1356 int &i2, int &j2, int &k2,
1357 double fraction) const
1358 {
1359 return BoundingBox(image->Attributes(), dof, i1, j1, k1, i2, j2, k2, fraction);
1360 }
1361
1362 // =============================================================================
1363 // Transformation parameters (DoFs)
1364 // =============================================================================
1365
1366 // -----------------------------------------------------------------------------
DOFGradientNorm(const double * gradient)1367 inline double FreeFormTransformation::DOFGradientNorm(const double *gradient) const
1368 {
1369 double norm, max = .0;
1370 int x, y, z;
1371 const int ncps = this->NumberOfCPs();
1372 for (int cp = 0; cp < ncps; ++cp) {
1373 this->IndexToDOFs(cp, x, y, z);
1374 norm = sqrt(gradient[x] * gradient[x] + gradient[y] * gradient[y] + gradient[z] * gradient[z]);
1375 if (norm > max) max = norm;
1376 }
1377 return max;
1378 }
1379
1380 // -----------------------------------------------------------------------------
CopyFrom(const Transformation * other)1381 inline bool FreeFormTransformation::CopyFrom(const Transformation *other)
1382 {
1383 const FreeFormTransformation *ffd;
1384 ffd = dynamic_cast<const FreeFormTransformation *>(other);
1385 if (ffd && typeid(*ffd) == typeid(*this) && ffd->Attributes() == this->Attributes()) {
1386 return Transformation::CopyFrom(other);
1387 } else {
1388 return false;
1389 }
1390 }
1391
1392 // -----------------------------------------------------------------------------
Put(int cp,const Vector & x)1393 inline void FreeFormTransformation::Put(int cp, const Vector &x)
1394 {
1395 _CPImage(cp) = x;
1396 }
1397
1398 // -----------------------------------------------------------------------------
Put(int cp,double x,double y,double z)1399 inline void FreeFormTransformation::Put(int cp, double x, double y, double z)
1400 {
1401 _CPImage(cp) = Vector(x, y, z);
1402 }
1403
1404 // -----------------------------------------------------------------------------
Put(int i,int j,int k,int l,double x,double y,double z)1405 inline void FreeFormTransformation::Put(int i, int j, int k, int l,
1406 double x, double y, double z)
1407 {
1408 _CPImage(i, j, k, l) = Vector(x, y, z);
1409 }
1410
1411 // -----------------------------------------------------------------------------
Put(int i,int j,int k,double x,double y,double z)1412 inline void FreeFormTransformation::Put(int i, int j, int k,
1413 double x, double y, double z)
1414 {
1415 Put(i, j, k, 0, x, y, z);
1416 }
1417
1418 // -----------------------------------------------------------------------------
Get(int cp,Vector & x)1419 inline void FreeFormTransformation::Get(int cp, Vector &x) const
1420 {
1421 x = _CPImage(cp);
1422 }
1423
1424 // -----------------------------------------------------------------------------
Get(int cp,double & x,double & y,double & z)1425 inline void FreeFormTransformation::Get(int cp, double &x, double &y, double &z) const
1426 {
1427 const Vector ¶m = _CPImage(cp);
1428 x = param._x, y = param._y, z = param._z;
1429 }
1430
1431 // -----------------------------------------------------------------------------
Get(int i,int j,int k,int l,double & x,double & y,double & z)1432 inline void FreeFormTransformation::Get(int i, int j, int k, int l,
1433 double &x, double &y, double &z) const
1434 {
1435 const Vector ¶m = _CPImage(i, j, k, l);
1436 x = param._x, y = param._y, z = param._z;
1437 }
1438
1439 // -----------------------------------------------------------------------------
Get(int i,int j,int k,double & x,double & y,double & z)1440 inline void FreeFormTransformation::Get(int i, int j, int k,
1441 double &x, double &y, double &z) const
1442 {
1443 Get(i, j, k, 0, x, y, z);
1444 }
1445
1446 // -----------------------------------------------------------------------------
PutStatus(int cp,const CPStatus & status)1447 inline void FreeFormTransformation::PutStatus(int cp, const CPStatus &status)
1448 {
1449 _CPStatus[0][0][0][cp] = status;
1450 }
1451
1452 // -----------------------------------------------------------------------------
PutStatus(int i,int j,int k,int l,DOFStatus sx,DOFStatus sy,DOFStatus sz)1453 inline void FreeFormTransformation::PutStatus(int i, int j, int k, int l,
1454 DOFStatus sx, DOFStatus sy, DOFStatus sz)
1455 {
1456 CPStatus &status = _CPStatus[l][k][j][i];
1457 status._x = sx;
1458 status._y = sy;
1459 status._z = sz;
1460 }
1461
1462 // -----------------------------------------------------------------------------
PutStatus(int i,int j,int k,DOFStatus sx,DOFStatus sy,DOFStatus sz)1463 inline void FreeFormTransformation::PutStatus(int i, int j, int k,
1464 DOFStatus sx, DOFStatus sy, DOFStatus sz)
1465 {
1466 PutStatus(i, j, k, 0, sx, sy, sz);
1467 }
1468
1469 // -----------------------------------------------------------------------------
GetStatus(int cp,CPStatus & status)1470 inline void FreeFormTransformation::GetStatus(int cp, CPStatus &status) const
1471 {
1472 const CPStatus &s = _CPStatus[0][0][0][cp];
1473 status._x = s._x;
1474 status._y = s._y;
1475 status._z = s._z;
1476 }
1477
1478 // -----------------------------------------------------------------------------
GetStatus(int i,int j,int k,int l,DOFStatus & sx,DOFStatus & sy,DOFStatus & sz)1479 inline void FreeFormTransformation::GetStatus(int i, int j, int k, int l,
1480 DOFStatus &sx, DOFStatus &sy, DOFStatus &sz) const
1481 {
1482 const CPStatus &status = _CPStatus[l][k][j][i];
1483 sx = status._x;
1484 sy = status._y;
1485 sz = status._z;
1486 }
1487
1488 // -----------------------------------------------------------------------------
GetStatus(int i,int j,int k,DOFStatus & sx,DOFStatus & sy,DOFStatus & sz)1489 inline void FreeFormTransformation::GetStatus(int i, int j, int k,
1490 DOFStatus &sx, DOFStatus &sy, DOFStatus &sz) const
1491 {
1492 GetStatus(i, j, k, 0, sx, sy, sz);
1493 }
1494
1495 // -----------------------------------------------------------------------------
IsActive(int cp)1496 inline bool FreeFormTransformation::IsActive(int cp) const
1497 {
1498 const CPStatus &status = _CPStatus[0][0][0][cp];
1499 return status._x == Active || status._y == Active || status._z == Active;
1500 }
1501
1502 // -----------------------------------------------------------------------------
IsActive(int i,int j,int k,int l)1503 inline bool FreeFormTransformation::IsActive(int i, int j, int k, int l) const
1504 {
1505 const CPStatus &status = _CPStatus[l][k][j][i];
1506 return status._x == Active || status._y == Active || status._z == Active;
1507 }
1508
1509 // =============================================================================
1510 // Point transformation
1511 // =============================================================================
1512
1513 // -----------------------------------------------------------------------------
GlobalTransform(double &,double &,double &,double,double)1514 inline void FreeFormTransformation::GlobalTransform(double &, double &, double &, double, double) const
1515 {
1516 // T_global(x) = x
1517 }
1518
1519 // -----------------------------------------------------------------------------
GlobalInverse(double &,double &,double &,double,double)1520 inline void FreeFormTransformation::GlobalInverse(double &, double &, double &, double, double) const
1521 {
1522 // T_global(x) = x
1523 }
1524
1525 // -----------------------------------------------------------------------------
Transform(double & x,double & y,double & z,double t,double t0)1526 inline void FreeFormTransformation::Transform(double &x, double &y, double &z, double t, double t0) const
1527 {
1528 this->LocalTransform(x, y, z, t, t0);
1529 }
1530
1531 // -----------------------------------------------------------------------------
Inverse(double & x,double & y,double & z,double t,double t0)1532 inline bool FreeFormTransformation::Inverse(double &x, double &y, double &z, double t, double t0) const
1533 {
1534 return this->LocalInverse(x, y, z, t, t0);
1535 }
1536
1537 // =============================================================================
1538 // Derivatives
1539 // =============================================================================
1540
1541 // -----------------------------------------------------------------------------
JacobianToWorld(double & du,double & dv)1542 inline void FreeFormTransformation::JacobianToWorld(double &du, double &dv) const
1543 {
1544 double dx = du * _matW2L(0, 0) + dv * _matW2L(1, 0);
1545 double dy = du * _matW2L(0, 1) + dv * _matW2L(1, 1);
1546 du = dx, dv = dy;
1547 }
1548
1549 // -----------------------------------------------------------------------------
JacobianToWorld(double & du,double & dv,double & dw)1550 inline void FreeFormTransformation::JacobianToWorld(double &du, double &dv, double &dw) const
1551 {
1552 double dx = du * _matW2L(0, 0) + dv * _matW2L(1, 0) + dw * _matW2L(2, 0);
1553 double dy = du * _matW2L(0, 1) + dv * _matW2L(1, 1) + dw * _matW2L(2, 1);
1554 double dz = du * _matW2L(0, 2) + dv * _matW2L(1, 2) + dw * _matW2L(2, 2);
1555 du = dx, dv = dy, dw = dz;
1556 }
1557
1558 // -----------------------------------------------------------------------------
JacobianToWorld(Matrix & jac)1559 inline void FreeFormTransformation::JacobianToWorld(Matrix &jac) const
1560 {
1561 if (jac.Rows() == 2) {
1562 JacobianToWorld(jac(0, 0), jac(0, 1));
1563 JacobianToWorld(jac(1, 0), jac(1, 1));
1564 } else {
1565 JacobianToWorld(jac(0, 0), jac(0, 1), jac(0, 2));
1566 JacobianToWorld(jac(1, 0), jac(1, 1), jac(1, 2));
1567 JacobianToWorld(jac(2, 0), jac(2, 1), jac(2, 2));
1568 }
1569 }
1570
1571 // -----------------------------------------------------------------------------
JacobianToWorldOrientation(double & du,double & dv)1572 inline void FreeFormTransformation::JacobianToWorldOrientation(double &du, double &dv) const
1573 {
1574 double dx = du * _attr._xaxis[0] + dv * _attr._yaxis[0];
1575 double dy = du * _attr._xaxis[1] + dv * _attr._yaxis[1];
1576 du = dx, dv = dy;
1577 }
1578
1579 // -----------------------------------------------------------------------------
JacobianToWorldOrientation(double & du,double & dv,double & dw)1580 inline void FreeFormTransformation::JacobianToWorldOrientation(double &du, double &dv, double &dw) const
1581 {
1582 double dx = du * _attr._xaxis[0] + dv * _attr._yaxis[0] + dw * _attr._zaxis[0];
1583 double dy = du * _attr._xaxis[1] + dv * _attr._yaxis[1] + dw * _attr._zaxis[1];
1584 double dz = du * _attr._xaxis[2] + dv * _attr._yaxis[2] + dw * _attr._zaxis[2];
1585 du = dx, dv = dy, dw = dz;
1586 }
1587
1588 // -----------------------------------------------------------------------------
1589 inline void
HessianToWorld(double & duu,double & duv,double & dvv)1590 FreeFormTransformation::HessianToWorld(double &duu, double &duv, double &dvv) const
1591 {
1592 // The derivatives of the world to lattice coordinate transformation
1593 // w.r.t the world coordinates which are needed for the chain rule below
1594 const double &dudx = _matW2L(0, 0);
1595 const double &dudy = _matW2L(0, 1);
1596 const double &dvdx = _matW2L(1, 0);
1597 const double &dvdy = _matW2L(1, 1);
1598 // Expression computed here is transpose(R) * Hessian * R = transpose(Hessian * R) * R
1599 // where R is the 2x2 world to lattice reorientation and scaling matrix
1600 double du, dv, dxx, dxy, dyy;
1601 du = duu * dudx + duv * dvdx;
1602 dv = duv * dudx + dvv * dvdx;
1603 dxx = du * dudx + dv * dvdx;
1604 dxy = du * dudy + dv * dvdy;
1605 du = duu * dudy + duv * dvdy;
1606 dv = duv * dudy + dvv * dvdy;
1607 dyy = du * dudy + dv * dvdy;
1608 // Return computed derivatives
1609 duu = dxx, duv = dxy, dvv = dyy;
1610 }
1611
1612 // -----------------------------------------------------------------------------
1613 inline void
HessianToWorld(double & duu,double & duv,double & duw,double & dvv,double & dvw,double & dww)1614 FreeFormTransformation::HessianToWorld(double &duu, double &duv, double &duw,
1615 double &dvv, double &dvw,
1616 double &dww) const
1617 {
1618 // The derivatives of the world to lattice coordinate transformation
1619 // w.r.t the world coordinates which are needed for the chain rule below
1620 const double &dudx = _matW2L(0, 0);
1621 const double &dudy = _matW2L(0, 1);
1622 const double &dudz = _matW2L(0, 2);
1623 const double &dvdx = _matW2L(1, 0);
1624 const double &dvdy = _matW2L(1, 1);
1625 const double &dvdz = _matW2L(1, 2);
1626 const double &dwdx = _matW2L(2, 0);
1627 const double &dwdy = _matW2L(2, 1);
1628 const double &dwdz = _matW2L(2, 2);
1629 // Expression computed here is transpose(R) * Hessian * R = transpose(Hessian * R) * R
1630 // where R is the 3x3 world to lattice reorientation and scaling matrix
1631 double du, dv, dw, dxx, dxy, dxz, dyy, dyz, dzz;
1632 du = duu * dudx + duv * dvdx + duw * dwdx;
1633 dv = duv * dudx + dvv * dvdx + dvw * dwdx;
1634 dw = duw * dudx + dvw * dvdx + dww * dwdx;
1635 dxx = du * dudx + dv * dvdx + dw * dwdx;
1636 dxy = du * dudy + dv * dvdy + dw * dwdy;
1637 dxz = du * dudz + dv * dvdz + dw * dwdz;
1638 du = duu * dudy + duv * dvdy + duw * dwdy;
1639 dv = duv * dudy + dvv * dvdy + dvw * dwdy;
1640 dw = duw * dudy + dvw * dvdy + dww * dwdy;
1641 dyy = du * dudy + dv * dvdy + dw * dwdy;
1642 dyz = du * dudz + dv * dvdz + dw * dwdz;
1643 du = duu * dudz + duv * dvdz + duw * dwdz;
1644 dv = duv * dudz + dvv * dvdz + dvw * dwdz;
1645 dw = duw * dudz + dvw * dvdz + dww * dwdz;
1646 dzz = du * dudz + dv * dvdz + dw * dwdz;
1647 // Return computed derivatives
1648 duu = dxx, duv = dxy, duw = dxz, dvv = dyy, dvw = dyz, dww = dzz;
1649 }
1650
1651 // -----------------------------------------------------------------------------
HessianToWorld(Matrix & hessian)1652 inline void FreeFormTransformation::HessianToWorld(Matrix &hessian) const
1653 {
1654 if (hessian.Rows() == 2) {
1655 HessianToWorld(hessian(0, 0), hessian(0, 1), hessian(1, 1));
1656 hessian(1, 0) = hessian(0, 1);
1657 } else {
1658 HessianToWorld(hessian(0, 0), hessian(0, 1), hessian(0, 2),
1659 hessian(1, 1), hessian(1, 2),
1660 hessian(2, 2));
1661 hessian(1, 0) = hessian(0, 1);
1662 hessian(2, 0) = hessian(0, 2);
1663 hessian(2, 1) = hessian(1, 2);
1664 }
1665 }
1666
1667 // -----------------------------------------------------------------------------
HessianToWorld(Matrix hessian[3])1668 inline void FreeFormTransformation::HessianToWorld(Matrix hessian[3]) const
1669 {
1670 HessianToWorld(hessian[0]);
1671 HessianToWorld(hessian[1]);
1672 HessianToWorld(hessian[2]);
1673 }
1674
1675 // -----------------------------------------------------------------------------
FFDJacobianWorld(Matrix &,double,double,double,double,double)1676 inline void FreeFormTransformation::FFDJacobianWorld(Matrix &, double, double, double, double, double) const
1677 {
1678 Throw(ERR_NotImplemented, __FUNCTION__, "Not implemented");
1679 }
1680
1681 // -----------------------------------------------------------------------------
GlobalJacobian(Matrix & jac,double,double,double,double,double)1682 inline void FreeFormTransformation::GlobalJacobian(Matrix &jac, double, double, double, double, double) const
1683 {
1684 // T_global(x) = x
1685 jac.Initialize(3, 3);
1686 jac(0, 0) = 1;
1687 jac(1, 1) = 1;
1688 jac(2, 2) = 1;
1689 }
1690
1691 // -----------------------------------------------------------------------------
Jacobian(Matrix & jac,double x,double y,double z,double t,double t0)1692 inline void FreeFormTransformation::Jacobian(Matrix &jac, double x, double y, double z, double t, double t0) const
1693 {
1694 this->LocalJacobian(jac, x, y, z, t, t0);
1695 }
1696
1697 // -----------------------------------------------------------------------------
GlobalHessian(Matrix hessian[3],double,double,double,double,double)1698 inline void FreeFormTransformation::GlobalHessian(Matrix hessian[3], double, double, double, double, double) const
1699 {
1700 // T_global(x) = x
1701 hessian[0].Initialize(3, 3);
1702 hessian[1].Initialize(3, 3);
1703 hessian[2].Initialize(3, 3);
1704 }
1705
1706 // -----------------------------------------------------------------------------
Hessian(Matrix hessian[3],double x,double y,double z,double t,double t0)1707 inline void FreeFormTransformation::Hessian(Matrix hessian[3], double x, double y, double z, double t, double t0) const
1708 {
1709 this->LocalHessian(hessian, x, y, z, t, t0);
1710 }
1711
1712 // -----------------------------------------------------------------------------
JacobianDOFs(Matrix & jac,int cp,double x,double y,double z,double t,double t0)1713 inline void FreeFormTransformation::JacobianDOFs(Matrix &jac, int cp, double x, double y, double z, double t, double t0) const
1714 {
1715 int xdof, ydof, zdof;
1716 double dofgrad[3];
1717
1718 // initialize 3x3 Jacobian matrix
1719 jac.Initialize(3, 3);
1720 // compute indices of control point DoFs
1721 this->IndexToDOFs(cp, xdof, ydof, zdof);
1722 // derivatives w.r.t. x component of control point
1723 if (this->GetStatus(xdof) == Active) {
1724 this->JacobianDOFs(dofgrad, xdof, x, y, z, t, t0);
1725 jac(0, 0) = dofgrad[0];
1726 jac(1, 0) = dofgrad[1];
1727 jac(2, 0) = dofgrad[2];
1728 }
1729 // derivatives w.r.t. y component of control point
1730 if (this->GetStatus(ydof) == Active) {
1731 this->JacobianDOFs(dofgrad, ydof, x, y, z, t, t0);
1732 jac(0, 1) = dofgrad[0];
1733 jac(1, 1) = dofgrad[1];
1734 jac(2, 1) = dofgrad[2];
1735 }
1736 // derivatives w.r.t. z component of control point
1737 if (this->GetStatus(zdof) == Active) {
1738 this->JacobianDOFs(dofgrad, zdof, x, y, z, t, t0);
1739 jac(0, 2) = dofgrad[0];
1740 jac(1, 2) = dofgrad[1];
1741 jac(2, 2) = dofgrad[2];
1742 }
1743 }
1744
1745 // -----------------------------------------------------------------------------
JacobianDOFs(TransformationJacobian & jac,double x,double y,double z,double t,double t0)1746 inline void FreeFormTransformation::JacobianDOFs(TransformationJacobian &jac, double x, double y, double z, double t, double t0) const
1747 {
1748 int xdof, ydof, zdof;
1749 double dofgrad[3];
1750
1751 for (int cp = 0; cp < NumberOfCPs(); ++cp) {
1752 // compute indices of control point DoFs
1753 this->IndexToDOFs(cp, xdof, ydof, zdof);
1754 // derivatives w.r.t. x component of control point
1755 if (this->GetStatus(xdof) == Active) {
1756 this->JacobianDOFs(dofgrad, xdof, x, y, z, t, t0);
1757 if (dofgrad[0] != .0 && dofgrad[1] != .0 && dofgrad[2] != .0) {
1758 TransformationJacobian::ColumnType &xdofgrad = jac(xdof);
1759 xdofgrad._x = dofgrad[0];
1760 xdofgrad._y = dofgrad[1];
1761 xdofgrad._z = dofgrad[2];
1762 }
1763 }
1764 // derivatives w.r.t. y component of control point
1765 if (this->GetStatus(ydof) == Active) {
1766 this->JacobianDOFs(dofgrad, ydof, x, y, z, t, t0);
1767 if (dofgrad[0] != .0 && dofgrad[1] != .0 && dofgrad[2] != .0) {
1768 TransformationJacobian::ColumnType &ydofgrad = jac(ydof);
1769 ydofgrad._x = dofgrad[0];
1770 ydofgrad._y = dofgrad[1];
1771 ydofgrad._z = dofgrad[2];
1772 }
1773 }
1774 // derivatives w.r.t. z component of control point
1775 if (this->GetStatus(zdof) == Active) {
1776 this->JacobianDOFs(dofgrad, zdof, x, y, z, t, t0);
1777 if (dofgrad[0] != .0 && dofgrad[1] != .0 && dofgrad[2] != .0) {
1778 TransformationJacobian::ColumnType &zdofgrad = jac(zdof);
1779 zdofgrad._x = dofgrad[0];
1780 zdofgrad._y = dofgrad[1];
1781 zdofgrad._z = dofgrad[2];
1782 }
1783 }
1784 }
1785 }
1786
1787 // -----------------------------------------------------------------------------
1788 inline void FreeFormTransformation
JacobianDetDerivative(double dJ[3],const Matrix & adj,int cp,double x,double y,double z,double t,double t0,bool wrt_world,bool use_spacing)1789 ::JacobianDetDerivative(double dJ[3], const Matrix &adj,
1790 int cp, double x, double y, double z, double t, double t0,
1791 bool wrt_world, bool use_spacing) const
1792 {
1793 this->FFDJacobianDetDerivative(dJ, adj, cp, x, y, z, t, t0, wrt_world, use_spacing);
1794 }
1795
1796 // =============================================================================
1797 // Properties
1798 // =============================================================================
1799
1800 // -----------------------------------------------------------------------------
Bending3D(const Matrix hessian[3])1801 inline double FreeFormTransformation::Bending3D(const Matrix hessian[3])
1802 {
1803 const Matrix &hx = hessian[0];
1804 const Matrix &hy = hessian[1];
1805 const Matrix &hz = hessian[2];
1806
1807 const double &x_ii = hx(0, 0);
1808 const double &x_ij = hx(0, 1);
1809 const double &x_ik = hx(0, 2);
1810 const double &x_jj = hx(1, 1);
1811 const double &x_jk = hx(1, 2);
1812 const double &x_kk = hx(2, 2);
1813
1814 const double &y_ii = hy(0, 0);
1815 const double &y_ij = hy(0, 1);
1816 const double &y_ik = hy(0, 2);
1817 const double &y_jj = hy(1, 1);
1818 const double &y_jk = hy(1, 2);
1819 const double &y_kk = hy(2, 2);
1820
1821 const double &z_ii = hz(0, 0);
1822 const double &z_ij = hz(0, 1);
1823 const double &z_ik = hz(0, 2);
1824 const double &z_jj = hz(1, 1);
1825 const double &z_jk = hz(1, 2);
1826 const double &z_kk = hz(2, 2);
1827
1828 return ( x_ii * x_ii + x_jj * x_jj + x_kk * x_kk
1829 + y_ii * y_ii + y_jj * y_jj + y_kk * y_kk
1830 + z_ii * z_ii + z_jj * z_jj + z_kk * z_kk)
1831 + 2.0 * ( x_ij * x_ij + x_ik * x_ik + x_jk * x_jk
1832 + y_ij * y_ij + y_ik * y_ik + y_jk * y_jk
1833 + z_ij * z_ij + z_ik * z_ik + z_jk * z_jk);
1834 }
1835
1836
1837 } // namespace mirtk
1838
1839 #endif // MIRTK_FreeFormTransformation_H
1840