itk/itkLabelOverlapMeasuresImageFilter.txx

   1 /*=========================================================================
   2
   3   Program:   Insight Segmentation & Registration Toolkit
   4   Module:    $RCSfile: itkLabelOverlapMeasuresImageFilter.txx,v $
   5   Language:  C++
   6   Date:      $Date: $
   7   Version:   $Revision: $
   8
   9   Copyright (c) Insight Software Consortium. All rights reserved.
  10   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
  11
  12      This software is distributed WITHOUT ANY WARRANTY; without even
  13      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14      PURPOSE.  See the above copyright notices for more information.
  15
  16 =========================================================================*/
  17 #ifndef _itkLabelOverlapMeasuresImageFilter_txx
  18 #define _itkLabelOverlapMeasuresImageFilter_txx
  19
  20 #include "itkLabelOverlapMeasuresImageFilter.h"
  21
  22 #include "itkImageRegionConstIterator.h"
  23 #include "itkProgressReporter.h"
  24
  25 namespace itk {
  26
  27 #if defined(__GNUC__) && (__GNUC__ <= 2) //NOTE: This class needs a mutex for gnu 2.95
  28 /** Used for mutex locking */
  29 #define LOCK_HASHMAP this->m_Mutex.Lock()
  30 #define UNLOCK_HASHMAP this->m_Mutex.Unlock()
  31 #else
  32 #define LOCK_HASHMAP
  33 #define UNLOCK_HASHMAP
  34 #endif
  35
  36 template<class TLabelImage>
  37 LabelOverlapMeasuresImageFilter<TLabelImage>
  38 ::LabelOverlapMeasuresImageFilter()
  39 {
  40   // this filter requires two input images
  41   this->SetNumberOfRequiredInputs( 2 );
  42 }
  43
  44 template<class TLabelImage>
  45 void
  46 LabelOverlapMeasuresImageFilter<TLabelImage>
  47 ::GenerateInputRequestedRegion()
  48 {
  49   Superclass::GenerateInputRequestedRegion();
  50   if( this->GetSourceImage() )
  51     {
  52     LabelImagePointer source = const_cast
  53       <LabelImageType *>( this->GetSourceImage() );
  54     source->SetRequestedRegionToLargestPossibleRegion();
  55     }
  56   if( this->GetTargetImage() )
  57     {
  58     LabelImagePointer target = const_cast
  59       <LabelImageType *>( this->GetTargetImage() );
  60     target->SetRequestedRegionToLargestPossibleRegion();
  61     }
  62 }
  63
  64 template<class TLabelImage>
  65 void
  66 LabelOverlapMeasuresImageFilter<TLabelImage>
  67 ::EnlargeOutputRequestedRegion( DataObject *data )
  68 {
  69   Superclass::EnlargeOutputRequestedRegion( data );
  70   data->SetRequestedRegionToLargestPossibleRegion();
  71 }
  72
  73
  74 template<class TLabelImage>
  75 void
  76 LabelOverlapMeasuresImageFilter<TLabelImage>
  77 ::AllocateOutputs()
  78 {
  79   // Pass the source through as the output
  80   LabelImagePointer image =
  81     const_cast<TLabelImage *>( this->GetSourceImage() );
  82   this->SetNthOutput( 0, image );
  83
  84   // Nothing that needs to be allocated for the remaining outputs
  85 }
  86
  87 template<class TLabelImage>
  88 void
  89 LabelOverlapMeasuresImageFilter<TLabelImage>
  90 ::BeforeThreadedGenerateData()
  91 {
  92   int numberOfThreads = this->GetNumberOfThreads();
  93
  94   // Resize the thread temporaries
  95   this->m_LabelSetMeasuresPerThread.resize( numberOfThreads );
  96
  97   // Initialize the temporaries
  98   for( int n = 0; n < numberOfThreads; n++ )
  99     {
 100     this->m_LabelSetMeasuresPerThread[n].clear();
 101     }
 102
 103   // Initialize the final map
 104   this->m_LabelSetMeasures.clear();
 105 }
 106
 107 template<class TLabelImage>
 108 void
 109 LabelOverlapMeasuresImageFilter<TLabelImage>
 110 ::AfterThreadedGenerateData()
 111 {
 112   // Run through the map for each thread and accumulate the set measures.
 113   for( int n = 0; n < this->GetNumberOfThreads(); n++ )
 114     {
 115     // iterate over the map for this thread
 116     for( MapConstIterator threadIt = this->m_LabelSetMeasuresPerThread[n].begin();
 117       threadIt != this->m_LabelSetMeasuresPerThread[n].end();
 118       ++threadIt )
 119       {
 120       // does this label exist in the cumulative stucture yet?
 121       MapIterator mapIt = this->m_LabelSetMeasures.find( ( *threadIt ).first );
 122       if( mapIt == this->m_LabelSetMeasures.end() )
 123         {
 124         // create a new entry
 125         typedef typename MapType::value_type MapValueType;
 126         mapIt = this->m_LabelSetMeasures.insert( MapValueType(
 127           (*threadIt).first, LabelSetMeasures() ) ).first;
 128         }
 129
 130       // accumulate the information from this thread
 131       (*mapIt).second.m_Source += (*threadIt).second.m_Source;
 132       (*mapIt).second.m_Target += (*threadIt).second.m_Target;
 133       (*mapIt).second.m_Union += (*threadIt).second.m_Union;
 134       (*mapIt).second.m_Intersection +=
 135         (*threadIt).second.m_Intersection;
 136       (*mapIt).second.m_SourceComplement +=
 137         (*threadIt).second.m_SourceComplement;
 138       (*mapIt).second.m_TargetComplement +=
 139         (*threadIt).second.m_TargetComplement;
 140       } // end of thread map iterator loop
 141     } // end of thread loop
 142 }
 143
 144 template<class TLabelImage>
 145 void
 146 LabelOverlapMeasuresImageFilter<TLabelImage>
 147 ::ThreadedGenerateData( const RegionType& outputRegionForThread,
 148   int threadId )
 149 {
 150   ImageRegionConstIterator<LabelImageType> ItS( this->GetSourceImage(),
 151     outputRegionForThread );
 152   ImageRegionConstIterator<LabelImageType> ItT( this->GetTargetImage(),
 153     outputRegionForThread );
 154
 155   // support progress methods/callbacks
 156   ProgressReporter progress( this, threadId,
 157     2*outputRegionForThread.GetNumberOfPixels() );
 158
 159   for( ItS.GoToBegin(), ItT.GoToBegin(); !ItS.IsAtEnd(); ++ItS, ++ItT )
 160     {
 161     LabelType sourceLabel = ItS.Get();
 162     LabelType targetLabel = ItT.Get();
 163
 164     // is the label already in this thread?
 165     MapIterator mapItS =
 166       this->m_LabelSetMeasuresPerThread[threadId].find( sourceLabel );
 167     MapIterator mapItT =
 168       this->m_LabelSetMeasuresPerThread[threadId].find( targetLabel );
 169
 170     if( mapItS == this->m_LabelSetMeasuresPerThread[threadId].end() )
 171       {
 172       // create a new label set measures object
 173       typedef typename MapType::value_type MapValueType;
 174       LOCK_HASHMAP;
 175       mapItS = this->m_LabelSetMeasuresPerThread[threadId].insert(
 176         MapValueType( sourceLabel, LabelSetMeasures() ) ).first;
 177       UNLOCK_HASHMAP;
 178       }
 179
 180     if( mapItT == this->m_LabelSetMeasuresPerThread[threadId].end() )
 181       {
 182       // create a new label set measures object
 183       typedef typename MapType::value_type MapValueType;
 184       LOCK_HASHMAP;
 185       mapItT = this->m_LabelSetMeasuresPerThread[threadId].insert(
 186         MapValueType( targetLabel, LabelSetMeasures() ) ).first;
 187       UNLOCK_HASHMAP;
 188       }
 189
 190     (*mapItS).second.m_Source++;
 191     (*mapItT).second.m_Target++;
 192
 193     if( sourceLabel == targetLabel )
 194       {
 195       (*mapItS).second.m_Intersection++;
 196       (*mapItS).second.m_Union++;
 197       }
 198     else
 199       {
 200       (*mapItS).second.m_Union++;
 201       (*mapItT).second.m_Union++;
 202
 203       (*mapItS).second.m_SourceComplement++;
 204       (*mapItT).second.m_TargetComplement++;
 205       }
 206
 207     progress.CompletedPixel();
 208     }
 209 }
 210
 211 /**
 212  *  measures
 213  */
 214 template<class TLabelImage>
 215 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 216 LabelOverlapMeasuresImageFilter<TLabelImage>
 217 ::GetTotalOverlap()
 218 {
 219   RealType numerator = 0.0;
 220   RealType denominator = 0.0;
 221   for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
 222     mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
 223     {
 224     // Do not include the background in the final value.
 225     if( (*mapIt).first == NumericTraits<LabelType>::Zero )
 226       {
 227       continue;
 228       }
 229     numerator += static_cast<RealType>( (*mapIt).second.m_Intersection );
 230     denominator += static_cast<RealType>( (*mapIt).second.m_Target );
 231     }
 232   return ( numerator / denominator );
 233 }
 234
 235 template<class TLabelImage>
 236 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 237 LabelOverlapMeasuresImageFilter<TLabelImage>
 238 ::GetTargetOverlap( LabelType label )
 239 {
 240   MapIterator mapIt = this->m_LabelSetMeasures.find( label );
 241   if( mapIt == this->m_LabelSetMeasures.end() )
 242     {
 243     itkWarningMacro( "Label " << label << " not found." );
 244     return 0.0;
 245     }
 246   RealType value =
 247     static_cast<RealType>( (*mapIt).second.m_Intersection ) /
 248     static_cast<RealType>( (*mapIt).second.m_Target );
 249   return value;
 250 }
 251
 252 template<class TLabelImage>
 253 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 254 LabelOverlapMeasuresImageFilter<TLabelImage>
 255 ::GetUnionOverlap()
 256 {
 257   RealType numerator = 0.0;
 258   RealType denominator = 0.0;
 259   for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
 260     mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
 261     {
 262     // Do not include the background in the final value.
 263     if( (*mapIt).first == NumericTraits<LabelType>::Zero )
 264       {
 265       continue;
 266       }
 267     numerator += static_cast<RealType>( (*mapIt).second.m_Intersection );
 268     denominator += static_cast<RealType>( (*mapIt).second.m_Union );
 269     }
 270   return ( numerator / denominator );
 271 }
 272
 273 template<class TLabelImage>
 274 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 275 LabelOverlapMeasuresImageFilter<TLabelImage>
 276 ::GetUnionOverlap( LabelType label )
 277 {
 278   MapIterator mapIt = this->m_LabelSetMeasures.find( label );
 279   if( mapIt == this->m_LabelSetMeasures.end() )
 280     {
 281     itkWarningMacro( "Label " << label << " not found." );
 282     return 0.0;
 283     }
 284   RealType value =
 285     static_cast<RealType>( (*mapIt).second.m_Intersection ) /
 286     static_cast<RealType>( (*mapIt).second.m_Union );
 287   return value;
 288 }
 289
 290 template<class TLabelImage>
 291 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 292 LabelOverlapMeasuresImageFilter<TLabelImage>
 293 ::GetMeanOverlap()
 294 {
 295   RealType uo = this->GetUnionOverlap();
 296   return ( 2.0 * uo / ( 1.0 + uo ) );
 297 }
 298
 299 template<class TLabelImage>
 300 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 301 LabelOverlapMeasuresImageFilter<TLabelImage>
 302 ::GetMeanOverlap( LabelType label )
 303 {
 304   RealType uo = this->GetUnionOverlap( label );
 305   return ( 2.0 * uo / ( 1.0 + uo ) );
 306 }
 307
 308 template<class TLabelImage>
 309 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 310 LabelOverlapMeasuresImageFilter<TLabelImage>
 311 ::GetVolumeSimilarity()
 312 {
 313   RealType numerator = 0.0;
 314   RealType denominator = 0.0;
 315   for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
 316     mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
 317     {
 318     // Do not include the background in the final value.
 319     if( (*mapIt).first == NumericTraits<LabelType>::Zero )
 320       {
 321       continue;
 322       }
 323     numerator += ( ( static_cast<RealType>( (*mapIt).second.m_Source ) -
 324       static_cast<RealType>( (*mapIt).second.m_Target ) ) );
 325     denominator += ( ( static_cast<RealType>( (*mapIt).second.m_Source ) +
 326       static_cast<RealType>( (*mapIt).second.m_Target ) ) );
 327     }
 328   return ( 2.0 * numerator / denominator );
 329 }
 330
 331 template<class TLabelImage>
 332 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 333 LabelOverlapMeasuresImageFilter<TLabelImage>
 334 ::GetVolumeSimilarity( LabelType label )
 335 {
 336   MapIterator mapIt = this->m_LabelSetMeasures.find( label );
 337   if( mapIt == this->m_LabelSetMeasures.end() )
 338     {
 339     itkWarningMacro( "Label " << label << " not found." );
 340     return 0.0;
 341     }
 342   RealType value = 2.0 *
 343     ( static_cast<RealType>( (*mapIt).second.m_Source ) -
 344       static_cast<RealType>( (*mapIt).second.m_Target ) ) /
 345     ( static_cast<RealType>( (*mapIt).second.m_Source ) +
 346       static_cast<RealType>( (*mapIt).second.m_Target ) );
 347   return value;
 348 }
 349
 350 template<class TLabelImage>
 351 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 352 LabelOverlapMeasuresImageFilter<TLabelImage>
 353 ::GetFalseNegativeError()
 354 {
 355   RealType numerator = 0.0;
 356   RealType denominator = 0.0;
 357   for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
 358     mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
 359     {
 360     // Do not include the background in the final value.
 361     if( (*mapIt).first == NumericTraits<LabelType>::Zero )
 362       {
 363       continue;
 364       }
 365     numerator += static_cast<RealType>( (*mapIt).second.m_TargetComplement );
 366     denominator += static_cast<RealType>( (*mapIt).second.m_Target );
 367     }
 368   return ( numerator / denominator );
 369 }
 370
 371 template<class TLabelImage>
 372 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 373 LabelOverlapMeasuresImageFilter<TLabelImage>
 374 ::GetFalseNegativeError( LabelType label )
 375 {
 376   MapIterator mapIt = this->m_LabelSetMeasures.find( label );
 377   if( mapIt == this->m_LabelSetMeasures.end() )
 378     {
 379     itkWarningMacro( "Label " << label << " not found." );
 380     return 0.0;
 381     }
 382   RealType value =
 383     static_cast<RealType>( (*mapIt).second.m_TargetComplement ) /
 384     static_cast<RealType>( (*mapIt).second.m_Target );
 385   return value;
 386 }
 387
 388 template<class TLabelImage>
 389 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 390 LabelOverlapMeasuresImageFilter<TLabelImage>
 391 ::GetFalsePositiveError()
 392 {
 393   RealType numerator = 0.0;
 394   RealType denominator = 0.0;
 395   for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
 396     mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
 397     {
 398     // Do not include the background in the final value.
 399     if( (*mapIt).first == NumericTraits<LabelType>::Zero )
 400       {
 401       continue;
 402       }
 403     numerator += static_cast<RealType>( (*mapIt).second.m_SourceComplement );
 404     denominator += static_cast<RealType>( (*mapIt).second.m_Source );
 405     }
 406   return ( numerator / denominator );
 407 }
 408
 409 template<class TLabelImage>
 410 typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
 411 LabelOverlapMeasuresImageFilter<TLabelImage>
 412 ::GetFalsePositiveError( LabelType label )
 413 {
 414   MapIterator mapIt = this->m_LabelSetMeasures.find( label );
 415   if( mapIt == this->m_LabelSetMeasures.end() )
 416     {
 417     itkWarningMacro( "Label " << label << " not found." );
 418     return 0.0;
 419     }
 420   RealType value =
 421     static_cast<RealType>( (*mapIt).second.m_SourceComplement ) /
 422     static_cast<RealType>( (*mapIt).second.m_Source );
 423   return value;
 424 }
 425
 426 template<class TLabelImage>
 427 void
 428 LabelOverlapMeasuresImageFilter<TLabelImage>
 429 ::PrintSelf( std::ostream& os, Indent indent ) const
 430 {
 431   Superclass::PrintSelf( os, indent );
 432
 433 }
 434
 435
 436 }// end namespace itk
 437 #endif