1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
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"
32 #include <itkStatisticsLabelMapFilter.h>
33 #include <itkLabelImageToStatisticsLabelMapFilter.h>
34 #include <itkRegionOfInterestImageFilter.h>
35 #include <itkBinaryThresholdImageFilter.h>
36 #include <itkImageSliceConstIteratorWithIndex.h>
37 #include <itkBinaryThinningImageFilter.h>
40 #include "RelativePositionPropImageFilter.h"
42 //--------------------------------------------------------------------
43 template <class TImageType>
44 clitk::ExtractLymphStationsFilter<TImageType>::
45 ExtractLymphStationsFilter():
47 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
48 itk::ImageToImageFilter<TImageType, MaskImageType>()
50 this->SetNumberOfRequiredInputs(1);
51 SetBackgroundValue(0);
52 SetForegroundValue(1);
54 //--------------------------------------------------------------------
57 //--------------------------------------------------------------------
58 template <class TImageType>
60 clitk::ExtractLymphStationsFilter<TImageType>::
61 GenerateOutputInformation() {
62 // Superclass::GenerateOutputInformation();
66 m_mediastinum = GetAFDB()->template GetImage <MaskImageType>("mediastinum");
67 m_trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
68 ImagePointType carina;
69 GetAFDB()->GetPoint3D("carina", carina);
71 m_CarinaZPositionInMM = carina[2];
72 DD(m_CarinaZPositionInMM);
74 GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
75 m_MiddleLobeBronchusZPositionInMM = mlb[2];
76 DD(m_MiddleLobeBronchusZPositionInMM);
78 // ----------------------------------------------------------------
79 // ----------------------------------------------------------------
80 // Superior limit = carina
81 // Inferior limit = origin right middle lobe bronchus
82 StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
83 ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
84 ImageSizeType size = region.GetSize();
85 ImageIndexType index = region.GetIndex();
86 DD(m_CarinaZPositionInMM);
87 DD(m_MiddleLobeBronchusZPositionInMM);
88 index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
89 size[2] = 1+ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]); // +1 to
91 region.SetIndex(index);
93 typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> CropFilterType;
94 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
95 cropFilter->SetInput(m_mediastinum);
96 cropFilter->SetRegionOfInterest(region);
98 m_working_image = cropFilter->GetOutput();
99 // Auto Crop (because following rel pos is faster)
100 m_working_image = clitk::AutoCrop<MaskImageType>(m_working_image, 0);
101 StopCurrentStep<MaskImageType>(m_working_image);
102 m_output = m_working_image;
104 // ----------------------------------------------------------------
105 // ----------------------------------------------------------------
106 // Separate trachea in two CCL
107 StartNewStep("Separate trachea under carina");
109 ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
110 for(int i=0; i<3; i++) {
111 index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i]
112 -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]);
114 size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
120 if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) {
121 size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i];
126 region.SetIndex(index);
127 region.SetSize(size);
128 // typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> CropFilterType;
129 // typename CropFilterType::Pointer
130 cropFilter = CropFilterType::New();
131 // m_trachea.Print(std::cout);
132 cropFilter->SetInput(m_trachea);
133 cropFilter->SetRegionOfInterest(region);
134 cropFilter->Update();
135 m_working_trachea = cropFilter->GetOutput();
137 // Labelize and consider two main labels
138 m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
140 // Detect wich label is at Left
141 typedef itk::ImageSliceConstIteratorWithIndex<MaskImageType> SliceIteratorType;
142 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
143 iter.SetFirstDirection(0);
144 iter.SetSecondDirection(1);
147 ImagePixelType leftLabel;
148 ImagePixelType rightLabel;
150 if (iter.Get() != GetBackgroundValue()) {
151 // DD(iter.GetIndex());
152 // DD((int)iter.Get());
153 leftLabel = iter.Get();
158 if (leftLabel == 1) rightLabel = 2;
164 StopCurrentStep<MaskImageType>(m_working_trachea);
166 //-----------------------------------------------------
167 /* DD("TEST Skeleton");
168 typedef itk::BinaryThinningImageFilter<MaskImageType, MaskImageType> SkeletonFilterType;
169 typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New();
170 skeletonFilter->SetInput(m_working_trachea);
171 skeletonFilter->Update();
172 writeImage<MaskImageType>(skeletonFilter->GetOutput(), "skel.mhd");
173 writeImage<MaskImageType>(skeletonFilter->GetThinning(), "skel2.mhd");
176 //-----------------------------------------------------
177 StartNewStep("Left limits with bronchus (slice by slice)");
178 // Select LeftLabel (set right label to 0)
179 MaskImagePointer temp = SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, rightLabel, 0);
180 writeImage<MaskImageType>(temp, "temp1.mhd");
183 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_image,
185 2, GetFuzzyThreshold(), "RightTo");
187 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
188 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
189 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
190 sliceRelPosFilter->VerboseStepOn();
191 sliceRelPosFilter->WriteStepOn();
192 sliceRelPosFilter->SetInput(m_working_image);
193 sliceRelPosFilter->SetInputObject(temp);
194 sliceRelPosFilter->SetDirection(2);
195 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
196 sliceRelPosFilter->SetOrientationTypeString("RightTo");
197 sliceRelPosFilter->Update();
198 m_working_image = sliceRelPosFilter->GetOutput();*/
199 writeImage<MaskImageType>(m_working_image, "afterslicebyslice.mhd");
202 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
203 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
206 //-----------------------------------------------------
207 StartNewStep("Right limits with bronchus (slice by slice)");
208 // Select LeftLabel (set right label to 0)
209 temp = SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, leftLabel, 0);
210 writeImage<MaskImageType>(temp, "temp2.mhd");
212 sliceRelPosFilter = SliceRelPosFilterType::New();
213 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
214 sliceRelPosFilter->VerboseStepOff();
215 sliceRelPosFilter->WriteStepOff();
216 sliceRelPosFilter->SetInput(m_working_image);
217 sliceRelPosFilter->SetInputObject(temp);
218 sliceRelPosFilter->SetDirection(2);
219 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
220 sliceRelPosFilter->SetOrientationTypeString("LeftTo");
221 sliceRelPosFilter->Update();
222 m_working_image = sliceRelPosFilter->GetOutput();
223 writeImage<MaskImageType>(m_working_image, "afterslicebyslice.mhd");
225 m_output = m_working_image;
226 StopCurrentStep<MaskImageType>(m_output);
228 //-----------------------------------------------------
229 StartNewStep("Trial Post position according to trachea");
230 // Post: do not extend past the post wall of the 2 main bronchi
231 sliceRelPosFilter = SliceRelPosFilterType::New();
232 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
233 sliceRelPosFilter->VerboseStepOn();
234 sliceRelPosFilter->WriteStepOn();
235 sliceRelPosFilter->SetInput(m_output);
236 sliceRelPosFilter->SetInputObject(m_trachea);
237 sliceRelPosFilter->SetDirection(2);
238 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold());
239 // It means "not PostTo" (different from AntTo)
240 sliceRelPosFilter->NotFlagOn();
241 //sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::PostTo);
242 sliceRelPosFilter->SetOrientationTypeString("PostTo");
243 sliceRelPosFilter->Update();
244 m_output = sliceRelPosFilter->GetOutput();
245 writeImage<MaskImageType>(m_output, "postrelpos.mhd");
248 // Set output image information (required)
249 MaskImagePointer outputImage = this->GetOutput(0);
250 outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
251 outputImage->SetOrigin(m_working_image->GetOrigin());
252 outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
255 //--------------------------------------------------------------------
258 //--------------------------------------------------------------------
259 template <class TImageType>
261 clitk::ExtractLymphStationsFilter<TImageType>::
262 GenerateInputRequestedRegion() {
263 DD("GenerateInputRequestedRegion");
265 // Superclass::GenerateInputRequestedRegion();
266 // Following needed because output region can be greater than input (trachea)
267 //ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
268 //ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
270 m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
271 m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
274 //--------------------------------------------------------------------
277 //--------------------------------------------------------------------
278 template <class TImageType>
280 clitk::ExtractLymphStationsFilter<TImageType>::
283 // Final Step -> graft output (if SetNthOutput => redo)
284 this->GraftOutput(m_output);
286 //--------------------------------------------------------------------
289 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX