--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+ ======================================================================-====*/
+
+#ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
+#define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
+
+// clitk
+#include "clitkCommon.h"
+#include "clitkExtractLymphStationsFilter.h"
+#include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
+#include "clitkSegmentationUtils.h"
+#include "clitkAutoCropFilter.h"
+#include "clitkSegmentationUtils.h"
+#include "clitkSliceBySliceRelativePositionFilter.h"
+
+// itk
+#include <itkStatisticsLabelMapFilter.h>
+#include <itkLabelImageToStatisticsLabelMapFilter.h>
+#include <itkRegionOfInterestImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
+#include <itkImageSliceConstIteratorWithIndex.h>
+#include <itkBinaryThinningImageFilter.h>
+
+// itk ENST
+#include "RelativePositionPropImageFilter.h"
+
+//--------------------------------------------------------------------
+template <class TImageType>
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractLymphStationsFilter():
+ clitk::FilterBase(),
+ clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
+ itk::ImageToImageFilter<TImageType, MaskImageType>()
+{
+ this->SetNumberOfRequiredInputs(1);
+ SetBackgroundValue(0);
+ SetForegroundValue(1);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+GenerateOutputInformation() {
+ // Superclass::GenerateOutputInformation();
+
+ // Get input
+ LoadAFDB();
+ m_mediastinum = GetAFDB()->template GetImage <MaskImageType>("mediastinum");
+ m_trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
+ ImagePointType carina;
+ GetAFDB()->GetPoint3D("carina", carina);
+ DD(carina);
+ m_CarinaZPositionInMM = carina[2];
+ DD(m_CarinaZPositionInMM);
+ ImagePointType mlb;
+ GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
+ m_MiddleLobeBronchusZPositionInMM = mlb[2];
+ DD(m_MiddleLobeBronchusZPositionInMM);
+
+ // ----------------------------------------------------------------
+ // ----------------------------------------------------------------
+ // Superior limit = carina
+ // Inferior limit = origin right middle lobe bronchus
+ StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
+ ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
+ ImageSizeType size = region.GetSize();
+ ImageIndexType index = region.GetIndex();
+ DD(m_CarinaZPositionInMM);
+ DD(m_MiddleLobeBronchusZPositionInMM);
+ index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
+ size[2] = 1+ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]); // +1 to
+ region.SetSize(size);
+ region.SetIndex(index);
+ DD(region);
+ typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> 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<MaskImageType>(m_working_image, 0);
+ StopCurrentStep<MaskImageType>(m_working_image);
+ m_output = m_working_image;
+
+ // ----------------------------------------------------------------
+ // ----------------------------------------------------------------
+ // Separate trachea in two CCL
+ StartNewStep("Separate trachea under carina");
+ // 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];
+ }
+ }
+ // DD(index);
+ // DD(size);
+ region.SetIndex(index);
+ region.SetSize(size);
+ // typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> 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<MaskImageType>(m_working_trachea, 0, true, 1);
+
+ // Detect wich label is at Left
+ typedef itk::ImageSliceConstIteratorWithIndex<MaskImageType> SliceIteratorType;
+ SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
+ iter.SetFirstDirection(0);
+ iter.SetSecondDirection(1);
+ iter.GoToBegin();
+ bool stop = false;
+ ImagePixelType leftLabel;
+ ImagePixelType rightLabel;
+ while (!stop) {
+ if (iter.Get() != GetBackgroundValue()) {
+ // DD(iter.GetIndex());
+ // DD((int)iter.Get());
+ leftLabel = iter.Get();
+ stop = true;
+ }
+ ++iter;
+ }
+ if (leftLabel == 1) rightLabel = 2;
+ else rightLabel = 1;
+ DD((int)leftLabel);
+ DD((int)rightLabel);
+
+ // End step
+ StopCurrentStep<MaskImageType>(m_working_trachea);
+
+ //-----------------------------------------------------
+ /* DD("TEST Skeleton");
+ typedef itk::BinaryThinningImageFilter<MaskImageType, MaskImageType> SkeletonFilterType;
+ typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New();
+ skeletonFilter->SetInput(m_working_trachea);
+ skeletonFilter->Update();
+ writeImage<MaskImageType>(skeletonFilter->GetOutput(), "skel.mhd");
+ writeImage<MaskImageType>(skeletonFilter->GetThinning(), "skel2.mhd");
+ */
+
+ //-----------------------------------------------------
+ StartNewStep("Left limits with bronchus (slice by slice)");
+ // Select LeftLabel (set right label to 0)
+ MaskImagePointer temp = SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, rightLabel, 0);
+ writeImage<MaskImageType>(temp, "temp1.mhd");
+
+ m_working_image =
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_image,
+ temp,
+ 2, GetFuzzyThreshold(), "RightTo");
+ /*
+ typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> 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->SetFuzzyThreshold(GetFuzzyThreshold());
+ sliceRelPosFilter->SetOrientationTypeString("RightTo");
+ sliceRelPosFilter->Update();
+ m_working_image = sliceRelPosFilter->GetOutput();*/
+ writeImage<MaskImageType>(m_working_image, "afterslicebyslice.mhd");
+
+
+ typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
+ typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
+
+
+ //-----------------------------------------------------
+ StartNewStep("Right limits with bronchus (slice by slice)");
+ // Select LeftLabel (set right label to 0)
+ temp = SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, leftLabel, 0);
+ writeImage<MaskImageType>(temp, "temp2.mhd");
+
+ sliceRelPosFilter = SliceRelPosFilterType::New();
+ sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
+ sliceRelPosFilter->VerboseStepOff();
+ sliceRelPosFilter->WriteStepOff();
+ sliceRelPosFilter->SetInput(m_working_image);
+ sliceRelPosFilter->SetInputObject(temp);
+ sliceRelPosFilter->SetDirection(2);
+ sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
+ sliceRelPosFilter->SetOrientationTypeString("LeftTo");
+ sliceRelPosFilter->Update();
+ m_working_image = sliceRelPosFilter->GetOutput();
+ writeImage<MaskImageType>(m_working_image, "afterslicebyslice.mhd");
+ DD("end");
+ m_output = m_working_image;
+ StopCurrentStep<MaskImageType>(m_output);
+
+ //-----------------------------------------------------
+ StartNewStep("Trial Post position according to trachea");
+ // Post: do not extend past the post wall of the 2 main bronchi
+ sliceRelPosFilter = SliceRelPosFilterType::New();
+ sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
+ sliceRelPosFilter->VerboseStepOn();
+ sliceRelPosFilter->WriteStepOn();
+ sliceRelPosFilter->SetInput(m_output);
+ sliceRelPosFilter->SetInputObject(m_trachea);
+ sliceRelPosFilter->SetDirection(2);
+ sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
+ // It means "not PostTo" (different from AntTo)
+ sliceRelPosFilter->NotFlagOn();
+ //sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::PostTo);
+ sliceRelPosFilter->SetOrientationTypeString("PostTo");
+ sliceRelPosFilter->Update();
+ m_output = sliceRelPosFilter->GetOutput();
+ writeImage<MaskImageType>(m_output, "postrelpos.mhd");
+
+
+ // Set output image information (required)
+ MaskImagePointer 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>::
+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));
+ DD("set reg");
+ m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
+ m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
+ DD("end");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+GenerateData() {
+ DD("GenerateData");
+ // Final Step -> graft output (if SetNthOutput => redo)
+ this->GraftOutput(m_output);
+}
+//--------------------------------------------------------------------
+
+
+#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff