mitk: mitkCalculateGrayValueStatisticsTool.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date$
00006 Version:   $Revision: 1.12 $
00007  
00008 Copyright (c) German Cancer Research Center, Division of Medical and
00009 Biological Informatics. All rights reserved.
00010 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details.
00011  
00012 This software is distributed WITHOUT ANY WARRANTY; without even
00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014 PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #include "mitkCalculateGrayValueStatisticsTool.h"
00019 
00020 #include "mitkCalculateGrayValueStatisticsTool.xpm"
00021 
00022 #include "mitkProgressBar.h"
00023 #include "mitkStatusBar.h"
00024 #include "mitkImageTimeSelector.h"
00025 #include "mitkImageCast.h"
00026 #include "mitkToolManager.h"
00027 
00028 #include <itkCommand.h>
00029 
00030 #include <itkHistogram.h>
00031 
00032 namespace mitk {
00033   MITK_TOOL_MACRO(MitkExt_EXPORT, CalculateGrayValueStatisticsTool, "Statistics tool");
00034 }
00035 
00036 mitk::CalculateGrayValueStatisticsTool::CalculateGrayValueStatisticsTool()
00037 {
00038 }
00039 
00040 mitk::CalculateGrayValueStatisticsTool::~CalculateGrayValueStatisticsTool()
00041 {
00042 }
00043 
00044 const char** mitk::CalculateGrayValueStatisticsTool::GetXPM() const
00045 {
00046   return mitkCalculateGrayValueStatisticsTool_xpm;
00047 }
00048 
00049 const char* mitk::CalculateGrayValueStatisticsTool::GetName() const
00050 {
00051   return "Statistics";
00052 }
00053 
00054 std::string mitk::CalculateGrayValueStatisticsTool::GetErrorMessage()
00055 {
00056   return "No statistics generated for these segmentations:";
00057 }
00058 
00059 void mitk::CalculateGrayValueStatisticsTool::StartProcessingAllData()
00060 {
00061   // clear/prepare report
00062   m_CompleteReport.str("");
00063 }
00064 
00065 bool mitk::CalculateGrayValueStatisticsTool::ProcessOneWorkingData( DataNode* node )
00066 {
00067   if (node)
00068   {
00069     Image::Pointer image = dynamic_cast<Image*>( node->GetData() );
00070     if (image.IsNull()) return false;
00071 
00072     DataNode* referencenode = m_ToolManager->GetReferenceData(0);
00073     if (!referencenode) return false;
00074 
00075     try
00076     {
00077       ProgressBar::GetInstance()->AddStepsToDo(1);
00078 
00079       // add to report
00080       std::string nodename("structure");
00081       node->GetName(nodename);
00082 
00083       std::string message = "Calculating statistics for ";
00084       message += nodename;
00085 
00086       StatusBar::GetInstance()->DisplayText(message.c_str());
00087 
00088       Image::Pointer refImage = dynamic_cast<Image*>( referencenode->GetData() );
00089       Image::Pointer image = dynamic_cast<Image*>( node->GetData() );
00090 
00091       m_CompleteReport << "======== Gray value analysis of " << nodename << " ========\n";
00092 
00093       if (image.IsNotNull() && refImage.IsNotNull() )
00094       {
00095         for (unsigned int timeStep = 0; timeStep < image->GetTimeSteps(); ++timeStep)
00096         {
00097           ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New();
00098           timeSelector->SetInput( refImage );
00099           timeSelector->SetTimeNr( timeStep );
00100           timeSelector->UpdateLargestPossibleRegion();
00101           Image::Pointer refImage3D = timeSelector->GetOutput();
00102 
00103           ImageTimeSelector::Pointer timeSelector2 = ImageTimeSelector::New();
00104           timeSelector2->SetInput( image );
00105           timeSelector2->SetTimeNr( timeStep );
00106           timeSelector2->UpdateLargestPossibleRegion();
00107           Image::Pointer image3D = timeSelector2->GetOutput();
00108 
00109           if (image3D.IsNotNull() && refImage3D.IsNotNull() )
00110           {
00111             m_CompleteReport << "=== " << nodename << ", time step " << timeStep << " ===\n";
00112             AccessFixedDimensionByItk_2( refImage3D, ITKHistogramming, 3, image3D, m_CompleteReport );
00113           }
00114         }
00115       }
00116 
00117       m_CompleteReport << "======== End of analysis for " << nodename << " ===========\n\n\n";
00118 
00119       ProgressBar::GetInstance()->Progress();
00120     }
00121     catch(...)
00122     {
00123       return false;
00124     }
00125   }
00126 
00127   return true;
00128 }
00129 
00130 void mitk::CalculateGrayValueStatisticsTool::FinishProcessingAllData()
00131 {
00132   SegmentationsProcessingTool::FinishProcessingAllData();
00133 
00134   // show/send results
00135   StatisticsCompleted.Send();
00136   //MITK_INFO << m_CompleteReport.str() << std::endl;
00137 }
00138 
00139 #define ROUND_P(x) ((int)((x)+0.5))
00140 
00141 template<typename TPixel, unsigned int VImageDimension>
00142 void mitk::CalculateGrayValueStatisticsTool::CalculateMinMax(
00143     itk::Image<TPixel, VImageDimension>* referenceImage, Image* segmentation, TPixel& minimum,
00144     TPixel& maximum)
00145 {
00146   typedef itk::Image<TPixel, VImageDimension> ImageType;
00147   typedef itk::Image<Tool::DefaultSegmentationDataType, VImageDimension> SegmentationType;
00148 
00149   typename SegmentationType::Pointer segmentationItk;
00150   CastToItkImage(segmentation, segmentationItk);
00151 
00152   typename SegmentationType::RegionType segmentationRegion =
00153       segmentationItk->GetLargestPossibleRegion();
00154 
00155   segmentationRegion.Crop(referenceImage->GetLargestPossibleRegion());
00156 
00157   itk::ImageRegionConstIteratorWithIndex<SegmentationType> segmentationIterator(segmentationItk,
00158       segmentationRegion);
00159   itk::ImageRegionConstIteratorWithIndex<ImageType> referenceIterator(referenceImage,
00160       segmentationRegion);
00161 
00162   segmentationIterator.GoToBegin();
00163   referenceIterator.GoToBegin();
00164 
00165   minimum = std::numeric_limits<TPixel>::max();
00166   maximum = std::numeric_limits<TPixel>::min();
00167 
00168   while (!segmentationIterator.IsAtEnd())
00169   {
00170     itk::Point<float, 3> pt;
00171     segmentationItk->TransformIndexToPhysicalPoint(segmentationIterator.GetIndex(), pt);
00172 
00173     typename ImageType::IndexType ind;
00174     itk::ContinuousIndex<float, 3> contInd;
00175     if (referenceImage->TransformPhysicalPointToContinuousIndex(pt, contInd))
00176     {
00177       for (unsigned int i = 0; i < 3; ++i)
00178         ind[i] = ROUND_P(contInd[i]);
00179 
00180       referenceIterator.SetIndex(ind);
00181 
00182       if (segmentationIterator.Get() > 0)
00183       {
00184 
00185         if (referenceIterator.Get() < minimum)
00186           minimum = referenceIterator.Get();
00187         if (referenceIterator.Get() > maximum)
00188           maximum = referenceIterator.Get();
00189       }
00190     }
00191 
00192     ++segmentationIterator;
00193   }
00194 
00195 }
00196 
00197 template<typename TPixel, unsigned int VImageDimension>
00198 void mitk::CalculateGrayValueStatisticsTool::ITKHistogramming(
00199     itk::Image<TPixel, VImageDimension>* referenceImage, Image* segmentation,
00200     std::stringstream& report)
00201 {
00202   typedef itk::Image<TPixel, VImageDimension> ImageType;
00203   typedef itk::Image<Tool::DefaultSegmentationDataType, VImageDimension> SegmentationType;
00204 
00205   typename SegmentationType::Pointer segmentationItk;
00206   CastToItkImage( segmentation, segmentationItk );
00207 
00208   // generate histogram
00209   typename SegmentationType::RegionType segmentationRegion = segmentationItk->GetLargestPossibleRegion();
00210 
00211   segmentationRegion.Crop( referenceImage->GetLargestPossibleRegion() );
00212 
00213   itk::ImageRegionConstIteratorWithIndex< SegmentationType > segmentationIterator( segmentationItk, segmentationRegion);
00214   itk::ImageRegionConstIteratorWithIndex< ImageType >        referenceIterator( referenceImage, segmentationRegion);
00215 
00216   segmentationIterator.GoToBegin();
00217   referenceIterator.GoToBegin();
00218 
00219   m_ITKHistogram = HistogramType::New();
00220 
00221   TPixel minimum = std::numeric_limits<TPixel>::max();
00222   TPixel maximum = std::numeric_limits<TPixel>::min();
00223     
00224   CalculateMinMax(referenceImage, segmentation, minimum, maximum);
00225 
00226   //initialize the histogram to the range of the cropped region
00227   HistogramType::SizeType size;
00228 #if defined(ITK_USE_REVIEW_STATISTICS)
00229   typedef typename HistogramType::SizeType::ValueType HSizeValueType;
00230 #else
00231   typedef typename HistogramType::SizeType::SizeValueType HSizeValueType;
00232 #endif
00233   size.Fill(static_cast<HSizeValueType> (maximum - minimum + 1));
00234   HistogramType::MeasurementVectorType lowerBound;
00235   HistogramType::MeasurementVectorType upperBound;
00236   lowerBound[0] = minimum;
00237   upperBound[0] = maximum;
00238   m_ITKHistogram->Initialize(size, lowerBound, upperBound);
00239 
00240   double mean(0.0);
00241   double sd(0.0);
00242   double voxelCount(0.0);
00243 
00244   //iterate through the cropped region add the values to the histogram
00245   while (!segmentationIterator.IsAtEnd())
00246   {
00247     itk::Point< float, 3 > pt;
00248     segmentationItk->TransformIndexToPhysicalPoint( segmentationIterator.GetIndex(), pt );
00249 
00250     typename ImageType::IndexType ind;
00251     itk::ContinuousIndex<float, 3> contInd;
00252     if (referenceImage->TransformPhysicalPointToContinuousIndex(pt, contInd))
00253     {
00254       for (unsigned int i = 0; i < 3; ++i) ind[i] = ROUND_P(contInd[i]);
00255 
00256       referenceIterator.SetIndex( ind );
00257 
00258       if ( segmentationIterator.Get() > 0 )
00259       {
00260         HistogramType::MeasurementVectorType currentMeasurementVector;
00261 
00262         currentMeasurementVector[0]
00263             = static_cast<HistogramType::MeasurementType> (referenceIterator.Get());
00264         m_ITKHistogram->IncreaseFrequency(currentMeasurementVector, 1);
00265 
00266         mean = (mean * (static_cast<double> (voxelCount) / static_cast<double> (voxelCount + 1))) // 3 points:   old center * 2/3 + currentPoint * 1/3;
00267             + static_cast<double> (referenceIterator.Get()) / static_cast<double> (voxelCount + 1);
00268 
00269         voxelCount += 1.0;
00270       }
00271     }
00272 
00273     ++segmentationIterator;
00274   }
00275 
00276   // second pass for SD
00277   segmentationIterator.GoToBegin();
00278   referenceIterator.GoToBegin();
00279   while ( !segmentationIterator.IsAtEnd() )
00280   {
00281     itk::Point< float, 3 > pt;
00282     segmentationItk->TransformIndexToPhysicalPoint( segmentationIterator.GetIndex(), pt );
00283 
00284     typename ImageType::IndexType ind;
00285     itk::ContinuousIndex<float, 3> contInd;
00286     if (referenceImage->TransformPhysicalPointToContinuousIndex(pt, contInd))
00287     {
00288       for (unsigned int i = 0; i < 3; ++i) ind[i] = ROUND_P(contInd[i]);
00289 
00290       referenceIterator.SetIndex( ind );
00291 
00292       if ( segmentationIterator.Get() > 0 )
00293       {
00294         sd += ((static_cast<double>(referenceIterator.Get() ) - mean)*(static_cast<double>(referenceIterator.Get() ) - mean));
00295       }
00296     }
00297 
00298     ++segmentationIterator;
00299   }
00300 
00301   sd /= static_cast<double>(voxelCount - 1);
00302   sd = sqrt( sd );
00303 
00304   // generate quantiles
00305   TPixel histogramQuantileValues[5];
00306   histogramQuantileValues[0] = m_ITKHistogram->Quantile(0, 0.05);
00307   histogramQuantileValues[1] = m_ITKHistogram->Quantile(0, 0.25);
00308   histogramQuantileValues[2] = m_ITKHistogram->Quantile(0, 0.50);
00309   histogramQuantileValues[3] = m_ITKHistogram->Quantile(0, 0.75);
00310   histogramQuantileValues[4] = m_ITKHistogram->Quantile(0, 0.95);
00311 
00312   // report histogram values
00313 report     << "         Minimum:" << minimum
00314            << "\n  5% quantile: " << histogramQuantileValues[0]
00315            << "\n 25% quantile: " << histogramQuantileValues[1]
00316            << "\n 50% quantile: " << histogramQuantileValues[2]
00317            << "\n 75% quantile: " << histogramQuantileValues[3]
00318            << "\n 95% quantile: " << histogramQuantileValues[4]
00319            << "\n      Maximum: " << maximum
00320            << "\n         Mean: " << mean
00321            << "\n           SD: " << sd
00322            << "\n";
00323 }
00324 
00325 std::string mitk::CalculateGrayValueStatisticsTool::GetReport() const
00326 {
00327   return m_CompleteReport.str();
00328 }
00329 
00330 mitk::CalculateGrayValueStatisticsTool::HistogramType::ConstPointer mitk::CalculateGrayValueStatisticsTool::GetHistogram()
00331 {
00332   return m_ITKHistogram.GetPointer();
00333 }
00334