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, TImageType>()
50 this->SetNumberOfRequiredInputs(1);
51 SetBackgroundValueMediastinum(0);
52 SetBackgroundValueTrachea(0);
53 SetBackgroundValue(0);
54 SetForegroundValue(1);
55 SetIntermediateSpacing(6);
56 SetFuzzyThreshold1(0.6);
57 m_CarinaZPositionInMMIsSet = false;
59 //--------------------------------------------------------------------
62 //--------------------------------------------------------------------
63 template <class TImageType>
65 clitk::ExtractLymphStationsFilter<TImageType>::
66 SetInputMediastinumLabelImage(const TImageType * image, ImagePixelType bg) {
67 this->SetNthInput(0, const_cast<TImageType *>(image));
68 SetBackgroundValueMediastinum(bg);
70 //--------------------------------------------------------------------
73 //--------------------------------------------------------------------
74 template <class TImageType>
76 clitk::ExtractLymphStationsFilter<TImageType>::
77 SetInputTracheaLabelImage(const TImageType * image, ImagePixelType bg) {
78 this->SetNthInput(1, const_cast<TImageType *>(image));
79 SetBackgroundValueTrachea(bg);
81 //--------------------------------------------------------------------
84 //--------------------------------------------------------------------
85 template <class TImageType>
86 template<class ArgsInfoType>
88 clitk::ExtractLymphStationsFilter<TImageType>::
89 SetArgsInfo(ArgsInfoType mArgsInfo)
91 SetVerboseOption_GGO(mArgsInfo);
92 SetVerboseStep_GGO(mArgsInfo);
93 SetWriteStep_GGO(mArgsInfo);
94 SetVerboseWarningOff_GGO(mArgsInfo);
95 SetAFDBFilename_GGO(mArgsInfo);
96 SetCarinaZPositionInMM_GGO(mArgsInfo);
97 m_CarinaZPositionInMMIsSet = false;
98 SetMiddleLobeBronchusZPositionInMM_GGO(mArgsInfo);
99 SetIntermediateSpacing_GGO(mArgsInfo);
100 SetFuzzyThreshold1_GGO(mArgsInfo);
102 //--------------------------------------------------------------------
105 //--------------------------------------------------------------------
106 template <class TImageType>
108 clitk::ExtractLymphStationsFilter<TImageType>::
109 GenerateOutputInformation() {
110 // Superclass::GenerateOutputInformation();
113 m_mediastinum = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
114 m_trachea = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
116 // Get landmarks information
117 if (!m_CarinaZPositionInMMIsSet) {
118 ImagePointType carina;
120 GetAFDB()->GetPoint3D("Carina", carina);
122 m_CarinaZPositionInMM = carina[2];
124 DD(m_CarinaZPositionInMM);
126 // ----------------------------------------------------------------
127 // ----------------------------------------------------------------
128 // Superior limit = carina
129 // Inferior limit = origin right middle lobe bronchus
130 StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
131 ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region);
132 ImageSizeType size = region.GetSize();
133 ImageIndexType index = region.GetIndex();
134 DD(m_CarinaZPositionInMM);
135 DD(m_MiddleLobeBronchusZPositionInMM);
136 index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]);
137 size[2] = 1+ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]); // +1 to
138 region.SetSize(size);
139 region.SetIndex(index);
141 typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
142 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
143 cropFilter->SetInput(m_mediastinum);
144 cropFilter->SetRegionOfInterest(region);
145 cropFilter->Update();
146 m_working_image = cropFilter->GetOutput();
147 // Auto Crop (because following rel pos is faster)
148 m_working_image = clitk::AutoCrop<ImageType>(m_working_image, 0);
149 StopCurrentStep<ImageType>(m_working_image);
150 m_output = m_working_image;
152 // ----------------------------------------------------------------
153 // ----------------------------------------------------------------
154 // Separate trachea in two CCL
155 StartNewStep("Separate trachea under carina");
157 ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion();
158 for(int i=0; i<3; i++) {
159 index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i]
160 -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]);
162 size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]);
168 if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) {
169 size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i];
174 region.SetIndex(index);
175 region.SetSize(size);
176 // typedef itk::RegionOfInterestImageFilter<ImageType, ImageType> CropFilterType;
177 // typename CropFilterType::Pointer
178 cropFilter = CropFilterType::New();
179 // m_trachea.Print(std::cout);
180 cropFilter->SetInput(m_trachea);
181 cropFilter->SetRegionOfInterest(region);
182 cropFilter->Update();
183 m_working_trachea = cropFilter->GetOutput();
185 // Labelize and consider two main labels
186 m_working_trachea = Labelize<ImageType>(m_working_trachea, 0, true, 1);
188 // Detect wich label is at Left
189 typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
190 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
191 iter.SetFirstDirection(0);
192 iter.SetSecondDirection(1);
195 ImagePixelType leftLabel;
196 ImagePixelType rightLabel;
198 if (iter.Get() != m_BackgroundValueTrachea) {
199 // DD(iter.GetIndex());
200 // DD((int)iter.Get());
201 leftLabel = iter.Get();
206 if (leftLabel == 1) rightLabel = 2;
212 StopCurrentStep<ImageType>(m_working_trachea);
214 //-----------------------------------------------------
215 /* DD("TEST Skeleton");
216 typedef itk::BinaryThinningImageFilter<ImageType, ImageType> SkeletonFilterType;
217 typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New();
218 skeletonFilter->SetInput(m_working_trachea);
219 skeletonFilter->Update();
220 writeImage<ImageType>(skeletonFilter->GetOutput(), "skel.mhd");
221 writeImage<ImageType>(skeletonFilter->GetThinning(), "skel2.mhd");
224 //-----------------------------------------------------
225 StartNewStep("Left limits with bronchus (slice by slice)");
226 // Select LeftLabel (set right label to 0)
227 ImagePointer temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, rightLabel, 0);
228 writeImage<ImageType>(temp, "temp1.mhd");
230 typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
231 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
232 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
233 sliceRelPosFilter->VerboseStepOn();
234 sliceRelPosFilter->WriteStepOn();
235 sliceRelPosFilter->SetInput(m_working_image);
236 sliceRelPosFilter->SetInputObject(temp);
237 sliceRelPosFilter->SetDirection(2);
238 sliceRelPosFilter->SetFuzzyThreshold(0.6);
239 sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::RightTo);
240 sliceRelPosFilter->Update();
241 m_working_image = sliceRelPosFilter->GetOutput();
242 writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
245 //-----------------------------------------------------
246 StartNewStep("Right limits with bronchus (slice by slice)");
247 // Select LeftLabel (set right label to 0)
248 temp = SetBackground<ImageType, ImageType>(m_working_trachea, m_working_trachea, leftLabel, 0);
249 writeImage<ImageType>(temp, "temp2.mhd");
251 sliceRelPosFilter = SliceRelPosFilterType::New();
252 sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
253 sliceRelPosFilter->VerboseStepOn();
254 sliceRelPosFilter->WriteStepOn();
255 sliceRelPosFilter->SetInput(m_working_image);
256 sliceRelPosFilter->SetInputObject(temp);
257 sliceRelPosFilter->SetDirection(2);
258 sliceRelPosFilter->SetFuzzyThreshold(0.6);
259 sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::LeftTo);
260 sliceRelPosFilter->Update();
261 m_working_image = sliceRelPosFilter->GetOutput();
262 writeImage<ImageType>(m_working_image, "afterslicebyslice.mhd");
266 m_output = m_working_image;
267 StopCurrentStep<ImageType>(m_output);
269 // Set output image information (required)
270 ImagePointer outputImage = this->GetOutput(0);
271 outputImage->SetRegions(m_working_image->GetLargestPossibleRegion());
272 outputImage->SetOrigin(m_working_image->GetOrigin());
273 outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion());
276 //--------------------------------------------------------------------
279 //--------------------------------------------------------------------
280 template <class TImageType>
282 clitk::ExtractLymphStationsFilter<TImageType>::
283 GenerateInputRequestedRegion() {
284 DD("GenerateInputRequestedRegion");
286 Superclass::GenerateInputRequestedRegion();
287 // Following needed because output region can be greater than input (trachea)
288 ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
289 ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
290 mediastinum->SetRequestedRegion(mediastinum->GetLargestPossibleRegion());
291 trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
293 //--------------------------------------------------------------------
296 //--------------------------------------------------------------------
297 template <class TImageType>
299 clitk::ExtractLymphStationsFilter<TImageType>::
302 // Final Step -> graft output (if SetNthOutput => redo)
303 this->GraftOutput(m_output);
305 //--------------------------------------------------------------------
308 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX