go home Home | Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Data Structures | File List | Namespace Members | Data Fields | Globals | Related Pages
itkRecursiveBSplineTransformImplementation.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright UMC Utrecht and contributors
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkRecursiveBSplineTransformImplementation_h
19 #define __itkRecursiveBSplineTransformImplementation_h
20 
22 
23 namespace itk
24 {
25 
40 template< unsigned int OutputDimension, unsigned int SpaceDimension, unsigned int SplineOrder, class TScalar >
42 {
43 public:
44 
48  typedef TScalar ScalarType;
49  typedef double InternalFloatType;
50 
52  itkStaticConstMacro( HelperConstVariable, unsigned int,
53  ( SpaceDimension - 1 ) * ( SplineOrder + 1 ) );
54 
57  ScalarType, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType;
58  itkStaticConstMacro( BSplineNumberOfIndices, unsigned int,
59  RecursiveBSplineWeightFunctionType::NumberOfIndices );
60 
63 
65  static inline void TransformPoint(
67  const OffsetValueType * gridOffsetTable,
68  const double * weights1D )
69  {
71  ScalarType * tmp_mu[ OutputDimension ];
72  for( unsigned int j = 0; j < OutputDimension; ++j )
73  {
74  tmp_mu[ j ] = mu[ j ];
75  }
76 
78  ScalarType tmp_opp[ OutputDimension ];
79  for( unsigned int j = 0; j < OutputDimension; ++j )
80  {
81  opp[ j ] = 0.0;
82  }
83 
84  OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
85  for( unsigned int k = 0; k <= SplineOrder; ++k )
86  {
89  ::TransformPoint( tmp_opp, tmp_mu, gridOffsetTable, weights1D );
90 
92  for( unsigned int j = 0; j < OutputDimension; ++j )
93  {
94  opp[ j ] += tmp_opp[ j ] * weights1D[ k + HelperConstVariable ];
95 
96  // move to the next mu
97  tmp_mu[ j ] += bot;
98  }
99  }
100  } // end TransformPoint()
101 
102 
104  static inline void GetJacobian(
105  ScalarType * & jacobians, const double * weights1D, double value )
106  {
107  for( unsigned int k = 0; k <= SplineOrder; ++k )
108  {
111  ::GetJacobian( jacobians, weights1D, value * weights1D[ k + HelperConstVariable ] );
112  }
113  } // end GetJacobian()
114 
115 
118  ScalarType * & imageJacobian, const InternalFloatType * movingImageGradient,
119  const double * weights1D, double value )
120  {
121  for( unsigned int k = 0; k <= SplineOrder; ++k )
122  {
125  ::EvaluateJacobianWithImageGradientProduct( imageJacobian, movingImageGradient, weights1D,
126  value * weights1D[ k + HelperConstVariable ] );
127  }
128  } // end EvaluateJacobianWithImageGradientProduct()
129 
130 
132  static inline void ComputeNonZeroJacobianIndices(
133  unsigned long * & nzji,
134  const unsigned long parametersPerDim,
135  unsigned long currentIndex,
136  const OffsetValueType * gridOffsetTable )
137  {
138  const OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
139  for( unsigned int k = 0; k <= SplineOrder; ++k )
140  {
143  ::ComputeNonZeroJacobianIndices( nzji, parametersPerDim, currentIndex, gridOffsetTable );
144 
145  currentIndex += bot;
146  }
147  } // end ComputeNonZeroJacobianIndices()
148 
149 
154  static inline void GetSpatialJacobian(
155  InternalFloatType * sj,
157  const OffsetValueType * gridOffsetTable,
158  const double * weights1D, // normal B-spline weights
159  const double * derivativeWeights1D ) // 1st derivative of B-spline
160  {
162  ScalarType * tmp_mu[ OutputDimension ];
163  for( unsigned int j = 0; j < OutputDimension; ++j )
164  {
165  tmp_mu[ j ] = mu[ j ];
166  }
167 
169  InternalFloatType tmp_sj[ OutputDimension * SpaceDimension ];
170  for( unsigned int n = 0; n < OutputDimension * ( SpaceDimension + 1 ); ++n )
171  {
172  sj[ n ] = 0.0;
173  }
174 
175  OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
176  for( unsigned int k = 0; k <= SplineOrder; ++k )
177  {
180  ::GetSpatialJacobian( tmp_sj, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D );
181 
183  for( unsigned int n = 0; n < OutputDimension * SpaceDimension; ++n )
184  {
185  sj[ n ] += tmp_sj[ n ] * weights1D[ k + HelperConstVariable ];
186  }
187 
189  for( unsigned int j = 0; j < OutputDimension; ++j )
190  {
191  sj[ OutputDimension * SpaceDimension + j ]
192  += tmp_sj[ j ] * derivativeWeights1D[ k + HelperConstVariable ];
193 
194  // move to the next mu
195  tmp_mu[ j ] += bot;
196  }
197  }
198  } // end GetSpatialJacobian()
199 
200 
220  static inline void GetSpatialHessian(
221  InternalFloatType * sh,
223  const OffsetValueType * gridOffsetTable,
224  const double * weights1D, // normal B-spline weights
225  const double * derivativeWeights1D, // 1st derivative of B-spline
226  const double * hessianWeights1D ) // 2nd derivative of B-spline
227  {
228  const unsigned int helperDim1 = OutputDimension * SpaceDimension * ( SpaceDimension + 1 ) / 2;
229  const unsigned int helperDim2 = OutputDimension * ( SpaceDimension + 1 ) * ( SpaceDimension + 2 ) / 2;
230 
232  ScalarType * tmp_mu[ OutputDimension ];
233  for( unsigned int j = 0; j < OutputDimension; ++j )
234  {
235  tmp_mu[ j ] = mu[ j ];
236  }
237 
239  InternalFloatType tmp_sh[ helperDim1 ];
240  for( unsigned int n = 0; n < helperDim2; ++n )
241  {
242  sh[ n ] = 0.0;
243  }
244 
245  OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
246  for( unsigned int k = 0; k <= SplineOrder; ++k )
247  {
250  ::GetSpatialHessian( tmp_sh, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D, hessianWeights1D );
251 
253  for( unsigned int n = 0; n < helperDim1; ++n )
254  {
255  sh[ n ] += tmp_sh[ n ] * weights1D[ k + HelperConstVariable ];
256  }
257 
259  for( unsigned int n = 0; n < SpaceDimension; ++n )
260  {
261  for( unsigned int j = 0; j < OutputDimension; ++j )
262  {
263  sh[ OutputDimension * n + helperDim1 + j ]
264  += tmp_sh[ OutputDimension * n * ( n + 1 ) / 2 + j ] * derivativeWeights1D[ k + HelperConstVariable ];
265  }
266  }
267 
269  for( unsigned int j = 0; j < OutputDimension; ++j )
270  {
271  sh[ helperDim2 - OutputDimension + j ]
272  += tmp_sh[ j ] * hessianWeights1D[ k + HelperConstVariable ];
273 
274  // move to the next mu
275  tmp_mu[ j ] += bot;
276  }
277  }
278  } // end GetSpatialHessian()
279 
280 
284  static inline void GetJacobianOfSpatialJacobian(
285  InternalFloatType * & jsj_out,
286  const double * weights1D, // normal B-spline weights
287  const double * derivativeWeights1D, // 1st derivative of B-spline
288  const double * directionCosines,
289  InternalFloatType * jsj )
290  {
291  const unsigned int helperDim = OutputDimension - SpaceDimension + 1;
292 
294  InternalFloatType tmp_jsj[ helperDim + 1 ];
295 
296  for( unsigned int k = 0; k <= SplineOrder; ++k )
297  {
298  const double w = weights1D[ k + HelperConstVariable ];
299  const double dw = derivativeWeights1D[ k + HelperConstVariable ];
300 
302  for( unsigned int n = 0; n < helperDim; ++n )
303  {
304  tmp_jsj[ n ] = jsj[ n ] * w;
305  }
306 
308  tmp_jsj[ helperDim ] = jsj[ 0 ] * dw;
309 
312  ::GetJacobianOfSpatialJacobian( jsj_out, weights1D, derivativeWeights1D, directionCosines, tmp_jsj );
313  }
314  } // end GetJacobianOfSpatialJacobian()
315 
316 
320  static inline void GetJacobianOfSpatialHessian(
321  InternalFloatType * & jsh_out,
322  const double * weights1D, // normal B-spline weights
323  const double * derivativeWeights1D, // 1st derivative of B-spline
324  const double * hessianWeights1D, // 2nd derivative of B-spline
325  const double * directionCosines,
326  InternalFloatType * jsh )
327  {
328  const unsigned int helperDim = OutputDimension - SpaceDimension;
329  const unsigned int helperDimW = ( helperDim + 1 ) * ( helperDim + 2 ) / 2;
330  const unsigned int helperDimDW = helperDim + 1;
331 
333  InternalFloatType tmp_jsh[ helperDimW + helperDimDW + 1 ];
334 
335  for( unsigned int k = 0; k <= SplineOrder; ++k )
336  {
338  const double w = weights1D[ k + HelperConstVariable ];
339  const double dw = derivativeWeights1D[ k + HelperConstVariable ];
340  const double hw = hessianWeights1D[ k + HelperConstVariable ];
341 
343  for( unsigned int n = 0; n < helperDimW; ++n )
344  {
345  tmp_jsh[ n ] = jsh[ n ] * w;
346  }
347 
349  for( unsigned int n = 0; n < helperDimDW; ++n )
350  {
351  unsigned int nn = n * ( n + 1 ) / 2;
352  tmp_jsh[ n + helperDimW ] = jsh[ nn ] * dw;
353  }
354 
356  tmp_jsh[ helperDimW + helperDimDW ] = jsh[ 0 ] * hw;
357 
360  ::GetJacobianOfSpatialHessian( jsh_out, weights1D, derivativeWeights1D, hessianWeights1D, directionCosines, tmp_jsh );
361  }
362  } // end GetJacobianOfSpatialHessian()
363 
364 
365 };
366 
367 
373 template< unsigned int OutputDimension, unsigned int SplineOrder, class TScalar >
374 class RecursiveBSplineTransformImplementation< OutputDimension, 0, SplineOrder, TScalar >
375 {
376 public:
377 
381  typedef TScalar ScalarType;
382  typedef double InternalFloatType;
383 
386  TScalar, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType;
387  itkStaticConstMacro( BSplineNumberOfIndices, unsigned int,
388  RecursiveBSplineWeightFunctionType::NumberOfIndices );
389 
392 
394  static inline void TransformPoint(
396  const OffsetValueType * gridOffsetTable,
397  const double * weights1D )
398  {
399  for( unsigned int j = 0; j < OutputDimension; ++j )
400  {
401  opp[ j ] = *( mu[ j ] );
402  }
403  } // end TransformPoint()
404 
405 
407  static inline void GetJacobian(
408  ScalarType * & jacobians, const double * weights1D, double value )
409  {
410  unsigned long offset = 0;
411  for( unsigned int j = 0; j < OutputDimension; ++j )
412  {
413  offset = j * BSplineNumberOfIndices * ( OutputDimension + 1 );
414  jacobians[ offset ] = value;
415  }
416  ++jacobians;
417  } // end GetJacobian()
418 
419 
422  ScalarType * & imageJacobian, const InternalFloatType * movingImageGradient,
423  const double * weights1D, double value )
424  {
425  for( unsigned int j = 0; j < OutputDimension; ++j )
426  {
427  *( imageJacobian + j * BSplineNumberOfIndices ) = value * movingImageGradient[ j ];
428  }
429  ++imageJacobian;
430  } // end EvaluateJacobianWithImageGradientProduct()
431 
432 
434  static inline void ComputeNonZeroJacobianIndices(
435  unsigned long * & nzji,
436  const unsigned long parametersPerDim,
437  unsigned long currentIndex,
438  const OffsetValueType * gridOffsetTable )
439  {
440  for( unsigned int j = 0; j < OutputDimension; ++j )
441  {
442  nzji[ j * BSplineNumberOfIndices ] = currentIndex + j * parametersPerDim;
443  }
444  ++nzji;
445  } // end ComputeNonZeroJacobianIndices()
446 
447 
449  static inline void GetSpatialJacobian(
450  InternalFloatType * sj,
452  const OffsetValueType * gridOffsetTable,
453  const double * weights1D, // normal B-spline weights
454  const double * derivativeWeights1D ) // 1st derivative of B-spline
455  {
456  for( unsigned int j = 0; j < OutputDimension; ++j )
457  {
458  sj[ j ] = *( mu[ j ] );
459  }
460  } // end GetSpatialJacobian()
461 
462 
464  static inline void GetSpatialHessian(
465  InternalFloatType * sh,
467  const OffsetValueType * gridOffsetTable,
468  const double * weights1D, // normal B-spline weights
469  const double * derivativeWeights1D, // 1st derivative of B-spline
470  const double * hessianWeights1D ) // 2nd derivative of B-spline
471  {
472  for( unsigned int j = 0; j < OutputDimension; ++j )
473  {
474  sh[ j ] = *( mu[ j ] );
475  }
476  } // end GetSpatialHessian()
477 
478 
480  static inline void GetJacobianOfSpatialJacobian(
481  InternalFloatType * & jsj_out,
482  const double * weights1D, // normal B-spline weights
483  const double * derivativeWeights1D, // 1st derivative of B-spline
484  const double * directionCosines,
485  InternalFloatType * jsj )
486  {
492  for( unsigned int j = 0; j < OutputDimension; ++j )
493  {
494  jsj_out[ j ] = jsj[ OutputDimension ] * directionCosines[ j ];
495  for( unsigned int k = 1; k < OutputDimension; ++k )
496  {
497  jsj_out[ k ] += jsj[ OutputDimension - k ] * directionCosines[ k * OutputDimension + j ];
498  }
499  }
500 
502  unsigned int offset = 0;
503  for( unsigned int i = 0; i < OutputDimension; ++i )
504  {
505  offset = i * ( OutputDimension * ( BSplineNumberOfIndices * OutputDimension + 1 ) );
506  for( unsigned int j = 0; j < OutputDimension; ++j )
507  {
508  jsj_out[ j + offset ] = jsj_out[ j ];
509  }
510  }
511 
513  jsj_out += OutputDimension * OutputDimension;
514 
515  } // end GetJacobianOfSpatialJacobian()
516 
517 
519  static inline void GetJacobianOfSpatialHessian(
520  InternalFloatType * & jsh_out,
521  const double * weights1D, // normal B-spline weights
522  const double * derivativeWeights1D, // 1st derivative of B-spline
523  const double * hessianWeights1D, // 2nd derivative of B-spline
524  const double * directionCosines,
525  InternalFloatType * jsh )
526  {
527  double jsh_tmp[ OutputDimension * OutputDimension ];
528  double matrixProduct[ OutputDimension * OutputDimension ];
529 
536  if( OutputDimension == 3 )
537  {
538  jsh_tmp[ 0 ] = jsh[ 9 ]; jsh_tmp[ 1 ] = jsh[ 8 ]; jsh_tmp[ 2 ] = jsh[ 7 ];
539  jsh_tmp[ 3 ] = jsh_tmp[ 1 ]; jsh_tmp[ 4 ] = jsh[ 5 ]; jsh_tmp[ 5 ] = jsh[ 4 ];
540  jsh_tmp[ 6 ] = jsh_tmp[ 2 ]; jsh_tmp[ 7 ] = jsh_tmp[ 5 ]; jsh_tmp[ 8 ] = jsh[ 2 ];
541  }
542  else if( OutputDimension == 2 )
543  {
544  jsh_tmp[ 0 ] = jsh[ 5 ]; jsh_tmp[ 1 ] = jsh[ 4 ];
545  jsh_tmp[ 2 ] = jsh_tmp[ 1 ]; jsh_tmp[ 3 ] = jsh[ 2 ];
546  }
547  else // the general case
548  {
549  for( unsigned int j = 0; j < OutputDimension; ++j )
550  {
551  for( unsigned int i = 0; i <= j; ++i )
552  {
553  jsh_tmp[ j * OutputDimension + i ] = jsh[ ( OutputDimension - j ) + ( OutputDimension - i ) * ( OutputDimension - i + 1 ) / 2 ];
554  if( i != j )
555  {
556  jsh_tmp[ i * OutputDimension + j ] = jsh_tmp[ j * OutputDimension + i ];
557  }
558  }
559  }
560  }
561 
563  for( unsigned int i = 0; i < OutputDimension; ++i ) // row
564  {
565  for( unsigned int j = 0; j < OutputDimension; ++j ) // column
566  {
567  double accum = directionCosines[ i ] * jsh_tmp[ j ];
568  for( unsigned int k = 1; k < OutputDimension; ++k )
569  {
570  accum += directionCosines[ k * OutputDimension + i ] * jsh_tmp[ k * OutputDimension + j ];
571  }
572  matrixProduct[ i * OutputDimension + j ] = accum;
573  }
574  }
575 
577  for( unsigned int i = 0; i < OutputDimension; ++i ) // row
578  {
579  for( unsigned int j = 0; j < OutputDimension; ++j ) // column
580  {
581  double accum = matrixProduct[ i * OutputDimension ] * directionCosines[ j ];
582  for( unsigned int k = 1; k < OutputDimension; ++k )
583  {
584  accum += matrixProduct[ i * OutputDimension + k ] * directionCosines[ k * OutputDimension + j ];
585  }
586  jsh_out[ i * OutputDimension + j ] = accum;
587  }
588  }
589 
591  unsigned long offset = 0;
592  for( unsigned int i = 0; i < OutputDimension; ++i )
593  {
594  offset = i * ( OutputDimension * OutputDimension * ( BSplineNumberOfIndices * OutputDimension + 1 ) );
595  for( unsigned int j = 0; j < OutputDimension * OutputDimension; ++j )
596  {
597  jsh_out[ j + offset ] = jsh_out[ j ];
598  }
599  }
600 
602  jsh_out += OutputDimension * OutputDimension * OutputDimension;
603 
604  } // end GetJacobianOfSpatialHessian()
605 
606 
607 };
608 
609 
610 } // end namespace itk
611 
612 #endif /* __itkRecursiveBSplineTransformImplementation_h */
Returns the weights over the support region used for B-spline interpolation/reconstruction.
itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices)
static void GetJacobianOfSpatialHessian(InternalFloatType *&jsh_out, const double *weights1D, const double *derivativeWeights1D, const double *hessianWeights1D, const double *directionCosines, InternalFloatType *jsh)
static void GetSpatialHessian(InternalFloatType *sh, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D, const double *derivativeWeights1D, const double *hessianWeights1D)
static void EvaluateJacobianWithImageGradientProduct(ScalarType *&imageJacobian, const InternalFloatType *movingImageGradient, const double *weights1D, double value)
static void TransformPoint(OutputPointType opp, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D)
itk::RecursiveBSplineInterpolationWeightFunction< TScalar, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType
static void GetSpatialJacobian(InternalFloatType *sj, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D, const double *derivativeWeights1D)
static void ComputeNonZeroJacobianIndices(unsigned long *&nzji, const unsigned long parametersPerDim, unsigned long currentIndex, const OffsetValueType *gridOffsetTable)
static void GetJacobianOfSpatialJacobian(InternalFloatType *&jsj_out, const double *weights1D, const double *derivativeWeights1D, const double *directionCosines, InternalFloatType *jsj)
This helper class contains the actual implementation of the recursive B-spline transform.
static void GetSpatialHessian(InternalFloatType *sh, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D, const double *derivativeWeights1D, const double *hessianWeights1D)
itkStaticConstMacro(HelperConstVariable, unsigned int,(SpaceDimension - 1) *(SplineOrder+1))
static void EvaluateJacobianWithImageGradientProduct(ScalarType *&imageJacobian, const InternalFloatType *movingImageGradient, const double *weights1D, double value)
static void ComputeNonZeroJacobianIndices(unsigned long *&nzji, const unsigned long parametersPerDim, unsigned long currentIndex, const OffsetValueType *gridOffsetTable)
static void GetJacobianOfSpatialHessian(InternalFloatType *&jsh_out, const double *weights1D, const double *derivativeWeights1D, const double *hessianWeights1D, const double *directionCosines, InternalFloatType *jsh)
static void TransformPoint(OutputPointType opp, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D)
static void GetJacobian(ScalarType *&jacobians, const double *weights1D, double value)
static void GetSpatialJacobian(InternalFloatType *sj, const CoefficientPointerVectorType mu, const OffsetValueType *gridOffsetTable, const double *weights1D, const double *derivativeWeights1D)
itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices)
static void GetJacobianOfSpatialJacobian(InternalFloatType *&jsj_out, const double *weights1D, const double *derivativeWeights1D, const double *directionCosines, InternalFloatType *jsj)
itk::RecursiveBSplineInterpolationWeightFunction< ScalarType, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType


Generated on 1667476801 for elastix by doxygen 1.9.1 elastix logo