From 80b59acdc1ed4926224d27eb9acb13d5c1116b95 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 4 Nov 2011 10:31:35 +0100 Subject: [PATCH] add first version of LabelOverlapMeasuresImageFilter (dice computation) --- itk/itkLabelOverlapMeasuresImageFilter.h | 212 ++++++++++ itk/itkLabelOverlapMeasuresImageFilter.txx | 437 +++++++++++++++++++++ 2 files changed, 649 insertions(+) create mode 100644 itk/itkLabelOverlapMeasuresImageFilter.h create mode 100644 itk/itkLabelOverlapMeasuresImageFilter.txx diff --git a/itk/itkLabelOverlapMeasuresImageFilter.h b/itk/itkLabelOverlapMeasuresImageFilter.h new file mode 100644 index 0000000..978fed3 --- /dev/null +++ b/itk/itkLabelOverlapMeasuresImageFilter.h @@ -0,0 +1,212 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkLabelOverlapMeasuresImageFilter.h,v $ + Language: C++ + Date: $Date: $ + Version: $Revision: $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkLabelOverlapMeasuresImageFilter_h +#define __itkLabelOverlapMeasuresImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkFastMutexLock.h" +#include "itkNumericTraits.h" + +//#include "itk_hash_map.h" +#include "itksys/hash_map.hxx" + +namespace itk { + +/** \class LabelOverlapMeasuresImageFilter + * \brief Computes overlap measures between the set same set of labels of + * pixels of two images. Background is assumed to be 0. + * + * \sa LabelOverlapMeasuresImageFilter + * + * \ingroup MultiThreaded + */ +template +class ITK_EXPORT LabelOverlapMeasuresImageFilter : + public ImageToImageFilter +{ +public: + /** Standard Self typedef */ + typedef LabelOverlapMeasuresImageFilter Self; + typedef ImageToImageFilter Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro( Self ); + + /** Runtime information support. */ + itkTypeMacro( LabelOverlapMeasuresImageFilter, ImageToImageFilter ); + + virtual void VerifyInputInformation() { } + + + /** Image related typedefs. */ + typedef TLabelImage LabelImageType; + typedef typename TLabelImage::Pointer LabelImagePointer; + typedef typename TLabelImage::ConstPointer LabelImageConstPointer; + + typedef typename TLabelImage::RegionType RegionType; + typedef typename TLabelImage::SizeType SizeType; + typedef typename TLabelImage::IndexType IndexType; + + typedef typename TLabelImage::PixelType LabelType; + + /** Type to use form computations. */ + typedef typename NumericTraits::RealType RealType; + + /** \class LabelLabelOverlapMeasuress + * \brief Metrics stored per label */ + class LabelSetMeasures + { + public: + // default constructor + LabelSetMeasures() + { + m_Source = 0; + m_Target = 0; + m_Union = 0; + m_Intersection = 0; + m_SourceComplement = 0; + m_TargetComplement = 0; + } + + // added for completeness + LabelSetMeasures& operator=( const LabelSetMeasures& l ) + { + m_Source = l.m_Source; + m_Target = l.m_Target; + m_Union = l.m_Union; + m_Intersection = l.m_Intersection; + m_SourceComplement = l.m_SourceComplement; + m_TargetComplement = l.m_TargetComplement; + } + + unsigned long m_Source; + unsigned long m_Target; + unsigned long m_Union; + unsigned long m_Intersection; + unsigned long m_SourceComplement; + unsigned long m_TargetComplement; + }; + + /** Type of the map used to store data per label */ + typedef itksys::hash_map MapType; + typedef typename MapType::iterator MapIterator; + typedef typename MapType::const_iterator MapConstIterator; + + /** Image related typedefs. */ + itkStaticConstMacro( ImageDimension, unsigned int, + TLabelImage::ImageDimension ); + + /** Set the source image. */ + void SetSourceImage( const LabelImageType * image ) + { this->SetNthInput( 0, const_cast( image ) ); } + + /** Set the target image. */ + void SetTargetImage( const LabelImageType * image ) + { this->SetNthInput( 1, const_cast( image ) ); } + + /** Get the source image. */ + const LabelImageType * GetSourceImage( void ) + { return this->GetInput( 0 ); } + + /** Get the target image. */ + const LabelImageType * GetTargetImage( void ) + { return this->GetInput( 1 ); } + + /** Get the label set measures */ + MapType GetLabelSetMeasures() + { return this->m_LabelSetMeasures; } + + /** + * tric overlap measures + */ + /** measures over all labels */ + RealType GetTotalOverlap(); + RealType GetUnionOverlap(); + RealType GetMeanOverlap(); + RealType GetVolumeSimilarity(); + RealType GetFalseNegativeError(); + RealType GetFalsePositiveError(); + /** measures over individual labels */ + RealType GetTargetOverlap( LabelType ); + RealType GetUnionOverlap( LabelType ); + RealType GetMeanOverlap( LabelType ); + RealType GetVolumeSimilarity( LabelType ); + RealType GetFalseNegativeError( LabelType ); + RealType GetFalsePositiveError( LabelType ); + /** alternative names */ + RealType GetJaccardCoefficient() + { return this->GetUnionOverlap(); } + RealType GetJaccardCoefficient( LabelType label ) + { return this->GetUnionOverlap( label ); } + RealType GetDiceCoefficient() + { return this->GetMeanOverlap(); } + RealType GetDiceCoefficient( LabelType label ) + { return this->GetMeanOverlap( label ); } + + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro( Input1HasNumericTraitsCheck, + ( Concept::HasNumericTraits ) ); + /** End concept checking */ +#endif + +protected: + LabelOverlapMeasuresImageFilter(); + ~LabelOverlapMeasuresImageFilter(){}; + void PrintSelf( std::ostream& os, Indent indent ) const; + + /** + * Pass the input through unmodified. Do this by setting the output to the + * source this by setting the output to the source image in the + * AllocateOutputs() method. + */ + void AllocateOutputs(); + + void BeforeThreadedGenerateData(); + + void AfterThreadedGenerateData(); + + /** Multi-thread version GenerateData. */ + void ThreadedGenerateData( const RegionType&, int ); + + // Override since the filter needs all the data for the algorithm + void GenerateInputRequestedRegion(); + + // Override since the filter produces all of its output + void EnlargeOutputRequestedRegion( DataObject *data ); + +private: + LabelOverlapMeasuresImageFilter( const Self& ); //purposely not implemented + void operator=( const Self& ); //purposely not implemented + + std::vector m_LabelSetMeasuresPerThread; + MapType m_LabelSetMeasures; + + SimpleFastMutexLock m_Mutex; + +}; // end of class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkLabelOverlapMeasuresImageFilter.txx" +#endif + +#endif diff --git a/itk/itkLabelOverlapMeasuresImageFilter.txx b/itk/itkLabelOverlapMeasuresImageFilter.txx new file mode 100644 index 0000000..ced277b --- /dev/null +++ b/itk/itkLabelOverlapMeasuresImageFilter.txx @@ -0,0 +1,437 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkLabelOverlapMeasuresImageFilter.txx,v $ + Language: C++ + Date: $Date: $ + Version: $Revision: $ + + Copyright (c) Insight Software Consortium. All rights reserved. + See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef _itkLabelOverlapMeasuresImageFilter_txx +#define _itkLabelOverlapMeasuresImageFilter_txx + +#include "itkLabelOverlapMeasuresImageFilter.h" + +#include "itkImageRegionConstIterator.h" +#include "itkProgressReporter.h" + +namespace itk { + +#if defined(__GNUC__) && (__GNUC__ <= 2) //NOTE: This class needs a mutex for gnu 2.95 +/** Used for mutex locking */ +#define LOCK_HASHMAP this->m_Mutex.Lock() +#define UNLOCK_HASHMAP this->m_Mutex.Unlock() +#else +#define LOCK_HASHMAP +#define UNLOCK_HASHMAP +#endif + +template +LabelOverlapMeasuresImageFilter +::LabelOverlapMeasuresImageFilter() +{ + // this filter requires two input images + this->SetNumberOfRequiredInputs( 2 ); +} + +template +void +LabelOverlapMeasuresImageFilter +::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); + if( this->GetSourceImage() ) + { + LabelImagePointer source = const_cast + ( this->GetSourceImage() ); + source->SetRequestedRegionToLargestPossibleRegion(); + } + if( this->GetTargetImage() ) + { + LabelImagePointer target = const_cast + ( this->GetTargetImage() ); + target->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template +void +LabelOverlapMeasuresImageFilter +::EnlargeOutputRequestedRegion( DataObject *data ) +{ + Superclass::EnlargeOutputRequestedRegion( data ); + data->SetRequestedRegionToLargestPossibleRegion(); +} + + +template +void +LabelOverlapMeasuresImageFilter +::AllocateOutputs() +{ + // Pass the source through as the output + LabelImagePointer image = + const_cast( this->GetSourceImage() ); + this->SetNthOutput( 0, image ); + + // Nothing that needs to be allocated for the remaining outputs +} + +template +void +LabelOverlapMeasuresImageFilter +::BeforeThreadedGenerateData() +{ + int numberOfThreads = this->GetNumberOfThreads(); + + // Resize the thread temporaries + this->m_LabelSetMeasuresPerThread.resize( numberOfThreads ); + + // Initialize the temporaries + for( int n = 0; n < numberOfThreads; n++ ) + { + this->m_LabelSetMeasuresPerThread[n].clear(); + } + + // Initialize the final map + this->m_LabelSetMeasures.clear(); +} + +template +void +LabelOverlapMeasuresImageFilter +::AfterThreadedGenerateData() +{ + // Run through the map for each thread and accumulate the set measures. + for( int n = 0; n < this->GetNumberOfThreads(); n++ ) + { + // iterate over the map for this thread + for( MapConstIterator threadIt = this->m_LabelSetMeasuresPerThread[n].begin(); + threadIt != this->m_LabelSetMeasuresPerThread[n].end(); + ++threadIt ) + { + // does this label exist in the cumulative stucture yet? + MapIterator mapIt = this->m_LabelSetMeasures.find( ( *threadIt ).first ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + // create a new entry + typedef typename MapType::value_type MapValueType; + mapIt = this->m_LabelSetMeasures.insert( MapValueType( + (*threadIt).first, LabelSetMeasures() ) ).first; + } + + // accumulate the information from this thread + (*mapIt).second.m_Source += (*threadIt).second.m_Source; + (*mapIt).second.m_Target += (*threadIt).second.m_Target; + (*mapIt).second.m_Union += (*threadIt).second.m_Union; + (*mapIt).second.m_Intersection += + (*threadIt).second.m_Intersection; + (*mapIt).second.m_SourceComplement += + (*threadIt).second.m_SourceComplement; + (*mapIt).second.m_TargetComplement += + (*threadIt).second.m_TargetComplement; + } // end of thread map iterator loop + } // end of thread loop +} + +template +void +LabelOverlapMeasuresImageFilter +::ThreadedGenerateData( const RegionType& outputRegionForThread, + int threadId ) +{ + ImageRegionConstIterator ItS( this->GetSourceImage(), + outputRegionForThread ); + ImageRegionConstIterator ItT( this->GetTargetImage(), + outputRegionForThread ); + + // support progress methods/callbacks + ProgressReporter progress( this, threadId, + 2*outputRegionForThread.GetNumberOfPixels() ); + + for( ItS.GoToBegin(), ItT.GoToBegin(); !ItS.IsAtEnd(); ++ItS, ++ItT ) + { + LabelType sourceLabel = ItS.Get(); + LabelType targetLabel = ItT.Get(); + + // is the label already in this thread? + MapIterator mapItS = + this->m_LabelSetMeasuresPerThread[threadId].find( sourceLabel ); + MapIterator mapItT = + this->m_LabelSetMeasuresPerThread[threadId].find( targetLabel ); + + if( mapItS == this->m_LabelSetMeasuresPerThread[threadId].end() ) + { + // create a new label set measures object + typedef typename MapType::value_type MapValueType; + LOCK_HASHMAP; + mapItS = this->m_LabelSetMeasuresPerThread[threadId].insert( + MapValueType( sourceLabel, LabelSetMeasures() ) ).first; + UNLOCK_HASHMAP; + } + + if( mapItT == this->m_LabelSetMeasuresPerThread[threadId].end() ) + { + // create a new label set measures object + typedef typename MapType::value_type MapValueType; + LOCK_HASHMAP; + mapItT = this->m_LabelSetMeasuresPerThread[threadId].insert( + MapValueType( targetLabel, LabelSetMeasures() ) ).first; + UNLOCK_HASHMAP; + } + + (*mapItS).second.m_Source++; + (*mapItT).second.m_Target++; + + if( sourceLabel == targetLabel ) + { + (*mapItS).second.m_Intersection++; + (*mapItS).second.m_Union++; + } + else + { + (*mapItS).second.m_Union++; + (*mapItT).second.m_Union++; + + (*mapItS).second.m_SourceComplement++; + (*mapItT).second.m_TargetComplement++; + } + + progress.CompletedPixel(); + } +} + +/** + * measures + */ +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetTotalOverlap() +{ + RealType numerator = 0.0; + RealType denominator = 0.0; + for( MapIterator mapIt = this->m_LabelSetMeasures.begin(); + mapIt != this->m_LabelSetMeasures.end(); ++mapIt ) + { + // Do not include the background in the final value. + if( (*mapIt).first == NumericTraits::Zero ) + { + continue; + } + numerator += static_cast( (*mapIt).second.m_Intersection ); + denominator += static_cast( (*mapIt).second.m_Target ); + } + return ( numerator / denominator ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetTargetOverlap( LabelType label ) +{ + MapIterator mapIt = this->m_LabelSetMeasures.find( label ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + itkWarningMacro( "Label " << label << " not found." ); + return 0.0; + } + RealType value = + static_cast( (*mapIt).second.m_Intersection ) / + static_cast( (*mapIt).second.m_Target ); + return value; +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetUnionOverlap() +{ + RealType numerator = 0.0; + RealType denominator = 0.0; + for( MapIterator mapIt = this->m_LabelSetMeasures.begin(); + mapIt != this->m_LabelSetMeasures.end(); ++mapIt ) + { + // Do not include the background in the final value. + if( (*mapIt).first == NumericTraits::Zero ) + { + continue; + } + numerator += static_cast( (*mapIt).second.m_Intersection ); + denominator += static_cast( (*mapIt).second.m_Union ); + } + return ( numerator / denominator ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetUnionOverlap( LabelType label ) +{ + MapIterator mapIt = this->m_LabelSetMeasures.find( label ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + itkWarningMacro( "Label " << label << " not found." ); + return 0.0; + } + RealType value = + static_cast( (*mapIt).second.m_Intersection ) / + static_cast( (*mapIt).second.m_Union ); + return value; +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetMeanOverlap() +{ + RealType uo = this->GetUnionOverlap(); + return ( 2.0 * uo / ( 1.0 + uo ) ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetMeanOverlap( LabelType label ) +{ + RealType uo = this->GetUnionOverlap( label ); + return ( 2.0 * uo / ( 1.0 + uo ) ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetVolumeSimilarity() +{ + RealType numerator = 0.0; + RealType denominator = 0.0; + for( MapIterator mapIt = this->m_LabelSetMeasures.begin(); + mapIt != this->m_LabelSetMeasures.end(); ++mapIt ) + { + // Do not include the background in the final value. + if( (*mapIt).first == NumericTraits::Zero ) + { + continue; + } + numerator += ( ( static_cast( (*mapIt).second.m_Source ) - + static_cast( (*mapIt).second.m_Target ) ) ); + denominator += ( ( static_cast( (*mapIt).second.m_Source ) + + static_cast( (*mapIt).second.m_Target ) ) ); + } + return ( 2.0 * numerator / denominator ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetVolumeSimilarity( LabelType label ) +{ + MapIterator mapIt = this->m_LabelSetMeasures.find( label ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + itkWarningMacro( "Label " << label << " not found." ); + return 0.0; + } + RealType value = 2.0 * + ( static_cast( (*mapIt).second.m_Source ) - + static_cast( (*mapIt).second.m_Target ) ) / + ( static_cast( (*mapIt).second.m_Source ) + + static_cast( (*mapIt).second.m_Target ) ); + return value; +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetFalseNegativeError() +{ + RealType numerator = 0.0; + RealType denominator = 0.0; + for( MapIterator mapIt = this->m_LabelSetMeasures.begin(); + mapIt != this->m_LabelSetMeasures.end(); ++mapIt ) + { + // Do not include the background in the final value. + if( (*mapIt).first == NumericTraits::Zero ) + { + continue; + } + numerator += static_cast( (*mapIt).second.m_TargetComplement ); + denominator += static_cast( (*mapIt).second.m_Target ); + } + return ( numerator / denominator ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetFalseNegativeError( LabelType label ) +{ + MapIterator mapIt = this->m_LabelSetMeasures.find( label ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + itkWarningMacro( "Label " << label << " not found." ); + return 0.0; + } + RealType value = + static_cast( (*mapIt).second.m_TargetComplement ) / + static_cast( (*mapIt).second.m_Target ); + return value; +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetFalsePositiveError() +{ + RealType numerator = 0.0; + RealType denominator = 0.0; + for( MapIterator mapIt = this->m_LabelSetMeasures.begin(); + mapIt != this->m_LabelSetMeasures.end(); ++mapIt ) + { + // Do not include the background in the final value. + if( (*mapIt).first == NumericTraits::Zero ) + { + continue; + } + numerator += static_cast( (*mapIt).second.m_SourceComplement ); + denominator += static_cast( (*mapIt).second.m_Source ); + } + return ( numerator / denominator ); +} + +template +typename LabelOverlapMeasuresImageFilter::RealType +LabelOverlapMeasuresImageFilter +::GetFalsePositiveError( LabelType label ) +{ + MapIterator mapIt = this->m_LabelSetMeasures.find( label ); + if( mapIt == this->m_LabelSetMeasures.end() ) + { + itkWarningMacro( "Label " << label << " not found." ); + return 0.0; + } + RealType value = + static_cast( (*mapIt).second.m_SourceComplement ) / + static_cast( (*mapIt).second.m_Source ); + return value; +} + +template +void +LabelOverlapMeasuresImageFilter +::PrintSelf( std::ostream& os, Indent indent ) const +{ + Superclass::PrintSelf( os, indent ); + +} + + +}// end namespace itk +#endif -- 2.45.1