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 }
  32 //--------------------------------------------------------------------
  33
  34
  35 //--------------------------------------------------------------------
  36 template <class ImageType>
  37 void
  38 clitk::RelativePositionAnalyzerFilter<ImageType>::
  39 SetInputSupport(const ImageType * image)
  40 {
  41   // Process object is not const-correct so the const casting is required.
  42   this->SetNthInput(0, const_cast<ImageType *>(image));
  43 }
  44 //--------------------------------------------------------------------
  45
  46
  47 //--------------------------------------------------------------------
  48 template <class ImageType>
  49 void
  50 clitk::RelativePositionAnalyzerFilter<ImageType>::
  51 SetInputObject(const ImageType * image)
  52 {
  53   // Process object is not const-correct so the const casting is required.
  54   this->SetNthInput(1, const_cast<ImageType *>(image));
  55 }
  56 //--------------------------------------------------------------------
  57
  58
  59 //--------------------------------------------------------------------
  60 template <class ImageType>
  61 void
  62 clitk::RelativePositionAnalyzerFilter<ImageType>::
  63 SetInputTarget(const ImageType * image)
  64 {
  65   // Process object is not const-correct so the const casting is required.
  66   this->SetNthInput(2, const_cast<ImageType *>(image));
  67 }
  68 //--------------------------------------------------------------------
  69
  70
  71 //--------------------------------------------------------------------
  72 template <class ImageType>
  73 void
  74 clitk::RelativePositionAnalyzerFilter<ImageType>::
  75 PrintOptions()
  76 {
  77   DD("TODO");
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 template <class ImageType>
  84 void
  85 clitk::RelativePositionAnalyzerFilter<ImageType>::
  86 GenerateOutputInformation()
  87 {
  88   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  89   ImagePointer outputImage = this->GetOutput(0);
  90   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
  91 }
  92 //--------------------------------------------------------------------
  93
  94
  95 //--------------------------------------------------------------------
  96 template <class ImageType>
  97 void
  98 clitk::RelativePositionAnalyzerFilter<ImageType>::
  99 GenerateInputRequestedRegion()
 100 {
 101   // Call default
 102   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
 103   // Get input pointers and set requested region to common region
 104   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 105   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 106   ImagePointer input3 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
 107   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 108   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 109   input3->SetRequestedRegion(input3->GetLargestPossibleRegion());
 110 }
 111 //--------------------------------------------------------------------
 112
 113 //--------------------------------------------------------------------
 114 template <class ImageType>
 115 void
 116 clitk::RelativePositionAnalyzerFilter<ImageType>::
 117 GenerateData()
 118 {
 119   this->LoadAFDB();
 120
 121   // Get input pointer
 122   m_Support = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 123   m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 124   m_Target = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
 125   static const unsigned int dim = ImageType::ImageDimension;
 126
 127   // Remove object from support
 128   clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
 129
 130   // Resize object like target (to enable substraction later)
 131   ImagePointer objectLikeTarget = clitk::ResizeImageLike<ImageType>(m_Object, m_Target, GetBackgroundValue());
 132
 133   // Define filter to compute statics on mask image
 134   typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
 135   typename StatFilterType::Pointer statFilter = StatFilterType::New();
 136
 137   // Compute the initial support size
 138   statFilter->SetInput(m_Support);
 139   statFilter->SetLabelInput(m_Support);
 140   statFilter->Update();
 141   int m_SupportSize = statFilter->GetCount(GetForegroundValue());
 142   // DD(m_SupportSize);
 143
 144   // Compute the initial target size
 145   ImagePointer s = clitk::ResizeImageLike<ImageType>(m_Support, m_Target, GetBackgroundValue());
 146   statFilter->SetInput(s);
 147   statFilter->SetLabelInput(m_Target);
 148   statFilter->Update();
 149   int m_TargetSize = statFilter->GetCount(GetForegroundValue());
 150   // DD(m_TargetSize);
 151
 152   // Build the list of tested orientations
 153   std::vector<double> m_ListOfAngles;
 154   int m_NumberOfAngles = this->GetAFDB()->GetDouble("NumberOfAngles");
 155   for(uint i=0; i<m_NumberOfAngles; i++) {
 156     double a = i*360.0/m_NumberOfAngles;
 157     if (a>180) a = 180-a;
 158     m_ListOfAngles.push_back(clitk::deg2rad(a));
 159   }
 160
 161   // Loop on all orientations
 162   for(int i=0; i<m_ListOfAngles.size(); i++) {
 163
 164     // Compute Fuzzy map
 165     typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]);
 166     writeImage<FloatImageType>(map, "fuzzy_"+toString(i)+".mha");
 167
 168     // Compute the optimal thresholds (direct and inverse)
 169     double mThreshold;
 170     double mReverseThreshold;
 171     int bins = this->GetAFDB()->GetDouble("bins");
 172     double tolerance = this->GetAFDB()->GetDouble("TargetAreaLossTolerance");
 173     ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold);
 174
 175     // Use the threshold to compute new support
 176     int s1;
 177     if (mThreshold > 0.0) {
 178       ImagePointer support1 =
 179         clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 180                                                        mThreshold,
 181                                                        m_ListOfAngles[i],false,
 182                                                        false, -1, true, false);
 183       writeImage<ImageType>(support1, "sup_"+toString(i)+".mha");
 184       // Compute the new support size
 185       statFilter->SetInput(support1);
 186       statFilter->SetLabelInput(support1);
 187       statFilter->Update();
 188       s1 = statFilter->GetCount(GetForegroundValue());
 189     }
 190     else s1 = m_SupportSize;
 191
 192     int s2;
 193     if (mReverseThreshold < 1.0) {
 194       ImagePointer support2 =
 195         clitk::SliceBySliceRelativePosition<ImageType>(m_Support, m_Object, 2,
 196                                                        mReverseThreshold,
 197                                                        m_ListOfAngles[i],true,
 198                                                        false, -1, true, false);
 199       writeImage<ImageType>(support2, "sup_rev_"+toString(i)+".mha");
 200       // Compute the new support size
 201       statFilter = StatFilterType::New();
 202       statFilter->SetInput(support2);
 203       statFilter->SetLabelInput(support2);
 204       statFilter->Update();
 205       s2 = statFilter->GetCount(GetForegroundValue());
 206     }
 207     else s2 =m_SupportSize;
 208
 209     // Print results
 210     std::cout << i << " " << clitk::rad2deg(m_ListOfAngles[i]) << "\t"
 211               << m_SupportSize << " " << m_TargetSize << "\t"
 212               << s1/(double)m_SupportSize << " " << s2/(double)m_SupportSize << "\t"
 213               << mThreshold << " " << mReverseThreshold << std::endl;
 214
 215   } // end loop on orientations
 216
 217
 218   // Final Step -> set output TODO
 219   // this->SetNthOutput(0, working_image);
 220   //  this->GraftOutput(working_image);
 221 }
 222 //--------------------------------------------------------------------
 223
 224
 225 //--------------------------------------------------------------------
 226 template <class ImageType>
 227 typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
 228 clitk::RelativePositionAnalyzerFilter<ImageType>::
 229 ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
 230 {
 231   typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
 232   typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
 233   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 234   sliceRelPosFilter->VerboseStepFlagOff();
 235   sliceRelPosFilter->WriteStepFlagOff();
 236   sliceRelPosFilter->SetInput(target);
 237   sliceRelPosFilter->SetInputObject(object);
 238   sliceRelPosFilter->SetDirection(2);
 239   sliceRelPosFilter->SetIntermediateSpacingFlag(false);
 240   //sliceRelPosFilter->AddOrientationTypeString(orientation);
 241   sliceRelPosFilter->AddAngles(angle, 0.0);
 242   sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
 243   // sliceRelPosFilter->PrintOptions();
 244   sliceRelPosFilter->Update();
 245   typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
 246
 247   // Remove initial object from the fuzzy map
 248   map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
 249
 250   // Resize the fuzzy map like the target, put 2.0 when outside
 251   map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0);  // Put 2.0 when out of initial map
 252
 253   // end
 254   return map;
 255 }
 256 //--------------------------------------------------------------------
 257
 258
 259 //--------------------------------------------------------------------
 260 template <class ImageType>
 261 void
 262 clitk::RelativePositionAnalyzerFilter<ImageType>::
 263 ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
 264                          double & threshold, double & reverseThreshold)
 265 {
 266   // Get the histogram of fuzzy values inside the target image
 267   typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
 268   typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
 269   f->SetInput(map);
 270   f->SetLabelInput(target);
 271   f->UseHistogramsOn();
 272   f->SetHistogramParameters(bins, 0.0, 1.1);
 273   f->Update();
 274   int count = f->GetCount(GetForegroundValue());
 275   typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
 276
 277   // Debug : dump histogram
 278   static int i=0;
 279   std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
 280   for(int j=0; j<bins; j++) {
 281     histogramFile << h->GetMeasurement(j,0)
 282                   << "\t" << h->GetFrequency(j)
 283                   << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
 284   }
 285   histogramFile.close();
 286   i++;
 287
 288   // Analyze the histogram (direct)
 289   double sum = 0.0;
 290   bool found = false;
 291   threshold = 0.0;
 292   for(int j=0; j<bins; j++) {
 293     sum += ((double)h->GetFrequency(j)/(double)count);
 294     if ((!found) && (sum > tolerance)) {
 295       threshold = h->GetBinMin(0,j);
 296       found = true;
 297     }
 298   }
 299
 300   // Analyze the histogram (reverse)
 301   sum = 0.0;
 302   found = false;
 303   reverseThreshold = 1.0;
 304   for(int j=bins-1; j>=0; j--) {
 305     sum += ((double)h->GetFrequency(j)/(double)count);
 306     if ((!found) && (sum > tolerance)) {
 307       reverseThreshold = h->GetBinMax(0,j);
 308       found = true;
 309     }
 310   }
 311 }
 312 //--------------------------------------------------------------------
 313