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 itkJensenHavrdaCharvatTsallisPointSetToPointSetMetricv4_hxx
19 #define itkJensenHavrdaCharvatTsallisPointSetToPointSetMetricv4_hxx
20
21 #include "itkMath.h"
22 #include "itkJensenHavrdaCharvatTsallisPointSetToPointSetMetricv4.h"
23
24 namespace itk {
25
26 /** Constructor */
27 template<typename TPointSet, class TInternalComputationValueType>
28 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4()29 ::JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4() :
30 m_PointSetSigma( static_cast<RealType>( 1.0 ) ),
31 m_KernelSigma( static_cast<RealType>( 10.0 ) ),
32 m_CovarianceKNeighborhood( static_cast<unsigned int>( 5 ) ),
33 m_EvaluationKNeighborhood( static_cast<unsigned int>( 50 ) ),
34 m_Alpha( static_cast<RealType>( 1.0 ) ),
35 m_TotalNumberOfPoints( 0 ),
36 m_Prefactor0( 0.0 ),
37 m_Prefactor1( 0.0 )
38 {
39 }
40
41 /** Initialize the metric */
42 template<typename TPointSet, class TInternalComputationValueType>
43 void
44 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
Initialize()45 ::Initialize()
46 {
47 Superclass::Initialize();
48
49 // Initialize the moving density function
50 this->m_MovingDensityFunction = DensityFunctionType::New();
51 this->m_MovingDensityFunction->SetKernelSigma( this->m_KernelSigma );
52 this->m_MovingDensityFunction->SetRegularizationSigma( this->m_PointSetSigma );
53 this->m_MovingDensityFunction->SetNormalize( true );
54 this->m_MovingDensityFunction->SetUseAnisotropicCovariances( this->m_UseAnisotropicCovariances );
55 this->m_MovingDensityFunction->SetCovarianceKNeighborhood( this->m_CovarianceKNeighborhood );
56 this->m_MovingDensityFunction->SetEvaluationKNeighborhood( this->m_EvaluationKNeighborhood );
57 this->m_MovingDensityFunction->SetInputPointSet( this->m_MovingTransformedPointSet );
58
59 // Pre-calc some values for efficiency
60 this->m_TotalNumberOfPoints = static_cast<RealType>( this->m_NumberOfValidPoints + this->m_MovingDensityFunction->GetInputPointSet()->GetNumberOfPoints() );
61 this->m_Prefactor0 = -1.0 / static_cast<RealType>( this->m_TotalNumberOfPoints );
62 if( this->m_Alpha != 1.0 )
63 {
64 this->m_Prefactor0 /= ( this->m_Alpha - 1.0 );
65 }
66 this->m_Prefactor1 = 1.0 / ( this->m_TotalNumberOfPoints * this->m_TotalNumberOfPoints );
67 }
68
69 template<typename TPointSet, class TInternalComputationValueType>
70 typename JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
71 ::MeasureType
72 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
GetLocalNeighborhoodValue(const PointType & point,const PixelType & itkNotUsed (pixel)) const73 ::GetLocalNeighborhoodValue( const PointType & point, const PixelType & itkNotUsed( pixel ) ) const
74 {
75 MeasureType value;
76 LocalDerivativeType derivative;
77 this->ComputeValueAndDerivative( point, value, derivative, true, false );
78 return value;
79 }
80
81 /** Get both the match Measure and the Derivative Measure */
82 template<typename TPointSet, class TInternalComputationValueType>
83 void
84 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
GetLocalNeighborhoodValueAndDerivative(const PointType & point,MeasureType & value,LocalDerivativeType & derivative,const PixelType & itkNotUsed (pixel)) const85 ::GetLocalNeighborhoodValueAndDerivative( const PointType & point, MeasureType & value, LocalDerivativeType & derivative, const PixelType & itkNotUsed( pixel ) ) const
86 {
87 this->ComputeValueAndDerivative( point, value, derivative, true, true );
88 }
89
90 template<typename TPointSet, class TInternalComputationValueType>
91 void
92 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
ComputeValueAndDerivative(const PointType & samplePoint,MeasureType & value,LocalDerivativeType & derivativeReturn,bool calcValue,bool calcDerivative) const93 ::ComputeValueAndDerivative( const PointType & samplePoint, MeasureType & value, LocalDerivativeType & derivativeReturn, bool calcValue, bool calcDerivative ) const
94 {
95 if( calcDerivative )
96 {
97 derivativeReturn.Fill( NumericTraits<DerivativeValueType>::ZeroValue() );
98 }
99 value = NumericTraits<MeasureType>::ZeroValue();
100
101 /**
102 * first term only
103 */
104 typename PointSetType::PointIdentifier numberOfMovingPoints = this->m_MovingDensityFunction->GetInputPointSet()->GetNumberOfPoints();
105 RealType probabilityStar = this->m_MovingDensityFunction->Evaluate( samplePoint ) * static_cast<RealType>( numberOfMovingPoints );
106
107 probabilityStar /= this->m_TotalNumberOfPoints;
108
109 if( Math::AlmostEquals( probabilityStar, NumericTraits<RealType>::ZeroValue() ) )
110 {
111 return;
112 }
113
114 if( calcValue )
115 {
116 RealType realOne = NumericTraits<RealType>::OneValue();
117 if( Math::AlmostEquals( this->m_Alpha, realOne ) )
118 {
119 value = ( std::log( probabilityStar ) );
120 }
121 else
122 {
123 value = realOne * ( std::pow( probabilityStar, static_cast<RealType>( this->m_Alpha - realOne ) ) );
124 }
125 value *= this->m_Prefactor0;
126 }
127
128 if( calcDerivative )
129 {
130 RealType probabilityStarFactor = std::pow( probabilityStar, static_cast<RealType>( 2.0 - this->m_Alpha ) );
131
132 typename DensityFunctionType::NeighborsIdentifierType neighbors;
133 this->m_MovingDensityFunction->GetPointsLocator()->FindClosestNPoints( samplePoint, this->m_EvaluationKNeighborhood, neighbors );
134
135 for( SizeValueType n = 0; n < neighbors.size(); n++ )
136 {
137 RealType gaussian = this->m_MovingDensityFunction->GetGaussian( neighbors[n] )->Evaluate( samplePoint );
138
139 if( Math::AlmostEquals( gaussian, NumericTraits<RealType>::ZeroValue() ) )
140 {
141 continue;
142 }
143
144 typename GaussianType::MeanVectorType mean = this->m_MovingDensityFunction->GetGaussian( neighbors[n] )->GetMean();
145
146 Array<CoordRepType> diffMean( PointDimension );
147 for( unsigned int i = 0; i < PointDimension; i++ )
148 {
149 diffMean[i] = mean[i] - samplePoint[i];
150 }
151
152 if( this->m_UseAnisotropicCovariances )
153 {
154 typename GaussianType::CovarianceMatrixType Ci = this->m_MovingDensityFunction->GetGaussian( neighbors[n] )->GetInverseCovariance();
155 diffMean = Ci * diffMean;
156 }
157 else
158 {
159 diffMean /= this->m_MovingDensityFunction->GetGaussian( neighbors[n] )->GetCovariance()( 0, 0 );
160 }
161
162 DerivativeValueType factor = this->m_Prefactor1 * gaussian / probabilityStarFactor;
163 for( unsigned int i = 0; i < PointDimension; i++ )
164 {
165 derivativeReturn[i] += diffMean[i] * factor;
166 }
167 }
168 }
169 }
170
171 template<typename TPointSet, class TInternalComputationValueType>
172 typename LightObject::Pointer
173 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
InternalClone() const174 ::InternalClone() const
175 {
176 typename Self::Pointer rval = Self::New();
177 rval->SetMovingPointSet( this->m_MovingPointSet );
178 rval->SetFixedPointSet( this->m_FixedPointSet );
179 rval->SetPointSetSigma( this->m_PointSetSigma );
180 rval->SetEvaluationKNeighborhood( this->m_EvaluationKNeighborhood );
181 rval->SetAlpha( this->m_Alpha );
182 rval->SetKernelSigma( this->m_KernelSigma );
183 rval->SetCovarianceKNeighborhood( this->m_CovarianceKNeighborhood );
184 rval->SetUseAnisotropicCovariances( this->m_UseAnisotropicCovariances );
185
186 return rval.GetPointer();
187 }
188
189 template<typename TPointSet, class TInternalComputationValueType>
190 void
191 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4<TPointSet, TInternalComputationValueType>
PrintSelf(std::ostream & os,Indent indent) const192 ::PrintSelf( std::ostream& os, Indent indent ) const
193 {
194 Superclass::PrintSelf( os, indent );
195
196 os << indent << "Alpha: " << this->m_Alpha << std::endl;
197 os << indent << "Point set sigma: " << this->m_PointSetSigma << std::endl;
198 os << indent << "Evaluation k-neighborhood: " << this->m_EvaluationKNeighborhood << std::endl;
199
200 if( this->m_UseAnisotropicCovariances )
201 {
202 os << indent << "Kernel sigma: " << this->m_KernelSigma << std::endl;
203 os << indent << "Covariance k-neighborhood: " << this->m_CovarianceKNeighborhood << std::endl;
204 }
205 else
206 {
207 os << indent << "Isotropic covariances are used." << std::endl;
208 }
209
210 }
211 } // end namespace itk
212
213 #endif
214