]> Creatis software - clitk.git/blob - segmentation/clitkExtractLymphStationsFilter.txx.save
Add one slice for the S1RL support
[clitk.git] / segmentation / clitkExtractLymphStationsFilter.txx.save
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 <itkBinaryThinningImageFilter.h>
38
39 // itk ENST
40 #include "RelativePositionPropImageFilter.h"
41
42 //--------------------------------------------------------------------
43 template <class TImageType>
44 clitk::ExtractLymphStationsFilter<TImageType>::
45 ExtractLymphStationsFilter():
46   clitk::FilterBase(),
47   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
48   itk::ImageToImageFilter<TImageType, MaskImageType>()
49 {
50   this->SetNumberOfRequiredInputs(1);
51   SetBackgroundValue(0);
52   SetForegroundValue(1);
53 }
54 //--------------------------------------------------------------------
55
56
57 //--------------------------------------------------------------------
58 template <class TImageType>
59 void 
60 clitk::ExtractLymphStationsFilter<TImageType>::
61 GenerateOutputInformation() { 
62   //  Superclass::GenerateOutputInformation();
63   
64   // Get input
65   LoadAFDB();
66   m_mediastinum = GetAFDB()->template GetImage <MaskImageType>("mediastinum");  
67   m_trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");  
68   ImagePointType carina;
69   GetAFDB()->GetPoint3D("carina", carina);
70   DD(carina);
71   m_CarinaZPositionInMM = carina[2];
72   DD(m_CarinaZPositionInMM);
73   ImagePointType mlb;
74   GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
75   m_MiddleLobeBronchusZPositionInMM = mlb[2];
76   DD(m_MiddleLobeBronchusZPositionInMM);
77
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 
90   region.SetSize(size);
91   region.SetIndex(index);
92   DD(region);
93   typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> CropFilterType;
94   typename CropFilterType::Pointer cropFilter = CropFilterType::New();
95   cropFilter->SetInput(m_mediastinum);
96   cropFilter->SetRegionOfInterest(region);
97   cropFilter->Update();
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;
103
104   // ----------------------------------------------------------------
105   // ----------------------------------------------------------------
106   // Separate trachea in two CCL
107   StartNewStep("Separate trachea under carina");
108   // DD(region);
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]);
113     // DD(index[i]);
114     size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
115     //  DD(size[i]);
116     if (index[i] < 0) { 
117       size[i] += index[i];
118       index[i] = 0;       
119     }
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];
122     }
123   }
124   // DD(index);
125   //   DD(size);
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();
136
137   // Labelize and consider two main labels
138   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
139
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);
145   iter.GoToBegin();
146   bool stop = false;
147   ImagePixelType leftLabel;
148   ImagePixelType rightLabel;
149   while (!stop) {
150     if (iter.Get() != GetBackgroundValue()) {
151       //     DD(iter.GetIndex());
152       //       DD((int)iter.Get());
153       leftLabel = iter.Get();
154       stop = true;
155     }
156     ++iter;
157   }
158   if (leftLabel == 1) rightLabel = 2;
159   else rightLabel = 1;
160   DD((int)leftLabel);
161   DD((int)rightLabel);  
162
163   // End step
164   StopCurrentStep<MaskImageType>(m_working_trachea);
165   
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");  
174   */
175
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");
181
182   m_working_image = 
183     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_image, 
184                                                        temp, 
185                                                        2, GetFuzzyThreshold(), "RightTo");
186   /*
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");
200
201
202   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
203   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
204
205
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");
211
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");
224   DD("end");
225   m_output = m_working_image;
226   StopCurrentStep<MaskImageType>(m_output);
227
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");
246
247
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());
253   DD("end2");
254 }
255 //--------------------------------------------------------------------
256
257
258 //--------------------------------------------------------------------
259 template <class TImageType>
260 void 
261 clitk::ExtractLymphStationsFilter<TImageType>::
262 GenerateInputRequestedRegion() {
263   DD("GenerateInputRequestedRegion");
264   // Call default
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));
269   DD("set reg");
270   m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
271   m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
272   DD("end");
273 }
274 //--------------------------------------------------------------------
275
276
277 //--------------------------------------------------------------------
278 template <class TImageType>
279 void 
280 clitk::ExtractLymphStationsFilter<TImageType>::
281 GenerateData() {
282   DD("GenerateData");
283   // Final Step -> graft output (if SetNthOutput => redo)
284   this->GraftOutput(m_output);
285 }
286 //--------------------------------------------------------------------
287   
288
289 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX