separate airway tracking from extract lung

[clitk.git] / itk / clitkSliceBySliceRelativePositionFilter.txx

diff --git a/itk/clitkSliceBySliceRelativePositionFilter.txx b/itk/clitkSliceBySliceRelativePositionFilter.txx
index af9c4c5c35f3378e5d0c6e90c1da399c224acfab..a8a79ec85599d681518edfe869da79598fcbec60 100644 (file)
--- a/itk/clitkSliceBySliceRelativePositionFilter.txx
+++ b/itk/clitkSliceBySliceRelativePositionFilter.txx
@@ -18,24 +18,8 @@
 
 // clitk
 #include "clitkSegmentationUtils.h"
-// #include "clitkBooleanOperatorLabelImageFilter.h"
-// #include "clitkAutoCropFilter.h"
-// #include "clitkResampleImageWithOptionsFilter.h"
-// #include "clitkBooleanOperatorLabelImageFilter.h"
 #include "clitkExtractSliceFilter.h"
 
-// // itk
-// #include <deque>
-// #include "itkStatisticsLabelMapFilter.h"
-// #include "itkLabelImageToStatisticsLabelMapFilter.h"
-// #include "itkRegionOfInterestImageFilter.h"
-// #include "itkBinaryThresholdImageFilter.h"
-// #include "itkBinaryErodeImageFilter.h"
-// #include "itkBinaryBallStructuringElement.h"
-
-// // itk [Bloch et al] 
-// #include "RelativePositionPropImageFilter.h"
-
 //--------------------------------------------------------------------
 template <class ImageType>
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
@@ -45,7 +29,11 @@ SliceBySliceRelativePositionFilter():
 {
   this->SetNumberOfRequiredInputs(2);
   SetDirection(2);
-  SetObjectBackgroundValue(0);
+  SetObjectBackgroundValue(0);  
+  SetFuzzyThreshold(0.6);
+  SetOrientationType(RelPosFilterType::LeftTo);
+  SetIntermediateSpacing(10);
+  ResampleBeforeRelativePositionFilterOff();
 }
 //--------------------------------------------------------------------
 
@@ -54,7 +42,8 @@ SliceBySliceRelativePositionFilter():
 template <class ImageType>
 void 
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
-SetInput(const ImageType * image) {
+SetInput(const ImageType * image) 
+{
   // Process object is not const-correct so the const casting is required.
   this->SetNthInput(0, const_cast<ImageType *>(image));
 }
@@ -65,7 +54,8 @@ SetInput(const ImageType * image) {
 template <class ImageType>
 void 
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
-SetInputObject(const ImageType * image) {
+SetInputObject(const ImageType * image) 
+{
   // Process object is not const-correct so the const casting is required.
   this->SetNthInput(1, const_cast<ImageType *>(image));
 }
@@ -76,7 +66,8 @@ SetInputObject(const ImageType * image) {
 template <class ImageType>
 void 
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
-GenerateOutputInformation() { 
+GenerateOutputInformation() 
+{ 
   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   ImagePointer outputImage = this->GetOutput(0);
   outputImage->SetRegions(input->GetLargestPossibleRegion());
@@ -88,7 +79,8 @@ GenerateOutputInformation() {
 template <class ImageType>
 void 
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
-GenerateInputRequestedRegion() {
+GenerateInputRequestedRegion() 
+{
   // Call default
   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
   // Get input pointers and set requested region to common region
@@ -104,7 +96,8 @@ GenerateInputRequestedRegion() {
 template <class ImageType>
 void 
 clitk::SliceBySliceRelativePositionFilter<ImageType>::
-GenerateData() {
+GenerateData() 
+{
   // Get input pointer
   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
@@ -133,66 +126,102 @@ GenerateData() {
   else {
     DD("no pad");
   }
-  
+
   /*
+    - extract vector of slices in input, in object
+    - slice by slice rel position
+    - joint result
+    - post process
+  */
 
 
   //--------------------------------------------------------------------
   // Extract input slices
   StartNewStep("Extract input slices");
   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
-  typename ExtractSliceFilterType::Pointer extractSliceFilter = Extractslicefilter::New();
+  typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
   extractSliceFilter->SetInput(input);
   extractSliceFilter->SetDirection(GetDirection());
   extractSliceFilter->Update();
   typedef typename ExtractSliceFilterType::SliceType SliceType;
-  std::vector<typename SliceType::Pointer> & mInputSlices = extractSliceFilter->GetOutput();
-  DD(mInputSlices.size());
-  StopCurrentStep<SliceType>(mInputSlices[5]);
+  std::vector<typename SliceType::Pointer> mInputSlices;
+  extractSliceFilter->GetOutputSlices(mInputSlices);
+  StopCurrentStep<SliceType>(mInputSlices[0]);
   
   //--------------------------------------------------------------------
   // Extract object slices
-
   StartNewStep("Extract object slices");
-  extractSliceFilter = Extractslicefilter::New();
-  extractSliceFilter->SetInput(input);
+  extractSliceFilter = ExtractSliceFilterType::New();
+  extractSliceFilter->SetInput(object);
   extractSliceFilter->SetDirection(GetDirection());
   extractSliceFilter->Update();
-  std::vector<typename SliceType::Pointer> & mObjectSlices = extractSliceFilter->GetOutput();
-  DD(mObjectSlices.size());
-  StopCurrentStep<SliceType>(mInputSlices[5]);
+  std::vector<typename SliceType::Pointer> mObjectSlices;
+  extractSliceFilter->GetOutputSlices(mObjectSlices);
+  StopCurrentStep<SliceType>(mObjectSlices[0]);
 
-  
   //--------------------------------------------------------------------
   // Perform slice by slice relative position
   StartNewStep("Perform slice by slice relative position");
-  for(int i=0; i<mInputSlices.size(); i++) {
-    DD(i);
+  for(unsigned int i=0; i<mInputSlices.size(); i++) {
+    // DD(i);
+    //     DD(mInputSlices[i]->GetOrigin());
+    //     writeImage<SliceType>(mInputSlices[i], "inp"+clitk::toString(i)+".mhd");
+
+    // Select main CC in each object slice : this should be the main bronchus
+    mObjectSlices[i] = Labelize<SliceType>(mObjectSlices[i], 0, true, 1);
+    mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
+
+    // Relative position
     typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
     typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
     relPosFilter->VerboseStepOff();
     relPosFilter->WriteStepOff();
+    relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
     relPosFilter->SetInput(mInputSlices[i]); 
     relPosFilter->SetInputObject(mObjectSlices[i]); 
-    relPosFilter->SetOrientationType(GetOrientationType());
-    relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
-    relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
+    relPosFilter->SetOrientationType(this->GetOrientationType());
+    relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
+    relPosFilter->SetResampleBeforeRelativePositionFilter(this->GetResampleBeforeRelativePositionFilter());
+    relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
+    relPosFilter->AutoCropOff(); // important ! because we join the slices after this loop
     relPosFilter->Update();
+    // writeImage<SliceType>(relPosFilter->GetOutput(), "inp-after"+clitk::toString(i)+".mhd");
     mInputSlices[i] = relPosFilter->GetOutput();
   }
-  typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JointSeriesFilterType;
-  typename JointSeriesFilterType::Pointer jointFilter = JointSeriesFilterType::New();
-  for(int i=0; i<mInputSlices.size(); i++) {
-    jointFilter->SetInput(i, mInputSlices[i]);
+  DD(this->GetIntermediateSpacing());
+  DD(this->GetResampleBeforeRelativePositionFilter());
+  DD("End slice");
+
+  typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JoinSeriesFilterType;
+  typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
+  joinFilter->SetOrigin(input->GetOrigin()[GetDirection()]);
+  joinFilter->SetSpacing(input->GetSpacing()[GetDirection()]);
+  for(unsigned int i=0; i<mInputSlices.size(); i++) {
+  // DD(mInputSlices[i]->GetLargestPossibleRegion().GetIndex());
+//   DD(mInputSlices[i]->GetLargestPossibleRegion().GetSize());
+//   DD(mInputSlices[i]->GetRequestedRegion().GetIndex());
+//   DD(mInputSlices[i]->GetRequestedRegion().GetSize());
+    joinFilter->PushBackInput(mInputSlices[i]);
+    //SetInput(i, mInputSlices[i]);
   }
-  m_working_input = jointFilter->GetOutput();
+  DD("before update");
+  joinFilter->Update();
+  DD("after update");
+  m_working_input = joinFilter->GetOutput();
+  
+  // Update the origin
+  DD(input->GetSpacing());
+  DD(input->GetOrigin());
+  DD(mInputSlices[0]->GetSpacing());
+  DD(mInputSlices[0]->GetOrigin());
   DD(m_working_input->GetSpacing());
   DD(m_working_input->GetOrigin());
+  // typename ImageType::PointType origin = m_working_input->GetOrigin();
+//   origin[GetDirection()] = input->GetOrigin()[GetDirection()];
+//   m_working_input->SetOrigin(origin);
+//   DD(m_working_input->GetOrigin());
   StopCurrentStep<ImageType>(m_working_input);
 
-  */
-
-
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------  
   // Final Step -> set output