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();
66 ExtractStation_5_SetDefaultValues();
67 ExtractStation_6_SetDefaultValues();
70 ExtractStation_7_SetDefaultValues();
71 ExtractStation_8_SetDefaultValues();
73 //--------------------------------------------------------------------
76 //--------------------------------------------------------------------
77 template <class TImageType>
79 clitk::ExtractLymphStationsFilter<TImageType>::
80 GenerateOutputInformation() {
83 m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
84 m_Mediastinum = this->GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
86 // Clean some computer landmarks to force the recomputation
87 this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
89 // Global supports for stations
90 bool supportsExist = true;
92 m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
93 m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
94 m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
95 m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
96 m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
97 m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
99 m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
100 m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
101 m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
102 m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
103 m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
104 m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
105 m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
106 m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
107 m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
108 } catch(clitk::ExceptionObject o) {
109 supportsExist = false;
112 if (!supportsExist || GetComputeStationsSupportsFlag()) {
113 this->StartNewStep("Supports for stations");
114 this->StartSubStep();
115 this->GetAFDB()->RemoveTag("CarinaZ");
116 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
117 this->GetAFDB()->RemoveTag("ApexOfTheChest");
118 this->GetAFDB()->RemoveTag("RightBronchus");
119 this->GetAFDB()->RemoveTag("LeftBronchus");
120 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
121 this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
122 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
123 this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
124 ExtractStationSupports();
129 ExtractStation_1RL();
130 ExtractStation_2RL();
133 ExtractStation_4RL();
137 // ---------- TODO -----------------------
140 // ExtractStation_8();
143 //this->StartNewStep("Station 7");
144 //this->StartSubStep();
145 //ExtractStation_7();
146 //this->StopSubStep();
149 //--------------------------------------------------------------------
152 //--------------------------------------------------------------------
153 template <class TImageType>
155 clitk::ExtractLymphStationsFilter<TImageType>::
156 GenerateInputRequestedRegion() {
157 //DD("GenerateInputRequestedRegion (nothing?)");
159 //--------------------------------------------------------------------
162 //--------------------------------------------------------------------
163 template <class TImageType>
165 clitk::ExtractLymphStationsFilter<TImageType>::
167 // Final Step -> graft output (if SetNthOutput => redo)
168 this->GraftOutput(m_ListOfStations["8"]);
170 //--------------------------------------------------------------------
173 //--------------------------------------------------------------------
174 template <class TImageType>
176 clitk::ExtractLymphStationsFilter<TImageType>::
177 CheckForStation(std::string station)
179 // Compute Station name
180 std::string s = "Station"+station;
183 // Check if station already exist in DB
185 if (this->GetAFDB()->TagExist(s)) {
186 m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
190 // Define the starting support
191 if (found && GetComputeStation(station)) {
192 std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
194 if (!found || GetComputeStation(station)) {
195 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
199 std::cout << "Station " << station << " found. I used it" << std::endl;
203 //--------------------------------------------------------------------
206 //--------------------------------------------------------------------
207 template <class TImageType>
209 clitk::ExtractLymphStationsFilter<TImageType>::
210 GetComputeStation(std::string station)
212 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
214 //--------------------------------------------------------------------
217 //--------------------------------------------------------------------
218 template <class TImageType>
220 clitk::ExtractLymphStationsFilter<TImageType>::
221 AddComputeStation(std::string station)
223 m_ComputeStationMap[station] = true;
225 //--------------------------------------------------------------------
229 //--------------------------------------------------------------------
230 template<class PointType>
231 class comparePointsX {
233 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
235 //--------------------------------------------------------------------
238 //--------------------------------------------------------------------
239 template<class PairType>
240 class comparePointsWithAngle {
242 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
244 //--------------------------------------------------------------------
247 //--------------------------------------------------------------------
249 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
250 std::vector<itk::Point<double, Dim-1> > previous;
251 HypercubeCorners<Dim-1>(previous);
252 out.resize(previous.size()*2);
253 for(unsigned int i=0; i<out.size(); i++) {
254 itk::Point<double, Dim> p;
255 if (i<previous.size()) p[0] = 0;
257 for(int j=0; j<Dim-1; j++)
259 p[j+1] = previous[i%previous.size()][j];
264 //--------------------------------------------------------------------
267 //--------------------------------------------------------------------
269 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
274 //--------------------------------------------------------------------
277 //--------------------------------------------------------------------
278 template<class ImageType>
279 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
280 std::vector<typename ImageType::PointType> & bounds)
282 // Get image max/min coordinates
283 const unsigned int dim=ImageType::ImageDimension;
284 typedef typename ImageType::PointType PointType;
285 typedef typename ImageType::IndexType IndexType;
286 PointType min_c, max_c;
287 IndexType min_i, max_i;
288 min_i = image->GetLargestPossibleRegion().GetIndex();
289 for(unsigned int i=0; i<dim; i++)
290 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
291 image->TransformIndexToPhysicalPoint(min_i, min_c);
292 image->TransformIndexToPhysicalPoint(max_i, max_c);
294 // Get corners coordinates
295 HypercubeCorners<ImageType::ImageDimension>(bounds);
296 for(unsigned int i=0; i<bounds.size(); i++) {
297 for(unsigned int j=0; j<dim; j++) {
298 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
299 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
303 //--------------------------------------------------------------------
306 //--------------------------------------------------------------------
307 template <class ImageType>
309 clitk::ExtractLymphStationsFilter<ImageType>::
310 Remove_Structures(std::string station, std::string s)
313 this->StartNewStep("[Station"+station+"] Remove "+s);
314 MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
315 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
317 catch(clitk::ExceptionObject e) {
318 std::cout << s << " not found, skip." << std::endl;
321 //--------------------------------------------------------------------
324 //--------------------------------------------------------------------
325 template <class TImageType>
327 clitk::ExtractLymphStationsFilter<TImageType>::
328 SetFuzzyThreshold(std::string station, std::string tag, double value)
330 m_FuzzyThreshold[station][tag] = value;
332 //--------------------------------------------------------------------
335 //--------------------------------------------------------------------
336 template <class TImageType>
338 clitk::ExtractLymphStationsFilter<TImageType>::
339 SetThreshold(std::string station, std::string tag, double value)
341 m_Threshold[station][tag] = value;
343 //--------------------------------------------------------------------
346 //--------------------------------------------------------------------
347 template <class TImageType>
349 clitk::ExtractLymphStationsFilter<TImageType>::
350 GetFuzzyThreshold(std::string station, std::string tag)
352 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
353 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
357 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
358 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
362 return m_FuzzyThreshold[station][tag];
364 //--------------------------------------------------------------------
367 //--------------------------------------------------------------------
368 template <class TImageType>
370 clitk::ExtractLymphStationsFilter<TImageType>::
371 GetThreshold(std::string station, std::string tag)
373 if (m_Threshold.find(station) == m_Threshold.end()) {
374 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
378 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
379 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
383 return m_Threshold[station][tag];
385 //--------------------------------------------------------------------
388 //--------------------------------------------------------------------
389 template <class TImageType>
391 clitk::ExtractLymphStationsFilter<TImageType>::
392 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
394 // Create line from A to B with
395 // A = upper border of LLL at left
396 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
399 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
400 this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
402 catch(clitk::ExceptionObject & o) {
404 DD("FindLineForS7S8Separation");
405 // Load LeftLowerLobeBronchus and get centroid point
406 MaskImagePointer LeftLowerLobeBronchus =
407 this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
408 std::vector<MaskImagePointType> c;
409 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
412 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
413 MaskImagePointer RightMiddleLobeBronchus =
414 this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
415 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
416 GetBackgroundValue(),
419 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
422 // Insert into the DB
423 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
424 this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
427 //--------------------------------------------------------------------
430 //--------------------------------------------------------------------
431 template <class TImageType>
433 clitk::ExtractLymphStationsFilter<TImageType>::
438 z = this->GetAFDB()->GetDouble("CarinaZ");
440 catch(clitk::ExceptionObject e) {
441 DD("FindCarinaSlicePosition");
443 MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
445 // Get Centroid and Z value
446 std::vector<MaskImagePointType> centroids;
447 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
449 // We dont need Carina structure from now
450 this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
452 // Put inside the AFDB
453 this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
454 this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
460 //--------------------------------------------------------------------
463 //--------------------------------------------------------------------
464 template <class TImageType>
466 clitk::ExtractLymphStationsFilter<TImageType>::
471 z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
473 catch(clitk::ExceptionObject e) {
474 DD("FindApexOfTheChestPosition");
475 MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
476 MaskImagePointType p;
477 p[0] = p[1] = p[2] = 0.0; // to avoid warning
478 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
480 // We dont need Lungs structure from now
481 this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
483 // Put inside the AFDB
484 this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
485 p[2] -= 5; // We consider 5 mm lower
486 this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
492 //--------------------------------------------------------------------
495 //--------------------------------------------------------------------
496 template <class TImageType>
498 clitk::ExtractLymphStationsFilter<TImageType>::
499 FindLeftAndRightBronchi()
502 m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
503 m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
505 catch(clitk::ExceptionObject & o) {
507 DD("FindLeftAndRightBronchi");
508 // The goal is to separate the trachea inferiorly to the carina into
509 // a Left and Right bronchus.
512 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
514 // Get the Carina position
515 double m_CarinaZ = FindCarina();
517 // Consider only inferiorly to the Carina
518 MaskImagePointer m_Working_Trachea =
519 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
520 GetBackgroundValue());
522 // Labelize the trachea
523 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
525 // Carina position must at the first slice that separate the two
526 // main bronchus (not superiorly). We thus first check that the
527 // upper slice is composed of at least two labels
528 MaskImagePointer RightBronchus;
529 MaskImagePointer LeftBronchus;
530 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
531 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
532 iter.SetFirstDirection(0);
533 iter.SetSecondDirection(1);
534 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
536 while (!iter.IsAtReverseEndOfSlice()) {
537 while (!iter.IsAtReverseEndOfLine()) {
538 if (iter.Get() > maxLabel) maxLabel = iter.Get();
544 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
547 // Compute 3D centroids of both parts to identify the left from the
549 std::vector<ImagePointType> c;
550 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
551 ImagePointType C1 = c[1];
552 ImagePointType C2 = c[2];
554 ImagePixelType rightLabel;
555 ImagePixelType leftLabel;
556 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
557 else { rightLabel = 2; leftLabel = 1; }
559 // Select LeftLabel (set one label to Backgroundvalue)
561 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
562 GetBackgroundValue(), GetForegroundValue());
564 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
565 leftLabel, GetBackgroundValue(), false);
567 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
568 GetBackgroundValue(), GetForegroundValue());
570 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
571 rightLabel, GetBackgroundValue(), false);
575 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
576 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
578 // Insert int AFDB if need after
579 this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
580 RightBronchus, true);
581 this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
585 //--------------------------------------------------------------------
588 //--------------------------------------------------------------------
589 template <class TImageType>
591 clitk::ExtractLymphStationsFilter<TImageType>::
592 FindSuperiorBorderOfAorticArch()
596 z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
598 catch(clitk::ExceptionObject e) {
599 DD("FindSuperiorBorderOfAorticArch");
600 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
601 MaskImagePointType p;
602 p[0] = p[1] = p[2] = 0.0; // to avoid warning
603 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
604 p[2] += Aorta->GetSpacing()[2]; // the slice above
606 // We dont need Lungs structure from now
607 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
609 // Put inside the AFDB
610 this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
611 this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
617 //--------------------------------------------------------------------
620 //--------------------------------------------------------------------
621 template <class TImageType>
623 clitk::ExtractLymphStationsFilter<TImageType>::
624 FindInferiorBorderOfAorticArch()
628 z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
630 catch(clitk::ExceptionObject e) {
631 DD("FindInferiorBorderOfAorticArch");
632 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
633 std::vector<MaskSlicePointer> slices;
634 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
636 uint i = slices.size()-1;
638 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
639 std::vector<typename MaskSliceType::PointType> c;
640 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
648 MaskImageIndexType index;
649 index[0] = index[1] = 0.0;
651 MaskImagePointType lower;
652 Aorta->TransformIndexToPhysicalPoint(index, lower);
654 // We dont need Lungs structure from now
655 this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
657 // Put inside the AFDB
658 this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
659 this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
665 //--------------------------------------------------------------------
668 //--------------------------------------------------------------------
669 template <class ImageType>
670 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
671 clitk::ExtractLymphStationsFilter<ImageType>::
672 FindAntPostVesselsOLD()
674 // -----------------------------------------------------
675 /* Rod says: "The anterior border, as with the Atlas – UM, is
676 posterior to the vessels (right subclavian vein, left
677 brachiocephalic vein, right brachiocephalic vein, left subclavian
678 artery, left common carotid artery and brachiocephalic trunk).
679 These vessels are not included in the nodal station. The anterior
680 border is drawn to the midpoint of the vessel and an imaginary
681 line joins the middle of these vessels. Between the vessels,
682 station 2 is in contact with station 3a." */
684 // Check if no already done
686 MaskImagePointer binarizedContour;
688 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
690 catch(clitk::ExceptionObject e) {
694 return binarizedContour;
697 /* Here, we consider the vessels form a kind of anterior barrier. We
698 link all vessels centroids and remove what is post to it.
699 - select the list of structure
700 vessel1 = BrachioCephalicArtery
701 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
702 vessel3 = CommonCarotidArtery
703 vessel4 = SubclavianArtery
706 - crop images as needed
707 - slice by slice, choose the CCL and extract centroids
708 - slice by slice, sort according to polar angle wrt Trachea centroid.
709 - slice by slice, link centroids and close contour
710 - remove outside this contour
715 MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
716 MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
717 MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
718 MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
719 MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
720 MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
721 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
723 // Create a temporay support
724 // From first slice of BrachioCephalicVein to end of 3A
725 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
726 MaskImagePointType p;
727 p[0] = p[1] = p[2] = 0.0; // to avoid warning
728 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
730 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
731 GetBackgroundValue(), 2, false, p);
733 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
734 false, GetBackgroundValue());
736 // Resize all structures like support
737 BrachioCephalicArtery =
738 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
739 CommonCarotidArtery =
740 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
742 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
744 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
746 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
747 BrachioCephalicVein =
748 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
750 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
753 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
754 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
755 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
756 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
757 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
758 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
759 std::vector<MaskSlicePointer> slices_SubclavianArtery;
760 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
761 std::vector<MaskSlicePointer> slices_Thyroid;
762 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
763 std::vector<MaskSlicePointer> slices_Aorta;
764 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
765 std::vector<MaskSlicePointer> slices_Trachea;
766 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
767 unsigned int n= slices_BrachioCephalicArtery.size();
769 // Get the boundaries of one slice
770 std::vector<MaskSlicePointType> bounds;
771 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
773 // For all slices, for all structures, find the centroid and build the contour
774 // List of 3D points (for debug)
775 std::vector<MaskImagePointType> p3D;
777 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
778 for(unsigned int i=0; i<n; i++) {
779 // Labelize the slices
780 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
781 GetBackgroundValue(), true, 1);
782 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
783 GetBackgroundValue(), true, 1);
784 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
785 GetBackgroundValue(), true, 1);
786 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
787 GetBackgroundValue(), true, 1);
788 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
789 GetBackgroundValue(), true, 1);
790 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
791 GetBackgroundValue(), true, 1);
794 std::vector<MaskSlicePointType> points2D;
795 std::vector<MaskSlicePointType> centroids1;
796 std::vector<MaskSlicePointType> centroids2;
797 std::vector<MaskSlicePointType> centroids3;
798 std::vector<MaskSlicePointType> centroids4;
799 std::vector<MaskSlicePointType> centroids5;
800 std::vector<MaskSlicePointType> centroids6;
801 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
802 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
803 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
804 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
805 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
806 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
808 // BrachioCephalicVein -> when it is separated into two CCL, we
809 // only consider the most at Right one
810 if (centroids6.size() > 2) {
811 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
812 else centroids6.erase(centroids6.begin()+1);
815 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
816 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
817 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
821 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
822 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
823 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
824 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
825 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
826 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
828 // Sort by angle according to trachea centroid and vertical line,
829 // in polar coordinates :
830 // http://en.wikipedia.org/wiki/Polar_coordinate_system
831 std::vector<MaskSlicePointType> centroids_trachea;
832 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
833 typedef std::pair<MaskSlicePointType, double> PointAngleType;
834 std::vector<PointAngleType> angles;
835 for(unsigned int j=0; j<points2D.size(); j++) {
836 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
837 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
838 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
840 if (x>0) angle = atan(y/x);
841 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
842 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
844 if (y>0) angle = M_PI/2.0;
845 if (y<0) angle = -M_PI/2.0;
848 angle = clitk::rad2deg(angle);
849 // Angle is [-180;180] wrt the X axis. We change the X axis to
850 // be the vertical line, because we want to sort from Right to
851 // Left from Post to Ant.
856 if (angle<0) angle = 360-angle;
859 a.first = points2D[j];
864 // Do nothing if less than 2 points --> n
865 if (points2D.size() < 3) { //continue;
869 // Sort points2D according to polar angles
870 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
871 for(unsigned int j=0; j<angles.size(); j++) {
872 points2D[j] = angles[j].first;
874 // DDV(points2D, points2D.size());
876 /* When vessels are far away, we try to replace the line segment
877 with a curved line that join the two vessels but stay
878 approximately at the same distance from the trachea centroids
882 - let a and c be two successive vessels centroids
883 - id distance(a,c) < threshold, next point
887 - compute mid position between two successive points -
888 compute dist to trachea centroid for the 3 pts - if middle too
891 std::vector<MaskSlicePointType> toadd;
892 unsigned int index = 0;
894 while (index<points2D.size()-1) {
895 MaskSlicePointType a = points2D[index];
896 MaskSlicePointType c = points2D[index+1];
898 double dac = a.EuclideanDistanceTo(c);
901 MaskSlicePointType b;
902 b[0] = a[0]+(c[0]-a[0])/2.0;
903 b[1] = a[1]+(c[1]-a[1])/2.0;
905 // Compute distance to trachea centroid
906 MaskSlicePointType m = centroids_trachea[1];
907 double da = m.EuclideanDistanceTo(a);
908 double db = m.EuclideanDistanceTo(b);
909 //double dc = m.EuclideanDistanceTo(c);
911 // Mean distance, find point on the line from trachea centroid
913 double alpha = (da+db)/2.0;
914 MaskSlicePointType v;
915 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
916 v[0] = (b[0]-m[0])/n;
917 v[1] = (b[1]-m[1])/n;
918 MaskSlicePointType r;
919 r[0] = m[0]+alpha*v[0];
920 r[1] = m[1]+alpha*v[1];
921 points2D.insert(points2D.begin()+index+1, r);
927 // DDV(points2D, points2D.size());
929 // Add some points to close the contour
930 // - H line towards Right
931 MaskSlicePointType p = points2D[0];
933 points2D.insert(points2D.begin(), p);
934 // - corner Right/Post
936 points2D.insert(points2D.begin(), p);
937 // - H line towards Left
940 points2D.push_back(p);
941 // - corner Right/Post
943 points2D.push_back(p);
944 // Close contour with the first point
945 points2D.push_back(points2D[0]);
946 // DDV(points2D, points2D.size());
948 // Build 3D points from the 2D points
949 std::vector<ImagePointType> points3D;
950 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
951 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
953 // Build the mesh from the contour's points
954 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
955 append->AddInput(mesh);
958 // DEBUG: write points3D
959 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
961 // Build the final 3D mesh form the list 2D mesh
963 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
965 // Debug, write the mesh
967 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
969 w->SetFileName("bidon.vtk");
973 // Compute a single binary 3D image from the list of contours
974 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
975 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
976 filter->SetMesh(mesh);
977 filter->SetLikeImage(support);
979 binarizedContour = filter->GetOutput();
982 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
985 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
988 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
989 false, GetBackgroundValue());
991 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
992 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
994 return binarizedContour;
997 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
998 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
999 p[1] += 50; // 50 mm Post from this point must be kept
1000 ImageIndexType index;
1001 binarizedContour->TransformPhysicalPointToIndex(p, index);
1002 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1004 // remove from support
1005 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1006 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1007 boolFilter->InPlaceOn();
1008 boolFilter->SetInput1(m_Working_Support);
1009 boolFilter->SetInput2(binarizedContour);
1010 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1011 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1013 boolFilter->SetOperationType(BoolFilterType::And);
1015 boolFilter->SetOperationType(BoolFilterType::AndNot);
1016 boolFilter->Update();
1017 m_Working_Support = boolFilter->GetOutput();
1021 //StopCurrentStep<MaskImageType>(m_Working_Support);
1022 //m_ListOfStations["2R"] = m_Working_Support;
1023 //m_ListOfStations["2L"] = m_Working_Support;
1025 //--------------------------------------------------------------------
1028 //--------------------------------------------------------------------
1029 template <class ImageType>
1030 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1031 clitk::ExtractLymphStationsFilter<ImageType>::
1032 FindAntPostVessels2()
1034 // -----------------------------------------------------
1035 /* Rod says: "The anterior border, as with the Atlas – UM, is
1036 posterior to the vessels (right subclavian vein, left
1037 brachiocephalic vein, right brachiocephalic vein, left subclavian
1038 artery, left common carotid artery and brachiocephalic trunk).
1039 These vessels are not included in the nodal station. The anterior
1040 border is drawn to the midpoint of the vessel and an imaginary
1041 line joins the middle of these vessels. Between the vessels,
1042 station 2 is in contact with station 3a." */
1044 // Check if no already done
1046 MaskImagePointer binarizedContour;
1048 binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1050 catch(clitk::ExceptionObject e) {
1054 return binarizedContour;
1057 /* Here, we consider the vessels form a kind of anterior barrier. We
1058 link all vessels centroids and remove what is post to it.
1059 - select the list of structure
1060 vessel1 = BrachioCephalicArtery
1061 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1062 vessel3 = CommonCarotidArtery
1063 vessel4 = SubclavianArtery
1066 - crop images as needed
1067 - slice by slice, choose the CCL and extract centroids
1068 - slice by slice, sort according to polar angle wrt Trachea centroid.
1069 - slice by slice, link centroids and close contour
1070 - remove outside this contour
1071 - merge with support
1075 std::map<std::string, MaskImagePointer> MapOfStructures;
1076 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1077 MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1078 MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1079 MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftCommonCarotidArtery");
1080 MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightCommonCarotidArtery");
1081 MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftSubclavianArtery");
1082 MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightSubclavianArtery");
1083 MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1084 MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1085 MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1087 std::vector<std::string> ListOfStructuresNames;
1089 // Create a temporay support
1090 // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1091 MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1092 MaskImagePointType p;
1093 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1094 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1095 GetBackgroundValue(), 2, true, p);
1097 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1098 GetBackgroundValue(), 2, false, p);
1099 MaskImagePointType p2;
1100 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1101 GetBackgroundValue(), 2, false, p2);
1102 if (p2[2] > p[2]) p = p2;
1104 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1105 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1106 false, GetBackgroundValue());
1108 // Resize all structures like support
1109 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1110 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1114 typedef std::vector<MaskSlicePointer> SlicesType;
1115 std::map<std::string, SlicesType> MapOfSlices;
1116 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1118 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1119 MapOfSlices[it->first] = s;
1122 unsigned int n= MapOfSlices["Trachea"].size();
1124 // Get the boundaries of one slice
1125 std::vector<MaskSlicePointType> bounds;
1126 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1129 // For all slices, for all structures, find the centroid and build the contour
1130 // List of 3D points (for debug)
1131 std::vector<MaskImagePointType> p3D;
1132 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1133 for(unsigned int i=0; i<n; i++) {
1135 // Labelize the slices
1136 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1137 MaskSlicePointer & s = MapOfSlices[it->first][i];
1138 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1142 std::vector<MaskSlicePointType> points2D;
1143 typedef std::vector<MaskSlicePointType> CentroidsType;
1144 std::map<std::string, CentroidsType> MapOfCentroids;
1145 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1146 std::string structure = it->first;
1147 MaskSlicePointer & s = MapOfSlices[structure][i];
1149 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1150 MapOfCentroids[structure] = c;
1154 // BrachioCephalicVein -> when it is separated into two CCL, we
1155 // only consider the most at Right one
1156 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1158 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1159 else c.erase(c.begin()+1);
1163 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1164 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1165 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1166 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1167 MapOfCentroids["BrachioCephalicVein"].clear();
1171 // Add all 2D points
1172 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1173 std::string structure = it->first;
1174 if (structure != "Trachea") { // do not add centroid of trachea
1175 CentroidsType & c = MapOfCentroids[structure];
1176 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1180 // Sort by angle according to trachea centroid and vertical line,
1181 // in polar coordinates :
1182 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1183 // std::vector<MaskSlicePointType> centroids_trachea;
1184 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1185 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1186 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1187 std::vector<PointAngleType> angles;
1188 for(unsigned int j=0; j<points2D.size(); j++) {
1189 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1190 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1191 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1193 if (x>0) angle = atan(y/x);
1194 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1195 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1197 if (y>0) angle = M_PI/2.0;
1198 if (y<0) angle = -M_PI/2.0;
1199 if (y==0) angle = 0;
1201 angle = clitk::rad2deg(angle);
1202 // Angle is [-180;180] wrt the X axis. We change the X axis to
1203 // be the vertical line, because we want to sort from Right to
1204 // Left from Post to Ant.
1206 angle = (270-angle);
1209 if (angle<0) angle = 360-angle;
1212 a.first = points2D[j];
1214 angles.push_back(a);
1217 // Do nothing if less than 2 points --> n
1218 if (points2D.size() < 3) { //continue;
1222 // Sort points2D according to polar angles
1223 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1224 for(unsigned int j=0; j<angles.size(); j++) {
1225 points2D[j] = angles[j].first;
1227 // DDV(points2D, points2D.size());
1229 /* When vessels are far away, we try to replace the line segment
1230 with a curved line that join the two vessels but stay
1231 approximately at the same distance from the trachea centroids
1235 - let a and c be two successive vessels centroids
1236 - id distance(a,c) < threshold, next point
1240 - compute mid position between two successive points -
1241 compute dist to trachea centroid for the 3 pts - if middle too
1244 std::vector<MaskSlicePointType> toadd;
1245 unsigned int index = 0;
1247 while (index<points2D.size()-1) {
1248 MaskSlicePointType a = points2D[index];
1249 MaskSlicePointType c = points2D[index+1];
1251 double dac = a.EuclideanDistanceTo(c);
1254 MaskSlicePointType b;
1255 b[0] = a[0]+(c[0]-a[0])/2.0;
1256 b[1] = a[1]+(c[1]-a[1])/2.0;
1258 // Compute distance to trachea centroid
1259 MaskSlicePointType m = centroids_trachea[1];
1260 double da = m.EuclideanDistanceTo(a);
1261 double db = m.EuclideanDistanceTo(b);
1262 //double dc = m.EuclideanDistanceTo(c);
1264 // Mean distance, find point on the line from trachea centroid
1266 double alpha = (da+db)/2.0;
1267 MaskSlicePointType v;
1268 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1269 v[0] = (b[0]-m[0])/n;
1270 v[1] = (b[1]-m[1])/n;
1271 MaskSlicePointType r;
1272 r[0] = m[0]+alpha*v[0];
1273 r[1] = m[1]+alpha*v[1];
1274 points2D.insert(points2D.begin()+index+1, r);
1280 // DDV(points2D, points2D.size());
1282 // Add some points to close the contour
1283 // - H line towards Right
1284 MaskSlicePointType p = points2D[0];
1285 p[0] = bounds[0][0];
1286 points2D.insert(points2D.begin(), p);
1287 // - corner Right/Post
1289 points2D.insert(points2D.begin(), p);
1290 // - H line towards Left
1291 p = points2D.back();
1292 p[0] = bounds[2][0];
1293 points2D.push_back(p);
1294 // - corner Right/Post
1296 points2D.push_back(p);
1297 // Close contour with the first point
1298 points2D.push_back(points2D[0]);
1299 // DDV(points2D, points2D.size());
1301 // Build 3D points from the 2D points
1302 std::vector<ImagePointType> points3D;
1303 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1304 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1306 // Build the mesh from the contour's points
1307 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1308 append->AddInput(mesh);
1309 // if (i ==n-1) { // last slice
1310 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1311 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1312 // append->AddInput(mesh);
1316 // DEBUG: write points3D
1317 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1319 // Build the final 3D mesh form the list 2D mesh
1321 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1323 // Debug, write the mesh
1325 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1327 w->SetFileName("bidon.vtk");
1331 // Compute a single binary 3D image from the list of contours
1332 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1333 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1334 filter->SetMesh(mesh);
1335 filter->SetLikeImage(support);
1337 binarizedContour = filter->GetOutput();
1340 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1341 GetForegroundValue(), 2, true, p);
1343 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1344 GetForegroundValue(), 2, false, p);
1345 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1346 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1347 false, GetBackgroundValue());
1350 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1351 this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1353 return binarizedContour;
1356 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1357 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1358 p[1] += 50; // 50 mm Post from this point must be kept
1359 ImageIndexType index;
1360 binarizedContour->TransformPhysicalPointToIndex(p, index);
1361 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1363 // remove from support
1364 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1365 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1366 boolFilter->InPlaceOn();
1367 boolFilter->SetInput1(m_Working_Support);
1368 boolFilter->SetInput2(binarizedContour);
1369 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1370 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1372 boolFilter->SetOperationType(BoolFilterType::And);
1374 boolFilter->SetOperationType(BoolFilterType::AndNot);
1375 boolFilter->Update();
1376 m_Working_Support = boolFilter->GetOutput();
1380 //StopCurrentStep<MaskImageType>(m_Working_Support);
1381 //m_ListOfStations["2R"] = m_Working_Support;
1382 //m_ListOfStations["2L"] = m_Working_Support;
1384 //--------------------------------------------------------------------
1388 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX