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 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 mediastinum 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 under carena");
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);
180 // Detect wich label is at Left
181 typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
182 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
183 iter.SetFirstDirection(0);
184 iter.SetSecondDirection(1);
187 ImagePixelType leftLabel;
188 ImagePixelType rightLabel;
190 if (iter.Get() != m_BackgroundValueTrachea) {
191 // DD(iter.GetIndex());
192 // DD((int)iter.Get());
193 leftLabel = iter.Get();
198 if (leftLabel == 1) rightLabel = 2;
204 StopCurrentStep<ImageType>(m_working_trachea);
206 //-----------------------------------------------------
207 /* DD("TEST Skeleton");
208 typedef itk::BinaryThinningImageFilter<ImageType, ImageType> SkeletonFilterType;
209 typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New();
210 skeletonFilter->SetInput(m_working_trachea);
211 skeletonFilter->Update();
212 writeImage<ImageType>(skeletonFilter->GetOutput(), "skel.mhd");
213 writeImage<ImageType>(skeletonFilter->GetThinning(), "skel2.mhd");
216 //-----------------------------------------------------
217 StartNewStep("Left limits with bronchus (slice by slice)");
218 // Select LeftLabel (set right label to 0)
219 ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, rightLabel, 0);
220 writeImage<ImageType>(temp, "temp1.mhd");
222 typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
223 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
224 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
225 sliceRelPosFilter->VerboseStepOn();
226 sliceRelPosFilter->WriteStepOn();
227 sliceRelPosFilter->SetInput(m_working_image);
228 sliceRelPosFilter->SetInputObject(temp);
229 sliceRelPosFilter->SetDirection(2);
230 sliceRelPosFilter->SetFuzzyThreshold(0.5);
231 sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::RightTo);
232 sliceRelPosFilter->Update();
233 m_working_image = sliceRelPosFilter->GetOutput();
234 writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
237 //-----------------------------------------------------
238 StartNewStep("Right limits with bronchus (slice by slice)");
239 // Select LeftLabel (set right label to 0)
240 temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, leftLabel, 0);
241 writeImage<ImageType>(temp, "temp2.mhd");
243 sliceRelPosFilter = SliceRelPosFilterType::New();
244 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
245 sliceRelPosFilter->VerboseStepOn();
246 sliceRelPosFilter->WriteStepOn();
247 sliceRelPosFilter->SetInput(m_working_image);
248 sliceRelPosFilter->SetInputObject(temp);
249 sliceRelPosFilter->SetDirection(2);
250 sliceRelPosFilter->SetFuzzyThreshold(0.5);
251 sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::LeftTo);
252 sliceRelPosFilter->Update();
253 m_working_image = sliceRelPosFilter->GetOutput();
254 writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
258 m_output = m_working_image;
259 StopCurrentStep<ImageType>(m_output);
261 // Set output image information (required)
262 ImagePointer outputImage = this->GetOutput(0);
263 outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
264 outputImage->SetOrigin(m_working_image->GetOrigin());
265 outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
268 //--------------------------------------------------------------------
271 //--------------------------------------------------------------------
272 template <class TImageType>
274 clitk::ExtractLymphStationsFilter<TImageType>::
275 GenerateInputRequestedRegion() {
276 DD("GenerateInputRequestedRegion");
278 Superclass::GenerateInputRequestedRegion();
279 // Following needed because output region can be greater than input (trachea)
280 ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
281 ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
282 mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
283 trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
285 //--------------------------------------------------------------------
288 //--------------------------------------------------------------------
289 template <class TImageType>
291 clitk::ExtractLymphStationsFilter<TImageType>::
294 // Final Step -> graft output (if SetNthOutput => redo)
295 this->GraftOutput(m_output);
297 //--------------------------------------------------------------------
300 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX