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 ComputeStationsSupportsFlagOn();
61 ExtractStation_3P_SetDefaultValues();
62 ExtractStation_2RL_SetDefaultValues();
63 ExtractStation_3A_SetDefaultValues();
64 ExtractStation_1RL_SetDefaultValues();
66 ExtractStation_7_SetDefaultValues();
67 ExtractStation_8_SetDefaultValues();
69 //--------------------------------------------------------------------
72 //--------------------------------------------------------------------
73 template <class TImageType>
75 clitk::ExtractLymphStationsFilter<TImageType>::
76 GenerateOutputInformation() {
79 m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
80 m_Mediastinum = this->GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
82 // Clean some computer landmarks to force the recomputation
83 this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
85 // Global supports for stations
86 bool supportsExist = true;
88 m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
89 m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
90 m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
91 m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
92 m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
93 m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
95 m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
96 m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
97 m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
98 m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
99 m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
100 m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
101 m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
102 m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
103 m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
104 } catch(clitk::ExceptionObject o) {
105 supportsExist = false;
108 if (!supportsExist || GetComputeStationsSupportsFlag()) {
109 this->StartNewStep("Supports for stations");
110 this->StartSubStep();
111 this->GetAFDB()->RemoveTag("CarinaZ");
112 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
113 this->GetAFDB()->RemoveTag("ApexOfTheChest");
114 this->GetAFDB()->RemoveTag("RightBronchus");
115 this->GetAFDB()->RemoveTag("LeftBronchus");
116 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
117 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
118 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
119 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
120 ExtractStationSupports();
130 // Extract Station2RL
131 ExtractStation_2RL();
133 // Extract Station1RL
134 ExtractStation_1RL();
136 // ---------- TODO -----------------------
142 this->StartNewStep("Station 7");
143 this->StartSubStep();
147 if (0) { // temporary suppress
148 // Extract Station4RL
149 this->StartNewStep("Station 4RL");
150 this->StartSubStep();
151 //ExtractStation_4RL();
155 //--------------------------------------------------------------------
158 //--------------------------------------------------------------------
159 template <class TImageType>
161 clitk::ExtractLymphStationsFilter<TImageType>::
162 GenerateInputRequestedRegion() {
163 //DD("GenerateInputRequestedRegion (nothing?)");
165 //--------------------------------------------------------------------
168 //--------------------------------------------------------------------
169 template <class TImageType>
171 clitk::ExtractLymphStationsFilter<TImageType>::
173 // Final Step -> graft output (if SetNthOutput => redo)
174 this->GraftOutput(m_ListOfStations["8"]);
176 //--------------------------------------------------------------------
179 //--------------------------------------------------------------------
180 template <class TImageType>
182 clitk::ExtractLymphStationsFilter<TImageType>::
183 CheckForStation(std::string station)
185 // Compute Station name
186 std::string s = "Station"+station;
189 // Check if station already exist in DB
191 if (this->GetAFDB()->TagExist(s)) {
192 m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
196 // Define the starting support
197 if (found && GetComputeStation(station)) {
198 std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
200 if (!found || GetComputeStation(station)) {
201 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
205 std::cout << "Station " << station << " found. I used it" << std::endl;
209 //--------------------------------------------------------------------
212 //--------------------------------------------------------------------
213 template <class TImageType>
215 clitk::ExtractLymphStationsFilter<TImageType>::
216 GetComputeStation(std::string station)
218 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
220 //--------------------------------------------------------------------
223 //--------------------------------------------------------------------
224 template <class TImageType>
226 clitk::ExtractLymphStationsFilter<TImageType>::
227 AddComputeStation(std::string station)
229 m_ComputeStationMap[station] = true;
231 //--------------------------------------------------------------------
235 //--------------------------------------------------------------------
236 template<class PointType>
237 class comparePointsX {
239 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
241 //--------------------------------------------------------------------
244 //--------------------------------------------------------------------
245 template<class PairType>
246 class comparePointsWithAngle {
248 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
250 //--------------------------------------------------------------------
253 //--------------------------------------------------------------------
255 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
256 std::vector<itk::Point<double, Dim-1> > previous;
257 HypercubeCorners<Dim-1>(previous);
258 out.resize(previous.size()*2);
259 for(unsigned int i=0; i<out.size(); i++) {
260 itk::Point<double, Dim> p;
261 if (i<previous.size()) p[0] = 0;
263 for(int j=0; j<Dim-1; j++)
265 p[j+1] = previous[i%previous.size()][j];
270 //--------------------------------------------------------------------
273 //--------------------------------------------------------------------
275 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
280 //--------------------------------------------------------------------
283 //--------------------------------------------------------------------
284 template<class ImageType>
285 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
286 std::vector<typename ImageType::PointType> & bounds)
288 // Get image max/min coordinates
289 const unsigned int dim=ImageType::ImageDimension;
290 typedef typename ImageType::PointType PointType;
291 typedef typename ImageType::IndexType IndexType;
292 PointType min_c, max_c;
293 IndexType min_i, max_i;
294 min_i = image->GetLargestPossibleRegion().GetIndex();
295 for(unsigned int i=0; i<dim; i++)
296 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
297 image->TransformIndexToPhysicalPoint(min_i, min_c);
298 image->TransformIndexToPhysicalPoint(max_i, max_c);
300 // Get corners coordinates
301 HypercubeCorners<ImageType::ImageDimension>(bounds);
302 for(unsigned int i=0; i<bounds.size(); i++) {
303 for(unsigned int j=0; j<dim; j++) {
304 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
305 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
309 //--------------------------------------------------------------------
312 //--------------------------------------------------------------------
313 template <class TImageType>
315 clitk::ExtractLymphStationsFilter<TImageType>::
316 ExtractStation_4RL() {
320 WARNING ONLY 4R FIRST !!! (not same inf limits)
322 ExtractStation_4RL_SI_Limits();
323 ExtractStation_4RL_LR_Limits();
324 ExtractStation_4RL_AP_Limits();
326 //--------------------------------------------------------------------
329 //--------------------------------------------------------------------
330 template <class ImageType>
332 clitk::ExtractLymphStationsFilter<ImageType>::
333 Remove_Structures(std::string station, std::string s)
336 this->StartNewStep("[Station"+station+"] Remove "+s);
337 MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
338 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
340 catch(clitk::ExceptionObject e) {
341 std::cout << s << " not found, skip." << std::endl;
344 //--------------------------------------------------------------------
347 //--------------------------------------------------------------------
348 template <class TImageType>
350 clitk::ExtractLymphStationsFilter<TImageType>::
351 SetFuzzyThreshold(std::string station, std::string tag, double value)
353 m_FuzzyThreshold[station][tag] = value;
355 //--------------------------------------------------------------------
358 //--------------------------------------------------------------------
359 template <class TImageType>
361 clitk::ExtractLymphStationsFilter<TImageType>::
362 SetThreshold(std::string station, std::string tag, double value)
364 m_Threshold[station][tag] = value;
366 //--------------------------------------------------------------------
369 //--------------------------------------------------------------------
370 template <class TImageType>
372 clitk::ExtractLymphStationsFilter<TImageType>::
373 GetFuzzyThreshold(std::string station, std::string tag)
375 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
376 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
380 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
381 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
385 return m_FuzzyThreshold[station][tag];
387 //--------------------------------------------------------------------
390 //--------------------------------------------------------------------
391 template <class TImageType>
393 clitk::ExtractLymphStationsFilter<TImageType>::
394 GetThreshold(std::string station, std::string tag)
396 if (m_Threshold.find(station) == m_Threshold.end()) {
397 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
401 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
402 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
406 return m_Threshold[station][tag];
408 //--------------------------------------------------------------------
411 //--------------------------------------------------------------------
412 template <class TImageType>
414 clitk::ExtractLymphStationsFilter<TImageType>::
415 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
417 // Create line from A to B with
418 // A = upper border of LLL at left
419 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
422 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
423 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
425 catch(clitk::ExceptionObject & o) {
427 DD("FindLineForS7S8Separation");
428 // Load LeftLowerLobeBronchus and get centroid point
429 MaskImagePointer LeftLowerLobeBronchus =
430 this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
431 std::vector<MaskImagePointType> c;
432 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
435 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
436 MaskImagePointer RightMiddleLobeBronchus =
437 this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
438 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
439 GetBackgroundValue(),
442 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
445 // Insert into the DB
446 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
447 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
450 //--------------------------------------------------------------------
453 //--------------------------------------------------------------------
454 template <class TImageType>
456 clitk::ExtractLymphStationsFilter<TImageType>::
461 z = this->GetAFDB()->GetDouble("CarinaZ");
463 catch(clitk::ExceptionObject e) {
464 DD("FindCarinaSlicePosition");
466 MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
468 // Get Centroid and Z value
469 std::vector<MaskImagePointType> centroids;
470 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
472 // We dont need Carina structure from now
473 this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
475 // Put inside the AFDB
476 this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
477 this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
483 //--------------------------------------------------------------------
486 //--------------------------------------------------------------------
487 template <class TImageType>
489 clitk::ExtractLymphStationsFilter<TImageType>::
494 z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
496 catch(clitk::ExceptionObject e) {
497 DD("FindApexOfTheChestPosition");
498 MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
499 MaskImagePointType p;
500 p[0] = p[1] = p[2] = 0.0; // to avoid warning
501 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
503 // We dont need Lungs structure from now
504 this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
506 // Put inside the AFDB
507 this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
508 p[2] -= 5; // We consider 5 mm lower
509 this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
515 //--------------------------------------------------------------------
518 //--------------------------------------------------------------------
519 template <class TImageType>
521 clitk::ExtractLymphStationsFilter<TImageType>::
522 FindLeftAndRightBronchi()
525 m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
526 m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
528 catch(clitk::ExceptionObject & o) {
530 DD("FindLeftAndRightBronchi");
531 // The goal is to separate the trachea inferiorly to the carina into
532 // a Left and Right bronchus.
535 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
537 // Get the Carina position
538 double m_CarinaZ = FindCarina();
540 // Consider only inferiorly to the Carina
541 MaskImagePointer m_Working_Trachea =
542 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
543 GetBackgroundValue());
545 // Labelize the trachea
546 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
548 // Carina position must at the first slice that separate the two
549 // main bronchus (not superiorly). We thus first check that the
550 // upper slice is composed of at least two labels
551 MaskImagePointer RightBronchus;
552 MaskImagePointer LeftBronchus;
553 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
554 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
555 iter.SetFirstDirection(0);
556 iter.SetSecondDirection(1);
557 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
559 while (!iter.IsAtReverseEndOfSlice()) {
560 while (!iter.IsAtReverseEndOfLine()) {
561 if (iter.Get() > maxLabel) maxLabel = iter.Get();
567 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
570 // Compute 3D centroids of both parts to identify the left from the
572 std::vector<ImagePointType> c;
573 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
574 ImagePointType C1 = c[1];
575 ImagePointType C2 = c[2];
577 ImagePixelType rightLabel;
578 ImagePixelType leftLabel;
579 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
580 else { rightLabel = 2; leftLabel = 1; }
582 // Select LeftLabel (set one label to Backgroundvalue)
584 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
585 GetBackgroundValue(), GetForegroundValue());
587 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
588 leftLabel, GetBackgroundValue(), false);
590 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
591 GetBackgroundValue(), GetForegroundValue());
593 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
594 rightLabel, GetBackgroundValue(), false);
598 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
599 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
601 // Insert int AFDB if need after
602 this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
603 RightBronchus, true);
604 this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
608 //--------------------------------------------------------------------
611 //--------------------------------------------------------------------
612 template <class TImageType>
614 clitk::ExtractLymphStationsFilter<TImageType>::
615 FindSuperiorBorderOfAorticArch()
619 z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
621 catch(clitk::ExceptionObject e) {
622 DD("FindSuperiorBorderOfAorticArch");
623 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
624 MaskImagePointType p;
625 p[0] = p[1] = p[2] = 0.0; // to avoid warning
626 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
627 p[2] += Aorta->GetSpacing()[2]; // the slice above
629 // We dont need Lungs structure from now
630 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
632 // Put inside the AFDB
633 this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
634 this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
640 //--------------------------------------------------------------------
643 //--------------------------------------------------------------------
644 template <class TImageType>
646 clitk::ExtractLymphStationsFilter<TImageType>::
647 FindInferiorBorderOfAorticArch()
651 z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
653 catch(clitk::ExceptionObject e) {
654 DD("FindInferiorBorderOfAorticArch");
655 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
656 std::vector<MaskSlicePointer> slices;
657 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
659 uint i = slices.size()-1;
661 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
662 std::vector<typename MaskSliceType::PointType> c;
663 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
671 MaskImageIndexType index;
672 index[0] = index[1] = 0.0;
674 MaskImagePointType lower;
675 Aorta->TransformIndexToPhysicalPoint(index, lower);
677 // We dont need Lungs structure from now
678 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
680 // Put inside the AFDB
681 this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
682 this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
688 //--------------------------------------------------------------------
691 //--------------------------------------------------------------------
692 template <class ImageType>
693 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
694 clitk::ExtractLymphStationsFilter<ImageType>::
695 FindAntPostVesselsOLD()
697 // -----------------------------------------------------
698 /* Rod says: "The anterior border, as with the Atlas – UM, is
699 posterior to the vessels (right subclavian vein, left
700 brachiocephalic vein, right brachiocephalic vein, left subclavian
701 artery, left common carotid artery and brachiocephalic trunk).
702 These vessels are not included in the nodal station. The anterior
703 border is drawn to the midpoint of the vessel and an imaginary
704 line joins the middle of these vessels. Between the vessels,
705 station 2 is in contact with station 3a." */
707 // Check if no already done
709 MaskImagePointer binarizedContour;
711 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
713 catch(clitk::ExceptionObject e) {
717 return binarizedContour;
720 /* Here, we consider the vessels form a kind of anterior barrier. We
721 link all vessels centroids and remove what is post to it.
722 - select the list of structure
723 vessel1 = BrachioCephalicArtery
724 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
725 vessel3 = CommonCarotidArtery
726 vessel4 = SubclavianArtery
729 - crop images as needed
730 - slice by slice, choose the CCL and extract centroids
731 - slice by slice, sort according to polar angle wrt Trachea centroid.
732 - slice by slice, link centroids and close contour
733 - remove outside this contour
738 MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
739 MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
740 MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
741 MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
742 MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
743 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
744 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
746 // Create a temporay support
747 // From first slice of BrachioCephalicVein to end of 3A
748 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
749 MaskImagePointType p;
750 p[0] = p[1] = p[2] = 0.0; // to avoid warning
751 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
753 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
754 GetBackgroundValue(), 2, false, p);
756 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
757 false, GetBackgroundValue());
759 // Resize all structures like support
760 BrachioCephalicArtery =
761 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
762 CommonCarotidArtery =
763 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
765 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
767 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
769 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
770 BrachioCephalicVein =
771 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
773 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
776 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
777 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
778 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
779 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
780 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
781 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
782 std::vector<MaskSlicePointer> slices_SubclavianArtery;
783 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
784 std::vector<MaskSlicePointer> slices_Thyroid;
785 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
786 std::vector<MaskSlicePointer> slices_Aorta;
787 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
788 std::vector<MaskSlicePointer> slices_Trachea;
789 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
790 unsigned int n= slices_BrachioCephalicArtery.size();
792 // Get the boundaries of one slice
793 std::vector<MaskSlicePointType> bounds;
794 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
796 // For all slices, for all structures, find the centroid and build the contour
797 // List of 3D points (for debug)
798 std::vector<MaskImagePointType> p3D;
800 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
801 for(unsigned int i=0; i<n; i++) {
802 // Labelize the slices
803 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
804 GetBackgroundValue(), true, 1);
805 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
806 GetBackgroundValue(), true, 1);
807 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
808 GetBackgroundValue(), true, 1);
809 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
810 GetBackgroundValue(), true, 1);
811 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
812 GetBackgroundValue(), true, 1);
813 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
814 GetBackgroundValue(), true, 1);
817 std::vector<MaskSlicePointType> points2D;
818 std::vector<MaskSlicePointType> centroids1;
819 std::vector<MaskSlicePointType> centroids2;
820 std::vector<MaskSlicePointType> centroids3;
821 std::vector<MaskSlicePointType> centroids4;
822 std::vector<MaskSlicePointType> centroids5;
823 std::vector<MaskSlicePointType> centroids6;
824 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
825 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
826 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
827 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
828 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
829 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
831 // BrachioCephalicVein -> when it is separated into two CCL, we
832 // only consider the most at Right one
833 if (centroids6.size() > 2) {
834 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
835 else centroids6.erase(centroids6.begin()+1);
838 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
839 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
840 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
844 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
845 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
846 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
847 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
848 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
849 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
851 // Sort by angle according to trachea centroid and vertical line,
852 // in polar coordinates :
853 // http://en.wikipedia.org/wiki/Polar_coordinate_system
854 std::vector<MaskSlicePointType> centroids_trachea;
855 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
856 typedef std::pair<MaskSlicePointType, double> PointAngleType;
857 std::vector<PointAngleType> angles;
858 for(unsigned int j=0; j<points2D.size(); j++) {
859 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
860 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
861 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
863 if (x>0) angle = atan(y/x);
864 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
865 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
867 if (y>0) angle = M_PI/2.0;
868 if (y<0) angle = -M_PI/2.0;
871 angle = clitk::rad2deg(angle);
872 // Angle is [-180;180] wrt the X axis. We change the X axis to
873 // be the vertical line, because we want to sort from Right to
874 // Left from Post to Ant.
879 if (angle<0) angle = 360-angle;
882 a.first = points2D[j];
887 // Do nothing if less than 2 points --> n
888 if (points2D.size() < 3) { //continue;
892 // Sort points2D according to polar angles
893 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
894 for(unsigned int j=0; j<angles.size(); j++) {
895 points2D[j] = angles[j].first;
897 // DDV(points2D, points2D.size());
899 /* When vessels are far away, we try to replace the line segment
900 with a curved line that join the two vessels but stay
901 approximately at the same distance from the trachea centroids
905 - let a and c be two successive vessels centroids
906 - id distance(a,c) < threshold, next point
910 - compute mid position between two successive points -
911 compute dist to trachea centroid for the 3 pts - if middle too
914 std::vector<MaskSlicePointType> toadd;
915 unsigned int index = 0;
917 while (index<points2D.size()-1) {
918 MaskSlicePointType a = points2D[index];
919 MaskSlicePointType c = points2D[index+1];
921 double dac = a.EuclideanDistanceTo(c);
924 MaskSlicePointType b;
925 b[0] = a[0]+(c[0]-a[0])/2.0;
926 b[1] = a[1]+(c[1]-a[1])/2.0;
928 // Compute distance to trachea centroid
929 MaskSlicePointType m = centroids_trachea[1];
930 double da = m.EuclideanDistanceTo(a);
931 double db = m.EuclideanDistanceTo(b);
932 //double dc = m.EuclideanDistanceTo(c);
934 // Mean distance, find point on the line from trachea centroid
936 double alpha = (da+db)/2.0;
937 MaskSlicePointType v;
938 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
939 v[0] = (b[0]-m[0])/n;
940 v[1] = (b[1]-m[1])/n;
941 MaskSlicePointType r;
942 r[0] = m[0]+alpha*v[0];
943 r[1] = m[1]+alpha*v[1];
944 points2D.insert(points2D.begin()+index+1, r);
950 // DDV(points2D, points2D.size());
952 // Add some points to close the contour
953 // - H line towards Right
954 MaskSlicePointType p = points2D[0];
956 points2D.insert(points2D.begin(), p);
957 // - corner Right/Post
959 points2D.insert(points2D.begin(), p);
960 // - H line towards Left
963 points2D.push_back(p);
964 // - corner Right/Post
966 points2D.push_back(p);
967 // Close contour with the first point
968 points2D.push_back(points2D[0]);
969 // DDV(points2D, points2D.size());
971 // Build 3D points from the 2D points
972 std::vector<ImagePointType> points3D;
973 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
974 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
976 // Build the mesh from the contour's points
977 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
978 append->AddInput(mesh);
981 // DEBUG: write points3D
982 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
984 // Build the final 3D mesh form the list 2D mesh
986 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
988 // Debug, write the mesh
990 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
992 w->SetFileName("bidon.vtk");
996 // Compute a single binary 3D image from the list of contours
997 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
998 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
999 filter->SetMesh(mesh);
1000 filter->SetLikeImage(support);
1002 binarizedContour = filter->GetOutput();
1005 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1008 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1011 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1012 false, GetBackgroundValue());
1014 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1015 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1017 return binarizedContour;
1020 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1021 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1022 p[1] += 50; // 50 mm Post from this point must be kept
1023 ImageIndexType index;
1024 binarizedContour->TransformPhysicalPointToIndex(p, index);
1025 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1027 // remove from support
1028 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1029 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1030 boolFilter->InPlaceOn();
1031 boolFilter->SetInput1(m_Working_Support);
1032 boolFilter->SetInput2(binarizedContour);
1033 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1034 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1036 boolFilter->SetOperationType(BoolFilterType::And);
1038 boolFilter->SetOperationType(BoolFilterType::AndNot);
1039 boolFilter->Update();
1040 m_Working_Support = boolFilter->GetOutput();
1044 //StopCurrentStep<MaskImageType>(m_Working_Support);
1045 //m_ListOfStations["2R"] = m_Working_Support;
1046 //m_ListOfStations["2L"] = m_Working_Support;
1048 //--------------------------------------------------------------------
1051 //--------------------------------------------------------------------
1052 template <class ImageType>
1053 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1054 clitk::ExtractLymphStationsFilter<ImageType>::
1055 FindAntPostVessels2()
1057 // -----------------------------------------------------
1058 /* Rod says: "The anterior border, as with the Atlas – UM, is
1059 posterior to the vessels (right subclavian vein, left
1060 brachiocephalic vein, right brachiocephalic vein, left subclavian
1061 artery, left common carotid artery and brachiocephalic trunk).
1062 These vessels are not included in the nodal station. The anterior
1063 border is drawn to the midpoint of the vessel and an imaginary
1064 line joins the middle of these vessels. Between the vessels,
1065 station 2 is in contact with station 3a." */
1067 // Check if no already done
1069 MaskImagePointer binarizedContour;
1071 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1073 catch(clitk::ExceptionObject e) {
1077 return binarizedContour;
1080 /* Here, we consider the vessels form a kind of anterior barrier. We
1081 link all vessels centroids and remove what is post to it.
1082 - select the list of structure
1083 vessel1 = BrachioCephalicArtery
1084 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1085 vessel3 = CommonCarotidArtery
1086 vessel4 = SubclavianArtery
1089 - crop images as needed
1090 - slice by slice, choose the CCL and extract centroids
1091 - slice by slice, sort according to polar angle wrt Trachea centroid.
1092 - slice by slice, link centroids and close contour
1093 - remove outside this contour
1094 - merge with support
1098 std::map<std::string, MaskImagePointer> MapOfStructures;
1099 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1100 MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1101 MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1102 MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1103 MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1104 MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1105 MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1106 MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1107 MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1108 MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1110 std::vector<std::string> ListOfStructuresNames;
1112 // Create a temporay support
1113 // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1114 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1115 MaskImagePointType p;
1116 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1117 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1118 GetBackgroundValue(), 2, true, p);
1120 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1121 GetBackgroundValue(), 2, false, p);
1122 MaskImagePointType p2;
1123 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1124 GetBackgroundValue(), 2, false, p2);
1125 if (p2[2] > p[2]) p = p2;
1127 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1128 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1129 false, GetBackgroundValue());
1131 // Resize all structures like support
1132 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1133 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1137 typedef std::vector<MaskSlicePointer> SlicesType;
1138 std::map<std::string, SlicesType> MapOfSlices;
1139 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1141 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1142 MapOfSlices[it->first] = s;
1145 unsigned int n= MapOfSlices["Trachea"].size();
1147 // Get the boundaries of one slice
1148 std::vector<MaskSlicePointType> bounds;
1149 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1152 // For all slices, for all structures, find the centroid and build the contour
1153 // List of 3D points (for debug)
1154 std::vector<MaskImagePointType> p3D;
1155 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1156 for(unsigned int i=0; i<n; i++) {
1158 // Labelize the slices
1159 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1160 MaskSlicePointer & s = MapOfSlices[it->first][i];
1161 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1165 std::vector<MaskSlicePointType> points2D;
1166 typedef std::vector<MaskSlicePointType> CentroidsType;
1167 std::map<std::string, CentroidsType> MapOfCentroids;
1168 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1169 std::string structure = it->first;
1170 MaskSlicePointer & s = MapOfSlices[structure][i];
1172 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1173 MapOfCentroids[structure] = c;
1177 // BrachioCephalicVein -> when it is separated into two CCL, we
1178 // only consider the most at Right one
1179 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1181 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1182 else c.erase(c.begin()+1);
1186 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1187 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1188 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1189 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1190 MapOfCentroids["BrachioCephalicVein"].clear();
1194 // Add all 2D points
1195 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1196 std::string structure = it->first;
1197 if (structure != "Trachea") { // do not add centroid of trachea
1198 CentroidsType & c = MapOfCentroids[structure];
1199 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1203 // Sort by angle according to trachea centroid and vertical line,
1204 // in polar coordinates :
1205 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1206 // std::vector<MaskSlicePointType> centroids_trachea;
1207 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1208 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1209 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1210 std::vector<PointAngleType> angles;
1211 for(unsigned int j=0; j<points2D.size(); j++) {
1212 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1213 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1214 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1216 if (x>0) angle = atan(y/x);
1217 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1218 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1220 if (y>0) angle = M_PI/2.0;
1221 if (y<0) angle = -M_PI/2.0;
1222 if (y==0) angle = 0;
1224 angle = clitk::rad2deg(angle);
1225 // Angle is [-180;180] wrt the X axis. We change the X axis to
1226 // be the vertical line, because we want to sort from Right to
1227 // Left from Post to Ant.
1229 angle = (270-angle);
1232 if (angle<0) angle = 360-angle;
1235 a.first = points2D[j];
1237 angles.push_back(a);
1240 // Do nothing if less than 2 points --> n
1241 if (points2D.size() < 3) { //continue;
1245 // Sort points2D according to polar angles
1246 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1247 for(unsigned int j=0; j<angles.size(); j++) {
1248 points2D[j] = angles[j].first;
1250 // DDV(points2D, points2D.size());
1252 /* When vessels are far away, we try to replace the line segment
1253 with a curved line that join the two vessels but stay
1254 approximately at the same distance from the trachea centroids
1258 - let a and c be two successive vessels centroids
1259 - id distance(a,c) < threshold, next point
1263 - compute mid position between two successive points -
1264 compute dist to trachea centroid for the 3 pts - if middle too
1267 std::vector<MaskSlicePointType> toadd;
1268 unsigned int index = 0;
1270 while (index<points2D.size()-1) {
1271 MaskSlicePointType a = points2D[index];
1272 MaskSlicePointType c = points2D[index+1];
1274 double dac = a.EuclideanDistanceTo(c);
1277 MaskSlicePointType b;
1278 b[0] = a[0]+(c[0]-a[0])/2.0;
1279 b[1] = a[1]+(c[1]-a[1])/2.0;
1281 // Compute distance to trachea centroid
1282 MaskSlicePointType m = centroids_trachea[1];
1283 double da = m.EuclideanDistanceTo(a);
1284 double db = m.EuclideanDistanceTo(b);
1285 //double dc = m.EuclideanDistanceTo(c);
1287 // Mean distance, find point on the line from trachea centroid
1289 double alpha = (da+db)/2.0;
1290 MaskSlicePointType v;
1291 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1292 v[0] = (b[0]-m[0])/n;
1293 v[1] = (b[1]-m[1])/n;
1294 MaskSlicePointType r;
1295 r[0] = m[0]+alpha*v[0];
1296 r[1] = m[1]+alpha*v[1];
1297 points2D.insert(points2D.begin()+index+1, r);
1303 // DDV(points2D, points2D.size());
1305 // Add some points to close the contour
1306 // - H line towards Right
1307 MaskSlicePointType p = points2D[0];
1308 p[0] = bounds[0][0];
1309 points2D.insert(points2D.begin(), p);
1310 // - corner Right/Post
1312 points2D.insert(points2D.begin(), p);
1313 // - H line towards Left
1314 p = points2D.back();
1315 p[0] = bounds[2][0];
1316 points2D.push_back(p);
1317 // - corner Right/Post
1319 points2D.push_back(p);
1320 // Close contour with the first point
1321 points2D.push_back(points2D[0]);
1322 // DDV(points2D, points2D.size());
1324 // Build 3D points from the 2D points
1325 std::vector<ImagePointType> points3D;
1326 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1327 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1329 // Build the mesh from the contour's points
1330 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1331 append->AddInput(mesh);
1332 // if (i ==n-1) { // last slice
1333 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1334 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1335 // append->AddInput(mesh);
1339 // DEBUG: write points3D
1340 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1342 // Build the final 3D mesh form the list 2D mesh
1344 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1346 // Debug, write the mesh
1348 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1350 w->SetFileName("bidon.vtk");
1354 // Compute a single binary 3D image from the list of contours
1355 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1356 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1357 filter->SetMesh(mesh);
1358 filter->SetLikeImage(support);
1360 binarizedContour = filter->GetOutput();
1363 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1364 GetForegroundValue(), 2, true, p);
1366 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1367 GetForegroundValue(), 2, false, p);
1368 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1369 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1370 false, GetBackgroundValue());
1373 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1374 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1376 return binarizedContour;
1379 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1380 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1381 p[1] += 50; // 50 mm Post from this point must be kept
1382 ImageIndexType index;
1383 binarizedContour->TransformPhysicalPointToIndex(p, index);
1384 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1386 // remove from support
1387 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1388 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1389 boolFilter->InPlaceOn();
1390 boolFilter->SetInput1(m_Working_Support);
1391 boolFilter->SetInput2(binarizedContour);
1392 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1393 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1395 boolFilter->SetOperationType(BoolFilterType::And);
1397 boolFilter->SetOperationType(BoolFilterType::AndNot);
1398 boolFilter->Update();
1399 m_Working_Support = boolFilter->GetOutput();
1403 //StopCurrentStep<MaskImageType>(m_Working_Support);
1404 //m_ListOfStations["2R"] = m_Working_Support;
1405 //m_ListOfStations["2L"] = m_Working_Support;
1407 //--------------------------------------------------------------------
1411 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX