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():
49 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
50 itk::ImageToImageFilter<ImageType, MaskImageType>()
52 this->SetNumberOfRequiredInputs(1);
53 SetBackgroundValuePatient(0);
54 SetBackgroundValueLung(0);
55 SetBackgroundValueBones(0);
56 SetForegroundValueLeftLung(1);
57 SetForegroundValueRightLung(2);
58 SetBackgroundValue(0);
59 SetForegroundValue(1);
60 SetIntermediateSpacing(6);
61 SetFuzzyThreshold("LR_lungs", 0.3);
62 SetFuzzyThreshold("bones", 0.6);
63 SetFuzzyThreshold("inf_lungs", 0.05);
64 SetDistanceMaxToAnteriorPartOfTheVertebralBody(10);
65 SetOutputMediastinumFilename("mediastinum.mhd");
68 //--------------------------------------------------------------------
71 //--------------------------------------------------------------------
72 template <class ImageType>
74 clitk::ExtractMediastinumFilter<ImageType>::
75 SetInputPatientLabelImage(const MaskImageType * image, MaskImagePixelType bg)
77 this->SetNthInput(0, const_cast<MaskImageType *>(image));
78 m_BackgroundValuePatient = bg;
80 //--------------------------------------------------------------------
83 //--------------------------------------------------------------------
84 template <class ImageType>
86 clitk::ExtractMediastinumFilter<ImageType>::
87 SetInputLungLabelImage(const MaskImageType * image, MaskImagePixelType bg,
88 MaskImagePixelType fgLeft, MaskImagePixelType fgRight)
90 this->SetNthInput(1, const_cast<MaskImageType *>(image));
91 m_BackgroundValueLung = bg;
92 m_ForegroundValueLeftLung = fgLeft;
93 m_ForegroundValueRightLung = fgRight;
95 //--------------------------------------------------------------------
98 //--------------------------------------------------------------------
99 template <class ImageType>
101 clitk::ExtractMediastinumFilter<ImageType>::
102 SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg)
104 this->SetNthInput(2, const_cast<MaskImageType *>(image));
105 m_BackgroundValueBones = bg;
107 //--------------------------------------------------------------------
110 //--------------------------------------------------------------------
111 template <class ImageType>
113 clitk::ExtractMediastinumFilter<ImageType>::
114 SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg)
116 this->SetNthInput(3, const_cast<MaskImageType *>(image));
117 m_BackgroundValueTrachea = bg;
119 //--------------------------------------------------------------------
121 //--------------------------------------------------------------------
122 template <class ImageType>
124 clitk::ExtractMediastinumFilter<ImageType>::
125 SetFuzzyThreshold(std::string tag, double value)
127 m_FuzzyThreshold[tag] = value;
129 //--------------------------------------------------------------------
132 //--------------------------------------------------------------------
133 template <class ImageType>
135 clitk::ExtractMediastinumFilter<ImageType>::
136 GetFuzzyThreshold(std::string tag)
138 if (m_FuzzyThreshold.find(tag) == m_FuzzyThreshold.end()) {
139 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThresholds.");
143 return m_FuzzyThreshold[tag];
145 //--------------------------------------------------------------------
148 //--------------------------------------------------------------------
149 template <class ImageType>
151 clitk::ExtractMediastinumFilter<ImageType>::
152 GenerateInputRequestedRegion()
154 // DD("GenerateInputRequestedRegion");
155 // Do not call default
156 // Superclass::GenerateInputRequestedRegion();
157 // DD("End GenerateInputRequestedRegion");
160 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template <class ImageType>
166 clitk::ExtractMediastinumFilter<ImageType>::
167 SetInput(const ImageType * image)
169 this->SetNthInput(0, const_cast<ImageType *>(image));
171 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 template <class ImageType>
177 clitk::ExtractMediastinumFilter<ImageType>::
178 GenerateOutputInformation() {
179 // Do not call default
180 // Superclass::GenerateOutputInformation();
182 //--------------------------------------------------------------------
183 // Get input pointers
185 ImageConstPointer input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
186 MaskImagePointer patient, lung, bones, trachea;
187 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
188 lung = GetAFDB()->template GetImage <MaskImageType>("Lungs");
190 bones = GetAFDB()->template GetImage <MaskImageType>("Bones");
192 trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
194 //--------------------------------------------------------------------
195 // Step : Crop support (patient) to lung extend in RL
196 StartNewStep("Crop support like lungs along LR");
197 typedef clitk::CropLikeImageFilter<MaskImageType> CropFilterType;
198 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
199 cropFilter->SetInput(patient);
200 cropFilter->SetCropLikeImage(lung, 0);// Indicate that we only crop in X (Left-Right) axe
201 cropFilter->Update();
202 output = cropFilter->GetOutput();
203 this->template StopCurrentStep<MaskImageType>(output);
205 //--------------------------------------------------------------------
206 // Step : Crop support (previous) to bones extend in AP
208 StartNewStep("Crop support like bones along AP");
209 cropFilter = CropFilterType::New();
210 cropFilter->SetInput(output);
211 cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe
212 cropFilter->Update();
213 output = cropFilter->GetOutput();
214 this->template StopCurrentStep<MaskImageType>(output);
217 //--------------------------------------------------------------------
218 // Step : patient minus lungs, minus bones, minus trachea
219 StartNewStep("Patient contours minus lungs, trachea [and bones]");
220 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
221 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
222 boolFilter->InPlaceOn();
223 boolFilter->SetInput1(output);
224 boolFilter->SetInput2(lung);
225 boolFilter->SetOperationType(BoolFilterType::AndNot);
226 boolFilter->Update();
228 boolFilter->SetInput1(boolFilter->GetOutput());
229 boolFilter->SetInput2(bones);
230 boolFilter->SetOperationType(BoolFilterType::AndNot);
231 boolFilter->Update();
233 boolFilter->SetInput1(boolFilter->GetOutput());
234 boolFilter->SetInput2(trachea);
235 boolFilter->SetOperationType(BoolFilterType::AndNot);
236 boolFilter->Update();
237 output = boolFilter->GetOutput();
239 // Auto crop to gain support area
240 output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
241 this->template StopCurrentStep<MaskImageType>(output);
243 //--------------------------------------------------------------------
244 // Step : LR limits from lung (need separate lung ?)
245 // Get separate lung images to get only the right and left lung
246 // (because RelativePositionPropImageFilter only consider fg=1);
247 // (label must be '1' because right is greater than left). (WE DO
248 // NOT NEED TO SEPARATE ? )
249 StartNewStep("Left/Right limits with lungs");
251 // The following cannot be "inplace" because mask is the same than input ...
252 MaskImagePointer right_lung =
253 clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 2, 0, false);
254 MaskImagePointer left_lung =
255 clitk::SetBackground<MaskImageType, MaskImageType>(lung, lung, 1, 0, false);
256 right_lung = clitk::ResizeImageLike<MaskImageType>(right_lung, output, GetBackgroundValue());
257 left_lung = clitk::ResizeImageLike<MaskImageType>(left_lung, output, GetBackgroundValue());
258 // writeImage<MaskImageType>(right_lung, "right.mhd");
259 // writeImage<MaskImageType>(left_lung, "left.mhd");
261 typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
262 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
263 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
264 relPosFilter->VerboseStepFlagOff();
265 relPosFilter->WriteStepFlagOff();
266 relPosFilter->SetInput(output);
267 relPosFilter->SetInputObject(left_lung);
268 relPosFilter->AddOrientationType(RelPosFilterType::AtRightTo); // warning left lung is at right ;)
269 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
270 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs"));
271 relPosFilter->Update();
272 output = relPosFilter->GetOutput();
274 relPosFilter = RelPosFilterType::New();
275 relPosFilter->SetInput(output);
276 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
277 relPosFilter->VerboseStepFlagOff();
278 relPosFilter->WriteStepFlagOff();
279 relPosFilter->SetInput(output);
280 relPosFilter->SetInputObject(right_lung);
281 relPosFilter->AddOrientationType(RelPosFilterType::AtLeftTo);
282 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
283 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs"));
284 relPosFilter->Update();
285 output = relPosFilter->GetOutput();
286 this->template StopCurrentStep<MaskImageType>(output);
288 //--------------------------------------------------------------------
289 // Step : superior limits
290 StartNewStep("Keep inferior to CricoidCartilag");
291 // load Cricoid, get centroid, cut above (or below), lower bound
292 MaskImagePointer CricoidCartilag = GetAFDB()->template GetImage <MaskImageType>("CricoidCartilag");
293 MaskImagePointType p;
294 p[0] = p[1] = p[2] = 0.0; // to avoid warning
295 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(CricoidCartilag, GetBackgroundValue(), 2, true, p);
296 output = clitk::CropImageRemoveGreaterThan<MaskImageType>(output, 2, p[2], true, GetBackgroundValue());
297 this->template StopCurrentStep<MaskImageType>(output);
299 //--------------------------------------------------------------------
300 // Step : AP limits from bones
301 // Separate the bones in the ant-post middle
302 MaskImageConstPointer bones_ant;
303 MaskImageConstPointer bones_post;
304 MaskImagePointType middle_AntPost__position;
306 StartNewStep("Ant/Post limits with bones");
307 middle_AntPost__position[0] = middle_AntPost__position[2] = 0;
308 middle_AntPost__position[1] = bones->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0;
309 MaskImageIndexType index_bones_middle;
310 bones->TransformPhysicalPointToIndex(middle_AntPost__position, index_bones_middle);
312 // Split bone image first into two parts (ant and post), and crop
313 // lateraly to get vertebral
314 typedef itk::RegionOfInterestImageFilter<MaskImageType, MaskImageType> ROIFilterType;
315 // typedef itk::ExtractImageFilter<MaskImageType, MaskImageType> ROIFilterType;
316 typename ROIFilterType::Pointer roiFilter = ROIFilterType::New();
317 MaskImageRegionType region = bones->GetLargestPossibleRegion();
318 MaskImageSizeType size = region.GetSize();
319 MaskImageIndexType index = region.GetIndex();
321 // crop LR to keep 1/4 center part
322 index[0] = size[0]/4+size[0]/8;
324 // crop AP to keep first (ant) part
325 size[1] = index_bones_middle[1]; //size[1]/2.0;
326 region.SetSize(size);
327 region.SetIndex(index);
328 roiFilter->SetInput(bones);
329 roiFilter->SetRegionOfInterest(region);
330 roiFilter->ReleaseDataFlagOff();
332 bones_ant = roiFilter->GetOutput();
333 // writeImage<MaskImageType>(bones_ant, "b_ant.mhd");
335 roiFilter = ROIFilterType::New();
336 index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1;
337 size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1];
338 region.SetIndex(index);
339 region.SetSize(size);
340 roiFilter->SetInput(bones);
341 roiFilter->SetRegionOfInterest(region);
342 roiFilter->ReleaseDataFlagOff();
344 bones_post = roiFilter->GetOutput();
345 // writeImage<MaskImageType>(bones_post, "b_post.mhd");
348 relPosFilter->SetCurrentStepNumber(0);
349 relPosFilter->ResetPipeline();// = RelPosFilterType::New();
350 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
351 relPosFilter->VerboseStepFlagOff();
352 relPosFilter->WriteStepFlagOff();
353 relPosFilter->SetInput(output);
354 relPosFilter->SetInputObject(bones_post);
355 relPosFilter->AddOrientationType(RelPosFilterType::AntTo);
356 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
357 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones"));
358 relPosFilter->Update();
359 output = relPosFilter->GetOutput();
360 // writeImage<MaskImageType>(output, "post.mhd");
362 relPosFilter->SetInput(relPosFilter->GetOutput());
363 relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
364 relPosFilter->VerboseStepFlagOff();
365 relPosFilter->WriteStepFlagOff();
366 relPosFilter->SetInput(output);
367 relPosFilter->SetInputObject(bones_ant);
368 relPosFilter->AddOrientationType(RelPosFilterType::PostTo);
369 relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
370 relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones"));
371 relPosFilter->Update();
372 output = relPosFilter->GetOutput();
373 this->template StopCurrentStep<MaskImageType>(output);
376 //--------------------------------------------------------------------
378 StartNewStep("Keep main connected component");
379 output = clitk::Labelize<MaskImageType>(output, GetBackgroundValue(), false, 500);
380 // output = RemoveLabels<MaskImageType>(output, BG, param->GetLabelsToRemove());
381 output = clitk::KeepLabels<MaskImageType>(output, GetBackgroundValue(),
382 GetForegroundValue(), 1, 1, 0);
383 this->template StopCurrentStep<MaskImageType>(output);
385 //--------------------------------------------------------------------
386 // Step: Remove post part from VertebralBody
387 StartNewStep("Remove post part according to VertebralBody");
388 RemovePostPartOfVertebralBody();
389 this->template StopCurrentStep<MaskImageType>(output);
391 //--------------------------------------------------------------------
392 // Step: Remove ant part according to Sternum
393 StartNewStep("Remove ant part according to Sternum");
394 MaskImagePointer Sternum = GetAFDB()->template GetImage <MaskImageType>("Sternum");
395 output = clitk::SliceBySliceRelativePosition<MaskImageType>(output, Sternum, 2,
396 GetFuzzyThreshold("ant_sternum"),
397 "PostTo", false, 3, true, false);
398 this->template StopCurrentStep<MaskImageType>(output);
400 //--------------------------------------------------------------------
401 // Step: Slice by Slice CCL
402 StartNewStep("Slice by Slice keep only one component");
403 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
404 // typename ExtractSliceFilterType::Pointer
405 ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
406 extractSliceFilter->SetInput(output);
407 extractSliceFilter->SetDirection(2);
408 extractSliceFilter->Update();
409 typedef typename ExtractSliceFilterType::SliceType SliceType;
410 std::vector<typename SliceType::Pointer> mSlices;
411 extractSliceFilter->GetOutputSlices(mSlices);
412 for(unsigned int i=0; i<mSlices.size(); i++) {
413 mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 100);
414 mSlices[i] = KeepLabels<SliceType>(mSlices[i], 0, 1, 1, 1, true);
416 typedef itk::JoinSeriesImageFilter<SliceType, MaskImageType> JoinSeriesFilterType;
417 typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
418 joinFilter->SetOrigin(output->GetOrigin()[2]);
419 joinFilter->SetSpacing(output->GetSpacing()[2]);
420 for(unsigned int i=0; i<mSlices.size(); i++) {
421 joinFilter->PushBackInput(mSlices[i]);
423 joinFilter->Update();
424 output = joinFilter->GetOutput();
425 this->template StopCurrentStep<MaskImageType>(output);
427 //--------------------------------------------------------------------
428 // Step 10 : AutoCrop
429 StartNewStep("AutoCrop");
430 output = clitk::AutoCrop<MaskImageType>(output, GetBackgroundValue());
431 this->template StopCurrentStep<MaskImageType>(output);
433 // End, set the real size
434 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
435 this->GetOutput(0)->SetLargestPossibleRegion(output->GetLargestPossibleRegion());
436 this->GetOutput(0)->SetRequestedRegion(output->GetLargestPossibleRegion());
437 this->GetOutput(0)->SetBufferedRegion(output->GetLargestPossibleRegion());
439 //--------------------------------------------------------------------
442 //--------------------------------------------------------------------
443 template <class ImageType>
445 clitk::ExtractMediastinumFilter<ImageType>::
448 this->GraftOutput(output);
449 // Store image filenames into AFDB
450 GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename());
453 //--------------------------------------------------------------------
456 //--------------------------------------------------------------------
457 template <class ImageType>
459 clitk::ExtractMediastinumFilter<ImageType>::
460 RemovePostPartOfVertebralBody()
464 Posteriorly, Station 8 abuts the descending aorta and the anterior
465 aspect of the vertebral body until an imaginary horizontal line
466 running 1 cm posterior to the anterior border of the vertebral
469 => We use this definition for all the mediastinum
471 Find most Ant point in VertebralBody. Consider the horizontal line
472 which is 'DistanceMaxToAnteriorPartOfTheVertebralBody' away from
476 // Get VertebralBody mask image
477 MaskImagePointer VertebralBody =
478 GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
480 // Consider vertebral body slice by slice
481 std::vector<MaskSlicePointer> vertebralSlices;
482 clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vertebralSlices);
484 // For each slice, compute the most anterior point
485 std::map<int, MaskSlicePointType> vertebralAntPositionBySlice;
486 for(uint i=0; i<vertebralSlices.size(); i++) {
487 MaskSlicePointType p;
488 bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vertebralSlices[i],
489 GetBackgroundValue(),
492 vertebralAntPositionBySlice[i] = p;
495 // It should not happen ! But sometimes, a contour is missing or
496 // the VertebralBody is not delineated enough inferiorly ... in
497 // those cases, we consider the first found slice.
498 // std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
499 // [ Possible alternative -> consider previous limit ]
503 // Convert 2D points in slice into 3D points
504 std::vector<MaskImagePointType> vertebralAntPositions;
505 clitk::PointsUtils<MaskImageType>::Convert2DMapTo3DList(vertebralAntPositionBySlice,
507 vertebralAntPositions);
509 // DEBUG : write list of points
510 clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
511 "VertebralBodyMostAnteriorPoints.txt");
513 // Cut support posteriorly 1cm the most anterior point of the
514 // VertebralBody. Do nothing below and above the VertebralBody. To
515 // do that compute several region, slice by slice and fill.
516 MaskImageRegionType region;
517 MaskImageSizeType size;
518 MaskImageIndexType start;
519 size[2] = 1; // a single slice
520 start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
521 size[0] = output->GetLargestPossibleRegion().GetSize()[0];
522 for(uint i=0; i<vertebralAntPositions.size(); i++) {
523 typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
525 IteratorType(output, output->GetLargestPossibleRegion());
526 MaskImageIndexType index;
527 // Consider some cm posterior to most anterior positions (usually
529 vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheVertebralBody();
530 // Get index of this point
531 output->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
532 // Compute region (a single slice)
534 start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
535 size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
536 region.SetSize(size);
537 region.SetIndex(start);
539 if (output->GetLargestPossibleRegion().IsInside(start)) {
540 itk::ImageRegionIterator<MaskImageType> it(output, region);
542 while (!it.IsAtEnd()) {
543 it.Set(GetBackgroundValue());
549 //--------------------------------------------------------------------
552 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX