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://www.centreleonberard.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 CLITKEXTRACTMEDIASTINUMFILTER_TXX
20 #define CLITKEXTRACTMEDIASTINUMFILTER_TXX
23 #include "clitkCommon.h"
24 #include "clitkExtractMediastinumFilter.h"
25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
26 #include "clitkSliceBySliceRelativePositionFilter.h"
27 #include "clitkSegmentationUtils.h"
28 #include "clitkExtractAirwaysTreeInfoFilter.h"
29 #include "clitkCropLikeImageFilter.h"
35 #include "itkStatisticsLabelMapFilter.h"
36 #include "itkLabelImageToStatisticsLabelMapFilter.h"
37 #include "itkRegionOfInterestImageFilter.h"
38 #include "itkBinaryThresholdImageFilter.h"
39 #include "itkScalarImageKmeansImageFilter.h"
42 #include "RelativePositionPropImageFilter.h"
44 //--------------------------------------------------------------------
45 template <class ImageType>
46 clitk::ExtractMediastinumFilter<ImageType>::
47 ExtractMediastinumFilter():
48 clitk::StructuresExtractionFilter<ImageType>()
50 this->SetNumberOfRequiredInputs(1);
51 SetBackgroundValuePatient(0);
52 SetBackgroundValueLung(0);
53 SetBackgroundValueBones(0);
54 SetForegroundValueLeftLung(1);
55 SetForegroundValueRightLung(2);
56 SetDistanceMaxToAnteriorPartOfTheVertebralBody(10);
57 SetOutputMediastinumFilename("mediastinum.mhd");
60 //--------------------------------------------------------------------
63 //--------------------------------------------------------------------
64 template <class ImageType>
66 clitk::ExtractMediastinumFilter<ImageType>::
67 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
69 this->SetNthInput(0, const_cast<MaskImageType *>(image));
70 m_BackgroundValuePatient = bg;
72 //--------------------------------------------------------------------
75 //--------------------------------------------------------------------
76 template <class ImageType>
78 clitk::ExtractMediastinumFilter<ImageType>::
79 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
80 MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
82 this->SetNthInput(1, const_cast<MaskImageType *>(image));
83 m_BackgroundValueLung = bg;
84 m_ForegroundValueLeftLung = fgLeft;
85 m_ForegroundValueRightLung = fgRight;
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template <class ImageType>
93 clitk::ExtractMediastinumFilter<ImageType>::
94 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
96 this->SetNthInput(2, const_cast<MaskImageType *>(image));
97 m_BackgroundValueBones = bg;
99 //--------------------------------------------------------------------
102 //--------------------------------------------------------------------
103 template <class ImageType>
105 clitk::ExtractMediastinumFilter<ImageType>::
106 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
108 this->SetNthInput(3, const_cast<MaskImageType *>(image));
109 m_BackgroundValueTrachea = bg;
111 //--------------------------------------------------------------------
113 //--------------------------------------------------------------------
114 template <class ImageType>
116 clitk::ExtractMediastinumFilter<ImageType>::
117 SetFuzzyThreshold(std::string tag, double value)
119 m_FuzzyThreshold[tag] = value;
121 //--------------------------------------------------------------------
124 //--------------------------------------------------------------------
125 template <class ImageType>
127 clitk::ExtractMediastinumFilter<ImageType>::
128 GetFuzzyThreshold(std::string tag)
130 if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) {
131 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds.");
135 return m_FuzzyThreshold[tag];
137 //--------------------------------------------------------------------
140 //--------------------------------------------------------------------
141 template <class ImageType>
143 clitk::ExtractMediastinumFilter<ImageType>::
144 GenerateInputRequestedRegion()
146 // DD("GenerateInputRequestedRegion");
147 // Do not call default
148 // Superclass::GenerateInputRequestedRegion();
149 // DD("End GenerateInputRequestedRegion");
151 //--------------------------------------------------------------------
154 //--------------------------------------------------------------------
155 template <class ImageType>
157 clitk::ExtractMediastinumFilter<ImageType>::
158 SetInput(const ImageType * image)
160 this->SetNthInput(0, const_cast<ImageType *>(image));
162 //--------------------------------------------------------------------
166 //--------------------------------------------------------------------
167 template <class ImageType>
169 clitk::ExtractMediastinumFilter<ImageType>::
170 GenerateOutputInformation() {
171 // Do not call default
172 // Superclass::GenerateOutputInformation();
174 //--------------------------------------------------------------------
175 // Get input pointers
177 ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
178 MaskImagePointer patient, lung, bones, trachea;
179 patient = this->GetAFDB()->template GetImage <MaskImageType>("Patient");
180 lung = this->GetAFDB()->template GetImage <MaskImageType>("Lungs");
182 bones = this->GetAFDB()->template GetImage <MaskImageType>("Bones");
184 trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");
186 //this->VerboseImageSizeFlagOn();
188 //--------------------------------------------------------------------
189 // Step : Crop support (patient) to lung extend in RL
190 this->StartNewStep("[Mediastinum] Crop support like lungs along LR");
191 typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
192 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
193 cropFilter->SetInput(patient);
194 cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
195 cropFilter->Update();
196 output = cropFilter->GetOutput();
197 this->template StopCurrentStep<MaskImageType>(output);
199 //--------------------------------------------------------------------
200 // Step : Crop support (previous) to bones extend in AP
202 this->StartNewStep("[Mediastinum] Crop support like bones along AP");
203 cropFilter = CropFilterType::New();
204 cropFilter->SetInput(output);
205 cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
206 cropFilter->Update();
207 output = cropFilter->GetOutput();
208 this->template StopCurrentStep<MaskImageType>(output);
211 //--------------------------------------------------------------------
212 // Step : patient minus lungs, minus bones, minus trachea
213 this->StartNewStep("[Mediastinum] Patient contours minus lungs, trachea [and bones]");
214 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
215 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
216 boolFilter->InPlaceOn();
217 boolFilter->SetInput1(output);
218 boolFilter->SetInput2(lung);
219 boolFilter->SetOperationType(BoolFilterType::AndNot);
220 boolFilter->Update();
222 boolFilter->SetInput1(boolFilter->GetOutput());
223 boolFilter->SetInput2(bones);
224 boolFilter->SetOperationType(BoolFilterType::AndNot);
225 boolFilter->Update();
227 boolFilter->SetInput1(boolFilter->GetOutput());
228 boolFilter->SetInput2(trachea);
229 boolFilter->SetOperationType(BoolFilterType::AndNot);
230 boolFilter->Update();
231 output = boolFilter->GetOutput();
233 // Auto crop to gain support area
234 output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
235 this->template StopCurrentStep<MaskImageType>(output);
237 //--------------------------------------------------------------------
238 // Step : LR limits from lung (need separate lung ?)
239 // Get separate lung images to get only the right and left lung
240 // (because RelativePositionPropImageFilter only consider fg=1);
241 // (label must be '1' because right is greater than left). (WE DO
242 // NOT NEED TO SEPARATE ? )
243 this->StartNewStep("[Mediastinum] Left/Right limits with lungs");
245 // The following cannot be "inplace" because mask is the same than input ...
246 MaskImagePointer right_lung =
247 clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0, false);
248 MaskImagePointer left_lung =
249 clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0, false);
250 left_lung = clitk::SetBackground<MaskImageType, MaskImageType>(left_lung, left_lung, 2, 1, false);
251 right_lung = clitk::ResizeImageLike<MaskImageType>(right_lung, output, this->GetBackgroundValue());
252 left_lung = clitk::ResizeImageLike<MaskImageType>(left_lung, output, this->GetBackgroundValue());
253 this->GetAFDB()->template SetImage<MaskImageType>("RightLung", "seg/RightLung.mhd",
255 this->GetAFDB()->template SetImage<MaskImageType>("LeftLung", "seg/LeftLung.mhd",
257 this->GetAFDB()->Write();
258 this->template StopCurrentStep<MaskImageType>(output);
260 //--------------------------------------------------------------------
261 // Step : AP limits from bones
262 // Separate the bones in the ant-post middle
263 MaskImagePointer bones_ant;
264 MaskImagePointer bones_post;
265 MaskImagePointType middle_AntPost_position;
267 this->StartNewStep("[Mediastinum] Ant/Post limits with bones");
269 // To define ant and post part of the bones with a single horizontal line
270 MaskImageIndexType index_bones_middle;
272 // Method1: cut in the middle (not optimal)
274 middle_AntPost_position[0] = middle_AntPost_position[2] = 0;
275 middle_AntPost_position[1] = bones->GetOrigin()[1]+
276 (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
277 DD(middle_AntPost_position);
278 bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
281 // Method2: Use VertebralBody, take most ant point
282 MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
283 FindExtremaPointInAGivenDirection<MaskImageType>(VertebralBody, this->GetBackgroundValue(),
284 1, true, middle_AntPost_position);
285 bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle);
287 // Split bone image first into two parts (ant and post), and crop
288 // lateraly to get vertebral
289 typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
290 // typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
291 typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
292 MaskImageRegionType region = bones->GetLargestPossibleRegion();
293 MaskImageSizeType size = region.GetSize();
294 MaskImageIndexType index = region.GetIndex();
296 // crop LR to keep 1/4 center part
297 // index[0] = size[0]/4+size[0]/8;
298 //size[0] = size[0]/4;
299 // crop AP to keep first (ant) part
300 size[1] = index_bones_middle[1]; //size[1]/2.0;
301 region.SetSize(size);
302 region.SetIndex(index);
303 roiFilter->SetInput(bones);
304 roiFilter->SetRegionOfInterest(region);
305 roiFilter->ReleaseDataFlagOff();
307 bones_ant = roiFilter->GetOutput();
308 // writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
310 roiFilter = ROIFilterType::New();
311 index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
312 size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
313 region.SetIndex(index);
314 region.SetSize(size);
315 roiFilter->SetInput(bones);
316 roiFilter->SetRegionOfInterest(region);
317 roiFilter->ReleaseDataFlagOff();
319 bones_post = roiFilter->GetOutput();
320 // writeImage<MaskImageType>(bones_post, "b_post.mhd");
322 // Now, insert this image in the AFDB
323 this->GetAFDB()->template SetImage<MaskImageType>("Bones_Post", "seg/Bones_Post.mhd",
325 this->GetAFDB()->template SetImage<MaskImageType>("Bones_Ant", "seg/Bones_Ant.mhd",
327 this->GetAFDB()->Write();
329 this->template StopCurrentStep<MaskImageType>(output);
332 //--------------------------------------------------------------------
333 // Generic RelativePosition processes
334 output = this->ApplyRelativePositionList("Mediastinum", output);
336 //--------------------------------------------------------------------
337 // Step : SI limits It is better to do this limit *AFTER* the
338 // RelativePosition to avoid some issue due to superior boundaries.
339 this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag");
340 // load Cricoid, get centroid, cut above (or below), lower bound
341 MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
342 MaskImagePointType p;
343 p[0] = p[1] = p[2] = 0.0; // to avoid warning
344 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag,
345 this->GetBackgroundValue(), 2, true, p);
346 output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, this->GetBackgroundValue());
347 this->template StopCurrentStep<MaskImageType>(output);
349 //--------------------------------------------------------------------
351 this->StartNewStep("[Mediastinum] Keep main connected component");
352 output = clitk::Labelize<MaskImageType>(output, this->GetBackgroundValue(), false, 500);
353 // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
354 output = clitk::KeepLabels<MaskImageType>(output, this->GetBackgroundValue(),
355 this->GetForegroundValue(), 1, 1, 0);
356 this->template StopCurrentStep<MaskImageType>(output);
358 //--------------------------------------------------------------------
359 // Step: Remove post part from VertebralBody
360 this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody");
361 RemovePostPartOfVertebralBody();
362 this->template StopCurrentStep<MaskImageType>(output);
364 //--------------------------------------------------------------------
365 // Step: Slice by Slice CCL
366 this->StartNewStep("[Mediastinum] Slice by Slice keep only one component");
367 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
368 // typename ExtractSliceFilterType::Pointer
369 ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
370 extractSliceFilter->SetInput(output);
371 extractSliceFilter->SetDirection(2);
372 extractSliceFilter->Update();
373 typedef typename ExtractSliceFilterType::SliceType SliceType;
374 std::vector<typename SliceType::Pointer> mSlices;
375 extractSliceFilter->GetOutputSlices(mSlices);
376 for(unsigned int i=0; i<mSlices.size(); i++) {
377 mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
378 mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
380 typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
381 typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
382 joinFilter->SetOrigin(output->GetOrigin()[2]);
383 joinFilter->SetSpacing(output->GetSpacing()[2]);
384 for(unsigned int i=0; i<mSlices.size(); i++) {
385 joinFilter->PushBackInput(mSlices[i]);
387 joinFilter->Update();
388 output = joinFilter->GetOutput();
389 this->template StopCurrentStep<MaskImageType>(output);
391 //--------------------------------------------------------------------
392 // Step 10 : AutoCrop
393 this->StartNewStep("[Mediastinum] AutoCrop");
394 output = clitk::AutoCrop<MaskImageType>(output, this->GetBackgroundValue());
395 this->template StopCurrentStep<MaskImageType>(output);
397 // End, set the real size
398 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
399 this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
400 this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
401 this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
403 //--------------------------------------------------------------------
406 //--------------------------------------------------------------------
407 template <class ImageType>
409 clitk::ExtractMediastinumFilter<ImageType>::
412 this->GraftOutput(output);
413 // Store image filenames into AFDB
414 this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
417 //--------------------------------------------------------------------
420 //--------------------------------------------------------------------
421 template <class ImageType>
423 clitk::ExtractMediastinumFilter<ImageType>::
424 RemovePostPartOfVertebralBody()
428 Posteriorly, Station 8 abuts the descending aorta and the anterior
429 aspect of the vertebral body until an imaginary horizontal line
430 running 1 cm posterior to the anterior border of the vertebral
433 => We use this definition for all the mediastinum
435 Find most Ant point in VertebralBody. Consider the horizontal line
436 which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
440 // Get VertebralBody mask image
441 MaskImagePointer VertebralBody =
442 this->GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
444 // Consider vertebral body slice by slice
445 std::vector<MaskSlicePointer> vertebralSlices;
446 clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
448 // For each slice, compute the most anterior point
449 std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
450 for(uint i=0; i<vertebralSlices.size(); i++) {
451 MaskSlicePointType p;
452 bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i],
453 this->GetBackgroundValue(),
456 vertebralAntPositionBySlice[i] = p;
459 // It should not happen ! But sometimes, a contour is missing or
460 // the VertebralBody is not delineated enough inferiorly ... in
461 // those cases, we consider the first found slice.
462 // std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
463 // [ Possible alternative -> consider previous limit ]
467 // Convert 2D points in slice into 3D points
468 std::vector<MaskImagePointType> vertebralAntPositions;
469 clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice,
471 vertebralAntPositions);
473 // DEBUG : write list of points
474 clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
475 "VertebralBodyMostAnteriorPoints.txt");
477 // Cut support posteriorly 1cm the most anterior point of the
478 // VertebralBody. Do nothing below and above the VertebralBody. To
479 // do that compute several region, slice by slice and fill.
480 MaskImageRegionType region;
481 MaskImageSizeType size;
482 MaskImageIndexType start;
483 size[2] = 1; // a single slice
484 start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
485 size[0] = output->GetLargestPossibleRegion().GetSize()[0];
486 for(uint i=0; i<vertebralAntPositions.size(); i++) {
487 typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
489 IteratorType(output, output->GetLargestPossibleRegion());
490 MaskImageIndexType index;
491 // Consider some cm posterior to most anterior positions (usually
493 vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
494 // Get index of this point
495 output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
496 // Compute region (a single slice)
498 start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
499 size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
500 region.SetSize(size);
501 region.SetIndex(start);
503 if (output->GetLargestPossibleRegion().IsInside(start)) {
504 itk::ImageRegionIterator<MaskImageType> it(output, region);
506 while (!it.IsAtEnd()) {
507 it.Set(this->GetBackgroundValue());
513 //--------------------------------------------------------------------
516 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX