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 itkShapePriorMAPCostFunction_hxx
19 #define itkShapePriorMAPCostFunction_hxx
20
21 #include "itkShapePriorMAPCostFunction.h"
22
23 namespace itk
24 {
25 /**
26 * Constructor
27 */
28 template< typename TFeatureImage, typename TOutputPixel >
29 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
ShapePriorMAPCostFunction()30 ::ShapePriorMAPCostFunction()
31 {
32 m_GaussianFunction = GaussianKernelFunction<double>::New();
33 m_ShapeParameterMeans = ArrayType(1);
34 m_ShapeParameterMeans.Fill(0.0);
35 m_ShapeParameterStandardDeviations = ArrayType(1);
36 m_ShapeParameterStandardDeviations.Fill(0.0);
37 m_Weights.Fill(1.0);
38 }
39
40 /**
41 * PrintSelf
42 */
43 template< typename TFeatureImage, typename TOutputPixel >
44 void
45 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
PrintSelf(std::ostream & os,Indent indent) const46 ::PrintSelf(std::ostream & os, Indent indent) const
47 {
48 Superclass::PrintSelf(os, indent);
49 os << indent << "ShapeParameterMeans: " << m_ShapeParameterMeans << std::endl;
50 os << indent << "ShapeParameterStandardDeviations: ";
51 os << m_ShapeParameterStandardDeviations << std::endl;
52 os << indent << "Weights: " << m_Weights << std::endl;
53 }
54
55 /**
56 *
57 */
58 template< typename TFeatureImage, typename TOutputPixel >
59 typename ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
60 ::MeasureType
61 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
ComputeLogInsideTerm(const ParametersType & parameters) const62 ::ComputeLogInsideTerm(const ParametersType & parameters) const
63 {
64 this->m_ShapeFunction->SetParameters(parameters);
65
66 typename NodeContainerType::ConstIterator iter = this->GetActiveRegion()->Begin();
67 typename NodeContainerType::ConstIterator end = this->GetActiveRegion()->End();
68
69 MeasureType counter = 0.0;
70
71 // count the number of pixels inside the current contour but outside the
72 // current shape
73 while ( iter != end )
74 {
75 NodeType node = iter.Value();
76 typename ShapeFunctionType::PointType point;
77
78 this->GetFeatureImage()->TransformIndexToPhysicalPoint(node.GetIndex(), point);
79
80 if ( node.GetValue() <= 0.0 )
81 {
82 double value = this->m_ShapeFunction->Evaluate(point);
83 if ( value > 0.0 )
84 {
85 counter += 1.0;
86 }
87 else if ( value > -1.0 )
88 {
89 counter += ( 1.0 + value );
90 }
91 }
92
93 ++iter;
94 }
95
96 MeasureType output = counter * m_Weights[0];
97 return output;
98 }
99
100 /**
101 *
102 */
103 template< typename TFeatureImage, typename TOutputPixel >
104 typename ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
105 ::MeasureType
106 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
ComputeLogShapePriorTerm(const ParametersType & parameters) const107 ::ComputeLogShapePriorTerm(const ParametersType & parameters) const
108 {
109 // assume the shape parameters is from a independent gaussian distributions
110 MeasureType measure = 0.0;
111
112 for ( unsigned int j = 0; j < this->m_ShapeFunction->GetNumberOfShapeParameters(); j++ )
113 {
114 measure += itk::Math::sqr( ( parameters[j] - m_ShapeParameterMeans[j] )
115 / m_ShapeParameterStandardDeviations[j] );
116 }
117 measure *= m_Weights[2];
118 return measure;
119 }
120
121 /**
122 *
123 */
124 template< typename TFeatureImage, typename TOutputPixel >
125 typename ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
126 ::MeasureType
127 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
ComputeLogGradientTerm(const ParametersType & parameters) const128 ::ComputeLogGradientTerm(const ParametersType & parameters) const
129 {
130 this->m_ShapeFunction->SetParameters(parameters);
131
132 typename NodeContainerType::ConstIterator iter = this->GetActiveRegion()->Begin();
133 typename NodeContainerType::ConstIterator end = this->GetActiveRegion()->End();
134 MeasureType sum = 0.0;
135
136 // Assume that ( 1 - FeatureImage ) approximates a Gaussian (zero mean, unit
137 // variance)
138 // along the normal of the evolving contour.
139 // The GradientTerm is then given by a Laplacian of the goodness of fit of
140 // the Gaussian.
141 while ( iter != end )
142 {
143 NodeType node = iter.Value();
144 typename ShapeFunctionType::PointType point;
145
146 this->GetFeatureImage()->TransformIndexToPhysicalPoint(node.GetIndex(), point);
147
148 sum += itk::Math::sqr( m_GaussianFunction->Evaluate( this->m_ShapeFunction->Evaluate(point) )
149 - 1.0 + this->GetFeatureImage()->GetPixel( node.GetIndex() ) );
150
151 ++iter;
152 }
153
154 sum *= m_Weights[1];
155 // std::cout << sum << " ";
156 return sum;
157 }
158
159 /**
160 *
161 */
162 template< typename TFeatureImage, typename TOutputPixel >
163 typename ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
164 ::MeasureType
165 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
ComputeLogPosePriorTerm(const ParametersType & itkNotUsed (parameters)) const166 ::ComputeLogPosePriorTerm( const ParametersType & itkNotUsed(parameters) ) const
167 {
168 return 0.0;
169 }
170
171 /**
172 *
173 */
174 template< typename TFeatureImage, typename TOutputPixel >
175 void
176 ShapePriorMAPCostFunction< TFeatureImage, TOutputPixel >
Initialize()177 ::Initialize()
178 {
179 this->Superclass::Initialize();
180
181 // check if the mean and variances array are of the right size
182 if ( m_ShapeParameterMeans.Size() <
183 this->m_ShapeFunction->GetNumberOfShapeParameters() )
184 {
185 itkExceptionMacro(<< "ShapeParameterMeans does not have at least "
186 << this->m_ShapeFunction->GetNumberOfShapeParameters()
187 << " number of elements.");
188 }
189
190 if ( m_ShapeParameterStandardDeviations.Size() <
191 this->m_ShapeFunction->GetNumberOfShapeParameters() )
192 {
193 itkExceptionMacro(<< "ShapeParameterStandardDeviations does not have at least "
194 << this->m_ShapeFunction->GetNumberOfShapeParameters()
195 << " number of elements.");
196 }
197 }
198 } // end namespace itk
199
200 #endif
201