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 itkMutualInformationHistogramImageToImageMetric_hxx
19 #define itkMutualInformationHistogramImageToImageMetric_hxx
20 
21 #include "itkMutualInformationHistogramImageToImageMetric.h"
22 #include "itkHistogram.h"
23 
24 namespace itk
25 {
26 template< typename TFixedImage, typename TMovingImage >
27 typename MutualInformationHistogramImageToImageMetric< TFixedImage, TMovingImage >::MeasureType
28 MutualInformationHistogramImageToImageMetric< TFixedImage, TMovingImage >
EvaluateMeasure(HistogramType & histogram) const29 ::EvaluateMeasure(HistogramType & histogram) const
30 {
31   MeasureType entropyX = NumericTraits< MeasureType >::ZeroValue();
32   MeasureType entropyY = NumericTraits< MeasureType >::ZeroValue();
33   MeasureType jointEntropy = NumericTraits< MeasureType >::ZeroValue();
34 
35   using HistogramFrequencyRealType = typename NumericTraits< HistogramFrequencyType >::RealType;
36 
37   auto totalFreq = static_cast< HistogramFrequencyRealType >( histogram.GetTotalFrequency() );
38 
39   for ( unsigned int i = 0; i < this->GetHistogramSize()[0]; i++ )
40     {
41     auto freq = static_cast< HistogramFrequencyRealType >( histogram.GetFrequency(i, 0) );
42     if ( freq > 0 )
43       {
44       entropyX += freq * std::log(freq);
45       }
46     }
47 
48   entropyX = -entropyX / static_cast< MeasureType >( totalFreq ) + std::log(totalFreq);
49 
50   for ( unsigned int i = 0; i < this->GetHistogramSize()[1]; i++ )
51     {
52     auto freq = static_cast< HistogramFrequencyRealType >( histogram.GetFrequency(i, 1) );
53     if ( freq > 0 )
54       {
55       entropyY += freq * std::log(freq);
56       }
57     }
58 
59   entropyY = -entropyY / static_cast< MeasureType >( totalFreq ) + std::log(totalFreq);
60 
61   HistogramIteratorType it = histogram.Begin();
62   HistogramIteratorType end = histogram.End();
63   while ( it != end )
64     {
65     auto freq = static_cast< HistogramFrequencyRealType >( it.GetFrequency() );
66     if ( freq > 0 )
67       {
68       jointEntropy += freq * std::log(freq);
69       }
70     ++it;
71     }
72 
73   jointEntropy = -jointEntropy
74                  / static_cast< MeasureType >( totalFreq ) + std::log(totalFreq);
75 
76   return entropyX + entropyY - jointEntropy;
77 }
78 } // End namespace itk
79 
80 #endif // itkMutualInformationHistogramImageToImageMetric_hxx
81