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