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();
65 ExtractStation_4RL_SetDefaultValues();
67 ExtractStation_7_SetDefaultValues();
68 ExtractStation_8_SetDefaultValues();
70 //--------------------------------------------------------------------
73 //--------------------------------------------------------------------
74 template <class TImageType>
76 clitk::ExtractLymphStationsFilter<TImageType>::
77 GenerateOutputInformation() {
80 m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
81 m_Mediastinum = this->GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
83 // Clean some computer landmarks to force the recomputation
84 this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
86 // Global supports for stations
87 bool supportsExist = true;
89 m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
90 m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
91 m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
92 m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
93 m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
94 m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
96 m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
97 m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
98 m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
99 m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
100 m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
101 m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
102 m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
103 m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
104 m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
105 } catch(clitk::ExceptionObject o) {
106 supportsExist = false;
109 if (!supportsExist || GetComputeStationsSupportsFlag()) {
110 this->StartNewStep("Supports for stations");
111 this->StartSubStep();
112 this->GetAFDB()->RemoveTag("CarinaZ");
113 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
114 this->GetAFDB()->RemoveTag("ApexOfTheChest");
115 this->GetAFDB()->RemoveTag("RightBronchus");
116 this->GetAFDB()->RemoveTag("LeftBronchus");
117 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
118 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
119 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
120 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
121 ExtractStationSupports();
131 // Extract Station2RL
132 ExtractStation_2RL();
134 // Extract Station1RL
135 ExtractStation_1RL();
137 // Extract Station1RL
138 ExtractStation_4RL();
140 // ---------- TODO -----------------------
146 this->StartNewStep("Station 7");
147 this->StartSubStep();
152 //--------------------------------------------------------------------
155 //--------------------------------------------------------------------
156 template <class TImageType>
158 clitk::ExtractLymphStationsFilter<TImageType>::
159 GenerateInputRequestedRegion() {
160 //DD("GenerateInputRequestedRegion (nothing?)");
162 //--------------------------------------------------------------------
165 //--------------------------------------------------------------------
166 template <class TImageType>
168 clitk::ExtractLymphStationsFilter<TImageType>::
170 // Final Step -> graft output (if SetNthOutput => redo)
171 this->GraftOutput(m_ListOfStations["8"]);
173 //--------------------------------------------------------------------
176 //--------------------------------------------------------------------
177 template <class TImageType>
179 clitk::ExtractLymphStationsFilter<TImageType>::
180 CheckForStation(std::string station)
182 // Compute Station name
183 std::string s = "Station"+station;
186 // Check if station already exist in DB
188 if (this->GetAFDB()->TagExist(s)) {
189 m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
193 // Define the starting support
194 if (found && GetComputeStation(station)) {
195 std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
197 if (!found || GetComputeStation(station)) {
198 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
202 std::cout << "Station " << station << " found. I used it" << std::endl;
206 //--------------------------------------------------------------------
209 //--------------------------------------------------------------------
210 template <class TImageType>
212 clitk::ExtractLymphStationsFilter<TImageType>::
213 GetComputeStation(std::string station)
215 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
217 //--------------------------------------------------------------------
220 //--------------------------------------------------------------------
221 template <class TImageType>
223 clitk::ExtractLymphStationsFilter<TImageType>::
224 AddComputeStation(std::string station)
226 m_ComputeStationMap[station] = true;
228 //--------------------------------------------------------------------
232 //--------------------------------------------------------------------
233 template<class PointType>
234 class comparePointsX {
236 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
238 //--------------------------------------------------------------------
241 //--------------------------------------------------------------------
242 template<class PairType>
243 class comparePointsWithAngle {
245 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
247 //--------------------------------------------------------------------
250 //--------------------------------------------------------------------
252 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
253 std::vector<itk::Point<double, Dim-1> > previous;
254 HypercubeCorners<Dim-1>(previous);
255 out.resize(previous.size()*2);
256 for(unsigned int i=0; i<out.size(); i++) {
257 itk::Point<double, Dim> p;
258 if (i<previous.size()) p[0] = 0;
260 for(int j=0; j<Dim-1; j++)
262 p[j+1] = previous[i%previous.size()][j];
267 //--------------------------------------------------------------------
270 //--------------------------------------------------------------------
272 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
277 //--------------------------------------------------------------------
280 //--------------------------------------------------------------------
281 template<class ImageType>
282 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
283 std::vector<typename ImageType::PointType> & bounds)
285 // Get image max/min coordinates
286 const unsigned int dim=ImageType::ImageDimension;
287 typedef typename ImageType::PointType PointType;
288 typedef typename ImageType::IndexType IndexType;
289 PointType min_c, max_c;
290 IndexType min_i, max_i;
291 min_i = image->GetLargestPossibleRegion().GetIndex();
292 for(unsigned int i=0; i<dim; i++)
293 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
294 image->TransformIndexToPhysicalPoint(min_i, min_c);
295 image->TransformIndexToPhysicalPoint(max_i, max_c);
297 // Get corners coordinates
298 HypercubeCorners<ImageType::ImageDimension>(bounds);
299 for(unsigned int i=0; i<bounds.size(); i++) {
300 for(unsigned int j=0; j<dim; j++) {
301 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
302 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
306 //--------------------------------------------------------------------
309 //--------------------------------------------------------------------
310 template <class ImageType>
312 clitk::ExtractLymphStationsFilter<ImageType>::
313 Remove_Structures(std::string station, std::string s)
316 this->StartNewStep("[Station"+station+"] Remove "+s);
317 MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
318 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
320 catch(clitk::ExceptionObject e) {
321 std::cout << s << " not found, skip." << std::endl;
324 //--------------------------------------------------------------------
327 //--------------------------------------------------------------------
328 template <class TImageType>
330 clitk::ExtractLymphStationsFilter<TImageType>::
331 SetFuzzyThreshold(std::string station, std::string tag, double value)
333 m_FuzzyThreshold[station][tag] = value;
335 //--------------------------------------------------------------------
338 //--------------------------------------------------------------------
339 template <class TImageType>
341 clitk::ExtractLymphStationsFilter<TImageType>::
342 SetThreshold(std::string station, std::string tag, double value)
344 m_Threshold[station][tag] = value;
346 //--------------------------------------------------------------------
349 //--------------------------------------------------------------------
350 template <class TImageType>
352 clitk::ExtractLymphStationsFilter<TImageType>::
353 GetFuzzyThreshold(std::string station, std::string tag)
355 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
356 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
360 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
361 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
365 return m_FuzzyThreshold[station][tag];
367 //--------------------------------------------------------------------
370 //--------------------------------------------------------------------
371 template <class TImageType>
373 clitk::ExtractLymphStationsFilter<TImageType>::
374 GetThreshold(std::string station, std::string tag)
376 if (m_Threshold.find(station) == m_Threshold.end()) {
377 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
381 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
382 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
386 return m_Threshold[station][tag];
388 //--------------------------------------------------------------------
391 //--------------------------------------------------------------------
392 template <class TImageType>
394 clitk::ExtractLymphStationsFilter<TImageType>::
395 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
397 // Create line from A to B with
398 // A = upper border of LLL at left
399 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
402 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
403 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
405 catch(clitk::ExceptionObject & o) {
407 DD("FindLineForS7S8Separation");
408 // Load LeftLowerLobeBronchus and get centroid point
409 MaskImagePointer LeftLowerLobeBronchus =
410 this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
411 std::vector<MaskImagePointType> c;
412 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
415 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
416 MaskImagePointer RightMiddleLobeBronchus =
417 this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
418 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
419 GetBackgroundValue(),
422 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
425 // Insert into the DB
426 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
427 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
430 //--------------------------------------------------------------------
433 //--------------------------------------------------------------------
434 template <class TImageType>
436 clitk::ExtractLymphStationsFilter<TImageType>::
441 z = this->GetAFDB()->GetDouble("CarinaZ");
443 catch(clitk::ExceptionObject e) {
444 DD("FindCarinaSlicePosition");
446 MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
448 // Get Centroid and Z value
449 std::vector<MaskImagePointType> centroids;
450 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
452 // We dont need Carina structure from now
453 this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
455 // Put inside the AFDB
456 this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
457 this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
463 //--------------------------------------------------------------------
466 //--------------------------------------------------------------------
467 template <class TImageType>
469 clitk::ExtractLymphStationsFilter<TImageType>::
474 z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
476 catch(clitk::ExceptionObject e) {
477 DD("FindApexOfTheChestPosition");
478 MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
479 MaskImagePointType p;
480 p[0] = p[1] = p[2] = 0.0; // to avoid warning
481 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
483 // We dont need Lungs structure from now
484 this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
486 // Put inside the AFDB
487 this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
488 p[2] -= 5; // We consider 5 mm lower
489 this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
495 //--------------------------------------------------------------------
498 //--------------------------------------------------------------------
499 template <class TImageType>
501 clitk::ExtractLymphStationsFilter<TImageType>::
502 FindLeftAndRightBronchi()
505 m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
506 m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
508 catch(clitk::ExceptionObject & o) {
510 DD("FindLeftAndRightBronchi");
511 // The goal is to separate the trachea inferiorly to the carina into
512 // a Left and Right bronchus.
515 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
517 // Get the Carina position
518 double m_CarinaZ = FindCarina();
520 // Consider only inferiorly to the Carina
521 MaskImagePointer m_Working_Trachea =
522 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
523 GetBackgroundValue());
525 // Labelize the trachea
526 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
528 // Carina position must at the first slice that separate the two
529 // main bronchus (not superiorly). We thus first check that the
530 // upper slice is composed of at least two labels
531 MaskImagePointer RightBronchus;
532 MaskImagePointer LeftBronchus;
533 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
534 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
535 iter.SetFirstDirection(0);
536 iter.SetSecondDirection(1);
537 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
539 while (!iter.IsAtReverseEndOfSlice()) {
540 while (!iter.IsAtReverseEndOfLine()) {
541 if (iter.Get() > maxLabel) maxLabel = iter.Get();
547 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
550 // Compute 3D centroids of both parts to identify the left from the
552 std::vector<ImagePointType> c;
553 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
554 ImagePointType C1 = c[1];
555 ImagePointType C2 = c[2];
557 ImagePixelType rightLabel;
558 ImagePixelType leftLabel;
559 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
560 else { rightLabel = 2; leftLabel = 1; }
562 // Select LeftLabel (set one label to Backgroundvalue)
564 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
565 GetBackgroundValue(), GetForegroundValue());
567 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
568 leftLabel, GetBackgroundValue(), false);
570 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
571 GetBackgroundValue(), GetForegroundValue());
573 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
574 rightLabel, GetBackgroundValue(), false);
578 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
579 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
581 // Insert int AFDB if need after
582 this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
583 RightBronchus, true);
584 this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
588 //--------------------------------------------------------------------
591 //--------------------------------------------------------------------
592 template <class TImageType>
594 clitk::ExtractLymphStationsFilter<TImageType>::
595 FindSuperiorBorderOfAorticArch()
599 z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
601 catch(clitk::ExceptionObject e) {
602 DD("FindSuperiorBorderOfAorticArch");
603 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
604 MaskImagePointType p;
605 p[0] = p[1] = p[2] = 0.0; // to avoid warning
606 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
607 p[2] += Aorta->GetSpacing()[2]; // the slice above
609 // We dont need Lungs structure from now
610 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
612 // Put inside the AFDB
613 this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
614 this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
620 //--------------------------------------------------------------------
623 //--------------------------------------------------------------------
624 template <class TImageType>
626 clitk::ExtractLymphStationsFilter<TImageType>::
627 FindInferiorBorderOfAorticArch()
631 z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
633 catch(clitk::ExceptionObject e) {
634 DD("FindInferiorBorderOfAorticArch");
635 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
636 std::vector<MaskSlicePointer> slices;
637 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
639 uint i = slices.size()-1;
641 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
642 std::vector<typename MaskSliceType::PointType> c;
643 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
651 MaskImageIndexType index;
652 index[0] = index[1] = 0.0;
654 MaskImagePointType lower;
655 Aorta->TransformIndexToPhysicalPoint(index, lower);
657 // We dont need Lungs structure from now
658 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
660 // Put inside the AFDB
661 this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
662 this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
668 //--------------------------------------------------------------------
671 //--------------------------------------------------------------------
672 template <class ImageType>
673 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
674 clitk::ExtractLymphStationsFilter<ImageType>::
675 FindAntPostVesselsOLD()
677 // -----------------------------------------------------
678 /* Rod says: "The anterior border, as with the Atlas – UM, is
679 posterior to the vessels (right subclavian vein, left
680 brachiocephalic vein, right brachiocephalic vein, left subclavian
681 artery, left common carotid artery and brachiocephalic trunk).
682 These vessels are not included in the nodal station. The anterior
683 border is drawn to the midpoint of the vessel and an imaginary
684 line joins the middle of these vessels. Between the vessels,
685 station 2 is in contact with station 3a." */
687 // Check if no already done
689 MaskImagePointer binarizedContour;
691 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
693 catch(clitk::ExceptionObject e) {
697 return binarizedContour;
700 /* Here, we consider the vessels form a kind of anterior barrier. We
701 link all vessels centroids and remove what is post to it.
702 - select the list of structure
703 vessel1 = BrachioCephalicArtery
704 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
705 vessel3 = CommonCarotidArtery
706 vessel4 = SubclavianArtery
709 - crop images as needed
710 - slice by slice, choose the CCL and extract centroids
711 - slice by slice, sort according to polar angle wrt Trachea centroid.
712 - slice by slice, link centroids and close contour
713 - remove outside this contour
718 MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
719 MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
720 MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
721 MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
722 MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
723 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
724 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
726 // Create a temporay support
727 // From first slice of BrachioCephalicVein to end of 3A
728 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
729 MaskImagePointType p;
730 p[0] = p[1] = p[2] = 0.0; // to avoid warning
731 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
733 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
734 GetBackgroundValue(), 2, false, p);
736 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
737 false, GetBackgroundValue());
739 // Resize all structures like support
740 BrachioCephalicArtery =
741 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
742 CommonCarotidArtery =
743 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
745 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
747 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
749 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
750 BrachioCephalicVein =
751 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
753 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
756 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
757 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
758 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
759 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
760 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
761 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
762 std::vector<MaskSlicePointer> slices_SubclavianArtery;
763 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
764 std::vector<MaskSlicePointer> slices_Thyroid;
765 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
766 std::vector<MaskSlicePointer> slices_Aorta;
767 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
768 std::vector<MaskSlicePointer> slices_Trachea;
769 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
770 unsigned int n= slices_BrachioCephalicArtery.size();
772 // Get the boundaries of one slice
773 std::vector<MaskSlicePointType> bounds;
774 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
776 // For all slices, for all structures, find the centroid and build the contour
777 // List of 3D points (for debug)
778 std::vector<MaskImagePointType> p3D;
780 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
781 for(unsigned int i=0; i<n; i++) {
782 // Labelize the slices
783 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
784 GetBackgroundValue(), true, 1);
785 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
786 GetBackgroundValue(), true, 1);
787 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
788 GetBackgroundValue(), true, 1);
789 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
790 GetBackgroundValue(), true, 1);
791 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
792 GetBackgroundValue(), true, 1);
793 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
794 GetBackgroundValue(), true, 1);
797 std::vector<MaskSlicePointType> points2D;
798 std::vector<MaskSlicePointType> centroids1;
799 std::vector<MaskSlicePointType> centroids2;
800 std::vector<MaskSlicePointType> centroids3;
801 std::vector<MaskSlicePointType> centroids4;
802 std::vector<MaskSlicePointType> centroids5;
803 std::vector<MaskSlicePointType> centroids6;
804 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
805 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
806 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
807 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
808 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
809 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
811 // BrachioCephalicVein -> when it is separated into two CCL, we
812 // only consider the most at Right one
813 if (centroids6.size() > 2) {
814 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
815 else centroids6.erase(centroids6.begin()+1);
818 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
819 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
820 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
824 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
825 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
826 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
827 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
828 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
829 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
831 // Sort by angle according to trachea centroid and vertical line,
832 // in polar coordinates :
833 // http://en.wikipedia.org/wiki/Polar_coordinate_system
834 std::vector<MaskSlicePointType> centroids_trachea;
835 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
836 typedef std::pair<MaskSlicePointType, double> PointAngleType;
837 std::vector<PointAngleType> angles;
838 for(unsigned int j=0; j<points2D.size(); j++) {
839 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
840 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
841 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
843 if (x>0) angle = atan(y/x);
844 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
845 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
847 if (y>0) angle = M_PI/2.0;
848 if (y<0) angle = -M_PI/2.0;
851 angle = clitk::rad2deg(angle);
852 // Angle is [-180;180] wrt the X axis. We change the X axis to
853 // be the vertical line, because we want to sort from Right to
854 // Left from Post to Ant.
859 if (angle<0) angle = 360-angle;
862 a.first = points2D[j];
867 // Do nothing if less than 2 points --> n
868 if (points2D.size() < 3) { //continue;
872 // Sort points2D according to polar angles
873 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
874 for(unsigned int j=0; j<angles.size(); j++) {
875 points2D[j] = angles[j].first;
877 // DDV(points2D, points2D.size());
879 /* When vessels are far away, we try to replace the line segment
880 with a curved line that join the two vessels but stay
881 approximately at the same distance from the trachea centroids
885 - let a and c be two successive vessels centroids
886 - id distance(a,c) < threshold, next point
890 - compute mid position between two successive points -
891 compute dist to trachea centroid for the 3 pts - if middle too
894 std::vector<MaskSlicePointType> toadd;
895 unsigned int index = 0;
897 while (index<points2D.size()-1) {
898 MaskSlicePointType a = points2D[index];
899 MaskSlicePointType c = points2D[index+1];
901 double dac = a.EuclideanDistanceTo(c);
904 MaskSlicePointType b;
905 b[0] = a[0]+(c[0]-a[0])/2.0;
906 b[1] = a[1]+(c[1]-a[1])/2.0;
908 // Compute distance to trachea centroid
909 MaskSlicePointType m = centroids_trachea[1];
910 double da = m.EuclideanDistanceTo(a);
911 double db = m.EuclideanDistanceTo(b);
912 //double dc = m.EuclideanDistanceTo(c);
914 // Mean distance, find point on the line from trachea centroid
916 double alpha = (da+db)/2.0;
917 MaskSlicePointType v;
918 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
919 v[0] = (b[0]-m[0])/n;
920 v[1] = (b[1]-m[1])/n;
921 MaskSlicePointType r;
922 r[0] = m[0]+alpha*v[0];
923 r[1] = m[1]+alpha*v[1];
924 points2D.insert(points2D.begin()+index+1, r);
930 // DDV(points2D, points2D.size());
932 // Add some points to close the contour
933 // - H line towards Right
934 MaskSlicePointType p = points2D[0];
936 points2D.insert(points2D.begin(), p);
937 // - corner Right/Post
939 points2D.insert(points2D.begin(), p);
940 // - H line towards Left
943 points2D.push_back(p);
944 // - corner Right/Post
946 points2D.push_back(p);
947 // Close contour with the first point
948 points2D.push_back(points2D[0]);
949 // DDV(points2D, points2D.size());
951 // Build 3D points from the 2D points
952 std::vector<ImagePointType> points3D;
953 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
954 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
956 // Build the mesh from the contour's points
957 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
958 append->AddInput(mesh);
961 // DEBUG: write points3D
962 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
964 // Build the final 3D mesh form the list 2D mesh
966 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
968 // Debug, write the mesh
970 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
972 w->SetFileName("bidon.vtk");
976 // Compute a single binary 3D image from the list of contours
977 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
978 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
979 filter->SetMesh(mesh);
980 filter->SetLikeImage(support);
982 binarizedContour = filter->GetOutput();
985 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
988 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
991 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
992 false, GetBackgroundValue());
994 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
995 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
997 return binarizedContour;
1000 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1001 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1002 p[1] += 50; // 50 mm Post from this point must be kept
1003 ImageIndexType index;
1004 binarizedContour->TransformPhysicalPointToIndex(p, index);
1005 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1007 // remove from support
1008 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1009 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1010 boolFilter->InPlaceOn();
1011 boolFilter->SetInput1(m_Working_Support);
1012 boolFilter->SetInput2(binarizedContour);
1013 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1014 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1016 boolFilter->SetOperationType(BoolFilterType::And);
1018 boolFilter->SetOperationType(BoolFilterType::AndNot);
1019 boolFilter->Update();
1020 m_Working_Support = boolFilter->GetOutput();
1024 //StopCurrentStep<MaskImageType>(m_Working_Support);
1025 //m_ListOfStations["2R"] = m_Working_Support;
1026 //m_ListOfStations["2L"] = m_Working_Support;
1028 //--------------------------------------------------------------------
1031 //--------------------------------------------------------------------
1032 template <class ImageType>
1033 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1034 clitk::ExtractLymphStationsFilter<ImageType>::
1035 FindAntPostVessels2()
1037 // -----------------------------------------------------
1038 /* Rod says: "The anterior border, as with the Atlas – UM, is
1039 posterior to the vessels (right subclavian vein, left
1040 brachiocephalic vein, right brachiocephalic vein, left subclavian
1041 artery, left common carotid artery and brachiocephalic trunk).
1042 These vessels are not included in the nodal station. The anterior
1043 border is drawn to the midpoint of the vessel and an imaginary
1044 line joins the middle of these vessels. Between the vessels,
1045 station 2 is in contact with station 3a." */
1047 // Check if no already done
1049 MaskImagePointer binarizedContour;
1051 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1053 catch(clitk::ExceptionObject e) {
1057 return binarizedContour;
1060 /* Here, we consider the vessels form a kind of anterior barrier. We
1061 link all vessels centroids and remove what is post to it.
1062 - select the list of structure
1063 vessel1 = BrachioCephalicArtery
1064 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1065 vessel3 = CommonCarotidArtery
1066 vessel4 = SubclavianArtery
1069 - crop images as needed
1070 - slice by slice, choose the CCL and extract centroids
1071 - slice by slice, sort according to polar angle wrt Trachea centroid.
1072 - slice by slice, link centroids and close contour
1073 - remove outside this contour
1074 - merge with support
1078 std::map<std::string, MaskImagePointer> MapOfStructures;
1079 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1080 MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1081 MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1082 MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1083 MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1084 MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1085 MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1086 MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1087 MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1088 MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1090 std::vector<std::string> ListOfStructuresNames;
1092 // Create a temporay support
1093 // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1094 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1095 MaskImagePointType p;
1096 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1097 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1098 GetBackgroundValue(), 2, true, p);
1100 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1101 GetBackgroundValue(), 2, false, p);
1102 MaskImagePointType p2;
1103 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1104 GetBackgroundValue(), 2, false, p2);
1105 if (p2[2] > p[2]) p = p2;
1107 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1108 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1109 false, GetBackgroundValue());
1111 // Resize all structures like support
1112 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1113 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1117 typedef std::vector<MaskSlicePointer> SlicesType;
1118 std::map<std::string, SlicesType> MapOfSlices;
1119 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1121 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1122 MapOfSlices[it->first] = s;
1125 unsigned int n= MapOfSlices["Trachea"].size();
1127 // Get the boundaries of one slice
1128 std::vector<MaskSlicePointType> bounds;
1129 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1132 // For all slices, for all structures, find the centroid and build the contour
1133 // List of 3D points (for debug)
1134 std::vector<MaskImagePointType> p3D;
1135 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1136 for(unsigned int i=0; i<n; i++) {
1138 // Labelize the slices
1139 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1140 MaskSlicePointer & s = MapOfSlices[it->first][i];
1141 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1145 std::vector<MaskSlicePointType> points2D;
1146 typedef std::vector<MaskSlicePointType> CentroidsType;
1147 std::map<std::string, CentroidsType> MapOfCentroids;
1148 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1149 std::string structure = it->first;
1150 MaskSlicePointer & s = MapOfSlices[structure][i];
1152 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1153 MapOfCentroids[structure] = c;
1157 // BrachioCephalicVein -> when it is separated into two CCL, we
1158 // only consider the most at Right one
1159 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1161 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1162 else c.erase(c.begin()+1);
1166 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1167 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1168 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1169 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1170 MapOfCentroids["BrachioCephalicVein"].clear();
1174 // Add all 2D points
1175 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1176 std::string structure = it->first;
1177 if (structure != "Trachea") { // do not add centroid of trachea
1178 CentroidsType & c = MapOfCentroids[structure];
1179 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1183 // Sort by angle according to trachea centroid and vertical line,
1184 // in polar coordinates :
1185 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1186 // std::vector<MaskSlicePointType> centroids_trachea;
1187 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1188 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1189 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1190 std::vector<PointAngleType> angles;
1191 for(unsigned int j=0; j<points2D.size(); j++) {
1192 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1193 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1194 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1196 if (x>0) angle = atan(y/x);
1197 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1198 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1200 if (y>0) angle = M_PI/2.0;
1201 if (y<0) angle = -M_PI/2.0;
1202 if (y==0) angle = 0;
1204 angle = clitk::rad2deg(angle);
1205 // Angle is [-180;180] wrt the X axis. We change the X axis to
1206 // be the vertical line, because we want to sort from Right to
1207 // Left from Post to Ant.
1209 angle = (270-angle);
1212 if (angle<0) angle = 360-angle;
1215 a.first = points2D[j];
1217 angles.push_back(a);
1220 // Do nothing if less than 2 points --> n
1221 if (points2D.size() < 3) { //continue;
1225 // Sort points2D according to polar angles
1226 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1227 for(unsigned int j=0; j<angles.size(); j++) {
1228 points2D[j] = angles[j].first;
1230 // DDV(points2D, points2D.size());
1232 /* When vessels are far away, we try to replace the line segment
1233 with a curved line that join the two vessels but stay
1234 approximately at the same distance from the trachea centroids
1238 - let a and c be two successive vessels centroids
1239 - id distance(a,c) < threshold, next point
1243 - compute mid position between two successive points -
1244 compute dist to trachea centroid for the 3 pts - if middle too
1247 std::vector<MaskSlicePointType> toadd;
1248 unsigned int index = 0;
1250 while (index<points2D.size()-1) {
1251 MaskSlicePointType a = points2D[index];
1252 MaskSlicePointType c = points2D[index+1];
1254 double dac = a.EuclideanDistanceTo(c);
1257 MaskSlicePointType b;
1258 b[0] = a[0]+(c[0]-a[0])/2.0;
1259 b[1] = a[1]+(c[1]-a[1])/2.0;
1261 // Compute distance to trachea centroid
1262 MaskSlicePointType m = centroids_trachea[1];
1263 double da = m.EuclideanDistanceTo(a);
1264 double db = m.EuclideanDistanceTo(b);
1265 //double dc = m.EuclideanDistanceTo(c);
1267 // Mean distance, find point on the line from trachea centroid
1269 double alpha = (da+db)/2.0;
1270 MaskSlicePointType v;
1271 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1272 v[0] = (b[0]-m[0])/n;
1273 v[1] = (b[1]-m[1])/n;
1274 MaskSlicePointType r;
1275 r[0] = m[0]+alpha*v[0];
1276 r[1] = m[1]+alpha*v[1];
1277 points2D.insert(points2D.begin()+index+1, r);
1283 // DDV(points2D, points2D.size());
1285 // Add some points to close the contour
1286 // - H line towards Right
1287 MaskSlicePointType p = points2D[0];
1288 p[0] = bounds[0][0];
1289 points2D.insert(points2D.begin(), p);
1290 // - corner Right/Post
1292 points2D.insert(points2D.begin(), p);
1293 // - H line towards Left
1294 p = points2D.back();
1295 p[0] = bounds[2][0];
1296 points2D.push_back(p);
1297 // - corner Right/Post
1299 points2D.push_back(p);
1300 // Close contour with the first point
1301 points2D.push_back(points2D[0]);
1302 // DDV(points2D, points2D.size());
1304 // Build 3D points from the 2D points
1305 std::vector<ImagePointType> points3D;
1306 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1307 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1309 // Build the mesh from the contour's points
1310 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1311 append->AddInput(mesh);
1312 // if (i ==n-1) { // last slice
1313 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1314 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1315 // append->AddInput(mesh);
1319 // DEBUG: write points3D
1320 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1322 // Build the final 3D mesh form the list 2D mesh
1324 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1326 // Debug, write the mesh
1328 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1330 w->SetFileName("bidon.vtk");
1334 // Compute a single binary 3D image from the list of contours
1335 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1336 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1337 filter->SetMesh(mesh);
1338 filter->SetLikeImage(support);
1340 binarizedContour = filter->GetOutput();
1343 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1344 GetForegroundValue(), 2, true, p);
1346 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1347 GetForegroundValue(), 2, false, p);
1348 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1349 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1350 false, GetBackgroundValue());
1353 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1354 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1356 return binarizedContour;
1359 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1360 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1361 p[1] += 50; // 50 mm Post from this point must be kept
1362 ImageIndexType index;
1363 binarizedContour->TransformPhysicalPointToIndex(p, index);
1364 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1366 // remove from support
1367 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1368 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1369 boolFilter->InPlaceOn();
1370 boolFilter->SetInput1(m_Working_Support);
1371 boolFilter->SetInput2(binarizedContour);
1372 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1373 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1375 boolFilter->SetOperationType(BoolFilterType::And);
1377 boolFilter->SetOperationType(BoolFilterType::AndNot);
1378 boolFilter->Update();
1379 m_Working_Support = boolFilter->GetOutput();
1383 //StopCurrentStep<MaskImageType>(m_Working_Support);
1384 //m_ListOfStations["2R"] = m_Working_Support;
1385 //m_ListOfStations["2L"] = m_Working_Support;
1387 //--------------------------------------------------------------------
1391 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX