itk/clitkRelativePositionAnalyzerFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ===========================================================================**/
  18
  19 //--------------------------------------------------------------------
  20 template <class ImageType>
  21 clitk::RelativePositionAnalyzerFilter<ImageType>::
  22 RelativePositionAnalyzerFilter():
  23   // clitk::FilterBase(),
  24   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  25   itk::ImageToImageFilter<ImageType, ImageType>()
  26 {
  27   this->SetNumberOfRequiredInputs(3); // support, object, target
  28   VerboseFlagOff();
  29   SetBackgroundValue(0);
  30   SetForegroundValue(1);
  31   SetNumberOfBins(100);
  32   SetNumberOfAngles(4);
  33   SetAreaLossTolerance(0.01);
  34   m_ListOfAngles.clear();
  35   SetSupportSize(0);
  36   SetTargetSize(0);
  37   SetSizeWithThreshold(0);
  38   SetSizeWithReverseThreshold(0);
  39 }
  40 //--------------------------------------------------------------------
  41
  42
  43 //--------------------------------------------------------------------
  44 template <class ImageType>
  45 void
  46 clitk::RelativePositionAnalyzerFilter<ImageType>::
  47 SetInputSupport(const ImageType * image)
  48 {
  49   // Process object is not const-correct so the const casting is required.
  50   this->SetNthInput(0, const_cast<ImageType *>(image));
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 template <class ImageType>
  57 void
  58 clitk::RelativePositionAnalyzerFilter<ImageType>::
  59 SetInputObject(const ImageType * image)
  60 {
  61   // Process object is not const-correct so the const casting is required.
  62   this->SetNthInput(1, const_cast<ImageType *>(image));
  63 }
  64 //--------------------------------------------------------------------
  65
  66
  67 //--------------------------------------------------------------------
  68 template <class ImageType>
  69 void
  70 clitk::RelativePositionAnalyzerFilter<ImageType>::
  71 SetInputTarget(const ImageType * image)
  72 {
  73   // Process object is not const-correct so the const casting is required.
  74   this->SetNthInput(2, const_cast<ImageType *>(image));
  75 }
  76 //--------------------------------------------------------------------
  77
  78
  79 //--------------------------------------------------------------------
  80 template <class ImageType>
  81 void
  82 clitk::RelativePositionAnalyzerFilter<ImageType>::
  83 PrintOptions()
  84 {
  85   DD("TODO");
  86 }
  87 //--------------------------------------------------------------------
  88
  89
  90 //--------------------------------------------------------------------
  91 template <class ImageType>
  92 void
  93 clitk::RelativePositionAnalyzerFilter<ImageType>::
  94 GenerateOutputInformation()
  95 {
  96   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  97   ImagePointer outputImage = this->GetOutput(0);
  98   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
  99 }
 100 //--------------------------------------------------------------------
 101
 102
 103 //--------------------------------------------------------------------
 104 template <class ImageType>
 105 void
 106 clitk::RelativePositionAnalyzerFilter<ImageType>::
 107 GenerateData()
 108 {
 109   this->LoadAFDB();
 110
 111   // Get input pointer
 112   m_Support = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 113   m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 114   m_Target = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
 115   static const unsigned int dim = ImageType::ImageDimension;
 116
 117   // Remove object from support
 118   clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
 119
 120   // Resize object like target (to enable substraction later)
 121   ImagePointer objectLikeTarget = clitk::ResizeImageLike<ImageType>(m_Object, m_Target, GetBackgroundValue());
 122
 123   // Define filter to compute statics on mask image
 124   typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
 125   typename StatFilterType::Pointer statFilter = StatFilterType::New();
 126
 127   // Compute the initial support size
 128   statFilter->SetInput(m_Support);
 129   statFilter->SetLabelInput(m_Support);
 130   statFilter->Update();
 131   SetSupportSize(statFilter->GetCount(GetForegroundValue()));
 132   // DD(GetSupportSize());
 133
 134   // Compute the initial target size
 135   ImagePointer s = clitk::ResizeImageLike<ImageType>(m_Support, m_Target, GetBackgroundValue());
 136   statFilter->SetInput(s);
 137   statFilter->SetLabelInput(m_Target);
 138   statFilter->Update();
 139   SetTargetSize(statFilter->GetCount(GetForegroundValue()));
 140   // DD(GetTargetSize());
 141
 142   // Build the list of tested orientations
 143   m_ListOfAngles.clear();
 144   for(uint i=0; i<GetNumberOfAngles(); i++) {
 145     double a = i*360.0/GetNumberOfAngles();
 146     if (a>180) a = 180-a;
 147     m_ListOfAngles.push_back(clitk::deg2rad(a));
 148     RelativePositionOrientationType r;
 149     r.angle1 = clitk::deg2rad(a);
 150     r.angle2 = 0;
 151     r.notFlag = false;
 152     m_ListOfOrientation.push_back(r);
 153     r.notFlag = true;
 154     m_ListOfOrientation.push_back(r);
 155   }
 156
 157   // Loop on all orientations
 158   int bins = GetNumberOfBins();
 159   double tolerance = GetAreaLossTolerance();
 160   for(int i=0; i<m_ListOfAngles.size(); i++) {
 161     // Compute Fuzzy map
 162     typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]);
 163     writeImage<FloatImageType>(map, "fuzzy_"+toString(i)+".mha");
 164
 165     // Compute the optimal thresholds (direct and inverse)
 166     double mThreshold=0.0;
 167     double mReverseThreshold=1.0;
 168     ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold);
 169
 170     // Use the threshold to compute new support
 171     int s1 = GetSupportSize();
 172     // DD(mThreshold);
 173     // DD(mReverseThreshold);
 174     if (mThreshold > 0.0) {
 175       ImagePointer support1 =
 176         clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 177                                                        mThreshold,
 178                                                        m_ListOfAngles[i],false,
 179                                                        false, -1, true, false);
 180       // writeImage<ImageType>(support1, "sup_"+toString(i)+".mha");
 181       // Compute the new support size
 182       statFilter->SetInput(support1);
 183       statFilter->SetLabelInput(support1);
 184       statFilter->Update();
 185       s1 = statFilter->GetCount(GetForegroundValue());
 186     }
 187
 188     int s2 = GetSupportSize();
 189     if (mReverseThreshold < 1.0) {
 190       // DD(m_ListOfAngles[1]);
 191       ImagePointer support2 =
 192         clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 193                                                        mReverseThreshold,
 194                                                        m_ListOfAngles[i],true,
 195                                                        false, -1, true, false);
 196       writeImage<ImageType>(support2, "sup_rev_"+toString(i)+".mha");
 197       // Compute the new support size
 198       statFilter = StatFilterType::New();
 199       statFilter->SetInput(support2);
 200       statFilter->SetLabelInput(support2);
 201       statFilter->Update();
 202       s2 = statFilter->GetCount(GetForegroundValue());
 203     }
 204
 205     // Set results values
 206     RelativePositionInformationType r;
 207     r.threshold = mThreshold;
 208     r.sizeAfterThreshold = s1; // DD(s1);
 209     r.sizeBeforeThreshold = GetSupportSize();
 210     r.sizeReference = GetTargetSize();
 211     m_ListOfInformation.push_back(r);
 212
 213     r.threshold = mReverseThreshold;
 214     r.sizeAfterThreshold = s2;     // DD(s2);
 215     m_ListOfInformation.push_back(r);
 216     // Print();
 217   } // end loop on orientations
 218 }
 219 //--------------------------------------------------------------------
 220
 221
 222 //--------------------------------------------------------------------
 223 template <class ImageType>
 224 typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
 225 clitk::RelativePositionAnalyzerFilter<ImageType>::
 226 ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
 227 {
 228   typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
 229   typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
 230   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 231   sliceRelPosFilter->VerboseStepFlagOff();
 232   sliceRelPosFilter->WriteStepFlagOff();
 233   sliceRelPosFilter->SetInput(target);
 234   sliceRelPosFilter->SetInputObject(object);
 235   sliceRelPosFilter->SetDirection(2);
 236   sliceRelPosFilter->SetIntermediateSpacingFlag(false);
 237   //sliceRelPosFilter->AddOrientationTypeString(orientation);
 238   sliceRelPosFilter->AddAngles(angle, 0.0);
 239   sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
 240   // sliceRelPosFilter->PrintOptions();
 241   sliceRelPosFilter->Update();
 242   typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
 243
 244   // Remove initial object from the fuzzy map
 245   map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
 246
 247   // Resize the fuzzy map like the target, put 2.0 when outside
 248   map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0);  // Put 2.0 when out of initial map
 249
 250   // end
 251   return map;
 252 }
 253 //--------------------------------------------------------------------
 254
 255
 256 //--------------------------------------------------------------------
 257 template <class ImageType>
 258 void
 259 clitk::RelativePositionAnalyzerFilter<ImageType>::
 260 ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
 261                          double & threshold, double & reverseThreshold)
 262 {
 263   // Get the histogram of fuzzy values inside the target image
 264   typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
 265   typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
 266   f->SetInput(map);
 267   f->SetLabelInput(target);
 268   f->UseHistogramsOn();
 269   f->SetHistogramParameters(bins, 0.0, 1.1);
 270   f->Update();
 271   int count = f->GetCount(GetForegroundValue());
 272   typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
 273
 274   // Debug : dump histogram
 275   static int i=0;
 276   std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
 277   for(int j=0; j<bins; j++) {
 278     histogramFile << h->GetMeasurement(j,0)
 279                   << "\t" << h->GetFrequency(j)
 280                   << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
 281   }
 282   histogramFile.close();
 283   i++;
 284
 285   // Analyze the histogram (direct)
 286   double sum = 0.0;
 287   bool found = false;
 288   threshold = 0.0;
 289   for(int j=0; j<bins; j++) {
 290     sum += ((double)h->GetFrequency(j)/(double)count);
 291     if ((!found) && (sum > tolerance)) {
 292       if (j==0) threshold = h->GetBinMin(0,j);
 293       else threshold = h->GetBinMin(0,j-1); // the last before reaching the threshold
 294       found = true;
 295     }
 296   }
 297
 298   // Analyze the histogram (reverse)
 299   sum = 0.0;
 300   found = false;
 301   reverseThreshold = 1.0;
 302   for(int j=bins-1; j>=0; j--) {
 303     sum += ((double)h->GetFrequency(j)/(double)count);
 304     if ((!found) && (sum > tolerance)) {
 305       if (j==bins-1) reverseThreshold = h->GetBinMax(0,j);
 306       else reverseThreshold = h->GetBinMax(0,j+1);
 307       found = true;
 308     }
 309   }
 310 }
 311 //--------------------------------------------------------------------
 312
 313
 314 //--------------------------------------------------------------------
 315 template <class ImageType>
 316 void
 317 clitk::RelativePositionAnalyzerFilter<ImageType>::
 318 Print(std::string s, std::ostream & os)
 319 {
 320   for(int i=0; i<m_ListOfOrientation.size(); i++) {
 321     os << s << " ";
 322     m_ListOfOrientation[i].Print(os);
 323     m_ListOfInformation[i].Println(os);
 324   }
 325 }
 326 //--------------------------------------------------------------------