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 // Clean some computer landmarks to force the recomputation
89 // FIXME -> to put elsewhere ?
90 this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
92 // Must I compute the supports ?
93 bool supportsExist = true;
94 if (!GetForceSupportsFlag()) {
96 m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
97 m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
98 m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
99 m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
100 m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
101 m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
103 m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
104 m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
105 m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
106 m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
107 m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
108 m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
109 m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
110 m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
111 m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
112 } catch(clitk::ExceptionObject o) {
113 supportsExist = false;
117 if (!supportsExist || GetForceSupportsFlag()) {
118 this->StartNewStep("Supports for stations");
119 this->StartSubStep();
121 // FIXME : why should I remove theses tags ???
122 this->GetAFDB()->RemoveTag("CarinaZ");
123 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
124 this->GetAFDB()->RemoveTag("ApexOfTheChest");
125 this->GetAFDB()->RemoveTag("RightBronchus");
126 this->GetAFDB()->RemoveTag("LeftBronchus");
127 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
128 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
129 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
130 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
131 ExtractStationSupports();
136 ExtractStation_1RL();
137 ExtractStation_2RL();
145 // ---------- todo -----------------------
148 // ExtractStation_8();
151 //this->StartNewStep("Station 7");
152 //this->StartSubStep();
153 //ExtractStation_7();
154 //this->StopSubStep();
157 //--------------------------------------------------------------------
160 //--------------------------------------------------------------------
161 template <class TImageType>
163 clitk::ExtractLymphStationsFilter<TImageType>::
164 GenerateInputRequestedRegion() {
165 //DD("GenerateInputRequestedRegion (nothing?)");
167 //--------------------------------------------------------------------
170 //--------------------------------------------------------------------
171 template <class TImageType>
173 clitk::ExtractLymphStationsFilter<TImageType>::
175 // Final Step -> graft output (if SetNthOutput => redo)
176 // this->GraftOutput(m_ListOfStations["8"]);
178 //--------------------------------------------------------------------
181 //--------------------------------------------------------------------
182 template <class TImageType>
184 clitk::ExtractLymphStationsFilter<TImageType>::
185 CheckForStation(std::string station)
187 // Compute Station name
188 std::string s = "Station"+station;
191 // Define the starting support
192 // if (GetComputeStation(station)) {
193 // std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
195 if (GetComputeStation(station)) {
196 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
201 // std::cout << "Station " << station << " found. I used it" << std::endl;
205 // Check if station already exist in DB
207 // FIXME -> do nothing if not on the command line. Is it what I want ?
208 //bool found = false;
209 if (this->GetAFDB()->TagExist(s)) {
210 m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
215 //--------------------------------------------------------------------
218 //--------------------------------------------------------------------
219 template <class TImageType>
221 clitk::ExtractLymphStationsFilter<TImageType>::
222 GetComputeStation(std::string station)
224 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
226 //--------------------------------------------------------------------
229 //--------------------------------------------------------------------
230 template <class TImageType>
232 clitk::ExtractLymphStationsFilter<TImageType>::
233 AddComputeStation(std::string station)
235 m_ComputeStationMap[station] = true;
237 //--------------------------------------------------------------------
241 //--------------------------------------------------------------------
242 template<class PointType>
243 class comparePointsX {
245 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
247 //--------------------------------------------------------------------
250 //--------------------------------------------------------------------
251 template<class PairType>
252 class comparePointsWithAngle {
254 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
256 //--------------------------------------------------------------------
259 //--------------------------------------------------------------------
261 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
262 std::vector<itk::Point<double, Dim-1> > previous;
263 HypercubeCorners<Dim-1>(previous);
264 out.resize(previous.size()*2);
265 for(unsigned int i=0; i<out.size(); i++) {
266 itk::Point<double, Dim> p;
267 if (i<previous.size()) p[0] = 0;
269 for(int j=0; j<Dim-1; j++)
271 p[j+1] = previous[i%previous.size()][j];
276 //--------------------------------------------------------------------
279 //--------------------------------------------------------------------
281 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
286 //--------------------------------------------------------------------
289 //--------------------------------------------------------------------
290 template<class ImageType>
291 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
292 std::vector<typename ImageType::PointType> & bounds)
294 // Get image max/min coordinates
295 const unsigned int dim=ImageType::ImageDimension;
296 typedef typename ImageType::PointType PointType;
297 typedef typename ImageType::IndexType IndexType;
298 PointType min_c, max_c;
299 IndexType min_i, max_i;
300 min_i = image->GetLargestPossibleRegion().GetIndex();
301 for(unsigned int i=0; i<dim; i++)
302 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
303 image->TransformIndexToPhysicalPoint(min_i, min_c);
304 image->TransformIndexToPhysicalPoint(max_i, max_c);
306 // Get corners coordinates
307 HypercubeCorners<ImageType::ImageDimension>(bounds);
308 for(unsigned int i=0; i<bounds.size(); i++) {
309 for(unsigned int j=0; j<dim; j++) {
310 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
311 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
315 //--------------------------------------------------------------------
318 //--------------------------------------------------------------------
319 template <class ImageType>
321 clitk::ExtractLymphStationsFilter<ImageType>::
322 Remove_Structures(std::string station, std::string s)
325 this->StartNewStep("[Station"+station+"] Remove "+s);
326 MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
327 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
329 catch(clitk::ExceptionObject e) {
330 std::cout << s << " not found, skip." << std::endl;
333 //--------------------------------------------------------------------
336 //--------------------------------------------------------------------
337 template <class TImageType>
339 clitk::ExtractLymphStationsFilter<TImageType>::
340 SetFuzzyThreshold(std::string station, std::string tag, double value)
342 m_FuzzyThreshold[station][tag] = value;
344 //--------------------------------------------------------------------
347 //--------------------------------------------------------------------
348 template <class TImageType>
350 clitk::ExtractLymphStationsFilter<TImageType>::
351 SetThreshold(std::string station, std::string tag, double value)
353 m_Threshold[station][tag] = value;
355 //--------------------------------------------------------------------
358 //--------------------------------------------------------------------
359 template <class TImageType>
361 clitk::ExtractLymphStationsFilter<TImageType>::
362 GetFuzzyThreshold(std::string station, std::string tag)
364 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
365 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
369 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
370 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
374 return m_FuzzyThreshold[station][tag];
376 //--------------------------------------------------------------------
379 //--------------------------------------------------------------------
380 template <class TImageType>
382 clitk::ExtractLymphStationsFilter<TImageType>::
383 GetThreshold(std::string station, std::string tag)
385 if (m_Threshold.find(station) == m_Threshold.end()) {
386 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
390 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
391 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
395 return m_Threshold[station][tag];
397 //--------------------------------------------------------------------
400 //--------------------------------------------------------------------
401 template <class TImageType>
403 clitk::ExtractLymphStationsFilter<TImageType>::
404 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
406 // Create line from A to B with
407 // A = upper border of LLL at left
408 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
411 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
412 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
414 catch(clitk::ExceptionObject & o) {
416 DD("FindLineForS7S8Separation");
417 // Load LeftLowerLobeBronchus and get centroid point
418 MaskImagePointer LeftLowerLobeBronchus =
419 this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
420 std::vector<MaskImagePointType> c;
421 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
424 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
425 MaskImagePointer RightMiddleLobeBronchus =
426 this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
427 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
428 GetBackgroundValue(),
431 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
434 // Insert into the DB
435 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
436 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
439 //--------------------------------------------------------------------
442 //--------------------------------------------------------------------
443 template <class TImageType>
445 clitk::ExtractLymphStationsFilter<TImageType>::
450 z = this->GetAFDB()->GetDouble("CarinaZ");
452 catch(clitk::ExceptionObject e) {
453 //DD("FindCarinaSlicePosition");
455 MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
457 // Get Centroid and Z value
458 std::vector<MaskImagePointType> centroids;
459 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
461 // We dont need Carina structure from now
462 this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
464 // Put inside the AFDB
465 this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
466 this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
472 //--------------------------------------------------------------------
475 //--------------------------------------------------------------------
476 template <class TImageType>
478 clitk::ExtractLymphStationsFilter<TImageType>::
483 z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
485 catch(clitk::ExceptionObject e) {
488 //DD("FindApexOfTheChestPosition");
489 MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
490 MaskImagePointType p;
491 p[0] = p[1] = p[2] = 0.0; // to avoid warning
492 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
494 // We dont need Lungs structure from now
495 this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
497 // Put inside the AFDB
498 this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
499 p[2] -= 5; // We consider 5 mm lower
500 this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
505 // the superior border becomes the more inferior of the two apices
506 MaskImagePointer RightLung = this->GetAFDB()->template GetImage<MaskImageType>("RightLung");
507 MaskImagePointer LeftLung = this->GetAFDB()->template GetImage<MaskImageType>("LeftLung");
508 MaskImagePointType pr;
509 MaskImagePointType pl;
510 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(RightLung, GetBackgroundValue(), 2, false, pr);
511 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(LeftLung, GetBackgroundValue(), 2, false, pl);
512 // We dont need Lungs structure from now
513 this->GetAFDB()->template ReleaseImage<MaskImageType>("RightLung");
514 this->GetAFDB()->template ReleaseImage<MaskImageType>("LeftLung");
515 // Put inside the AFDB
516 if (pr[2] < pl[2]) z = pr[2];
518 this->GetAFDB()->SetDouble("ApexOfTheChestZ", z);
523 //--------------------------------------------------------------------
526 //--------------------------------------------------------------------
527 template <class TImageType>
529 clitk::ExtractLymphStationsFilter<TImageType>::
530 FindLeftAndRightBronchi()
533 m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
534 m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
536 catch(clitk::ExceptionObject & o) {
538 DD("FindLeftAndRightBronchi");
539 // The goal is to separate the trachea inferiorly to the carina into
540 // a Left and Right bronchus.
543 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
545 // Get the Carina position
546 double m_CarinaZ = FindCarina();
548 // Consider only inferiorly to the Carina
549 MaskImagePointer m_Working_Trachea =
550 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
551 GetBackgroundValue());
553 // Labelize the trachea
554 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
556 // Carina position must at the first slice that separate the two
557 // main bronchus (not superiorly). We thus first check that the
558 // upper slice is composed of at least two labels
559 MaskImagePointer RightBronchus;
560 MaskImagePointer LeftBronchus;
561 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
562 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
563 iter.SetFirstDirection(0);
564 iter.SetSecondDirection(1);
565 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
567 while (!iter.IsAtReverseEndOfSlice()) {
568 while (!iter.IsAtReverseEndOfLine()) {
569 if (iter.Get() > maxLabel) maxLabel = iter.Get();
575 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
578 // Compute 3D centroids of both parts to identify the left from the
580 std::vector<ImagePointType> c;
581 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
582 ImagePointType C1 = c[1];
583 ImagePointType C2 = c[2];
585 ImagePixelType rightLabel;
586 ImagePixelType leftLabel;
587 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
588 else { rightLabel = 2; leftLabel = 1; }
590 // Select LeftLabel (set one label to Backgroundvalue)
592 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
593 GetBackgroundValue(), GetForegroundValue());
595 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
596 leftLabel, GetBackgroundValue(), false);
598 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
599 GetBackgroundValue(), GetForegroundValue());
601 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
602 rightLabel, GetBackgroundValue(), false);
606 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
607 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
609 // Insert int AFDB if need after
610 this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
611 RightBronchus, true);
612 this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
616 //--------------------------------------------------------------------
619 //--------------------------------------------------------------------
620 template <class TImageType>
622 clitk::ExtractLymphStationsFilter<TImageType>::
623 FindSuperiorBorderOfAorticArch()
627 z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
629 catch(clitk::ExceptionObject e) {
630 // DD("FindSuperiorBorderOfAorticArch");
631 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
632 MaskImagePointType p;
633 p[0] = p[1] = p[2] = 0.0; // to avoid warning
634 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
635 p[2] += Aorta->GetSpacing()[2]; // the slice above
637 // We dont need Lungs structure from now
638 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
640 // Put inside the AFDB
641 this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
642 this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
648 //--------------------------------------------------------------------
651 //--------------------------------------------------------------------
652 template <class TImageType>
654 clitk::ExtractLymphStationsFilter<TImageType>::
655 FindInferiorBorderOfAorticArch()
659 z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
661 catch(clitk::ExceptionObject e) {
662 //DD("FindInferiorBorderOfAorticArch");
663 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
664 std::vector<MaskSlicePointer> slices;
665 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
667 uint i = slices.size()-1;
669 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
670 std::vector<typename MaskSliceType::PointType> c;
671 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
679 MaskImageIndexType index;
680 index[0] = index[1] = 0.0;
682 MaskImagePointType lower;
683 Aorta->TransformIndexToPhysicalPoint(index, lower);
685 // We dont need Lungs structure from now
686 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
688 // Put inside the AFDB
689 this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
690 this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
696 //--------------------------------------------------------------------
699 //--------------------------------------------------------------------
700 template <class ImageType>
701 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
702 clitk::ExtractLymphStationsFilter<ImageType>::
703 FindAntPostVesselsOLD()
705 // -----------------------------------------------------
706 /* Rod says: "The anterior border, as with the Atlas – UM, is
707 posterior to the vessels (right subclavian vein, left
708 brachiocephalic vein, right brachiocephalic vein, left subclavian
709 artery, left common carotid artery and brachiocephalic trunk).
710 These vessels are not included in the nodal station. The anterior
711 border is drawn to the midpoint of the vessel and an imaginary
712 line joins the middle of these vessels. Between the vessels,
713 station 2 is in contact with station 3a." */
715 // Check if no already done
717 MaskImagePointer binarizedContour;
719 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
721 catch(clitk::ExceptionObject e) {
725 return binarizedContour;
728 /* Here, we consider the vessels form a kind of anterior barrier. We
729 link all vessels centroids and remove what is post to it.
730 - select the list of structure
731 vessel1 = BrachioCephalicArtery
732 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
733 vessel3 = CommonCarotidArtery
734 vessel4 = SubclavianArtery
737 - crop images as needed
738 - slice by slice, choose the CCL and extract centroids
739 - slice by slice, sort according to polar angle wrt Trachea centroid.
740 - slice by slice, link centroids and close contour
741 - remove outside this contour
746 MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
747 MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
748 MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
749 MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
750 MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
751 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
752 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
754 // Create a temporay support
755 // From first slice of BrachioCephalicVein to end of 3A
756 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
757 MaskImagePointType p;
758 p[0] = p[1] = p[2] = 0.0; // to avoid warning
759 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
761 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
762 GetBackgroundValue(), 2, false, p);
764 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
765 false, GetBackgroundValue());
767 // Resize all structures like support
768 BrachioCephalicArtery =
769 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
770 CommonCarotidArtery =
771 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
773 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
775 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
777 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
778 BrachioCephalicVein =
779 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
781 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
784 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
785 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
786 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
787 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
788 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
789 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
790 std::vector<MaskSlicePointer> slices_SubclavianArtery;
791 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
792 std::vector<MaskSlicePointer> slices_Thyroid;
793 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
794 std::vector<MaskSlicePointer> slices_Aorta;
795 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
796 std::vector<MaskSlicePointer> slices_Trachea;
797 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
798 unsigned int n= slices_BrachioCephalicArtery.size();
800 // Get the boundaries of one slice
801 std::vector<MaskSlicePointType> bounds;
802 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
804 // For all slices, for all structures, find the centroid and build the contour
805 // List of 3D points (for debug)
806 std::vector<MaskImagePointType> p3D;
808 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
809 for(unsigned int i=0; i<n; i++) {
810 // Labelize the slices
811 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
812 GetBackgroundValue(), true, 1);
813 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
814 GetBackgroundValue(), true, 1);
815 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
816 GetBackgroundValue(), true, 1);
817 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
818 GetBackgroundValue(), true, 1);
819 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
820 GetBackgroundValue(), true, 1);
821 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
822 GetBackgroundValue(), true, 1);
825 std::vector<MaskSlicePointType> points2D;
826 std::vector<MaskSlicePointType> centroids1;
827 std::vector<MaskSlicePointType> centroids2;
828 std::vector<MaskSlicePointType> centroids3;
829 std::vector<MaskSlicePointType> centroids4;
830 std::vector<MaskSlicePointType> centroids5;
831 std::vector<MaskSlicePointType> centroids6;
832 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
833 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
834 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
835 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
836 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
837 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
839 // BrachioCephalicVein -> when it is separated into two CCL, we
840 // only consider the most at Right one
841 if (centroids6.size() > 2) {
842 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
843 else centroids6.erase(centroids6.begin()+1);
846 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
847 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
848 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
852 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
853 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
854 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
855 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
856 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
857 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
859 // Sort by angle according to trachea centroid and vertical line,
860 // in polar coordinates :
861 // http://en.wikipedia.org/wiki/Polar_coordinate_system
862 std::vector<MaskSlicePointType> centroids_trachea;
863 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
864 typedef std::pair<MaskSlicePointType, double> PointAngleType;
865 std::vector<PointAngleType> angles;
866 for(unsigned int j=0; j<points2D.size(); j++) {
867 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
868 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
869 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
871 if (x>0) angle = atan(y/x);
872 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
873 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
875 if (y>0) angle = M_PI/2.0;
876 if (y<0) angle = -M_PI/2.0;
879 angle = clitk::rad2deg(angle);
880 // Angle is [-180;180] wrt the X axis. We change the X axis to
881 // be the vertical line, because we want to sort from Right to
882 // Left from Post to Ant.
887 if (angle<0) angle = 360-angle;
890 a.first = points2D[j];
895 // Do nothing if less than 2 points --> n
896 if (points2D.size() < 3) { //continue;
900 // Sort points2D according to polar angles
901 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
902 for(unsigned int j=0; j<angles.size(); j++) {
903 points2D[j] = angles[j].first;
905 // DDV(points2D, points2D.size());
907 /* When vessels are far away, we try to replace the line segment
908 with a curved line that join the two vessels but stay
909 approximately at the same distance from the trachea centroids
913 - let a and c be two successive vessels centroids
914 - id distance(a,c) < threshold, next point
918 - compute mid position between two successive points -
919 compute dist to trachea centroid for the 3 pts - if middle too
922 std::vector<MaskSlicePointType> toadd;
923 unsigned int index = 0;
925 while (index<points2D.size()-1) {
926 MaskSlicePointType a = points2D[index];
927 MaskSlicePointType c = points2D[index+1];
929 double dac = a.EuclideanDistanceTo(c);
932 MaskSlicePointType b;
933 b[0] = a[0]+(c[0]-a[0])/2.0;
934 b[1] = a[1]+(c[1]-a[1])/2.0;
936 // Compute distance to trachea centroid
937 MaskSlicePointType m = centroids_trachea[1];
938 double da = m.EuclideanDistanceTo(a);
939 double db = m.EuclideanDistanceTo(b);
940 //double dc = m.EuclideanDistanceTo(c);
942 // Mean distance, find point on the line from trachea centroid
944 double alpha = (da+db)/2.0;
945 MaskSlicePointType v;
946 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
947 v[0] = (b[0]-m[0])/n;
948 v[1] = (b[1]-m[1])/n;
949 MaskSlicePointType r;
950 r[0] = m[0]+alpha*v[0];
951 r[1] = m[1]+alpha*v[1];
952 points2D.insert(points2D.begin()+index+1, r);
958 // DDV(points2D, points2D.size());
960 // Add some points to close the contour
961 // - H line towards Right
962 MaskSlicePointType p = points2D[0];
964 points2D.insert(points2D.begin(), p);
965 // - corner Right/Post
967 points2D.insert(points2D.begin(), p);
968 // - H line towards Left
971 points2D.push_back(p);
972 // - corner Right/Post
974 points2D.push_back(p);
975 // Close contour with the first point
976 points2D.push_back(points2D[0]);
977 // DDV(points2D, points2D.size());
979 // Build 3D points from the 2D points
980 std::vector<ImagePointType> points3D;
981 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
982 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
984 // Build the mesh from the contour's points
985 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
986 append->AddInput(mesh);
989 // DEBUG: write points3D
990 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
992 // Build the final 3D mesh form the list 2D mesh
994 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
996 // Debug, write the mesh
998 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1000 w->SetFileName("bidon.vtk");
1004 // Compute a single binary 3D image from the list of contours
1005 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1006 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1007 filter->SetMesh(mesh);
1008 filter->SetLikeImage(support);
1010 binarizedContour = filter->GetOutput();
1013 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1016 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1019 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1020 false, GetBackgroundValue());
1022 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1023 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1025 return binarizedContour;
1028 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1029 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1030 p[1] += 50; // 50 mm Post from this point must be kept
1031 ImageIndexType index;
1032 binarizedContour->TransformPhysicalPointToIndex(p, index);
1033 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1035 // remove from support
1036 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1037 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1038 boolFilter->InPlaceOn();
1039 boolFilter->SetInput1(m_Working_Support);
1040 boolFilter->SetInput2(binarizedContour);
1041 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1042 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1044 boolFilter->SetOperationType(BoolFilterType::And);
1046 boolFilter->SetOperationType(BoolFilterType::AndNot);
1047 boolFilter->Update();
1048 m_Working_Support = boolFilter->GetOutput();
1052 //StopCurrentStep<MaskImageType>(m_Working_Support);
1053 //m_ListOfStations["2R"] = m_Working_Support;
1054 //m_ListOfStations["2L"] = m_Working_Support;
1056 //--------------------------------------------------------------------
1059 //--------------------------------------------------------------------
1060 template <class ImageType>
1061 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1062 clitk::ExtractLymphStationsFilter<ImageType>::
1063 FindAntPostVessels2()
1065 // -----------------------------------------------------
1066 /* Rod says: "The anterior border, as with the Atlas – UM, is
1067 posterior to the vessels (right subclavian vein, left
1068 brachiocephalic vein, right brachiocephalic vein, left subclavian
1069 artery, left common carotid artery and brachiocephalic trunk).
1070 These vessels are not included in the nodal station. The anterior
1071 border is drawn to the midpoint of the vessel and an imaginary
1072 line joins the middle of these vessels. Between the vessels,
1073 station 2 is in contact with station 3a." */
1075 // Check if no already done
1077 MaskImagePointer binarizedContour;
1079 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1081 catch(clitk::ExceptionObject e) {
1085 return binarizedContour;
1088 /* Here, we consider the vessels form a kind of anterior barrier. We
1089 link all vessels centroids and remove what is post to it.
1090 - select the list of structure
1091 vessel1 = BrachioCephalicArtery
1092 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1093 vessel3 = CommonCarotidArtery
1094 vessel4 = SubclavianArtery
1097 - crop images as needed
1098 - slice by slice, choose the CCL and extract centroids
1099 - slice by slice, sort according to polar angle wrt Trachea centroid.
1100 - slice by slice, link centroids and close contour
1101 - remove outside this contour
1102 - merge with support
1106 std::map<std::string, MaskImagePointer> MapOfStructures;
1107 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1108 MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1109 MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1110 MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftCommonCarotidArtery");
1111 MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightCommonCarotidArtery");
1112 MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftSubclavianArtery");
1113 MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightSubclavianArtery");
1114 MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1115 MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1116 MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1118 std::vector<std::string> ListOfStructuresNames;
1120 // Create a temporay support
1121 // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1122 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1123 MaskImagePointType p;
1124 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1125 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1126 GetBackgroundValue(), 2, true, p);
1128 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1129 GetBackgroundValue(), 2, false, p);
1130 MaskImagePointType p2;
1131 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1132 GetBackgroundValue(), 2, false, p2);
1133 if (p2[2] > p[2]) p = p2;
1135 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1136 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1137 false, GetBackgroundValue());
1139 // Resize all structures like support
1140 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1141 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1145 typedef std::vector<MaskSlicePointer> SlicesType;
1146 std::map<std::string, SlicesType> MapOfSlices;
1147 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1149 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1150 MapOfSlices[it->first] = s;
1153 unsigned int n= MapOfSlices["Trachea"].size();
1155 // Get the boundaries of one slice
1156 std::vector<MaskSlicePointType> bounds;
1157 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1160 // For all slices, for all structures, find the centroid and build the contour
1161 // List of 3D points (for debug)
1162 std::vector<MaskImagePointType> p3D;
1163 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1164 for(unsigned int i=0; i<n; i++) {
1166 // Labelize the slices
1167 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1168 MaskSlicePointer & s = MapOfSlices[it->first][i];
1169 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1173 std::vector<MaskSlicePointType> points2D;
1174 typedef std::vector<MaskSlicePointType> CentroidsType;
1175 std::map<std::string, CentroidsType> MapOfCentroids;
1176 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1177 std::string structure = it->first;
1178 MaskSlicePointer & s = MapOfSlices[structure][i];
1180 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1181 MapOfCentroids[structure] = c;
1185 // BrachioCephalicVein -> when it is separated into two CCL, we
1186 // only consider the most at Right one
1187 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1189 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1190 else c.erase(c.begin()+1);
1194 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1195 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1196 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1197 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1198 MapOfCentroids["BrachioCephalicVein"].clear();
1202 // Add all 2D points
1203 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1204 std::string structure = it->first;
1205 if (structure != "Trachea") { // do not add centroid of trachea
1206 CentroidsType & c = MapOfCentroids[structure];
1207 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1211 // Sort by angle according to trachea centroid and vertical line,
1212 // in polar coordinates :
1213 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1214 // std::vector<MaskSlicePointType> centroids_trachea;
1215 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1216 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1217 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1218 std::vector<PointAngleType> angles;
1219 for(unsigned int j=0; j<points2D.size(); j++) {
1220 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1221 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1222 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1224 if (x>0) angle = atan(y/x);
1225 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1226 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1228 if (y>0) angle = M_PI/2.0;
1229 if (y<0) angle = -M_PI/2.0;
1230 if (y==0) angle = 0;
1232 angle = clitk::rad2deg(angle);
1233 // Angle is [-180;180] wrt the X axis. We change the X axis to
1234 // be the vertical line, because we want to sort from Right to
1235 // Left from Post to Ant.
1237 angle = (270-angle);
1240 if (angle<0) angle = 360-angle;
1243 a.first = points2D[j];
1245 angles.push_back(a);
1248 // Do nothing if less than 2 points --> n
1249 if (points2D.size() < 3) { //continue;
1253 // Sort points2D according to polar angles
1254 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1255 for(unsigned int j=0; j<angles.size(); j++) {
1256 points2D[j] = angles[j].first;
1258 // DDV(points2D, points2D.size());
1260 /* When vessels are far away, we try to replace the line segment
1261 with a curved line that join the two vessels but stay
1262 approximately at the same distance from the trachea centroids
1266 - let a and c be two successive vessels centroids
1267 - id distance(a,c) < threshold, next point
1271 - compute mid position between two successive points -
1272 compute dist to trachea centroid for the 3 pts - if middle too
1275 std::vector<MaskSlicePointType> toadd;
1276 unsigned int index = 0;
1278 while (index<points2D.size()-1) {
1279 MaskSlicePointType a = points2D[index];
1280 MaskSlicePointType c = points2D[index+1];
1282 double dac = a.EuclideanDistanceTo(c);
1285 MaskSlicePointType b;
1286 b[0] = a[0]+(c[0]-a[0])/2.0;
1287 b[1] = a[1]+(c[1]-a[1])/2.0;
1289 // Compute distance to trachea centroid
1290 MaskSlicePointType m = centroids_trachea[1];
1291 double da = m.EuclideanDistanceTo(a);
1292 double db = m.EuclideanDistanceTo(b);
1293 //double dc = m.EuclideanDistanceTo(c);
1295 // Mean distance, find point on the line from trachea centroid
1297 double alpha = (da+db)/2.0;
1298 MaskSlicePointType v;
1299 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1300 v[0] = (b[0]-m[0])/n;
1301 v[1] = (b[1]-m[1])/n;
1302 MaskSlicePointType r;
1303 r[0] = m[0]+alpha*v[0];
1304 r[1] = m[1]+alpha*v[1];
1305 points2D.insert(points2D.begin()+index+1, r);
1311 // DDV(points2D, points2D.size());
1313 // Add some points to close the contour
1314 // - H line towards Right
1315 MaskSlicePointType p = points2D[0];
1316 p[0] = bounds[0][0];
1317 points2D.insert(points2D.begin(), p);
1318 // - corner Right/Post
1320 points2D.insert(points2D.begin(), p);
1321 // - H line towards Left
1322 p = points2D.back();
1323 p[0] = bounds[2][0];
1324 points2D.push_back(p);
1325 // - corner Right/Post
1327 points2D.push_back(p);
1328 // Close contour with the first point
1329 points2D.push_back(points2D[0]);
1330 // DDV(points2D, points2D.size());
1332 // Build 3D points from the 2D points
1333 std::vector<ImagePointType> points3D;
1334 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1335 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1337 // Build the mesh from the contour's points
1338 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1339 append->AddInput(mesh);
1340 // if (i ==n-1) { // last slice
1341 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1342 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1343 // append->AddInput(mesh);
1347 // DEBUG: write points3D
1348 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1350 // Build the final 3D mesh form the list 2D mesh
1352 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1354 // Debug, write the mesh
1356 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1358 w->SetFileName("bidon.vtk");
1362 // Compute a single binary 3D image from the list of contours
1363 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1364 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1365 filter->SetMesh(mesh);
1366 filter->SetLikeImage(support);
1368 binarizedContour = filter->GetOutput();
1371 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1372 GetForegroundValue(), 2, true, p);
1374 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1375 GetForegroundValue(), 2, false, p);
1376 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1377 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1378 false, GetBackgroundValue());
1381 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1382 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1384 return binarizedContour;
1387 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1388 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1389 p[1] += 50; // 50 mm Post from this point must be kept
1390 ImageIndexType index;
1391 binarizedContour->TransformPhysicalPointToIndex(p, index);
1392 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1394 // remove from support
1395 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1396 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1397 boolFilter->InPlaceOn();
1398 boolFilter->SetInput1(m_Working_Support);
1399 boolFilter->SetInput2(binarizedContour);
1400 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1401 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1403 boolFilter->SetOperationType(BoolFilterType::And);
1405 boolFilter->SetOperationType(BoolFilterType::AndNot);
1406 boolFilter->Update();
1407 m_Working_Support = boolFilter->GetOutput();
1411 //StopCurrentStep<MaskImageType>(m_Working_Support);
1412 //m_ListOfStations["2R"] = m_Working_Support;
1413 //m_ListOfStations["2L"] = m_Working_Support;
1415 //--------------------------------------------------------------------
1417 //--------------------------------------------------------------------
1418 template <class ImageType>
1420 clitk::ExtractLymphStationsFilter<ImageType>::
1421 WriteImageSupport(std::string support)
1423 writeImage<MaskImageType>(m_ListOfSupports[support], this->GetAFDBPath()+"/"+"seg/Support_"+support+".mha");
1424 this->GetAFDB()->SetImageFilename("Support_"+support, "seg/Support_"+support+".mha");
1426 //--------------------------------------------------------------------
1429 //--------------------------------------------------------------------
1430 template <class ImageType>
1432 clitk::ExtractLymphStationsFilter<ImageType>::
1433 WriteImageStation(std::string station)
1435 writeImage<MaskImageType>(m_ListOfStations[station], GetAFDB()->GetPath()+"/seg/Station"+station+".mha");
1436 GetAFDB()->SetImageFilename("Station"+station, "seg/Station"+station+".mha");
1439 //--------------------------------------------------------------------
1442 //--------------------------------------------------------------------
1443 template <class ImageType>
1445 clitk::ExtractLymphStationsFilter<ImageType>::
1446 ComputeOverlapWithRef(std::string station)
1448 if (GetComputeStation(station)) {
1449 MaskImagePointer ref = this->GetAFDB()->template GetImage <MaskImageType>("Station"+station+"_Ref");
1450 typedef clitk::LabelImageOverlapMeasureFilter<MaskImageType> FilterType;
1451 typename FilterType::Pointer filter = FilterType::New();
1452 filter->SetInput(0, m_ListOfStations[station]);
1453 filter->SetInput(1, ref);
1457 //--------------------------------------------------------------------
1459 //--------------------------------------------------------------------
1460 template <class ImageType>
1462 clitk::ExtractLymphStationsFilter<ImageType>::
1463 ReadSupportLimits(std::string filename)
1465 m_ListOfSupportLimits.clear();
1467 openFileForReading(is, filename);
1470 SupportLimitsType s;
1472 if (is) m_ListOfSupportLimits.push_back(s);
1475 //--------------------------------------------------------------------
1477 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX