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 itkExpectationBasedPointSetToPointSetMetricv4_hxx
19 #define itkExpectationBasedPointSetToPointSetMetricv4_hxx
20
21 #include "itkExpectationBasedPointSetToPointSetMetricv4.h"
22 #include "itkArray.h"
23 #include "itkCompensatedSummation.h"
24
25 namespace itk {
26
27 /** Constructor */
28 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
29 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
ExpectationBasedPointSetToPointSetMetricv4()30 ::ExpectationBasedPointSetToPointSetMetricv4() :
31 m_PointSetSigma( 1.0 ),
32 m_PreFactor( 0.0 ),
33 m_Denominator( 0.0 )
34
35 {
36 }
37
38 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
39 void
40 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
Initialize()41 ::Initialize()
42 {
43 Superclass::Initialize();
44
45 if( this->m_PointSetSigma <= NumericTraits<CoordRepType>::epsilon() )
46 {
47 itkExceptionMacro("m_PointSetSigma is too small. <= epsilon");
48 }
49 this->m_PreFactor = 1.0 / ( std::sqrt( 2 * itk::Math::pi ) * this->m_PointSetSigma );
50 this->m_Denominator = 2.0 * itk::Math::sqr( this->m_PointSetSigma );
51 }
52
53 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
54 typename ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
55 ::MeasureType
56 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
GetLocalNeighborhoodValue(const PointType & point,const PixelType & itkNotUsed (pixel)) const57 ::GetLocalNeighborhoodValue( const PointType & point, const PixelType & itkNotUsed( pixel ) ) const
58 {
59 CompensatedSummation< MeasureType > localValue;
60
61 NeighborsIdentifierType neighborhood;
62 this->m_MovingTransformedPointsLocator->FindClosestNPoints( point, this->m_EvaluationKNeighborhood, neighborhood );
63
64 for( NeighborsIterator it = neighborhood.begin(); it != neighborhood.end(); ++it )
65 {
66 PointType neighbor = this->m_MovingTransformedPointSet->GetPoint( *it );
67 const MeasureType distance = point.SquaredEuclideanDistanceTo( neighbor );
68 localValue -= this->m_PreFactor * std::exp( -distance / this->m_Denominator );
69 }
70
71 return localValue.GetSum();
72 }
73
74 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
75 void
76 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
GetLocalNeighborhoodValueAndDerivative(const PointType & point,MeasureType & measure,LocalDerivativeType & localDerivative,const PixelType & itkNotUsed (pixel)) const77 ::GetLocalNeighborhoodValueAndDerivative( const PointType & point,
78 MeasureType &measure, LocalDerivativeType &localDerivative, const PixelType & itkNotUsed( pixel ) ) const
79 {
80 Array<MeasureType> measureValues;
81 measureValues.SetSize( this->m_EvaluationKNeighborhood );
82 measureValues.Fill( 0.0 );
83
84 CompensatedSummation< MeasureType > measureSum;
85
86 localDerivative.Fill( 0.0 );
87
88 PointType weightedPoint;
89 weightedPoint.Fill( 0.0 );
90
91 NeighborsIdentifierType neighborhood;
92
93 this->m_MovingTransformedPointsLocator->FindClosestNPoints( point, this->m_EvaluationKNeighborhood, neighborhood );
94
95 for( NeighborsIterator it = neighborhood.begin(); it != neighborhood.end(); ++it )
96 {
97 PointType neighbor = this->m_MovingTransformedPointSet->GetPoint( *it );
98 const MeasureType distance = point.SquaredEuclideanDistanceTo( neighbor );
99 measureValues[it - neighborhood.begin()] = -this->m_PreFactor * std::exp( -distance / this->m_Denominator );
100 measureSum += measureValues[it - neighborhood.begin()];
101 }
102
103 measure = measureSum.GetSum();
104 if ( std::fabs(measure) <= NumericTraits<MeasureType>::epsilon() )
105 {
106 return;
107 }
108
109 for( NeighborsIterator it = neighborhood.begin(); it != neighborhood.end(); ++it )
110 {
111 PointType neighbor = this->m_MovingTransformedPointSet->GetPoint( *it );
112 VectorType neighborVector = neighbor.GetVectorFromOrigin();
113 weightedPoint += ( neighborVector * measureValues[it - neighborhood.begin()] / measure );
114 }
115
116 const MeasureType distance = point.SquaredEuclideanDistanceTo( weightedPoint );
117
118 const MeasureType weight = this->m_PreFactor * std::exp( -distance / this->m_Denominator ) / -measure;
119
120 VectorType force = ( weightedPoint - point ) * weight;
121
122 for( unsigned int d = 0; d < localDerivative.Size(); d++ )
123 {
124 localDerivative[d] = force[d];
125 }
126 }
127
128 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
129 typename LightObject::Pointer
130 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
InternalClone() const131 ::InternalClone() const
132 {
133 typename Self::Pointer rval = Self::New();
134 rval->SetMovingPointSet( this->m_MovingPointSet );
135 rval->SetFixedPointSet( this->m_FixedPointSet );
136 rval->SetPointSetSigma( this->m_PointSetSigma );
137 rval->SetEvaluationKNeighborhood( this->m_EvaluationKNeighborhood );
138
139 return rval.GetPointer();
140 }
141
142 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
143 void
144 ExpectationBasedPointSetToPointSetMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
PrintSelf(std::ostream & os,Indent indent) const145 ::PrintSelf( std::ostream& os, Indent indent ) const
146 {
147 Superclass::PrintSelf( os, indent );
148
149 os << indent << "PointSetSigma: " << this->m_PointSetSigma << std::endl;
150 os << indent << "EvaluateKNeighborhood: " << this->m_EvaluationKNeighborhood << std::endl;
151 }
152 } // end namespace itk
153
154
155 #endif
156