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"
33 #include <itkStatisticsLabelMapFilter.h>
34 #include <itkLabelImageToStatisticsLabelMapFilter.h>
35 #include <itkRegionOfInterestImageFilter.h>
36 #include <itkBinaryThresholdImageFilter.h>
37 #include <itkImageSliceConstIteratorWithIndex.h>
40 #include "RelativePositionPropImageFilter.h"
42 //--------------------------------------------------------------------
43 template <class TImageType>
44 clitk::ExtractLymphStationsFilter<TImageType>::
45 ExtractLymphStationsFilter():
47 itk::ImageToImageFilter<TImageType, TImageType>()
49 this->SetNumberOfRequiredInputs(1);
50 SetBackgroundValueMediastinum(0);
51 SetBackgroundValueTrachea(0);
52 SetBackgroundValue(0);
53 SetForegroundValue(1);
55 SetIntermediateSpacing(6);
56 SetFuzzyThreshold1(0.6);
58 //--------------------------------------------------------------------
61 //--------------------------------------------------------------------
62 template <class TImageType>
64 clitk::ExtractLymphStationsFilter<TImageType>::
65 SetInputMediastinumLabelImage(const TImageType * image, ImagePixelType bg) {
66 this->SetNthInput(0, const_cast<TImageType *>(image));
67 m_BackgroundValueMediastinum = bg;
68 SetCarenaZPositionInMM(image->GetOrigin()[2]+image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]);
69 SetMiddleLobeBronchusZPositionInMM(image->GetOrigin()[2]);
70 // DD(m_CarenaZPositionInMM);
71 // DD(m_MiddleLobeBronchusZPositionInMM);
73 //--------------------------------------------------------------------
76 //--------------------------------------------------------------------
77 template <class TImageType>
79 clitk::ExtractLymphStationsFilter<TImageType>::
80 SetInputTracheaLabelImage(const TImageType * image, ImagePixelType bg) {
81 this->SetNthInput(1, const_cast<TImageType *>(image));
82 m_BackgroundValueTrachea = bg;
84 //--------------------------------------------------------------------
87 //--------------------------------------------------------------------
88 template <class TImageType>
89 template<class ArgsInfoType>
91 clitk::ExtractLymphStationsFilter<TImageType>::
92 SetArgsInfo(ArgsInfoType mArgsInfo)
94 SetVerboseOption_GGO(mArgsInfo);
95 SetVerboseStep_GGO(mArgsInfo);
96 SetWriteStep_GGO(mArgsInfo);
97 SetVerboseWarningOff_GGO(mArgsInfo);
98 SetCarenaZPositionInMM_GGO(mArgsInfo);
99 SetMiddleLobeBronchusZPositionInMM_GGO(mArgsInfo);
100 SetIntermediateSpacing_GGO(mArgsInfo);
101 SetFuzzyThreshold1_GGO(mArgsInfo);
102 //SetBackgroundValueMediastinum_GGO(mArgsInfo);
104 //--------------------------------------------------------------------
107 //--------------------------------------------------------------------
108 template <class TImageType>
110 clitk::ExtractLymphStationsFilter<TImageType>::
111 GenerateOutputInformation() {
112 // Superclass::GenerateOutputInformation();
115 m_mediastinum = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
116 m_trachea = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
118 // ----------------------------------------------------------------
119 // ----------------------------------------------------------------
120 // Superior limit = carena
121 // Inferior limit = origine middle lobe bronchus
122 StartNewStep("Inf/Sup limits with carena/bronchus");
123 ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
124 ImageSizeType size = region.GetSize();
125 ImageIndexType index = region.GetIndex();
126 DD(m_CarenaZPositionInMM);
127 DD(m_MiddleLobeBronchusZPositionInMM);
128 index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
129 size[2] = ceil((m_CarenaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]);
130 region.SetSize(size);
131 region.SetIndex(index);
133 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
134 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
135 cropFilter->SetInput(m_mediastinum);
136 cropFilter->SetRegionOfInterest(region);
137 cropFilter->Update();
138 m_working_image = cropFilter->GetOutput();
139 // Auto Crop (because following rel pos is faster)
140 m_working_image = clitk::AutoCrop<ImageType>(m_working_image, 0);
141 StopCurrentStep<ImageType>(m_working_image);
142 m_output = m_working_image;
144 // ----------------------------------------------------------------
145 // ----------------------------------------------------------------
146 // Separate trachea in two CCL
147 StartNewStep("Separate trachea");
149 ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
150 for(int i=0; i<3; i++) {
151 index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i]
152 -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]);
154 size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
160 if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) {
161 size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i];
166 region.SetIndex(index);
167 region.SetSize(size);
168 // typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
169 // typename CropFilterType::Pointer
170 cropFilter = CropFilterType::New();
171 // m_trachea.Print(std::cout);
172 cropFilter->SetInput(m_trachea);
173 cropFilter->SetRegionOfInterest(region);
174 cropFilter->Update();
175 m_working_trachea = cropFilter->GetOutput();
177 // Labelize and consider two main labels
178 m_working_trachea = Labelize<ImageType>(m_working_trachea, 0, true, 1);
179 StopCurrentStep<ImageType>(m_working_trachea);
181 // Detect wich label is at Left
182 typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
183 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
184 iter.SetFirstDirection(0);
185 iter.SetSecondDirection(1);
188 ImagePixelType leftLabel;
190 if (iter.Get() != m_BackgroundValueTrachea) {
191 // DD(iter.GetIndex());
192 // DD((int)iter.Get());
193 leftLabel = iter.Get();
202 StartNewStep("Left/Right limits with trachea");
204 // Select LeftLabel (set label 1 to 0)
205 ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
206 writeImage<ImageType>(temp, "temp1.mhd");
208 // Left relative position
209 typedef clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType> RelPosFilterType;
210 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
211 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
212 relPosFilter->VerboseStepOff();
213 relPosFilter->WriteStepOff();
214 relPosFilter->SetInput(m_working_image);
215 relPosFilter->SetInputObject(temp);
216 relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
217 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
218 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
219 relPosFilter->Update();
220 m_working_image = relPosFilter->GetOutput();
222 // Select RightLabel (set label 2 to 0)
223 temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 2, 0);
224 writeImage<ImageType>(temp, "temp2.mhd");
226 // Left relative position
227 relPosFilter = RelPosFilterType::New();
228 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
229 relPosFilter->VerboseStepOn();
230 relPosFilter->WriteStepOn();
231 relPosFilter->SetInput(m_working_image);
232 relPosFilter->SetInputObject(temp);
233 relPosFilter->SetOrientationType(RelPosFilterType::LeftTo);
234 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
235 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
236 relPosFilter->Update();
237 m_working_image = relPosFilter->GetOutput();
240 //-----------------------------------------------------
241 StartNewStep("Left/Right limits with trachea (slice by slice");
243 // Select LeftLabel (set label 1 to 0)
244 ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, 1, 0);
245 writeImage<ImageType>(temp, "temp1.mhd");
248 - écrire utilisation de filter SliceBySliceRelPosFilter (combine in a 3D)
249 - filter SliceBySliceRelPosFilter + combine in a 3D
250 - resampling, affine transfo to pair the slices
251 - extract list of images (ecrire utilisation de ExtractSliceFilter)
254 typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
255 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
256 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
257 sliceRelPosFilter->VerboseStepOn();
258 sliceRelPosFilter->WriteStepOn();
259 sliceRelPosFilter->SetInput(m_working_image);
260 sliceRelPosFilter->SetInputObject(temp);
261 sliceRelPosFilter->SetDirection(2);
262 sliceRelPosFilter->Update();
263 m_working_image = sliceRelPosFilter->GetOutput();
264 writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
269 itk::ImageSliceConstIteratorWithIndex<ImageType> it(temp, temp->GetRequestedRegion());
270 itk::ImageRegionIterator<SliceType> * ot;
271 typename SliceType::Pointer slice;
272 it.SetFirstDirection(0);
273 it.SetSecondDirection(1);
275 typename SliceType::RegionType mSliceRegion;
276 typename SliceType::IndexType mSliceIndex;
277 typename SliceType::SizeType mSliceSize;
278 typename SliceType::SpacingType mSliceSpacing;
279 typename SliceType::PointType mSliceOrigin;
280 mSliceIndex[0] = temp->GetLargestPossibleRegion().GetIndex()[0];
281 mSliceIndex[1] = temp->GetLargestPossibleRegion().GetIndex()[1];
282 mSliceSize[0] = temp->GetLargestPossibleRegion().GetSize()[0];
283 mSliceSize[1] = temp->GetLargestPossibleRegion().GetSize()[1];
284 mSliceSpacing[0] = temp->GetSpacing()[0];
285 mSliceSpacing[1] = temp->GetSpacing()[1];
286 mSliceOrigin[0] = temp->GetOrigin()[0];
287 mSliceOrigin[1] = temp->GetOrigin()[1];
288 mSliceRegion.SetIndex(mSliceIndex);
289 mSliceRegion.SetSize(mSliceSize);
291 while( !it.IsAtEnd() ) {
293 slice = SliceType::New();
294 slice->SetRegions(mSliceRegion);
295 slice->SetOrigin(mSliceOrigin);
296 slice->SetSpacing(mSliceSpacing);
298 ot = new itk::ImageRegionIterator<SliceType>(slice, slice->GetLargestPossibleRegion());
300 // DD(it.GetIndex());
301 while(!it.IsAtEndOfSlice()) {
302 // DD(ot->GetIndex());
303 while(!it.IsAtEndOfLine()) {
310 mImageSlices.push_back(slice);
314 writeImage<SliceType>(mImageSlices[10], "slice.mhd");
316 // Perform RelPos by slice
317 for(int i=0; i<mImageSlices.size(); i++) {
319 typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
320 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
321 // relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
322 relPosFilter->VerboseStepOff();
323 relPosFilter->WriteStepOff();
324 relPosFilter->SetInput(m_working_image);
325 relPosFilter->SetInputObject(temp);
326 relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
327 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
328 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1());
329 relPosFilter->Update();
330 m_working_image = relPosFilter->GetOutput();
336 m_output = m_working_image;
337 StopCurrentStep<ImageType>(m_output);
339 // Set output image information (required)
340 ImagePointer outputImage = this->GetOutput(0);
341 outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
342 outputImage->SetOrigin(m_working_image->GetOrigin());
343 outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
346 //--------------------------------------------------------------------
349 //--------------------------------------------------------------------
350 template <class TImageType>
352 clitk::ExtractLymphStationsFilter<TImageType>::
353 GenerateInputRequestedRegion() {
354 DD("GenerateInputRequestedRegion");
356 Superclass::GenerateInputRequestedRegion();
357 // Following needed because output region can be greater than input (trachea)
358 ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
359 ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
360 mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
361 trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
363 //--------------------------------------------------------------------
366 //--------------------------------------------------------------------
367 template <class TImageType>
369 clitk::ExtractLymphStationsFilter<TImageType>::
372 // Final Step -> graft output (if SetNthOutput => redo)
373 this->GraftOutput(m_output);
375 //--------------------------------------------------------------------
378 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX