1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
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 itkBSplineInterpolationWeightFunction_hxx
19 #define itkBSplineInterpolationWeightFunction_hxx
20
21 #include "itkBSplineInterpolationWeightFunction.h"
22 #include "itkImage.h"
23 #include "itkMatrix.h"
24 #include "itkMath.h"
25 #include "itkImageRegionConstIteratorWithIndex.h"
26
27 namespace itk
28 {
29 /** Constructor */
30 template< typename TCoordRep, unsigned int VSpaceDimension,
31 unsigned int VSplineOrder >
32 BSplineInterpolationWeightFunction< TCoordRep, VSpaceDimension, VSplineOrder >
BSplineInterpolationWeightFunction()33 ::BSplineInterpolationWeightFunction()
34 {
35 // Initialize the number of weights;
36 m_NumberOfWeights =
37 static_cast< unsigned int >( std::pow( static_cast< double >( SplineOrder + 1 ),
38 static_cast< double >( SpaceDimension ) ) );
39
40 // Initialize support region is a hypercube of length SplineOrder + 1
41 m_SupportSize.Fill(SplineOrder + 1);
42
43 // Initialize offset to index lookup table
44 m_OffsetToIndexTable.set_size(m_NumberOfWeights, SpaceDimension);
45
46 using CharImageType = Image< char, SpaceDimension >;
47 typename CharImageType::Pointer tempImage = CharImageType::New();
48 tempImage->SetRegions(m_SupportSize);
49 tempImage->Allocate(true); // initialize buffer to zero
50
51 using IteratorType = ImageRegionConstIteratorWithIndex< CharImageType >;
52 IteratorType iterator( tempImage, tempImage->GetBufferedRegion() );
53 unsigned int counter = 0;
54
55 while ( !iterator.IsAtEnd() )
56 {
57 for ( unsigned int j = 0; j < SpaceDimension; j++ )
58 {
59 m_OffsetToIndexTable[counter][j] = iterator.GetIndex()[j];
60 }
61 ++counter;
62 ++iterator;
63 }
64
65 // Initialize the interpolation kernel
66 m_Kernel = KernelType::New();
67 }
68
69 /**
70 * Standard "PrintSelf" method
71 */
72 template< typename TCoordRep, unsigned int VSpaceDimension,
73 unsigned int VSplineOrder >
74 void
75 BSplineInterpolationWeightFunction< TCoordRep, VSpaceDimension, VSplineOrder >
PrintSelf(std::ostream & os,Indent indent) const76 ::PrintSelf(
77 std::ostream & os,
78 Indent indent) const
79 {
80 Superclass::PrintSelf(os, indent);
81
82 os << indent << "NumberOfWeights: " << m_NumberOfWeights << std::endl;
83 os << indent << "SupportSize: " << m_SupportSize << std::endl;
84 }
85
86 /** Compute weights for interpolation at continuous index position */
87 template< typename TCoordRep, unsigned int VSpaceDimension,
88 unsigned int VSplineOrder >
89 typename BSplineInterpolationWeightFunction< TCoordRep, VSpaceDimension,
90 VSplineOrder >
91 ::WeightsType
92 BSplineInterpolationWeightFunction< TCoordRep, VSpaceDimension, VSplineOrder >
Evaluate(const ContinuousIndexType & index) const93 ::Evaluate(
94 const ContinuousIndexType & index) const
95 {
96 WeightsType weights(m_NumberOfWeights);
97 IndexType startIndex;
98
99 this->Evaluate(index, weights, startIndex);
100
101 return weights;
102 }
103
104 /** Compute weights for interpolation at continuous index position */
105 template< typename TCoordRep, unsigned int VSpaceDimension,
106 unsigned int VSplineOrder >
107 void BSplineInterpolationWeightFunction< TCoordRep, VSpaceDimension,
108 VSplineOrder >
Evaluate(const ContinuousIndexType & index,WeightsType & weights,IndexType & startIndex) const109 ::Evaluate(
110 const ContinuousIndexType & index,
111 WeightsType & weights,
112 IndexType & startIndex) const
113 {
114 unsigned int j, k;
115
116 // Find the starting index of the support region
117 for ( j = 0; j < SpaceDimension; j++ )
118 {
119 // Note that the expression passed to Math::Floor is adapted to work around
120 // a compiler bug which caused endless compilations (apparently), by
121 // Visual C++ 2015 Update 3, on 64-bit builds of Release configurations.
122 startIndex[j] = Math::Floor< IndexValueType >(index[j] + 0.5 - SplineOrder/2.0);
123 }
124
125 // Compute the weights
126 Matrix< double, SpaceDimension, SplineOrder + 1 > weights1D;
127 for ( j = 0; j < SpaceDimension; j++ )
128 {
129 double x = index[j] - static_cast< double >( startIndex[j] );
130
131 for ( k = 0; k <= SplineOrder; k++ )
132 {
133 weights1D[j][k] = m_Kernel->Evaluate(x);
134 x -= 1.0;
135 }
136 }
137
138 for ( k = 0; k < m_NumberOfWeights; k++ )
139 {
140 weights[k] = 1.0;
141
142 for ( j = 0; j < SpaceDimension; j++ )
143 {
144 weights[k] *= weights1D[j][m_OffsetToIndexTable[k][j]];
145 }
146 }
147 }
148 } // end namespace itk
149
150 #endif
151