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// Standard C++ header files:
24#include <algorithm> // For copy_n and fill_n.
25#include <cassert>
26#include <cstring> // For memcpy.
27
28namespace itk
29{
30
45template <unsigned int OutputDimension, unsigned int SpaceDimension, unsigned int SplineOrder, class TScalar>
47{
48public:
50
52 itkStaticConstMacro(HelperConstVariable, unsigned int, (SpaceDimension - 1) * (SplineOrder + 1));
53
57 itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices);
58
60 static void
61 TransformPoint(TScalar * const opp,
62 const TScalar * const * const mu,
63 const OffsetValueType * const gridOffsetTable,
64 const double * const weights1D)
65 {
67 const TScalar * tmp_mu[OutputDimension];
68 std::copy_n(mu, OutputDimension, tmp_mu);
69
71 TScalar tmp_opp[OutputDimension];
72 std::fill_n(opp, OutputDimension, 0.0);
73
74 const OffsetValueType bot = gridOffsetTable[SpaceDimension - 1];
75 for (unsigned int k = 0; k <= SplineOrder; ++k)
76 {
79 TransformPoint(tmp_opp, tmp_mu, gridOffsetTable, weights1D);
80
82 for (unsigned int j = 0; j < OutputDimension; ++j)
83 {
84 opp[j] += tmp_opp[j] * weights1D[k + HelperConstVariable];
85
86 // move to the next mu
87 tmp_mu[j] += bot;
88 }
89 }
90 } // end TransformPoint()
91
92
94 static void
95 GetJacobian(TScalar *& jacobians, const double * const weights1D, const double value)
96 {
97 for (unsigned int k = 0; k <= SplineOrder; ++k)
98 {
101 jacobians, weights1D, value * weights1D[k + HelperConstVariable]);
102 }
103 } // end GetJacobian()
104
105
107 static void
109 const InternalFloatType * const movingImageGradient,
110 const double * const weights1D,
111 const double value)
112 {
113 for (unsigned int k = 0; k <= SplineOrder; ++k)
114 {
118 imageJacobian, movingImageGradient, weights1D, value * weights1D[k + HelperConstVariable]);
119 }
120 } // end EvaluateJacobianWithImageGradientProduct()
121
122
124 static void
125 ComputeNonZeroJacobianIndices(unsigned long *& nzji,
126 const unsigned long parametersPerDim,
127 unsigned long currentIndex,
128 const OffsetValueType * const gridOffsetTable)
129 {
130 const OffsetValueType bot = gridOffsetTable[SpaceDimension - 1];
131 for (unsigned int k = 0; k <= SplineOrder; ++k)
132 {
135 ComputeNonZeroJacobianIndices(nzji, parametersPerDim, currentIndex, gridOffsetTable);
136
137 currentIndex += bot;
138 }
139 } // end ComputeNonZeroJacobianIndices()
140
141
146 static void
148 const TScalar * const * const mu,
149 const OffsetValueType * const gridOffsetTable,
150 const double * const weights1D, // normal B-spline weights
151 const double * const derivativeWeights1D) // 1st derivative of B-spline
152 {
154 const TScalar * tmp_mu[OutputDimension];
155 std::copy_n(mu, OutputDimension, tmp_mu);
156
158 InternalFloatType tmp_sj[OutputDimension * SpaceDimension];
159 for (unsigned int n = 0; n < OutputDimension * (SpaceDimension + 1); ++n)
160 {
161 sj[n] = 0.0;
162 }
163
164 const OffsetValueType bot = gridOffsetTable[SpaceDimension - 1];
165 for (unsigned int k = 0; k <= SplineOrder; ++k)
166 {
169 GetSpatialJacobian(tmp_sj, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D);
170
172 for (unsigned int n = 0; n < OutputDimension * SpaceDimension; ++n)
173 {
174 sj[n] += tmp_sj[n] * weights1D[k + HelperConstVariable];
175 }
176
178 for (unsigned int j = 0; j < OutputDimension; ++j)
179 {
180 sj[OutputDimension * SpaceDimension + j] += tmp_sj[j] * derivativeWeights1D[k + HelperConstVariable];
181
182 // move to the next mu
183 tmp_mu[j] += bot;
184 }
185 }
186 } // end GetSpatialJacobian()
187
188
208 static void
210 const TScalar * const * const mu,
211 const OffsetValueType * const gridOffsetTable,
212 const double * const weights1D, // normal B-spline weights
213 const double * const derivativeWeights1D, // 1st derivative of B-spline
214 const double * const hessianWeights1D) // 2nd derivative of B-spline
215 {
216 const unsigned int helperDim1 = OutputDimension * SpaceDimension * (SpaceDimension + 1) / 2;
217 const unsigned int helperDim2 = OutputDimension * (SpaceDimension + 1) * (SpaceDimension + 2) / 2;
218
220 const TScalar * tmp_mu[OutputDimension];
221 std::copy_n(mu, OutputDimension, tmp_mu);
222
224 InternalFloatType tmp_sh[helperDim1];
225 for (unsigned int n = 0; n < helperDim2; ++n)
226 {
227 sh[n] = 0.0;
228 }
229
230 const OffsetValueType bot = gridOffsetTable[SpaceDimension - 1];
231 for (unsigned int k = 0; k <= SplineOrder; ++k)
232 {
235 GetSpatialHessian(tmp_sh, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D, hessianWeights1D);
236
238 for (unsigned int n = 0; n < helperDim1; ++n)
239 {
240 sh[n] += tmp_sh[n] * weights1D[k + HelperConstVariable];
241 }
242
244 for (unsigned int n = 0; n < SpaceDimension; ++n)
245 {
246 for (unsigned int j = 0; j < OutputDimension; ++j)
247 {
248 sh[OutputDimension * n + helperDim1 + j] +=
249 tmp_sh[OutputDimension * n * (n + 1) / 2 + j] * derivativeWeights1D[k + HelperConstVariable];
250 }
251 }
252
254 for (unsigned int j = 0; j < OutputDimension; ++j)
255 {
256 sh[helperDim2 - OutputDimension + j] += tmp_sh[j] * hessianWeights1D[k + HelperConstVariable];
257
258 // move to the next mu
259 tmp_mu[j] += bot;
260 }
261 }
262 } // end GetSpatialHessian()
263
264
268 static void
270 const double * const weights1D, // normal B-spline weights
271 const double * const derivativeWeights1D, // 1st derivative of B-spline
272 const double * const directionCosines,
273 const InternalFloatType * const jsj)
274 {
275 const unsigned int helperDim = OutputDimension - SpaceDimension + 1;
276
278 InternalFloatType tmp_jsj[helperDim + 1];
279
280 for (unsigned int k = 0; k <= SplineOrder; ++k)
281 {
282 const double w = weights1D[k + HelperConstVariable];
283 const double dw = derivativeWeights1D[k + HelperConstVariable];
284
286 for (unsigned int n = 0; n < helperDim; ++n)
287 {
288 tmp_jsj[n] = jsj[n] * w;
289 }
290
292 tmp_jsj[helperDim] = jsj[0] * dw;
293
296 GetJacobianOfSpatialJacobian(jsj_out, weights1D, derivativeWeights1D, directionCosines, tmp_jsj);
297 }
298 } // end GetJacobianOfSpatialJacobian()
299
300
304 static void
306 const double * const weights1D, // normal B-spline weights
307 const double * const derivativeWeights1D, // 1st derivative of B-spline
308 const double * const hessianWeights1D, // 2nd derivative of B-spline
309 const double * const directionCosines,
310 const InternalFloatType * const jsh)
311 {
312 const unsigned int helperDim = OutputDimension - SpaceDimension;
313 const unsigned int helperDimW = (helperDim + 1) * (helperDim + 2) / 2;
314 const unsigned int helperDimDW = helperDim + 1;
315
317 InternalFloatType tmp_jsh[helperDimW + helperDimDW + 1];
318
319 for (unsigned int k = 0; k <= SplineOrder; ++k)
320 {
322 const double w = weights1D[k + HelperConstVariable];
323 const double dw = derivativeWeights1D[k + HelperConstVariable];
324 const double hw = hessianWeights1D[k + HelperConstVariable];
325
327 for (unsigned int n = 0; n < helperDimW; ++n)
328 {
329 tmp_jsh[n] = jsh[n] * w;
330 }
331
333 for (unsigned int n = 0; n < helperDimDW; ++n)
334 {
335 unsigned int nn = n * (n + 1) / 2;
336 tmp_jsh[n + helperDimW] = jsh[nn] * dw;
337 }
338
340 tmp_jsh[helperDimW + helperDimDW] = jsh[0] * hw;
341
345 jsh_out, weights1D, derivativeWeights1D, hessianWeights1D, directionCosines, tmp_jsh);
346 }
347 } // end GetJacobianOfSpatialHessian()
348};
349
350
356template <unsigned int OutputDimension, unsigned int SplineOrder, class TScalar>
357class ITK_TEMPLATE_EXPORT RecursiveBSplineTransformImplementation<OutputDimension, 0, SplineOrder, TScalar>
358{
359public:
361
365 itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices);
366
368 static void
369 TransformPoint(TScalar * const opp,
370 const TScalar * const * const mu,
371 const OffsetValueType * const gridOffsetTable,
372 const double * const weights1D)
373 {
374 for (unsigned int j = 0; j < OutputDimension; ++j)
375 {
376 opp[j] = *(mu[j]);
377 }
378 } // end TransformPoint()
379
380
382 static void
383 GetJacobian(TScalar *& jacobians, const double * const weights1D, const double value)
384 {
385 unsigned long offset = 0;
386 for (unsigned int j = 0; j < OutputDimension; ++j)
387 {
388 offset = j * BSplineNumberOfIndices * (OutputDimension + 1);
389 jacobians[offset] = value;
390 }
391 ++jacobians;
392 } // end GetJacobian()
393
394
396 static void
398 const InternalFloatType * const movingImageGradient,
399 const double * const weights1D,
400 const double value)
401 {
402 for (unsigned int j = 0; j < OutputDimension; ++j)
403 {
404 *(imageJacobian + j * BSplineNumberOfIndices) = value * movingImageGradient[j];
405 }
406 ++imageJacobian;
407 } // end EvaluateJacobianWithImageGradientProduct()
408
409
411 static void
412 ComputeNonZeroJacobianIndices(unsigned long *& nzji,
413 const unsigned long parametersPerDim,
414 const unsigned long currentIndex,
415 const OffsetValueType * const gridOffsetTable)
416 {
417 for (unsigned int j = 0; j < OutputDimension; ++j)
418 {
419 nzji[j * BSplineNumberOfIndices] = currentIndex + j * parametersPerDim;
420 }
421 ++nzji;
422 } // end ComputeNonZeroJacobianIndices()
423
424
426 static void
428 const TScalar * const * const mu,
429 const OffsetValueType * const gridOffsetTable,
430 const double * const weights1D, // normal B-spline weights
431 const double * const derivativeWeights1D) // 1st derivative of B-spline
432 {
433 for (unsigned int j = 0; j < OutputDimension; ++j)
434 {
435 sj[j] = *(mu[j]);
436 }
437 } // end GetSpatialJacobian()
438
439
441 static void
443 const TScalar * const * const mu,
444 const OffsetValueType * const gridOffsetTable,
445 const double * const weights1D, // normal B-spline weights
446 const double * const derivativeWeights1D, // 1st derivative of B-spline
447 const double * const hessianWeights1D) // 2nd derivative of B-spline
448 {
449 for (unsigned int j = 0; j < OutputDimension; ++j)
450 {
451 sh[j] = *(mu[j]);
452 }
453 } // end GetSpatialHessian()
454
455
457 static void
459 const double * const weights1D, // normal B-spline weights
460 const double * const derivativeWeights1D, // 1st derivative of B-spline
461 const double * const directionCosines,
462 const InternalFloatType * const jsj)
463 {
469 for (unsigned int j = 0; j < OutputDimension; ++j)
470 {
471 jsj_out[j] = jsj[OutputDimension] * directionCosines[j];
472 for (unsigned int k = 1; k < OutputDimension; ++k)
473 {
474 jsj_out[k] += jsj[OutputDimension - k] * directionCosines[k * OutputDimension + j];
475 }
476 }
477
479 unsigned int offset = 0;
480 for (unsigned int i = 0; i < OutputDimension; ++i)
481 {
482 offset = i * (OutputDimension * (BSplineNumberOfIndices * OutputDimension + 1));
483 for (unsigned int j = 0; j < OutputDimension; ++j)
484 {
485 jsj_out[j + offset] = jsj_out[j];
486 }
487 }
488
490 jsj_out += OutputDimension * OutputDimension;
491
492 } // end GetJacobianOfSpatialJacobian()
493
494
496 static void
498 const double * const weights1D, // normal B-spline weights
499 const double * const derivativeWeights1D, // 1st derivative of B-spline
500 const double * const hessianWeights1D, // 2nd derivative of B-spline
501 const double * const directionCosines,
502 const InternalFloatType * const jsh)
503 {
504 double jsh_tmp[OutputDimension * OutputDimension];
505 double matrixProduct[OutputDimension * OutputDimension];
506
513 if (OutputDimension == 3)
514 {
515 const double tmp[] = { jsh[9], jsh[8], jsh[7], jsh[8], jsh[5], jsh[4], jsh[7], jsh[4], jsh[2] };
516 FastBitwiseCopy(jsh_tmp, tmp);
517 }
518 else if (OutputDimension == 2)
519 {
520 const double tmp[] = { jsh[5], jsh[4], jsh[4], jsh[2] };
521 FastBitwiseCopy(jsh_tmp, tmp);
522 }
523 else // the general case
524 {
525 for (unsigned int j = 0; j < OutputDimension; ++j)
526 {
527 for (unsigned int i = 0; i <= j; ++i)
528 {
529 jsh_tmp[j * OutputDimension + i] =
530 jsh[(OutputDimension - j) + (OutputDimension - i) * (OutputDimension - i + 1) / 2];
531 if (i != j)
532 {
533 jsh_tmp[i * OutputDimension + j] = jsh_tmp[j * OutputDimension + i];
534 }
535 }
536 }
537 }
538
540 for (unsigned int i = 0; i < OutputDimension; ++i) // row
541 {
542 for (unsigned int j = 0; j < OutputDimension; ++j) // column
543 {
544 double accum = directionCosines[i] * jsh_tmp[j];
545 for (unsigned int k = 1; k < OutputDimension; ++k)
546 {
547 accum += directionCosines[k * OutputDimension + i] * jsh_tmp[k * OutputDimension + j];
548 }
549 matrixProduct[i * OutputDimension + j] = accum;
550 }
551 }
552
554 for (unsigned int i = 0; i < OutputDimension; ++i) // row
555 {
556 for (unsigned int j = 0; j < OutputDimension; ++j) // column
557 {
558 double accum = matrixProduct[i * OutputDimension] * directionCosines[j];
559 for (unsigned int k = 1; k < OutputDimension; ++k)
560 {
561 accum += matrixProduct[i * OutputDimension + k] * directionCosines[k * OutputDimension + j];
562 }
563 jsh_out[i * OutputDimension + j] = accum;
564 }
565 }
566
568 unsigned long offset = 0;
569 for (unsigned int i = 0; i < OutputDimension; ++i)
570 {
571 offset = i * (OutputDimension * OutputDimension * (BSplineNumberOfIndices * OutputDimension + 1));
572 for (unsigned int j = 0; j < OutputDimension * OutputDimension; ++j)
573 {
574 jsh_out[j + offset] = jsh_out[j];
575 }
576 }
577
579 jsh_out += OutputDimension * OutputDimension * OutputDimension;
580
581 } // end GetJacobianOfSpatialHessian()
582
583
584private:
585 template <typename T>
586 static void
587 FastBitwiseCopy(T & destination, const T & source)
588 {
589 std::memcpy(&destination, &source, sizeof(T));
590 }
591
592 template <typename T1, typename T2>
593 static void
594 FastBitwiseCopy(const T1 &, const T2 &)
595 {
596 assert(!"This FastBitwiseCopy overload should not be called!");
597 }
598};
599
600
601} // end namespace itk
602
603#endif /* itkRecursiveBSplineTransformImplementation_h */
Returns the weights over the support region used for B-spline interpolation/reconstruction.
itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices)
static void GetSpatialHessian(InternalFloatType *const sh, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D, const double *const derivativeWeights1D, const double *const hessianWeights1D)
static void GetJacobianOfSpatialHessian(InternalFloatType *&jsh_out, const double *const weights1D, const double *const derivativeWeights1D, const double *const hessianWeights1D, const double *const directionCosines, const InternalFloatType *const jsh)
static void GetSpatialJacobian(InternalFloatType *const sj, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D, const double *const derivativeWeights1D)
static void EvaluateJacobianWithImageGradientProduct(TScalar *&imageJacobian, const InternalFloatType *const movingImageGradient, const double *const weights1D, const double value)
static void TransformPoint(TScalar *const opp, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D)
static void GetJacobianOfSpatialJacobian(InternalFloatType *&jsj_out, const double *const weights1D, const double *const derivativeWeights1D, const double *const directionCosines, const InternalFloatType *const jsj)
static void GetJacobian(TScalar *&jacobians, const double *const weights1D, const double value)
static void ComputeNonZeroJacobianIndices(unsigned long *&nzji, const unsigned long parametersPerDim, const unsigned long currentIndex, const OffsetValueType *const gridOffsetTable)
This helper class contains the actual implementation of the recursive B-spline transform.
static void TransformPoint(TScalar *const opp, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D)
itkStaticConstMacro(HelperConstVariable, unsigned int,(SpaceDimension - 1) *(SplineOrder+1))
static void GetSpatialHessian(InternalFloatType *const sh, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D, const double *const derivativeWeights1D, const double *const hessianWeights1D)
static void ComputeNonZeroJacobianIndices(unsigned long *&nzji, const unsigned long parametersPerDim, unsigned long currentIndex, const OffsetValueType *const gridOffsetTable)
static void GetJacobianOfSpatialHessian(InternalFloatType *&jsh_out, const double *const weights1D, const double *const derivativeWeights1D, const double *const hessianWeights1D, const double *const directionCosines, const InternalFloatType *const jsh)
static void GetSpatialJacobian(InternalFloatType *const sj, const TScalar *const *const mu, const OffsetValueType *const gridOffsetTable, const double *const weights1D, const double *const derivativeWeights1D)
static void GetJacobian(TScalar *&jacobians, const double *const weights1D, const double value)
static void EvaluateJacobianWithImageGradientProduct(TScalar *&imageJacobian, const InternalFloatType *const movingImageGradient, const double *const weights1D, const double value)
static void GetJacobianOfSpatialJacobian(InternalFloatType *&jsj_out, const double *const weights1D, const double *const derivativeWeights1D, const double *const directionCosines, const InternalFloatType *const jsj)
itkStaticConstMacro(BSplineNumberOfIndices, unsigned int, RecursiveBSplineWeightFunctionType::NumberOfIndices)


Generated on 2023-01-13 for elastix by doxygen 1.9.6 elastix logo