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   itk::ImageToImageFilter<ImageType, ImageType>()
  24 {
  25   this->SetNumberOfRequiredInputs(3); // Input : support, object, target
  26   SetBackgroundValue(0);
  27   SetForegroundValue(1);
  28   SetNumberOfBins(100);
  29   SetAreaLossTolerance(0.01);
  30   SetSupportSize(0);
  31   SetTargetSize(0);
  32   SetSizeWithThreshold(0);
  33   SetSizeWithReverseThreshold(0);
  34 }
  35 //--------------------------------------------------------------------
  36
  37
  38 //--------------------------------------------------------------------
  39 template <class ImageType>
  40 void
  41 clitk::RelativePositionAnalyzerFilter<ImageType>::
  42 SetInputSupport(const ImageType * image)
  43 {
  44   // Process object is not const-correct so the const casting is required.
  45   this->SetNthInput(0, const_cast<ImageType *>(image));
  46 }
  47 //--------------------------------------------------------------------
  48
  49
  50 //--------------------------------------------------------------------
  51 template <class ImageType>
  52 void
  53 clitk::RelativePositionAnalyzerFilter<ImageType>::
  54 SetInputObject(const ImageType * image)
  55 {
  56   // Process object is not const-correct so the const casting is required.
  57   this->SetNthInput(1, const_cast<ImageType *>(image));
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 template <class ImageType>
  64 void
  65 clitk::RelativePositionAnalyzerFilter<ImageType>::
  66 SetInputTarget(const ImageType * image)
  67 {
  68   // Process object is not const-correct so the const casting is required.
  69   this->SetNthInput(2, const_cast<ImageType *>(image));
  70 }
  71 //--------------------------------------------------------------------
  72
  73
  74 //--------------------------------------------------------------------
  75 template <class ImageType>
  76 void
  77 clitk::RelativePositionAnalyzerFilter<ImageType>::
  78 PrintOptions()
  79 {
  80   DD("TODO");
  81 }
  82 //--------------------------------------------------------------------
  83
  84
  85 //--------------------------------------------------------------------
  86 template <class ImageType>
  87 void
  88 clitk::RelativePositionAnalyzerFilter<ImageType>::
  89 GenerateOutputInformation()
  90 {
  91   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  92   ImagePointer outputImage = this->GetOutput(0);
  93   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
  94 }
  95 //--------------------------------------------------------------------
  96
  97
  98 //--------------------------------------------------------------------
  99 template <class ImageType>
 100 void
 101 clitk::RelativePositionAnalyzerFilter<ImageType>::
 102 GenerateData()
 103 {
 104   static const unsigned int dim = ImageType::ImageDimension;
 105
 106   ImagePointer temp = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 107   m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 108   ImagePointer temp2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
 109
 110   // Remove object from support (keep initial image)
 111   m_Support = clitk::Clone<ImageType>(temp);
 112   clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
 113
 114   // Remove object from target. Important because sometimes, there is
 115   // overlap between target and object.
 116   m_Target = clitk::Clone<ImageType>(temp2);
 117   clitk::AndNot<ImageType>(m_Target, m_Object, GetBackgroundValue());
 118
 119   // Define filter to compute statics on mask image
 120   typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
 121   typename StatFilterType::Pointer statFilter = StatFilterType::New();
 122
 123   // Compute the initial support size
 124   statFilter->SetInput(m_Support);
 125   statFilter->SetLabelInput(m_Support);
 126   statFilter->Update();
 127   SetSupportSize(statFilter->GetCount(GetForegroundValue()));
 128   // DD(GetSupportSize());
 129
 130   // Compute the initial target size
 131   ImagePointer s = clitk::ResizeImageLike<ImageType>(m_Support, m_Target, GetBackgroundValue());
 132   statFilter->SetInput(s);
 133   statFilter->SetLabelInput(m_Target);
 134   statFilter->Update();
 135   SetTargetSize(statFilter->GetCount(GetForegroundValue()));
 136   // DD(GetTargetSize());
 137
 138   //
 139   int bins = GetNumberOfBins();
 140   double tolerance = GetAreaLossTolerance();
 141
 142   // Compute Fuzzy map
 143   double angle = GetDirection().angle1;
 144   typename FloatImageType::Pointer map = ComputeFuzzyMap(m_Object, m_Target, m_Support, angle);
 145   writeImage<FloatImageType>(map, "fuzzy_"+toString(clitk::rad2deg(angle))+".mha");
 146
 147   // Compute the optimal thresholds (direct and inverse)
 148   double mThreshold=0.0;
 149   double mReverseThreshold=1.0;
 150   ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold);
 151
 152   // DD(mThreshold);
 153   // DD(mReverseThreshold);
 154
 155   // Use the threshold to compute new support
 156   int s1 = GetSupportSize();
 157   if (mThreshold > 0.0) {
 158     ImagePointer support1 =
 159       clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 160                                                      mThreshold,
 161                                                      angle,false, // inverseFlag
 162                                                      false,  // uniqueConnectedComponent
 163                                                      -1, true,
 164                                                      false);//singleObjectCCL
 165     // Compute the new support size
 166     statFilter->SetInput(support1);
 167     statFilter->SetLabelInput(support1);
 168     statFilter->Update();
 169     s1 = statFilter->GetCount(GetForegroundValue());
 170   }
 171
 172   int s2 = GetSupportSize();
 173   if (mReverseThreshold < 1.0) {
 174     ImagePointer support2 =
 175       clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 176                                                      mReverseThreshold,
 177                                                      angle,true,// inverseFlag
 178                                                      false, // uniqueConnectedComponent
 179                                                      -1, true,
 180                                                      false); //singleObjectCCL
 181     // Compute the new support size
 182     statFilter = StatFilterType::New();
 183     statFilter->SetInput(support2);
 184     statFilter->SetLabelInput(support2);
 185     statFilter->Update();
 186     s2 = statFilter->GetCount(GetForegroundValue());
 187   }
 188
 189   // Check threshold, if we gain nothing, we force to max/min thresholds
 190   // DD(GetSupportSize());
 191   // DD(s1);
 192   // DD(s2);
 193   if (s1 >= GetSupportSize()) mThreshold = 0.0;
 194   if (s2 >= GetSupportSize()) mReverseThreshold = 1.0;
 195
 196   // Set results values
 197   m_Info.threshold = mThreshold;
 198   m_Info.sizeAfterThreshold = s1;
 199   m_Info.sizeBeforeThreshold = GetSupportSize();
 200   m_Info.sizeReference = GetTargetSize();
 201   m_InfoReverse.threshold = mReverseThreshold;
 202   m_InfoReverse.sizeAfterThreshold = s2;
 203   m_InfoReverse.sizeBeforeThreshold = GetSupportSize();
 204   m_InfoReverse.sizeReference = GetTargetSize();
 205 }
 206 //--------------------------------------------------------------------
 207
 208
 209 //--------------------------------------------------------------------
 210 template <class ImageType>
 211 typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
 212 clitk::RelativePositionAnalyzerFilter<ImageType>::
 213 ComputeFuzzyMap(ImageType * object, ImageType * target, ImageType * support, double angle)
 214 {
 215   typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
 216   typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
 217   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 218   sliceRelPosFilter->VerboseStepFlagOff();
 219   sliceRelPosFilter->WriteStepFlagOff();
 220   sliceRelPosFilter->SetInput(support);
 221   sliceRelPosFilter->SetInputObject(object);
 222   sliceRelPosFilter->SetDirection(2);
 223   sliceRelPosFilter->SetIntermediateSpacingFlag(false);
 224   //sliceRelPosFilter->AddOrientationTypeString(orientation);
 225   sliceRelPosFilter->AddAnglesInRad(angle, 0.0);
 226   sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
 227   // sliceRelPosFilter->PrintOptions();
 228   sliceRelPosFilter->Update();
 229   typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
 230   writeImage<FloatImageType>(map, "fuzzy_0_"+toString(clitk::rad2deg(angle))+".mha");
 231
 232   // Resize object like map to allow SetBackground
 233   ImagePointer temp = clitk::ResizeImageLike<ImageType>(object, map, GetBackgroundValue());
 234   //  writeImage<FloatImageType>(map, "fuzzy_1_"+toString(clitk::rad2deg(angle))+".mha");
 235
 236   // Remove initial object from the fuzzy map
 237   map = clitk::SetBackground<FloatImageType, ImageType>(map, temp, GetForegroundValue(), 0.0, true);
 238   writeImage<FloatImageType>(map, "fuzzy_2_"+toString(clitk::rad2deg(angle))+".mha");
 239
 240   // Resize the fuzzy map like the target, put 2.0 when outside
 241   map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0);  // Put 2.0 when out of initial map
 242   writeImage<FloatImageType>(map, "fuzzy_3_"+toString(clitk::rad2deg(angle))+".mha");
 243
 244   // end
 245   return map;
 246 }
 247 //--------------------------------------------------------------------
 248
 249
 250 //--------------------------------------------------------------------
 251 template <class ImageType>
 252 void
 253 clitk::RelativePositionAnalyzerFilter<ImageType>::
 254 ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
 255                          double & threshold, double & reverseThreshold)
 256 {
 257   // Get the histogram of fuzzy values inside the target image
 258   typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
 259   typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
 260   f->SetInput(map);
 261   f->SetLabelInput(target);
 262   f->UseHistogramsOn();
 263   f->SetHistogramParameters(bins, 0.0-(1.0/bins), 1.0+(1.0/bins));
 264   f->Update();
 265   int count = f->GetCount(GetForegroundValue());
 266   // DD(count);
 267   typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
 268
 269   // Debug : dump histogram
 270   static int i=0;
 271   std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
 272   for(int j=0; j<bins; j++) {
 273     histogramFile << h->GetMeasurement(j,0)
 274                   << "\t" << h->GetFrequency(j)
 275                   << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
 276   }
 277   histogramFile.close();
 278   std::ofstream histogramFile2(std::string("fuzzy_histo_R_"+toString(i)+".txt").c_str());
 279   for(int j=bins-1; j>=0; j--) {
 280     histogramFile2 << h->GetMeasurement(j,0)
 281                   << "\t" << h->GetFrequency(j)
 282                   << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
 283   }
 284   histogramFile2.close();
 285   i++;
 286
 287   // Analyze the histogram (direct)
 288   double sum = 0.0;
 289   bool found = false;
 290   threshold = 0.0;
 291   for(int j=0; j<bins-1; j++) {
 292     sum += ((double)h->GetFrequency(j)/(double)count);
 293      // DD(j);
 294      // DD(sum);
 295      // DD(threshold);
 296      // DD(h->GetBinMin(0,j));
 297      // DD(h->GetBinMax(0,j));
 298     if ((!found) && (sum > tolerance)) {
 299       // We consider as threshold the laste before current, because
 300       if (j==0)
 301         threshold = h->GetBinMin(0,j);
 302       else threshold = h->GetBinMin(0,j-1); // FIXME  ? the last before reaching the threshold
 303       // DD(threshold);
 304       found = true;
 305       j = bins;
 306     }
 307   }
 308
 309   // Analyze the histogram (reverse)
 310   sum = 0.0;
 311   found = false;
 312   reverseThreshold = 1.0;
 313   for(int j=bins-1; j>0; j--) {
 314     sum += ((double)h->GetFrequency(j)/(double)count);
 315      // DD(j);
 316      // DD(sum);
 317      // DD(reverseThreshold);
 318      // DD(h->GetBinMin(0,j));
 319      // DD(h->GetBinMax(0,j));
 320     if ((!found) && (sum > tolerance)) {
 321       if (j==bins-1)
 322       reverseThreshold = h->GetBinMax(0,j);
 323       else reverseThreshold = h->GetBinMax(0,j-1);// FIXME  ? the last before reaching the threshold
 324       // DD(reverseThreshold);
 325       found = true;
 326       j = -1;
 327     }
 328   }
 329
 330 }
 331 //--------------------------------------------------------------------
 332