segmentation/clitkExtractLymphStationsFilter.old.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://oncora1.lyon.fnclcc.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 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkCommon.h"
  24 #include "clitkExtractLymphStationsFilter.h"
  25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  26 #include "clitkSegmentationUtils.h"
  27 #include "clitkAutoCropFilter.h"
  28 #include "clitkSegmentationUtils.h"
  29 #include "clitkSliceBySliceRelativePositionFilter.h"
  30
  31 // itk
  32 #include <itkStatisticsLabelMapFilter.h>
  33 #include <itkLabelImageToStatisticsLabelMapFilter.h>
  34 #include <itkRegionOfInterestImageFilter.h>
  35 #include <itkBinaryThresholdImageFilter.h>
  36 #include <itkImageSliceConstIteratorWithIndex.h>
  37 #include <itkImageSliceIteratorWithIndex.h>
  38 #include <itkBinaryThinningImageFilter.h>
  39
  40 // itk ENST
  41 #include "RelativePositionPropImageFilter.h"
  42
  43 //--------------------------------------------------------------------
  44 template <class TImageType>
  45 clitk::ExtractLymphStationsFilter<TImageType>::
  46 ExtractLymphStationsFilter():
  47   clitk::FilterBase(),
  48   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  49   itk::ImageToImageFilter<TImageType, MaskImageType>()
  50 {
  51   this->SetNumberOfRequiredInputs(1);
  52   SetBackgroundValue(0);
  53   SetForegroundValue(1);
  54   SetFuzzyThreshold(0.5);
  55 }
  56 //--------------------------------------------------------------------
  57
  58
  59 //--------------------------------------------------------------------
  60 template <class TImageType>
  61 void
  62 clitk::ExtractLymphStationsFilter<TImageType>::
  63 GenerateOutputInformation() {
  64   //  Superclass::GenerateOutputInformation();
  65
  66   // Get input
  67   LoadAFDB();
  68   m_mediastinum = GetAFDB()->template GetImage <MaskImageType>("mediastinum");
  69   m_trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
  70   ImagePointType carina;
  71   GetAFDB()->GetPoint3D("carina", carina);
  72   m_CarinaZPositionInMM = carina[2];
  73   DD(m_CarinaZPositionInMM);
  74   ImagePointType mlb;
  75   GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
  76   m_MiddleLobeBronchusZPositionInMM = mlb[2];
  77   DD(m_MiddleLobeBronchusZPositionInMM);
  78
  79   // ----------------------------------------------------------------
  80   // ----------------------------------------------------------------
  81   // Superior limit = carina
  82   // Inferior limit = origin right middle lobe bronchus
  83   StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
  84   m_working_mediastinum =
  85     clitk::CropImageAlongOneAxis<MaskImageType>(m_mediastinum, 2,
  86                                                 m_MiddleLobeBronchusZPositionInMM,
  87                                                 m_CarinaZPositionInMM, true,
  88                                                 GetBackgroundValue());
  89   StopCurrentStep<MaskImageType>(m_working_mediastinum);
  90   m_output = m_working_mediastinum;
  91
  92   // ----------------------------------------------------------------
  93   // ----------------------------------------------------------------
  94   // Superior limit = carina
  95   // Inferior limit = origin right middle lobe bronchus
  96   StartNewStep("Inf/Sup trachea limits with carina/bronchus");
  97   m_working_trachea =
  98     clitk::CropImageAlongOneAxis<MaskImageType>(m_trachea, 2,
  99                                                 m_MiddleLobeBronchusZPositionInMM,
 100                                                 m_CarinaZPositionInMM, true,
 101                                                 GetBackgroundValue());
 102   StopCurrentStep<MaskImageType>(m_working_trachea);
 103
 104   // ----------------------------------------------------------------
 105   // ----------------------------------------------------------------
 106   // Separate trachea in two CCL
 107   StartNewStep("Separate trachea under carina");
 108
 109   // Labelize and consider two main labels
 110   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
 111
 112   // Carina position must at the first slice that separate the two main bronchus (not superiorly)
 113   // Check that upper slice is composed of at least two labels
 114   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 115   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
 116   iter.SetFirstDirection(0);
 117   iter.SetSecondDirection(1);
 118   iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 119   int maxLabel=0;
 120   while (!iter.IsAtReverseEndOfSlice()) {
 121     while (!iter.IsAtReverseEndOfLine()) {
 122       //  DD(iter.GetIndex());
 123       if (iter.Get() > maxLabel) maxLabel = iter.Get();
 124       --iter;
 125     }
 126     iter.PreviousLine();
 127   }
 128   if (maxLabel < 2) {
 129     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
 130   }
 131
 132   // Compute centroid of both parts to identify the left from the right bronchus
 133   typedef long LabelType;
 134   static const unsigned int Dim = ImageType::ImageDimension;
 135   typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 136   typedef itk::LabelMap< LabelObjectType > LabelMapType;
 137   typedef itk::LabelImageToLabelMapFilter<MaskImageType, LabelMapType> ImageToMapFilterType;
 138   typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 139   typedef itk::ShapeLabelMapFilter<LabelMapType, MaskImageType> ShapeFilterType;
 140   typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 141   imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
 142   imageToLabelFilter->SetInput(m_working_trachea);
 143   statFilter->SetInput(imageToLabelFilter->GetOutput());
 144   statFilter->Update();
 145   typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 146
 147   ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
 148   ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
 149
 150   ImagePixelType leftLabel;
 151   ImagePixelType rightLabel;
 152   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
 153   else { leftLabel = 2; rightLabel = 1; }
 154   DD(leftLabel);
 155   DD(rightLabel);
 156
 157   StopCurrentStep<MaskImageType>(m_working_trachea);
 158
 159   //-----------------------------------------------------
 160   StartNewStep("Left limits with bronchus (slice by slice)");
 161   // Select LeftLabel (set one label to Backgroundvalue)
 162   MaskImagePointer leftBronchus =
 163     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
 164                                                 rightLabel, GetBackgroundValue());
 165   MaskImagePointer rightBronchus  =
 166     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
 167                                                 leftLabel, GetBackgroundValue());
 168   writeImage<MaskImageType>(leftBronchus, "left.mhd");
 169   writeImage<MaskImageType>(rightBronchus, "right.mhd");
 170
 171   m_working_mediastinum =
 172     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 173                                                        leftBronchus, 2,
 174                                                        GetFuzzyThreshold(), "RightTo",
 175                                                        true, 4);
 176   m_working_mediastinum =
 177     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 178                                                        rightBronchus,
 179                                                        2, GetFuzzyThreshold(), "LeftTo",
 180                                                        true, 4);
 181   m_working_mediastinum =
 182     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 183                                                        leftBronchus,
 184                                                        2, GetFuzzyThreshold(), "AntTo",
 185                                                        true, 4, true); // NOT
 186   m_working_mediastinum =
 187     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 188                                                        rightBronchus,
 189                                                        2, GetFuzzyThreshold(), "AntTo",
 190                                                        true, 4, true);
 191   m_working_mediastinum =
 192     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 193                                                        leftBronchus,
 194                                                        2, GetFuzzyThreshold(), "PostTo",
 195                                                        true, 4, true);
 196   m_working_mediastinum =
 197     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
 198                                                        rightBronchus,
 199                                                        2, GetFuzzyThreshold(), "PostTo",
 200                                                        true, 4, true);
 201   m_output = m_working_mediastinum;
 202   StopCurrentStep<MaskImageType>(m_output);
 203
 204   //-----------------------------------------------------
 205   //-----------------------------------------------------
 206   StartNewStep("Posterior limits");
 207
 208   // Left Bronchus slices
 209   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 210   typedef typename ExtractSliceFilterType::SliceType SliceType;
 211   typename ExtractSliceFilterType::Pointer
 212     extractSliceFilter = ExtractSliceFilterType::New();
 213   extractSliceFilter->SetInput(leftBronchus);
 214   extractSliceFilter->SetDirection(2);
 215   extractSliceFilter->Update();
 216   std::vector<typename SliceType::Pointer> leftBronchusSlices;
 217   extractSliceFilter->GetOutputSlices(leftBronchusSlices);
 218
 219   // Right Bronchus slices
 220   extractSliceFilter = ExtractSliceFilterType::New();
 221   extractSliceFilter->SetInput(rightBronchus);
 222   extractSliceFilter->SetDirection(2);
 223   extractSliceFilter->Update();
 224   std::vector<typename SliceType::Pointer> rightBronchusSlices ;
 225   extractSliceFilter->GetOutputSlices(rightBronchusSlices);
 226
 227   assert(leftBronchusSlices.size() == rightBronchusSlices.size());
 228
 229   std::vector<MaskImageType::PointType> leftPoints;
 230   std::vector<MaskImageType::PointType> rightPoints;
 231   for(uint i=0; i<leftBronchusSlices.size(); i++) {
 232     // Keep main CCL
 233     leftBronchusSlices[i] = Labelize<SliceType>(leftBronchusSlices[i], 0, true, 10);
 234     leftBronchusSlices[i] = KeepLabels<SliceType>(leftBronchusSlices[i],
 235                                                   GetBackgroundValue(),
 236                                                   GetForegroundValue(), 1, 1, true);
 237     rightBronchusSlices[i] = Labelize<SliceType>(rightBronchusSlices[i], 0, true, 10);
 238     rightBronchusSlices[i] = KeepLabels<SliceType>(rightBronchusSlices[i],
 239                                                    GetBackgroundValue(),
 240                                                    GetForegroundValue(), 1, 1, true);
 241     double distance_max_from_center_point = 15;
 242
 243     // ------- Find point in left Bronchus -------
 244     // find right most point in left  = rightMost
 245     SliceType::PointType a;
 246     SliceType::PointType rightMost =
 247       clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
 248                                                           GetBackgroundValue(),
 249                                                           0, false, a, 0);
 250     // find post most point in left, not far away from rightMost
 251     SliceType::PointType p =
 252       clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
 253                                                           GetBackgroundValue(),
 254                                                           1, false, rightMost,
 255                                                           distance_max_from_center_point);
 256     MaskImageType::PointType pp;
 257     pp[0] = p[0]; pp[1] = p[1];
 258     pp[2] = i*leftBronchus->GetSpacing()[2] + leftBronchus->GetOrigin()[2];
 259     leftPoints.push_back(pp);
 260
 261     // ------- Find point in right Bronchus -------
 262     // find left most point in right  = leftMost
 263     SliceType::PointType leftMost =
 264       clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
 265                                                           GetBackgroundValue(),
 266                                                           0, true, a, 0);
 267     // find post most point in left, not far away from leftMost
 268     p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
 269                                                             GetBackgroundValue(),
 270                                                             1, false, leftMost,
 271                                                             distance_max_from_center_point);
 272     pp[0] = p[0]; pp[1] = p[1];
 273     pp[2] = i*rightBronchus->GetSpacing()[2] + rightBronchus->GetOrigin()[2];
 274     rightPoints.push_back(pp);
 275   }
 276
 277   // DEBUG
 278   std::ofstream osl;
 279   openFileForWriting(osl, "left.txt");
 280   osl << "LANDMARKS1" << std::endl;
 281   std::ofstream osr;
 282   openFileForWriting(osr, "right.txt");
 283   osr << "LANDMARKS1" << std::endl;
 284   for(uint i=0; i<leftBronchusSlices.size(); i++) {
 285     osl << i << " " << leftPoints[i][0] << " " << leftPoints[i][1]
 286         << " " << leftPoints[i][2] << " 0 0 " << std::endl;
 287     osr << i << " " << rightPoints[i][0] << " " << rightPoints[i][1]
 288         << " " << rightPoints[i][2] << " 0 0 " << std::endl;
 289   }
 290   osl.close();
 291   osr.close();
 292
 293   // Now uses these points to limit, slice by slice
 294   // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 295   /*
 296     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 297     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 298     This will equal zero if the point C is on the line formed by points A and B, and will have a different sign depending on the side. Which side this is depends on the orientation of your (x,y) coordinates, but you can plug test values for A,B and C into this formula to determine whether negative values are to the left or to the right.
 299
 300     => to accelerate, start with formula, when change sign -> stop and fill
 301   */
 302   //  typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 303   iter = SliceIteratorType(m_working_mediastinum, m_working_mediastinum->GetLargestPossibleRegion());
 304   iter.SetFirstDirection(0);
 305   iter.SetSecondDirection(1);
 306   iter.GoToBegin();
 307   int i=0;
 308   MaskImageType::PointType A;
 309   MaskImageType::PointType B;
 310   MaskImageType::PointType C;
 311   while (!iter.IsAtEnd()) {
 312     A = leftPoints[i];
 313     B = rightPoints[i];
 314     C = A;
 315     C[1] -= 10; // I know I must keep this point
 316     double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 317     bool isPositive = s<0;
 318     while (!iter.IsAtEndOfSlice()) {
 319       while (!iter.IsAtEndOfLine()) {
 320         // Very slow, I know ... but image should be very small
 321         m_working_mediastinum->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
 322         double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 323         if (s == 0) iter.Set(2);
 324         if (isPositive) {
 325           if (s > 0) iter.Set(GetBackgroundValue());
 326         }
 327         else {
 328           if (s < 0) iter.Set(GetBackgroundValue());
 329         }
 330         ++iter;
 331       }
 332       iter.NextLine();
 333     }
 334     iter.NextSlice();
 335     ++i;
 336   }
 337
 338   //-----------------------------------------------------
 339   //-----------------------------------------------------
 340   // StartNewStep("Anterior limits");
 341
 342
 343   // MISSING FROM NOW
 344
 345   // Station 4R, Station 4L, the right pulmonary artery, and/or the
 346   // left superior pulmonary vein
 347
 348
 349   //-----------------------------------------------------
 350   //-----------------------------------------------------
 351   // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
 352
 353
 354   // Set output image information (required)
 355   MaskImagePointer outputImage = this->GetOutput(0);
 356   outputImage->SetRegions(m_working_mediastinum->GetLargestPossibleRegion());
 357   outputImage->SetOrigin(m_working_mediastinum->GetOrigin());
 358   outputImage->SetRequestedRegion(m_working_mediastinum->GetLargestPossibleRegion());
 359
 360 }
 361 //--------------------------------------------------------------------
 362
 363
 364 //--------------------------------------------------------------------
 365 template <class TImageType>
 366 void
 367 clitk::ExtractLymphStationsFilter<TImageType>::
 368 GenerateInputRequestedRegion() {
 369   DD("GenerateInputRequestedRegion");
 370   // Call default
 371   //  Superclass::GenerateInputRequestedRegion();
 372   // Following needed because output region can be greater than input (trachea)
 373   //ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
 374   //ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
 375   DD("set reg");
 376   m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
 377   m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
 378   DD("end");
 379 }
 380 //--------------------------------------------------------------------
 381
 382
 383 //--------------------------------------------------------------------
 384 template <class TImageType>
 385 void
 386 clitk::ExtractLymphStationsFilter<TImageType>::
 387 GenerateData() {
 388   DD("GenerateData");
 389   // Final Step -> graft output (if SetNthOutput => redo)
 390   this->GraftOutput(m_output);
 391 }
 392 //--------------------------------------------------------------------
 393
 394
 395 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX