]> Creatis software - clitk.git/blobdiff - segmentation/clitkExtractLymphStationsFilter.txx
small improvement for S8
[clitk.git] / segmentation / clitkExtractLymphStationsFilter.txx
index 993d2888f52ee3dbc328d6c4bc4c1f5dd223f221..ea9e2b678c01fe40d75fffc8d2b009f3ada6d550 100644 (file)
 #include "clitkSliceBySliceRelativePositionFilter.h"
 
 // itk
-#include <deque>
 #include <itkStatisticsLabelMapFilter.h>
 #include <itkLabelImageToStatisticsLabelMapFilter.h>
 #include <itkRegionOfInterestImageFilter.h>
 #include <itkBinaryThresholdImageFilter.h>
 #include <itkImageSliceConstIteratorWithIndex.h>
+#include <itkImageSliceIteratorWithIndex.h>
+#include <itkBinaryThinningImageFilter.h>
 
 // itk ENST
 #include "RelativePositionPropImageFilter.h"
@@ -44,16 +45,20 @@ template <class TImageType>
 clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractLymphStationsFilter():
   clitk::FilterBase(),
-  itk::ImageToImageFilter<TImageType, TImageType>()
+  clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
+  itk::ImageToImageFilter<TImageType, MaskImageType>()
 {
   this->SetNumberOfRequiredInputs(1);
-  SetBackgroundValueMediastinum(0);
-  SetBackgroundValueTrachea(0);
   SetBackgroundValue(0);
   SetForegroundValue(1);
 
-  SetIntermediateSpacing(6);
-  SetFuzzyThreshold1(0.6);
+  // Default values
+  ExtractStation_8_SetDefaultValues();
+  ExtractStation_3P_SetDefaultValues();
+
+  // Station 7
+  SetFuzzyThreshold(0.5);
+  SetStation7Filename("station7.mhd");
 }
 //--------------------------------------------------------------------
 
@@ -62,13 +67,51 @@ ExtractLymphStationsFilter():
 template <class TImageType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
-SetInputMediastinumLabelImage(const TImageType * image, ImagePixelType bg) {
-  this->SetNthInput(0, const_cast<TImageType *>(image));
-  m_BackgroundValueMediastinum = bg;
-  SetCarenaZPositionInMM(image->GetOrigin()[2]+image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]);
-  SetMiddleLobeBronchusZPositionInMM(image->GetOrigin()[2]);
-  // DD(m_CarenaZPositionInMM);
-//   DD(m_MiddleLobeBronchusZPositionInMM);
+GenerateOutputInformation() { 
+  // Get inputs
+  LoadAFDB();
+  m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
+  m_Mediastinum = GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
+
+  // Extract Station8
+  StartNewStep("Station 8");
+  StartSubStep(); 
+  ExtractStation_8();
+  StopSubStep();
+
+  // Extract Station3P
+  StartNewStep("Station 3P");
+  StartSubStep(); 
+  ExtractStation_3P();
+  StopSubStep();
+
+  if (0) { // temporary suppress
+    // Extract Station7
+    StartNewStep("Station 7");
+    StartSubStep();
+    ExtractStation_7();
+    StopSubStep();
+
+    // Extract Station4RL
+    StartNewStep("Station 4RL");
+    StartSubStep();
+    //ExtractStation_4RL();
+    StopSubStep();
+  }
+
+
+  //
+  //  typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BFilter;
+  //BFilter::Pointer merge = BFilter::New();  
+  // writeImage<MaskImageType>(m_Output, "ouput.mhd");
+  //writeImage<MaskImageType>(m_Working_Support, "ws.mhd");
+  /*merge->SetInput1(m_Station7);
+    merge->SetInput2(m_Station4RL); // support
+    merge->SetOperationType(BFilter::AndNot); CHANGE OPERATOR
+    merge->SetForegroundValue(4);
+    merge->Update();
+    m_Output = merge->GetOutput();
+  */
 }
 //--------------------------------------------------------------------
 
@@ -77,271 +120,124 @@ SetInputMediastinumLabelImage(const TImageType * image, ImagePixelType bg) {
 template <class TImageType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
-SetInputTracheaLabelImage(const TImageType * image, ImagePixelType bg) {
-  this->SetNthInput(1, const_cast<TImageType *>(image));
-  m_BackgroundValueTrachea = bg;
+GenerateInputRequestedRegion() {
+  //DD("GenerateInputRequestedRegion (nothing?)");
 }
 //--------------------------------------------------------------------
 
 
 //--------------------------------------------------------------------
 template <class TImageType>
-template<class ArgsInfoType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
-SetArgsInfo(ArgsInfoType mArgsInfo)
-{
-  SetVerboseOption_GGO(mArgsInfo);
-  SetVerboseStep_GGO(mArgsInfo);
-  SetWriteStep_GGO(mArgsInfo);
-  SetVerboseWarningOff_GGO(mArgsInfo);
-  SetCarenaZPositionInMM_GGO(mArgsInfo);
-  SetMiddleLobeBronchusZPositionInMM_GGO(mArgsInfo);
-  SetIntermediateSpacing_GGO(mArgsInfo);
-  SetFuzzyThreshold1_GGO(mArgsInfo);
-  //SetBackgroundValueMediastinum_GGO(mArgsInfo);
+GenerateData() {
+  DD("GenerateData, graft output");
+
+  // Final Step -> graft output (if SetNthOutput => redo)
+  this->GraftOutput(m_ListOfStations["8"]);
 }
 //--------------------------------------------------------------------
-
+  
 
 //--------------------------------------------------------------------
 template <class TImageType>
-void 
+bool 
 clitk::ExtractLymphStationsFilter<TImageType>::
-GenerateOutputInformation() { 
-  //  Superclass::GenerateOutputInformation();
+CheckForStation(std::string station) 
+{
+  // Compute Station name
+  std::string s = "Station"+station;
   
-  // Get input
-  m_mediastinum = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
-  m_trachea = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
-    
-  // ----------------------------------------------------------------
-  // ----------------------------------------------------------------
-  // Superior limit = carena
-  // Inferior limit = origine middle lobe bronchus
-  StartNewStep("Inf/Sup limits with carena/bronchus");
-  ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
-  ImageSizeType size = region.GetSize();
-  ImageIndexType index = region.GetIndex();
-  DD(m_CarenaZPositionInMM);
-  DD(m_MiddleLobeBronchusZPositionInMM);
-  index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
-  size[2] = ceil((m_CarenaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]);
-  region.SetSize(size);
-  region.SetIndex(index);
-  DD(region);
-  typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
-  typename CropFilterType::Pointer cropFilter = CropFilterType::New();
-  cropFilter->SetInput(m_mediastinum);
-  cropFilter->SetRegionOfInterest(region);
-  cropFilter->Update();
-  m_working_image = cropFilter->GetOutput();
-  // Auto Crop (because following rel pos is faster)
-  m_working_image = clitk::AutoCrop<ImageType>(m_working_image, 0); 
-  StopCurrentStep<ImageType>(m_working_image);
-  m_output = m_working_image;
-
-  // ----------------------------------------------------------------
-  // ----------------------------------------------------------------
-  // Separate trachea in two CCL
-  StartNewStep("Separate trachea");
-  // DD(region);
-  ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
-  for(int i=0; i<3; i++) {
-    index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i]
-                      -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]);
-    // DD(index[i]);
-    size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
-    //  DD(size[i]);
-    if (index[i] < 0) { 
-      size[i] += index[i];
-      index[i] = 0;       
-    }
-    if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) {
-      size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i];
-    }
+
+  // Check if station already exist in DB
+  bool found = false;
+  if (GetAFDB()->TagExist(s)) {
+    m_ListOfStations[station] = GetAFDB()->template GetImage<MaskImageType>(s);
+    found = true;
   }
-  // DD(index);
-  //   DD(size);
-  region.SetIndex(index);
-  region.SetSize(size);  
-  //  typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
-  //  typename CropFilterType::Pointer 
-  cropFilter = CropFilterType::New();
- //  m_trachea.Print(std::cout);
-  cropFilter->SetInput(m_trachea);
-  cropFilter->SetRegionOfInterest(region);
-  cropFilter->Update();
-  m_working_trachea = cropFilter->GetOutput();
-
-  // Labelize and consider two main labels
-  m_working_trachea = Labelize<ImageType>(m_working_trachea, 0, true, 1);
-  StopCurrentStep<ImageType>(m_working_trachea);
-  
-  // Detect wich label is at Left
-  typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
-  SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
-  iter.SetFirstDirection(0);
-  iter.SetSecondDirection(1);
-  iter.GoToBegin();
-  bool stop = false;
-  ImagePixelType leftLabel;
-  while (!stop) {
-    if (iter.Get() != m_BackgroundValueTrachea) {
-      //     DD(iter.GetIndex());
-      //       DD((int)iter.Get());
-      leftLabel = iter.Get();
-      stop = true;
-    }
-    ++iter;
+
+  // Define the starting support 
+  if (found && GetComputeStation(station)) {
+    std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
   }
-  DD((int)leftLabel);
-  
-  /*
-  // Relative position 
-  StartNewStep("Left/Right limits with trachea");
-
-  // Select LeftLabel (set label 1 to 0)
-  ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
-  writeImage<ImageType>(temp, "temp1.mhd");
-
-  // Left relative position
-  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
-  relPosFilter->VerboseStepOff();
-  relPosFilter->WriteStepOff();
-  relPosFilter->SetInput(m_working_image); 
-  relPosFilter->SetInputObject(temp); 
-  relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
-  relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
-  relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
-  relPosFilter->Update();
-  m_working_image = relPosFilter->GetOutput();
-
-  // Select RightLabel (set label 2 to 0)
-  temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 2, 0);
-  writeImage<ImageType>(temp, "temp2.mhd");
-
-  // Left relative position
-  relPosFilter = RelPosFilterType::New();
-  relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
-  relPosFilter->VerboseStepOn();
-  relPosFilter->WriteStepOn();
-  relPosFilter->SetInput(m_working_image); 
-  relPosFilter->SetInputObject(temp); 
-  relPosFilter->SetOrientationType(RelPosFilterType::LeftTo);
-  relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
-  relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
-  relPosFilter->Update();
-  m_working_image = relPosFilter->GetOutput();
-  */
+  if (!found || GetComputeStation(station)) {
+    m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
+    return true;
+  }
+  else {
+    std::cout << "Station " << station << " found. I used it" << std::endl;
+    return false;
+  }
+}
+//--------------------------------------------------------------------
 
-  //-----------------------------------------------------
-  StartNewStep("Left/Right limits with trachea (slice by slice");
-  
-  // Select LeftLabel (set label 1 to 0)
-  ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
-  writeImage<ImageType>(temp, "temp1.mhd");
 
-  /*
-    - écrire utilisation de filter SliceBySliceRelPosFilter (combine in a 3D)
-    - filter SliceBySliceRelPosFilter + combine in a 3D 
-          - resampling, affine transfo to pair the slices
-          - extract list of images (ecrire utilisation de ExtractSliceFilter)
-   */
-
-  typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
-  typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
-  sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
-  sliceRelPosFilter->VerboseStepOn();
-  sliceRelPosFilter->WriteStepOn();
-  sliceRelPosFilter->SetInput(m_working_image);
-  sliceRelPosFilter->SetInputObject(temp);
-  sliceRelPosFilter->SetDirection(2);
-  sliceRelPosFilter->Update();
-  m_working_image = sliceRelPosFilter->GetOutput();
-  writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
+//--------------------------------------------------------------------
+template <class TImageType>
+bool
+clitk::ExtractLymphStationsFilter<TImageType>::
+GetComputeStation(std::string station) 
+{
+  return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
+}
+//--------------------------------------------------------------------
 
-  /*
-  
 
-  itk::ImageSliceConstIteratorWithIndex<ImageType> it(temp, temp->GetRequestedRegion());
-  itk::ImageRegionIterator<SliceType> * ot;
-  typename SliceType::Pointer slice;
-  it.SetFirstDirection(0);
-  it.SetSecondDirection(1);
-  it.GoToBegin();
-  typename SliceType::RegionType mSliceRegion;
-  typename SliceType::IndexType mSliceIndex;
-  typename SliceType::SizeType mSliceSize;
-  typename SliceType::SpacingType mSliceSpacing;
-  typename SliceType::PointType mSliceOrigin;
-  mSliceIndex[0] = temp->GetLargestPossibleRegion().GetIndex()[0];
-  mSliceIndex[1] = temp->GetLargestPossibleRegion().GetIndex()[1];
-  mSliceSize[0] = temp->GetLargestPossibleRegion().GetSize()[0];
-  mSliceSize[1] = temp->GetLargestPossibleRegion().GetSize()[1];
-  mSliceSpacing[0] = temp->GetSpacing()[0];
-  mSliceSpacing[1] = temp->GetSpacing()[1];
-  mSliceOrigin[0] = temp->GetOrigin()[0];
-  mSliceOrigin[1] = temp->GetOrigin()[1];
-  mSliceRegion.SetIndex(mSliceIndex);
-  mSliceRegion.SetSize(mSliceSize);
-  int i=0;
-  while( !it.IsAtEnd() ) {
-    // DD(i);
-    slice = SliceType::New();
-    slice->SetRegions(mSliceRegion);  
-    slice->SetOrigin(mSliceOrigin);
-    slice->SetSpacing(mSliceSpacing);
-    slice->Allocate();
-    ot = new itk::ImageRegionIterator<SliceType>(slice, slice->GetLargestPossibleRegion());
-    ot->GoToBegin();
-    // DD(it.GetIndex());
-    while(!it.IsAtEndOfSlice()) {
-      // DD(ot->GetIndex());
-      while(!it.IsAtEndOfLine()) {
-        ot->Set(it.Get());
-        ++it;
-        ++(*ot);
-      }
-      it.NextLine();
-    }
-    mImageSlices.push_back(slice);
-    it.NextSlice();
-    ++i;
-  } 
-  writeImage<SliceType>(mImageSlices[10], "slice.mhd");
-
-  // Perform RelPos by slice
-  for(int i=0; i<mImageSlices.size(); i++) {
-    DD(i);
-    typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
-    typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-    //    relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
-    relPosFilter->VerboseStepOff();
-    relPosFilter->WriteStepOff();
-    relPosFilter->SetInput(m_working_image); 
-    relPosFilter->SetInputObject(temp); 
-    relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
-    relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
-    relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
-    relPosFilter->Update();
-    m_working_image = relPosFilter->GetOutput();
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+AddComputeStation(std::string station) 
+{
+  m_ComputeStationMap[station] = true;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void 
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_8() 
+{
+  if (CheckForStation("8")) {
+    ExtractStation_8_SI_Limits();
+    ExtractStation_8_Post_Limits();
+    ExtractStation_8_Ant_Sup_Limits();
+    ExtractStation_8_Ant_Inf_Limits();
+    ExtractStation_8_LR_1_Limits();
+    ExtractStation_8_LR_2_Limits();
+    ExtractStation_8_LR_Limits();
+    ExtractStation_8_Ant_Injected_Limits();
+    ExtractStation_8_Single_CCL_Limits();
+    ExtractStation_8_Remove_Structures();
+    // Store image filenames into AFDB 
+    writeImage<MaskImageType>(m_ListOfStations["8"], "seg/Station8.mhd");
+    GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd");  
+    WriteAFDB();
   }
+}
+//--------------------------------------------------------------------
 
-  */
-  
-  DD("end");
-  m_output = m_working_image;
-  StopCurrentStep<ImageType>(m_output);
-
-  // Set output image information (required)
-  ImagePointer outputImage = this->GetOutput(0);
-  outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
-  outputImage->SetOrigin(m_working_image->GetOrigin());
-  outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
-  DD("end2");
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void 
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_3P()
+{
+  if (CheckForStation("3P")) {
+    ExtractStation_3P_SI_Limits();
+    ExtractStation_3P_Remove_Structures();
+    ExtractStation_3P_Ant_Limits();
+    ExtractStation_3P_Post_Limits();
+    ExtractStation_3P_LR_sup_Limits();
+    ExtractStation_3P_LR_inf_Limits();
+    // Store image filenames into AFDB 
+    writeImage<MaskImageType>(m_ListOfStations["3P"], "seg/Station3P.mhd");
+    GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); 
+    WriteAFDB(); 
+  }
 }
 //--------------------------------------------------------------------
 
@@ -350,15 +246,15 @@ GenerateOutputInformation() {
 template <class TImageType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
-GenerateInputRequestedRegion() {
-  DD("GenerateInputRequestedRegion");
-  // Call default
-  Superclass::GenerateInputRequestedRegion();
-  // Following needed because output region can be greater than input (trachea)
-  ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
-  ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
-  mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
-  trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
+ExtractStation_7() {
+  if (m_ListOfStations["7"]) {
+    DD("Station 7 support already exist -> use it");
+    m_Working_Support = m_ListOfStations["7"];
+  }
+  else m_Working_Support = m_Mediastinum;
+  ExtractStation_7_SI_Limits();
+  ExtractStation_7_RL_Limits();
+  ExtractStation_7_Posterior_Limits();
 }
 //--------------------------------------------------------------------
 
@@ -367,12 +263,19 @@ GenerateInputRequestedRegion() {
 template <class TImageType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
-GenerateData() {
-  DD("GenerateData");
-  // Final Step -> graft output (if SetNthOutput => redo)
-  this->GraftOutput(m_output);
+ExtractStation_4RL() {
+  /*
+    WARNING ONLY 4R FIRST !!! (not same inf limits)
+  */    
+  ExtractStation_4RL_SI_Limits();
+  ExtractStation_4RL_LR_Limits();
+  ExtractStation_4RL_AP_Limits();
 }
 //--------------------------------------------------------------------
-  
+
+
+//--------------------------------------------------------------------
+
+
 
 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX