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"
30 #include "clitkMeshToBinaryImageFilter.h"
33 #include <itkStatisticsLabelMapFilter.h>
34 #include <itkLabelImageToStatisticsLabelMapFilter.h>
35 #include <itkRegionOfInterestImageFilter.h>
36 #include <itkBinaryThresholdImageFilter.h>
37 #include <itkImageSliceConstIteratorWithIndex.h>
38 #include <itkImageSliceIteratorWithIndex.h>
39 #include <itkBinaryThinningImageFilter.h>
42 #include "RelativePositionPropImageFilter.h"
45 #include <vtkAppendPolyData.h>
46 #include <vtkPolyDataWriter.h>
47 #include <vtkCellArray.h>
49 //--------------------------------------------------------------------
50 template <class TImageType>
51 clitk::ExtractLymphStationsFilter<TImageType>::
52 ExtractLymphStationsFilter():
53 clitk::StructuresExtractionFilter<ImageType>()
55 this->SetNumberOfRequiredInputs(1);
56 SetBackgroundValue(0);
57 SetForegroundValue(1);
58 ForceSupportsFlagOn();
59 SetSupportLimitsFilename("none");
60 CheckSupportFlagOff();
63 ExtractStation_3P_SetDefaultValues();
64 ExtractStation_2RL_SetDefaultValues();
65 ExtractStation_3A_SetDefaultValues();
66 ExtractStation_1RL_SetDefaultValues();
67 ExtractStation_4RL_SetDefaultValues();
68 ExtractStation_5_SetDefaultValues();
69 ExtractStation_6_SetDefaultValues();
72 ExtractStation_7_SetDefaultValues();
73 ExtractStation_8_SetDefaultValues();
75 //--------------------------------------------------------------------
78 //--------------------------------------------------------------------
79 template <class TImageType>
81 clitk::ExtractLymphStationsFilter<TImageType>::
82 GenerateOutputInformation() {
85 m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
86 m_Mediastinum = this->GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
88 // DD(this->GetVerboseMemoryFlag());
89 // clitk::PrintMemory(this->GetVerboseMemoryFlag(), "Start");
91 // Clean some computer landmarks to force the recomputation
92 // FIXME -> to put elsewhere ?
93 this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
95 // Must I compute the supports ?
96 bool supportsExist = true;
97 if (!GetForceSupportsFlag()) {
99 m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
100 m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
101 m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
102 m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
103 m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
104 m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
106 m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
107 m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
108 m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
109 m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
110 m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
111 m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
112 m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
113 m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
114 m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
115 } catch(clitk::ExceptionObject o) {
116 supportsExist = false;
120 if (!supportsExist || GetForceSupportsFlag()) {
121 this->StartNewStep("Supports for stations");
122 this->StartSubStep();
124 // FIXME : why should I remove theses tags ???
125 this->GetAFDB()->RemoveTag("CarinaZ");
126 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
127 this->GetAFDB()->RemoveTag("ApexOfTheChest");
128 this->GetAFDB()->RemoveTag("RightBronchus");
129 this->GetAFDB()->RemoveTag("LeftBronchus");
130 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
131 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
132 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
133 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
134 ExtractStationSupports();
139 ExtractStation_1RL();
140 ExtractStation_2RL();
148 // ---------- todo -----------------------
151 // ExtractStation_8();
154 //this->StartNewStep("Station 7");
155 //this->StartSubStep();
156 //ExtractStation_7();
157 //this->StopSubStep();
160 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 template <class TImageType>
166 clitk::ExtractLymphStationsFilter<TImageType>::
167 GenerateInputRequestedRegion() {
168 //DD("GenerateInputRequestedRegion (nothing?)");
170 //--------------------------------------------------------------------
173 //--------------------------------------------------------------------
174 template <class TImageType>
176 clitk::ExtractLymphStationsFilter<TImageType>::
178 // Final Step -> graft output (if SetNthOutput => redo)
179 // this->GraftOutput(m_ListOfStations["8"]);
181 //--------------------------------------------------------------------
184 //--------------------------------------------------------------------
185 template <class TImageType>
187 clitk::ExtractLymphStationsFilter<TImageType>::
188 CheckForStation(std::string station)
190 // Compute Station name
191 std::string s = "Station"+station;
194 // Define the starting support
195 // if (GetComputeStation(station)) {
196 // std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
198 if (GetComputeStation(station)) {
199 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
204 // std::cout << "Station " << station << " found. I used it" << std::endl;
208 // Check if station already exist in DB
210 // FIXME -> do nothing if not on the command line. Is it what I want ?
211 //bool found = false;
212 if (this->GetAFDB()->TagExist(s)) {
213 m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
218 //--------------------------------------------------------------------
221 //--------------------------------------------------------------------
222 template <class TImageType>
224 clitk::ExtractLymphStationsFilter<TImageType>::
225 GetComputeStation(std::string station)
227 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
229 //--------------------------------------------------------------------
232 //--------------------------------------------------------------------
233 template <class TImageType>
235 clitk::ExtractLymphStationsFilter<TImageType>::
236 AddComputeStation(std::string station)
238 m_ComputeStationMap[station] = true;
240 //--------------------------------------------------------------------
244 //--------------------------------------------------------------------
245 template<class PointType>
246 class comparePointsX {
248 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
250 //--------------------------------------------------------------------
253 //--------------------------------------------------------------------
254 template<class PairType>
255 class comparePointsWithAngle {
257 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
259 //--------------------------------------------------------------------
262 //--------------------------------------------------------------------
264 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
265 std::vector<itk::Point<double, Dim-1> > previous;
266 HypercubeCorners<Dim-1>(previous);
267 out.resize(previous.size()*2);
268 for(unsigned int i=0; i<out.size(); i++) {
269 itk::Point<double, Dim> p;
270 if (i<previous.size()) p[0] = 0;
272 for(int j=0; j<Dim-1; j++)
274 p[j+1] = previous[i%previous.size()][j];
279 //--------------------------------------------------------------------
282 //--------------------------------------------------------------------
284 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
289 //--------------------------------------------------------------------
292 //--------------------------------------------------------------------
293 template<class ImageType>
294 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
295 std::vector<typename ImageType::PointType> & bounds)
297 // Get image max/min coordinates
298 const unsigned int dim=ImageType::ImageDimension;
299 typedef typename ImageType::PointType PointType;
300 typedef typename ImageType::IndexType IndexType;
301 PointType min_c, max_c;
302 IndexType min_i, max_i;
303 min_i = image->GetLargestPossibleRegion().GetIndex();
304 for(unsigned int i=0; i<dim; i++)
305 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
306 image->TransformIndexToPhysicalPoint(min_i, min_c);
307 image->TransformIndexToPhysicalPoint(max_i, max_c);
309 // Get corners coordinates
310 HypercubeCorners<ImageType::ImageDimension>(bounds);
311 for(unsigned int i=0; i<bounds.size(); i++) {
312 for(unsigned int j=0; j<dim; j++) {
313 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
314 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
318 //--------------------------------------------------------------------
321 //--------------------------------------------------------------------
322 template <class ImageType>
324 clitk::ExtractLymphStationsFilter<ImageType>::
325 Remove_Structures(std::string station, std::string s)
328 this->StartNewStep("[Station"+station+"] Remove "+s);
329 MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
330 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
332 catch(clitk::ExceptionObject e) {
333 std::cout << s << " not found, skip." << std::endl;
336 //--------------------------------------------------------------------
339 //--------------------------------------------------------------------
340 template <class TImageType>
342 clitk::ExtractLymphStationsFilter<TImageType>::
343 SetFuzzyThreshold(std::string station, std::string tag, double value)
345 m_FuzzyThreshold[station][tag] = value;
347 //--------------------------------------------------------------------
350 //--------------------------------------------------------------------
351 template <class TImageType>
353 clitk::ExtractLymphStationsFilter<TImageType>::
354 SetThreshold(std::string station, std::string tag, double value)
356 m_Threshold[station][tag] = value;
358 //--------------------------------------------------------------------
361 //--------------------------------------------------------------------
362 template <class TImageType>
364 clitk::ExtractLymphStationsFilter<TImageType>::
365 GetFuzzyThreshold(std::string station, std::string tag)
367 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
368 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
372 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
373 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
377 return m_FuzzyThreshold[station][tag];
379 //--------------------------------------------------------------------
382 //--------------------------------------------------------------------
383 template <class TImageType>
385 clitk::ExtractLymphStationsFilter<TImageType>::
386 GetThreshold(std::string station, std::string tag)
388 if (m_Threshold.find(station) == m_Threshold.end()) {
389 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
393 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
394 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
398 return m_Threshold[station][tag];
400 //--------------------------------------------------------------------
403 //--------------------------------------------------------------------
404 template <class TImageType>
406 clitk::ExtractLymphStationsFilter<TImageType>::
407 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
409 // Create line from A to B with
410 // A = upper border of LLL at left
411 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
414 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
415 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
417 catch(clitk::ExceptionObject & o) {
419 DD("FindLineForS7S8Separation");
420 // Load LeftLowerLobeBronchus and get centroid point
421 MaskImagePointer LeftLowerLobeBronchus =
422 this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
423 std::vector<MaskImagePointType> c;
424 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
427 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
428 MaskImagePointer RightMiddleLobeBronchus =
429 this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
430 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
431 GetBackgroundValue(),
434 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
437 // Insert into the DB
438 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
439 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
442 //--------------------------------------------------------------------
445 //--------------------------------------------------------------------
446 template <class TImageType>
448 clitk::ExtractLymphStationsFilter<TImageType>::
453 z = this->GetAFDB()->GetDouble("CarinaZ");
455 catch(clitk::ExceptionObject e) {
456 //DD("FindCarinaSlicePosition");
458 MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
460 // Get Centroid and Z value
461 std::vector<MaskImagePointType> centroids;
462 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
464 // We dont need Carina structure from now
465 this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
467 // Put inside the AFDB
468 this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
469 this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
475 //--------------------------------------------------------------------
478 //--------------------------------------------------------------------
479 template <class TImageType>
481 clitk::ExtractLymphStationsFilter<TImageType>::
486 z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
488 catch(clitk::ExceptionObject e) {
491 //DD("FindApexOfTheChestPosition");
492 MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
493 MaskImagePointType p;
494 p[0] = p[1] = p[2] = 0.0; // to avoid warning
495 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
497 // We dont need Lungs structure from now
498 this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
500 // Put inside the AFDB
501 this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
502 p[2] -= 5; // We consider 5 mm lower
503 this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
508 // the superior border becomes the more inferior of the two apices
509 MaskImagePointer RightLung = this->GetAFDB()->template GetImage<MaskImageType>("RightLung");
510 MaskImagePointer LeftLung = this->GetAFDB()->template GetImage<MaskImageType>("LeftLung");
511 MaskImagePointType pr;
512 MaskImagePointType pl;
513 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(RightLung, GetBackgroundValue(), 2, false, pr);
514 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(LeftLung, GetBackgroundValue(), 2, false, pl);
515 // We dont need Lungs structure from now
516 this->GetAFDB()->template ReleaseImage<MaskImageType>("RightLung");
517 this->GetAFDB()->template ReleaseImage<MaskImageType>("LeftLung");
518 // Put inside the AFDB
519 if (pr[2] < pl[2]) z = pr[2];
521 this->GetAFDB()->SetDouble("ApexOfTheChestZ", z);
526 //--------------------------------------------------------------------
529 //--------------------------------------------------------------------
530 template <class TImageType>
532 clitk::ExtractLymphStationsFilter<TImageType>::
533 FindLeftAndRightBronchi()
536 m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
537 m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
539 catch(clitk::ExceptionObject & o) {
541 DD("FindLeftAndRightBronchi");
542 // The goal is to separate the trachea inferiorly to the carina into
543 // a Left and Right bronchus.
546 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
548 // Get the Carina position
549 double m_CarinaZ = FindCarina();
551 // Consider only inferiorly to the Carina
552 MaskImagePointer m_Working_Trachea =
553 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
554 GetBackgroundValue());
556 // Labelize the trachea
557 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
559 // Carina position must at the first slice that separate the two
560 // main bronchus (not superiorly). We thus first check that the
561 // upper slice is composed of at least two labels
562 MaskImagePointer RightBronchus;
563 MaskImagePointer LeftBronchus;
564 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
565 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
566 iter.SetFirstDirection(0);
567 iter.SetSecondDirection(1);
568 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
570 while (!iter.IsAtReverseEndOfSlice()) {
571 while (!iter.IsAtReverseEndOfLine()) {
572 if (iter.Get() > maxLabel) maxLabel = iter.Get();
578 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
581 // Compute 3D centroids of both parts to identify the left from the
583 std::vector<ImagePointType> c;
584 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
585 ImagePointType C1 = c[1];
586 ImagePointType C2 = c[2];
588 ImagePixelType rightLabel;
589 ImagePixelType leftLabel;
590 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
591 else { rightLabel = 2; leftLabel = 1; }
593 // Select LeftLabel (set one label to Backgroundvalue)
595 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
596 GetBackgroundValue(), GetForegroundValue());
598 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
599 leftLabel, GetBackgroundValue(), false);
601 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
602 GetBackgroundValue(), GetForegroundValue());
604 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
605 rightLabel, GetBackgroundValue(), false);
609 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
610 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
612 // Insert int AFDB if need after
613 this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
614 RightBronchus, true);
615 this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
619 //--------------------------------------------------------------------
622 //--------------------------------------------------------------------
623 template <class TImageType>
625 clitk::ExtractLymphStationsFilter<TImageType>::
626 FindSuperiorBorderOfAorticArch()
630 z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
632 catch(clitk::ExceptionObject e) {
633 // DD("FindSuperiorBorderOfAorticArch");
634 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
635 MaskImagePointType p;
636 p[0] = p[1] = p[2] = 0.0; // to avoid warning
637 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
638 p[2] += Aorta->GetSpacing()[2]; // the slice above
640 // We dont need Lungs structure from now
641 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
643 // Put inside the AFDB
644 this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
645 this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
651 //--------------------------------------------------------------------
654 //--------------------------------------------------------------------
655 template <class TImageType>
657 clitk::ExtractLymphStationsFilter<TImageType>::
658 FindInferiorBorderOfAorticArch()
662 z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
664 catch(clitk::ExceptionObject e) {
665 //DD("FindInferiorBorderOfAorticArch");
666 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
667 std::vector<MaskSlicePointer> slices;
668 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
670 uint i = slices.size()-1;
672 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
673 std::vector<typename MaskSliceType::PointType> c;
674 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
682 MaskImageIndexType index;
683 index[0] = index[1] = 0.0;
685 MaskImagePointType lower;
686 Aorta->TransformIndexToPhysicalPoint(index, lower);
688 // We dont need Lungs structure from now
689 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
691 // Put inside the AFDB
692 this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
693 this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
699 //--------------------------------------------------------------------
702 //--------------------------------------------------------------------
703 template <class ImageType>
704 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
705 clitk::ExtractLymphStationsFilter<ImageType>::
706 FindAntPostVesselsOLD()
708 // -----------------------------------------------------
709 /* Rod says: "The anterior border, as with the Atlas – UM, is
710 posterior to the vessels (right subclavian vein, left
711 brachiocephalic vein, right brachiocephalic vein, left subclavian
712 artery, left common carotid artery and brachiocephalic trunk).
713 These vessels are not included in the nodal station. The anterior
714 border is drawn to the midpoint of the vessel and an imaginary
715 line joins the middle of these vessels. Between the vessels,
716 station 2 is in contact with station 3a." */
718 // Check if no already done
720 MaskImagePointer binarizedContour;
722 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
724 catch(clitk::ExceptionObject e) {
728 return binarizedContour;
731 /* Here, we consider the vessels form a kind of anterior barrier. We
732 link all vessels centroids and remove what is post to it.
733 - select the list of structure
734 vessel1 = BrachioCephalicArtery
735 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
736 vessel3 = CommonCarotidArtery
737 vessel4 = SubclavianArtery
740 - crop images as needed
741 - slice by slice, choose the CCL and extract centroids
742 - slice by slice, sort according to polar angle wrt Trachea centroid.
743 - slice by slice, link centroids and close contour
744 - remove outside this contour
749 MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
750 MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
751 MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
752 MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
753 MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
754 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
755 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
757 // Create a temporay support
758 // From first slice of BrachioCephalicVein to end of 3A
759 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
760 MaskImagePointType p;
761 p[0] = p[1] = p[2] = 0.0; // to avoid warning
762 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
764 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
765 GetBackgroundValue(), 2, false, p);
767 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
768 false, GetBackgroundValue());
770 // Resize all structures like support
771 BrachioCephalicArtery =
772 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
773 CommonCarotidArtery =
774 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
776 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
778 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
780 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
781 BrachioCephalicVein =
782 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
784 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
787 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
788 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
789 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
790 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
791 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
792 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
793 std::vector<MaskSlicePointer> slices_SubclavianArtery;
794 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
795 std::vector<MaskSlicePointer> slices_Thyroid;
796 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
797 std::vector<MaskSlicePointer> slices_Aorta;
798 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
799 std::vector<MaskSlicePointer> slices_Trachea;
800 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
801 unsigned int n= slices_BrachioCephalicArtery.size();
803 // Get the boundaries of one slice
804 std::vector<MaskSlicePointType> bounds;
805 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
807 // For all slices, for all structures, find the centroid and build the contour
808 // List of 3D points (for debug)
809 std::vector<MaskImagePointType> p3D;
811 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
812 for(unsigned int i=0; i<n; i++) {
813 // Labelize the slices
814 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
815 GetBackgroundValue(), true, 1);
816 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
817 GetBackgroundValue(), true, 1);
818 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
819 GetBackgroundValue(), true, 1);
820 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
821 GetBackgroundValue(), true, 1);
822 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
823 GetBackgroundValue(), true, 1);
824 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
825 GetBackgroundValue(), true, 1);
828 std::vector<MaskSlicePointType> points2D;
829 std::vector<MaskSlicePointType> centroids1;
830 std::vector<MaskSlicePointType> centroids2;
831 std::vector<MaskSlicePointType> centroids3;
832 std::vector<MaskSlicePointType> centroids4;
833 std::vector<MaskSlicePointType> centroids5;
834 std::vector<MaskSlicePointType> centroids6;
835 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
836 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
837 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
838 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
839 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
840 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
842 // BrachioCephalicVein -> when it is separated into two CCL, we
843 // only consider the most at Right one
844 if (centroids6.size() > 2) {
845 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
846 else centroids6.erase(centroids6.begin()+1);
849 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
850 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
851 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
855 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
856 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
857 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
858 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
859 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
860 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
862 // Sort by angle according to trachea centroid and vertical line,
863 // in polar coordinates :
864 // http://en.wikipedia.org/wiki/Polar_coordinate_system
865 std::vector<MaskSlicePointType> centroids_trachea;
866 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
867 typedef std::pair<MaskSlicePointType, double> PointAngleType;
868 std::vector<PointAngleType> angles;
869 for(unsigned int j=0; j<points2D.size(); j++) {
870 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
871 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
872 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
874 if (x>0) angle = atan(y/x);
875 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
876 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
878 if (y>0) angle = M_PI/2.0;
879 if (y<0) angle = -M_PI/2.0;
882 angle = clitk::rad2deg(angle);
883 // Angle is [-180;180] wrt the X axis. We change the X axis to
884 // be the vertical line, because we want to sort from Right to
885 // Left from Post to Ant.
890 if (angle<0) angle = 360-angle;
893 a.first = points2D[j];
898 // Do nothing if less than 2 points --> n
899 if (points2D.size() < 3) { //continue;
903 // Sort points2D according to polar angles
904 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
905 for(unsigned int j=0; j<angles.size(); j++) {
906 points2D[j] = angles[j].first;
908 // DDV(points2D, points2D.size());
910 /* When vessels are far away, we try to replace the line segment
911 with a curved line that join the two vessels but stay
912 approximately at the same distance from the trachea centroids
916 - let a and c be two successive vessels centroids
917 - id distance(a,c) < threshold, next point
921 - compute mid position between two successive points -
922 compute dist to trachea centroid for the 3 pts - if middle too
925 std::vector<MaskSlicePointType> toadd;
926 unsigned int index = 0;
928 while (index<points2D.size()-1) {
929 MaskSlicePointType a = points2D[index];
930 MaskSlicePointType c = points2D[index+1];
932 double dac = a.EuclideanDistanceTo(c);
935 MaskSlicePointType b;
936 b[0] = a[0]+(c[0]-a[0])/2.0;
937 b[1] = a[1]+(c[1]-a[1])/2.0;
939 // Compute distance to trachea centroid
940 MaskSlicePointType m = centroids_trachea[1];
941 double da = m.EuclideanDistanceTo(a);
942 double db = m.EuclideanDistanceTo(b);
943 //double dc = m.EuclideanDistanceTo(c);
945 // Mean distance, find point on the line from trachea centroid
947 double alpha = (da+db)/2.0;
948 MaskSlicePointType v;
949 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
950 v[0] = (b[0]-m[0])/n;
951 v[1] = (b[1]-m[1])/n;
952 MaskSlicePointType r;
953 r[0] = m[0]+alpha*v[0];
954 r[1] = m[1]+alpha*v[1];
955 points2D.insert(points2D.begin()+index+1, r);
961 // DDV(points2D, points2D.size());
963 // Add some points to close the contour
964 // - H line towards Right
965 MaskSlicePointType p = points2D[0];
967 points2D.insert(points2D.begin(), p);
968 // - corner Right/Post
970 points2D.insert(points2D.begin(), p);
971 // - H line towards Left
974 points2D.push_back(p);
975 // - corner Right/Post
977 points2D.push_back(p);
978 // Close contour with the first point
979 points2D.push_back(points2D[0]);
980 // DDV(points2D, points2D.size());
982 // Build 3D points from the 2D points
983 std::vector<ImagePointType> points3D;
984 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
985 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
987 // Build the mesh from the contour's points
988 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
989 append->AddInput(mesh);
992 // DEBUG: write points3D
993 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
995 // Build the final 3D mesh form the list 2D mesh
997 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
999 // Debug, write the mesh
1001 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1003 w->SetFileName("bidon.vtk");
1007 // Compute a single binary 3D image from the list of contours
1008 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1009 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1010 filter->SetMesh(mesh);
1011 filter->SetLikeImage(support);
1013 binarizedContour = filter->GetOutput();
1016 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1019 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1022 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1023 false, GetBackgroundValue());
1025 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1026 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1028 return binarizedContour;
1031 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1032 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1033 p[1] += 50; // 50 mm Post from this point must be kept
1034 ImageIndexType index;
1035 binarizedContour->TransformPhysicalPointToIndex(p, index);
1036 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1038 // remove from support
1039 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1040 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1041 boolFilter->InPlaceOn();
1042 boolFilter->SetInput1(m_Working_Support);
1043 boolFilter->SetInput2(binarizedContour);
1044 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1045 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1047 boolFilter->SetOperationType(BoolFilterType::And);
1049 boolFilter->SetOperationType(BoolFilterType::AndNot);
1050 boolFilter->Update();
1051 m_Working_Support = boolFilter->GetOutput();
1055 //StopCurrentStep<MaskImageType>(m_Working_Support);
1056 //m_ListOfStations["2R"] = m_Working_Support;
1057 //m_ListOfStations["2L"] = m_Working_Support;
1059 //--------------------------------------------------------------------
1062 //--------------------------------------------------------------------
1063 template <class ImageType>
1064 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1065 clitk::ExtractLymphStationsFilter<ImageType>::
1066 FindAntPostVessels2()
1068 // -----------------------------------------------------
1069 /* Rod says: "The anterior border, as with the Atlas – UM, is
1070 posterior to the vessels (right subclavian vein, left
1071 brachiocephalic vein, right brachiocephalic vein, left subclavian
1072 artery, left common carotid artery and brachiocephalic trunk).
1073 These vessels are not included in the nodal station. The anterior
1074 border is drawn to the midpoint of the vessel and an imaginary
1075 line joins the middle of these vessels. Between the vessels,
1076 station 2 is in contact with station 3a." */
1078 // Check if no already done
1080 MaskImagePointer binarizedContour;
1082 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1084 catch(clitk::ExceptionObject e) {
1088 return binarizedContour;
1091 /* Here, we consider the vessels form a kind of anterior barrier. We
1092 link all vessels centroids and remove what is post to it.
1093 - select the list of structure
1094 vessel1 = BrachioCephalicArtery
1095 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1096 vessel3 = CommonCarotidArtery
1097 vessel4 = SubclavianArtery
1100 - crop images as needed
1101 - slice by slice, choose the CCL and extract centroids
1102 - slice by slice, sort according to polar angle wrt Trachea centroid.
1103 - slice by slice, link centroids and close contour
1104 - remove outside this contour
1105 - merge with support
1109 std::map<std::string, MaskImagePointer> MapOfStructures;
1110 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1111 MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1112 MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1113 MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftCommonCarotidArtery");
1114 MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightCommonCarotidArtery");
1115 MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftSubclavianArtery");
1116 MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightSubclavianArtery");
1117 MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1118 MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1119 MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1121 std::vector<std::string> ListOfStructuresNames;
1123 // Create a temporay support
1124 // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1125 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1126 MaskImagePointType p;
1127 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1128 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1129 GetBackgroundValue(), 2, true, p);
1131 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1132 GetBackgroundValue(), 2, false, p);
1133 MaskImagePointType p2;
1134 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1135 GetBackgroundValue(), 2, false, p2);
1136 if (p2[2] > p[2]) p = p2;
1138 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1139 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1140 false, GetBackgroundValue());
1142 // Resize all structures like support
1143 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1144 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1148 typedef std::vector<MaskSlicePointer> SlicesType;
1149 std::map<std::string, SlicesType> MapOfSlices;
1150 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1152 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1153 MapOfSlices[it->first] = s;
1156 unsigned int n= MapOfSlices["Trachea"].size();
1158 // Get the boundaries of one slice
1159 std::vector<MaskSlicePointType> bounds;
1160 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1163 // For all slices, for all structures, find the centroid and build the contour
1164 // List of 3D points (for debug)
1165 std::vector<MaskImagePointType> p3D;
1166 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1167 for(unsigned int i=0; i<n; i++) {
1169 // Labelize the slices
1170 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1171 MaskSlicePointer & s = MapOfSlices[it->first][i];
1172 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1176 std::vector<MaskSlicePointType> points2D;
1177 typedef std::vector<MaskSlicePointType> CentroidsType;
1178 std::map<std::string, CentroidsType> MapOfCentroids;
1179 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1180 std::string structure = it->first;
1181 MaskSlicePointer & s = MapOfSlices[structure][i];
1183 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1184 MapOfCentroids[structure] = c;
1188 // BrachioCephalicVein -> when it is separated into two CCL, we
1189 // only consider the most at Right one
1190 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1192 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1193 else c.erase(c.begin()+1);
1197 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1198 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1199 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1200 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1201 MapOfCentroids["BrachioCephalicVein"].clear();
1205 // Add all 2D points
1206 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1207 std::string structure = it->first;
1208 if (structure != "Trachea") { // do not add centroid of trachea
1209 CentroidsType & c = MapOfCentroids[structure];
1210 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1214 // Sort by angle according to trachea centroid and vertical line,
1215 // in polar coordinates :
1216 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1217 // std::vector<MaskSlicePointType> centroids_trachea;
1218 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1219 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1220 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1221 std::vector<PointAngleType> angles;
1222 for(unsigned int j=0; j<points2D.size(); j++) {
1223 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1224 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1225 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1227 if (x>0) angle = atan(y/x);
1228 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1229 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1231 if (y>0) angle = M_PI/2.0;
1232 if (y<0) angle = -M_PI/2.0;
1233 if (y==0) angle = 0;
1235 angle = clitk::rad2deg(angle);
1236 // Angle is [-180;180] wrt the X axis. We change the X axis to
1237 // be the vertical line, because we want to sort from Right to
1238 // Left from Post to Ant.
1240 angle = (270-angle);
1243 if (angle<0) angle = 360-angle;
1246 a.first = points2D[j];
1248 angles.push_back(a);
1251 // Do nothing if less than 2 points --> n
1252 if (points2D.size() < 3) { //continue;
1256 // Sort points2D according to polar angles
1257 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1258 for(unsigned int j=0; j<angles.size(); j++) {
1259 points2D[j] = angles[j].first;
1261 // DDV(points2D, points2D.size());
1263 /* When vessels are far away, we try to replace the line segment
1264 with a curved line that join the two vessels but stay
1265 approximately at the same distance from the trachea centroids
1269 - let a and c be two successive vessels centroids
1270 - id distance(a,c) < threshold, next point
1274 - compute mid position between two successive points -
1275 compute dist to trachea centroid for the 3 pts - if middle too
1278 std::vector<MaskSlicePointType> toadd;
1279 unsigned int index = 0;
1281 while (index<points2D.size()-1) {
1282 MaskSlicePointType a = points2D[index];
1283 MaskSlicePointType c = points2D[index+1];
1285 double dac = a.EuclideanDistanceTo(c);
1288 MaskSlicePointType b;
1289 b[0] = a[0]+(c[0]-a[0])/2.0;
1290 b[1] = a[1]+(c[1]-a[1])/2.0;
1292 // Compute distance to trachea centroid
1293 MaskSlicePointType m = centroids_trachea[1];
1294 double da = m.EuclideanDistanceTo(a);
1295 double db = m.EuclideanDistanceTo(b);
1296 //double dc = m.EuclideanDistanceTo(c);
1298 // Mean distance, find point on the line from trachea centroid
1300 double alpha = (da+db)/2.0;
1301 MaskSlicePointType v;
1302 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1303 v[0] = (b[0]-m[0])/n;
1304 v[1] = (b[1]-m[1])/n;
1305 MaskSlicePointType r;
1306 r[0] = m[0]+alpha*v[0];
1307 r[1] = m[1]+alpha*v[1];
1308 points2D.insert(points2D.begin()+index+1, r);
1314 // DDV(points2D, points2D.size());
1316 // Add some points to close the contour
1317 // - H line towards Right
1318 MaskSlicePointType p = points2D[0];
1319 p[0] = bounds[0][0];
1320 points2D.insert(points2D.begin(), p);
1321 // - corner Right/Post
1323 points2D.insert(points2D.begin(), p);
1324 // - H line towards Left
1325 p = points2D.back();
1326 p[0] = bounds[2][0];
1327 points2D.push_back(p);
1328 // - corner Right/Post
1330 points2D.push_back(p);
1331 // Close contour with the first point
1332 points2D.push_back(points2D[0]);
1333 // DDV(points2D, points2D.size());
1335 // Build 3D points from the 2D points
1336 std::vector<ImagePointType> points3D;
1337 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1338 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1340 // Build the mesh from the contour's points
1341 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1342 append->AddInput(mesh);
1343 // if (i ==n-1) { // last slice
1344 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1345 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1346 // append->AddInput(mesh);
1350 // DEBUG: write points3D
1351 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1353 // Build the final 3D mesh form the list 2D mesh
1355 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1357 // Debug, write the mesh
1359 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1361 w->SetFileName("bidon.vtk");
1365 // Compute a single binary 3D image from the list of contours
1366 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1367 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1368 filter->SetMesh(mesh);
1369 filter->SetLikeImage(support);
1371 binarizedContour = filter->GetOutput();
1374 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1375 GetForegroundValue(), 2, true, p);
1377 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1378 GetForegroundValue(), 2, false, p);
1379 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1380 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1381 false, GetBackgroundValue());
1384 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1385 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1387 return binarizedContour;
1390 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1391 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1392 p[1] += 50; // 50 mm Post from this point must be kept
1393 ImageIndexType index;
1394 binarizedContour->TransformPhysicalPointToIndex(p, index);
1395 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1397 // remove from support
1398 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1399 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1400 boolFilter->InPlaceOn();
1401 boolFilter->SetInput1(m_Working_Support);
1402 boolFilter->SetInput2(binarizedContour);
1403 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1404 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1406 boolFilter->SetOperationType(BoolFilterType::And);
1408 boolFilter->SetOperationType(BoolFilterType::AndNot);
1409 boolFilter->Update();
1410 m_Working_Support = boolFilter->GetOutput();
1414 //StopCurrentStep<MaskImageType>(m_Working_Support);
1415 //m_ListOfStations["2R"] = m_Working_Support;
1416 //m_ListOfStations["2L"] = m_Working_Support;
1418 //--------------------------------------------------------------------
1420 //--------------------------------------------------------------------
1421 template <class ImageType>
1423 clitk::ExtractLymphStationsFilter<ImageType>::
1424 WriteImageSupport(std::string support)
1426 writeImage<MaskImageType>(m_ListOfSupports[support], this->GetAFDBPath()+"/"+"seg/Support_"+support+".mha");
1427 this->GetAFDB()->SetImageFilename("Support_"+support, "seg/Support_"+support+".mha");
1429 //--------------------------------------------------------------------
1432 //--------------------------------------------------------------------
1433 template <class ImageType>
1435 clitk::ExtractLymphStationsFilter<ImageType>::
1436 WriteImageStation(std::string station)
1438 writeImage<MaskImageType>(m_ListOfStations[station], GetAFDB()->GetPath()+"/seg/Station"+station+".mha");
1439 GetAFDB()->SetImageFilename("Station"+station, "seg/Station"+station+".mha");
1442 //--------------------------------------------------------------------
1445 //--------------------------------------------------------------------
1446 template <class ImageType>
1448 clitk::ExtractLymphStationsFilter<ImageType>::
1449 ComputeOverlapWithRef(std::string station)
1451 if (GetComputeStation(station)) {
1452 MaskImagePointer ref = this->GetAFDB()->template GetImage <MaskImageType>("Station"+station+"_Ref");
1453 typedef clitk::LabelImageOverlapMeasureFilter<MaskImageType> FilterType;
1454 typename FilterType::Pointer filter = FilterType::New();
1455 filter->SetInput(0, m_ListOfStations[station]);
1456 filter->SetInput(1, ref);
1460 //--------------------------------------------------------------------
1462 //--------------------------------------------------------------------
1463 template <class ImageType>
1465 clitk::ExtractLymphStationsFilter<ImageType>::
1466 ReadSupportLimits(std::string filename)
1468 m_ListOfSupportLimits.clear();
1470 openFileForReading(is, filename);
1473 SupportLimitsType s;
1475 if (is) m_ListOfSupportLimits.push_back(s);
1478 //--------------------------------------------------------------------
1480 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX