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 itkWarpHarmonicEnergyCalculator_hxx
19 #define itkWarpHarmonicEnergyCalculator_hxx
20
21 #include "itkWarpHarmonicEnergyCalculator.h"
22 #include "itkNeighborhoodAlgorithm.h"
23 #include "itkImageRegionIterator.h"
24 #include "itkZeroFluxNeumannBoundaryCondition.h"
25
26 #include "vnl/vnl_matrix.h"
27 #include "itkMath.h"
28
29 namespace itk
30 {
31
32 template< typename TInputImage >
33 WarpHarmonicEnergyCalculator< TInputImage >
WarpHarmonicEnergyCalculator()34 ::WarpHarmonicEnergyCalculator() :
35 m_HarmonicEnergy( 0.0 ),
36 m_RegionSetByUser( false ),
37 m_UseImageSpacing( true )
38 {
39 m_Image = TInputImage::New();
40 for ( unsigned int i = 0; i < ImageDimension; i++ )
41 {
42 m_NeighborhoodRadius[i] = 1; // radius of neighborhood we will use
43 m_DerivativeWeights[i] = 1.0;
44 }
45 }
46
47 template< typename TInputImage >
48 void
49 WarpHarmonicEnergyCalculator< TInputImage >
SetUseImageSpacing(bool f)50 ::SetUseImageSpacing(bool f)
51 {
52 if ( m_UseImageSpacing == f )
53 {
54 return;
55 }
56
57 // Only reset the weights if they were previously set to the image spacing,
58 // otherwise, the user may have provided their own weightings.
59 if ( f == false && m_UseImageSpacing == true )
60 {
61 for ( unsigned int i = 0; i < ImageDimension; ++i )
62 {
63 m_DerivativeWeights[i] = 1.0;
64 }
65 }
66
67 m_UseImageSpacing = f;
68 }
69
70 template< typename TInputImage >
71 void
72 WarpHarmonicEnergyCalculator< TInputImage >
Compute()73 ::Compute()
74 {
75 if ( !m_RegionSetByUser )
76 {
77 m_Region = m_Image->GetRequestedRegion();
78 }
79
80 // Set the weights on the derivatives.
81 // Are we using image spacing in the calculations? If so we must update now
82 // in case our input image has changed.
83 if ( m_UseImageSpacing == true )
84 {
85 for ( unsigned int i = 0; i < ImageDimension; i++ )
86 {
87 if ( m_Image->GetSpacing()[i] <= 0.0 )
88 {
89 itkExceptionMacro(<< "Image spacing in dimension " << i << " is zero.");
90 }
91 m_DerivativeWeights[i] = 1.0 / static_cast< double >( m_Image->GetSpacing()[i] );
92 }
93 }
94
95 m_HarmonicEnergy = 0.0;
96
97 ZeroFluxNeumannBoundaryCondition< ImageType > nBc;
98 ConstNeighborhoodIteratorType bIt;
99
100 // Find the data-set boundary "faces"
101 typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< ImageType >::
102 FaceListType faceList;
103 NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< ImageType > bC;
104 faceList = bC(m_Image, m_Region, m_NeighborhoodRadius);
105
106 typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< ImageType >::
107 FaceListType::iterator fIt;
108 fIt = faceList.begin();
109
110 // Process each of the data set faces. The iterator is reinitialized on each
111 // face so that it can determine whether or not to check for boundary
112 // conditions.
113 for ( fIt = faceList.begin(); fIt != faceList.end(); ++fIt )
114 {
115 bIt = ConstNeighborhoodIteratorType(m_NeighborhoodRadius,
116 m_Image,
117 *fIt);
118 bIt.OverrideBoundaryCondition(&nBc);
119 bIt.GoToBegin();
120
121 while ( !bIt.IsAtEnd() )
122 {
123 m_HarmonicEnergy += this->EvaluateAtNeighborhood(bIt);
124 ++bIt;
125 }
126 }
127
128 m_HarmonicEnergy /= m_Region.GetNumberOfPixels();
129 }
130
131 template< typename TInputImage >
132 double
133 WarpHarmonicEnergyCalculator< TInputImage >
EvaluateAtNeighborhood(ConstNeighborhoodIteratorType & it) const134 ::EvaluateAtNeighborhood(ConstNeighborhoodIteratorType & it) const
135 {
136 // Simple method using field derivatives
137
138 vnl_matrix_fixed< double, ImageDimension, VectorDimension > J;
139
140 PixelType next, prev;
141
142 double weight;
143
144 for ( unsigned int i = 0; i < ImageDimension; ++i )
145 {
146 next = it.GetNext(i);
147 prev = it.GetPrevious(i);
148
149 weight = 0.5 * m_DerivativeWeights[i];
150
151 for ( unsigned int j = 0; j < VectorDimension; ++j )
152 {
153 J[i][j] = weight * ( static_cast< double >( next[j] ) - static_cast< double >( prev[j] ) );
154 }
155
156 // Add one on the diagonal to consider the warping and not only the
157 // deformation field
158 //J[i][i] += 1.0;
159 }
160
161 const double norm = J.fro_norm();
162 return norm * norm;
163 }
164
165 template< typename TInputImage >
166 void
167 WarpHarmonicEnergyCalculator< TInputImage >
SetRegion(const RegionType & region)168 ::SetRegion(const RegionType & region)
169 {
170 m_Region = region;
171 m_RegionSetByUser = true;
172 }
173
174 template< typename TInputImage >
175 void
176 WarpHarmonicEnergyCalculator< TInputImage >
PrintSelf(std::ostream & os,Indent indent) const177 ::PrintSelf(std::ostream & os, Indent indent) const
178 {
179 Superclass::PrintSelf(os, indent);
180
181 os << indent << "HarmonicEnergy: " << m_HarmonicEnergy << std::endl;
182 itkPrintSelfObjectMacro( Image );
183 os << indent << "Region: " << std::endl;
184 m_Region.Print( os, indent.GetNextIndent() );
185 os << indent << "Region set by User: " << m_RegionSetByUser << std::endl;
186 os << indent << "Use image spacing: " << this->m_UseImageSpacing << std::endl;
187 os << indent << "Derivative Weights: " << this->m_DerivativeWeights << std::endl;
188 os << indent << "Neighborhood Radius: " << this->m_NeighborhoodRadius << std::endl;
189 }
190 } // end namespace itk
191
192 #endif
193