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():
54 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
55 itk::ImageToImageFilter<TImageType, MaskImageType>()
57 this->SetNumberOfRequiredInputs(1);
58 SetBackgroundValue(0);
59 SetForegroundValue(1);
60 ComputeStationsSupportsFlagOn();
63 ExtractStation_8_SetDefaultValues();
64 ExtractStation_3P_SetDefaultValues();
65 ExtractStation_2RL_SetDefaultValues();
66 ExtractStation_3A_SetDefaultValues();
67 ExtractStation_7_SetDefaultValues();
69 //--------------------------------------------------------------------
72 //--------------------------------------------------------------------
73 template <class TImageType>
75 clitk::ExtractLymphStationsFilter<TImageType>::
76 GenerateOutputInformation() {
79 m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
80 m_Mediastinum = GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
82 // Clean some computer landmarks to force the recomputation
83 GetAFDB()->RemoveTag("AntPostVesselsSeparation");
85 // Global supports for stations
86 if (GetComputeStationsSupportsFlag()) {
87 StartNewStep("Supports for stations");
89 GetAFDB()->RemoveTag("CarinaZ");
90 GetAFDB()->RemoveTag("ApexOfTheChestZ");
91 GetAFDB()->RemoveTag("ApexOfTheChest");
92 GetAFDB()->RemoveTag("RightBronchus");
93 GetAFDB()->RemoveTag("LeftBronchus");
94 GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
95 GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
96 GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
97 GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
98 ExtractStationSupports();
102 m_ListOfSupports["S1R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
103 m_ListOfSupports["S1L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
104 m_ListOfSupports["S2R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
105 m_ListOfSupports["S2L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
106 m_ListOfSupports["S4R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
107 m_ListOfSupports["S4L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
109 m_ListOfSupports["S3A"] = GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
110 m_ListOfSupports["S3P"] = GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
111 m_ListOfSupports["S5"] = GetAFDB()->template GetImage<MaskImageType>("Support_S5");
112 m_ListOfSupports["S6"] = GetAFDB()->template GetImage<MaskImageType>("Support_S6");
113 m_ListOfSupports["S7"] = GetAFDB()->template GetImage<MaskImageType>("Support_S7");
114 m_ListOfSupports["S8"] = GetAFDB()->template GetImage<MaskImageType>("Support_S8");
115 m_ListOfSupports["S9"] = GetAFDB()->template GetImage<MaskImageType>("Support_S9");
116 m_ListOfSupports["S10"] = GetAFDB()->template GetImage<MaskImageType>("Support_S10");
117 m_ListOfSupports["S11"] = GetAFDB()->template GetImage<MaskImageType>("Support_S11");
121 StartNewStep("Station 8");
127 StartNewStep("Station 3P");
133 StartNewStep("Station 3A");
141 // Extract Station2RL
142 StartNewStep("Station 2RL");
144 ExtractStation_2RL();
148 StartNewStep("Station 7");
153 if (0) { // temporary suppress
154 // Extract Station4RL
155 StartNewStep("Station 4RL");
157 //ExtractStation_4RL();
161 //--------------------------------------------------------------------
164 //--------------------------------------------------------------------
165 template <class TImageType>
167 clitk::ExtractLymphStationsFilter<TImageType>::
168 GenerateInputRequestedRegion() {
169 //DD("GenerateInputRequestedRegion (nothing?)");
171 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 template <class TImageType>
177 clitk::ExtractLymphStationsFilter<TImageType>::
179 // Final Step -> graft output (if SetNthOutput => redo)
180 this->GraftOutput(m_ListOfStations["8"]);
182 //--------------------------------------------------------------------
185 //--------------------------------------------------------------------
186 template <class TImageType>
188 clitk::ExtractLymphStationsFilter<TImageType>::
189 CheckForStation(std::string station)
191 // Compute Station name
192 std::string s = "Station"+station;
195 // Check if station already exist in DB
197 if (GetAFDB()->TagExist(s)) {
198 m_ListOfStations[station] = GetAFDB()->template GetImage<MaskImageType>(s);
202 // Define the starting support
203 if (found && GetComputeStation(station)) {
204 std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
206 if (!found || GetComputeStation(station)) {
207 m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
211 std::cout << "Station " << station << " found. I used it" << std::endl;
215 //--------------------------------------------------------------------
218 //--------------------------------------------------------------------
219 template <class TImageType>
221 clitk::ExtractLymphStationsFilter<TImageType>::
222 GetComputeStation(std::string station)
224 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
226 //--------------------------------------------------------------------
229 //--------------------------------------------------------------------
230 template <class TImageType>
232 clitk::ExtractLymphStationsFilter<TImageType>::
233 AddComputeStation(std::string station)
235 m_ComputeStationMap[station] = true;
237 //--------------------------------------------------------------------
241 //--------------------------------------------------------------------
242 template<class PointType>
243 class comparePointsX {
245 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
247 //--------------------------------------------------------------------
250 //--------------------------------------------------------------------
251 template<class PairType>
252 class comparePointsWithAngle {
254 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
256 //--------------------------------------------------------------------
259 //--------------------------------------------------------------------
261 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
262 std::vector<itk::Point<double, Dim-1> > previous;
263 HypercubeCorners<Dim-1>(previous);
264 out.resize(previous.size()*2);
265 for(unsigned int i=0; i<out.size(); i++) {
266 itk::Point<double, Dim> p;
267 if (i<previous.size()) p[0] = 0;
269 for(int j=0; j<Dim-1; j++)
271 p[j+1] = previous[i%previous.size()][j];
276 //--------------------------------------------------------------------
279 //--------------------------------------------------------------------
281 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
286 //--------------------------------------------------------------------
289 //--------------------------------------------------------------------
290 template<class ImageType>
291 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
292 std::vector<typename ImageType::PointType> & bounds)
294 // Get image max/min coordinates
295 const unsigned int dim=ImageType::ImageDimension;
296 typedef typename ImageType::PointType PointType;
297 typedef typename ImageType::IndexType IndexType;
298 PointType min_c, max_c;
299 IndexType min_i, max_i;
300 min_i = image->GetLargestPossibleRegion().GetIndex();
301 for(unsigned int i=0; i<dim; i++)
302 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
303 image->TransformIndexToPhysicalPoint(min_i, min_c);
304 image->TransformIndexToPhysicalPoint(max_i, max_c);
306 // Get corners coordinates
307 HypercubeCorners<ImageType::ImageDimension>(bounds);
308 for(unsigned int i=0; i<bounds.size(); i++) {
309 for(unsigned int j=0; j<dim; j++) {
310 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
311 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
315 //--------------------------------------------------------------------
318 //--------------------------------------------------------------------
319 template <class TImageType>
321 clitk::ExtractLymphStationsFilter<TImageType>::
322 ExtractStation_4RL() {
326 WARNING ONLY 4R FIRST !!! (not same inf limits)
328 ExtractStation_4RL_SI_Limits();
329 ExtractStation_4RL_LR_Limits();
330 ExtractStation_4RL_AP_Limits();
332 //--------------------------------------------------------------------
335 //--------------------------------------------------------------------
336 template <class ImageType>
338 clitk::ExtractLymphStationsFilter<ImageType>::
339 Remove_Structures(std::string station, std::string s)
342 StartNewStep("[Station"+station+"] Remove "+s);
343 MaskImagePointer Structure = GetAFDB()->template GetImage<MaskImageType>(s);
344 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
346 catch(clitk::ExceptionObject e) {
347 std::cout << s << " not found, skip." << std::endl;
350 //--------------------------------------------------------------------
353 //--------------------------------------------------------------------
354 template <class TImageType>
356 clitk::ExtractLymphStationsFilter<TImageType>::
357 SetFuzzyThreshold(std::string station, std::string tag, double value)
359 m_FuzzyThreshold[station][tag] = value;
361 //--------------------------------------------------------------------
364 //--------------------------------------------------------------------
365 template <class TImageType>
367 clitk::ExtractLymphStationsFilter<TImageType>::
368 GetFuzzyThreshold(std::string station, std::string tag)
370 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
371 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
375 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
376 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
380 return m_FuzzyThreshold[station][tag];
382 //--------------------------------------------------------------------
385 //--------------------------------------------------------------------
386 template <class TImageType>
388 clitk::ExtractLymphStationsFilter<TImageType>::
389 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
391 // Create line from A to B with
392 // A = upper border of LLL at left
393 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
396 GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
397 GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
399 catch(clitk::ExceptionObject & o) {
401 DD("FindLineForS7S8Separation");
402 // Load LeftLowerLobeBronchus and get centroid point
403 MaskImagePointer LeftLowerLobeBronchus =
404 GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
405 std::vector<MaskImagePointType> c;
406 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
409 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
410 MaskImagePointer RightMiddleLobeBronchus =
411 GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
412 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
413 GetBackgroundValue(),
416 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
419 // Insert into the DB
420 GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
421 GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
424 //--------------------------------------------------------------------
427 //--------------------------------------------------------------------
428 template <class TImageType>
430 clitk::ExtractLymphStationsFilter<TImageType>::
435 z = GetAFDB()->GetDouble("CarinaZ");
437 catch(clitk::ExceptionObject e) {
438 DD("FindCarinaSlicePosition");
440 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
442 // Get Centroid and Z value
443 std::vector<MaskImagePointType> centroids;
444 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
446 // We dont need Carina structure from now
447 GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
449 // Put inside the AFDB
450 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
451 GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
457 //--------------------------------------------------------------------
460 //--------------------------------------------------------------------
461 template <class TImageType>
463 clitk::ExtractLymphStationsFilter<TImageType>::
468 z = GetAFDB()->GetDouble("ApexOfTheChestZ");
470 catch(clitk::ExceptionObject e) {
471 DD("FindApexOfTheChestPosition");
472 MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
473 MaskImagePointType p;
474 p[0] = p[1] = p[2] = 0.0; // to avoid warning
475 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
477 // We dont need Lungs structure from now
478 GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
480 // Put inside the AFDB
481 GetAFDB()->SetPoint3D("ApexOfTheChest", p);
482 p[2] -= 5; // We consider 5 mm lower
483 GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
489 //--------------------------------------------------------------------
492 //--------------------------------------------------------------------
493 template <class TImageType>
495 clitk::ExtractLymphStationsFilter<TImageType>::
496 FindLeftAndRightBronchi()
499 m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
500 m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
502 catch(clitk::ExceptionObject & o) {
504 DD("FindLeftAndRightBronchi");
505 // The goal is to separate the trachea inferiorly to the carina into
506 // a Left and Right bronchus.
509 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
511 // Get the Carina position
512 double m_CarinaZ = FindCarina();
514 // Consider only inferiorly to the Carina
515 MaskImagePointer m_Working_Trachea =
516 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
517 GetBackgroundValue());
519 // Labelize the trachea
520 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
522 // Carina position must at the first slice that separate the two
523 // main bronchus (not superiorly). We thus first check that the
524 // upper slice is composed of at least two labels
525 MaskImagePointer RightBronchus;
526 MaskImagePointer LeftBronchus;
527 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
528 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
529 iter.SetFirstDirection(0);
530 iter.SetSecondDirection(1);
531 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
533 while (!iter.IsAtReverseEndOfSlice()) {
534 while (!iter.IsAtReverseEndOfLine()) {
535 if (iter.Get() > maxLabel) maxLabel = iter.Get();
541 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
544 // Compute 3D centroids of both parts to identify the left from the
546 std::vector<ImagePointType> c;
547 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
548 ImagePointType C1 = c[1];
549 ImagePointType C2 = c[2];
551 ImagePixelType rightLabel;
552 ImagePixelType leftLabel;
553 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
554 else { rightLabel = 2; leftLabel = 1; }
556 // Select LeftLabel (set one label to Backgroundvalue)
558 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
559 GetBackgroundValue(), GetForegroundValue());
561 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
562 leftLabel, GetBackgroundValue(), false);
564 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
565 GetBackgroundValue(), GetForegroundValue());
567 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
568 rightLabel, GetBackgroundValue(), false);
572 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
573 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
575 // Insert int AFDB if need after
576 GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
577 RightBronchus, true);
578 GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
582 //--------------------------------------------------------------------
585 //--------------------------------------------------------------------
586 template <class TImageType>
588 clitk::ExtractLymphStationsFilter<TImageType>::
589 FindSuperiorBorderOfAorticArch()
593 z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
595 catch(clitk::ExceptionObject e) {
596 DD("FindSuperiorBorderOfAorticArch");
597 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
598 MaskImagePointType p;
599 p[0] = p[1] = p[2] = 0.0; // to avoid warning
600 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
601 p[2] += Aorta->GetSpacing()[2]; // the slice above
603 // We dont need Lungs structure from now
604 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
606 // Put inside the AFDB
607 GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
608 GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
614 //--------------------------------------------------------------------
617 //--------------------------------------------------------------------
618 template <class TImageType>
620 clitk::ExtractLymphStationsFilter<TImageType>::
621 FindInferiorBorderOfAorticArch()
625 z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
627 catch(clitk::ExceptionObject e) {
628 DD("FindInferiorBorderOfAorticArch");
629 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
630 std::vector<MaskSlicePointer> slices;
631 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
633 uint i = slices.size()-1;
635 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
636 std::vector<typename MaskSliceType::PointType> c;
637 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
645 MaskImageIndexType index;
646 index[0] = index[1] = 0.0;
648 MaskImagePointType lower;
649 Aorta->TransformIndexToPhysicalPoint(index, lower);
651 // We dont need Lungs structure from now
652 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
654 // Put inside the AFDB
655 GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
656 GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
662 //--------------------------------------------------------------------
665 //--------------------------------------------------------------------
666 template <class ImageType>
667 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
668 clitk::ExtractLymphStationsFilter<ImageType>::
671 // -----------------------------------------------------
672 /* Rod says: "The anterior border, as with the Atlas – UM, is
673 posterior to the vessels (right subclavian vein, left
674 brachiocephalic vein, right brachiocephalic vein, left subclavian
675 artery, left common carotid artery and brachiocephalic trunk).
676 These vessels are not included in the nodal station. The anterior
677 border is drawn to the midpoint of the vessel and an imaginary
678 line joins the middle of these vessels. Between the vessels,
679 station 2 is in contact with station 3a." */
681 // Check if no already done
683 MaskImagePointer binarizedContour;
685 DD("FindAntPostVessels try to get");
686 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
688 catch(clitk::ExceptionObject e) {
693 return binarizedContour;
696 /* Here, we consider the vessels form a kind of anterior barrier. We
697 link all vessels centroids and remove what is post to it.
698 - select the list of structure
699 vessel1 = BrachioCephalicArtery
700 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
701 vessel3 = CommonCarotidArtery
702 vessel4 = SubclavianArtery
705 - crop images as needed
706 - slice by slice, choose the CCL and extract centroids
707 - slice by slice, sort according to polar angle wrt Trachea centroid.
708 - slice by slice, link centroids and close contour
709 - remove outside this contour
714 MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
715 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
716 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
717 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
718 MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
719 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
720 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
722 // Create a temporay support
723 // From first slice of BrachioCephalicVein to end of 3A
724 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
725 MaskImagePointType p;
726 p[0] = p[1] = p[2] = 0.0; // to avoid warning
727 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
729 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
730 GetBackgroundValue(), 2, false, p);
732 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
733 false, GetBackgroundValue());
735 // Resize all structures like support
736 BrachioCephalicArtery =
737 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
738 CommonCarotidArtery =
739 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
741 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
743 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
745 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
746 BrachioCephalicVein =
747 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
749 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
752 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
753 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
754 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
755 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
756 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
757 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
758 std::vector<MaskSlicePointer> slices_SubclavianArtery;
759 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
760 std::vector<MaskSlicePointer> slices_Thyroid;
761 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
762 std::vector<MaskSlicePointer> slices_Aorta;
763 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
764 std::vector<MaskSlicePointer> slices_Trachea;
765 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
766 unsigned int n= slices_BrachioCephalicArtery.size();
768 // Get the boundaries of one slice
769 std::vector<MaskSlicePointType> bounds;
770 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
772 // For all slices, for all structures, find the centroid and build the contour
773 // List of 3D points (for debug)
774 std::vector<MaskImagePointType> p3D;
776 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
777 for(unsigned int i=0; i<n; i++) {
778 // Labelize the slices
779 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
780 GetBackgroundValue(), true, 1);
781 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
782 GetBackgroundValue(), true, 1);
783 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
784 GetBackgroundValue(), true, 1);
785 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
786 GetBackgroundValue(), true, 1);
787 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
788 GetBackgroundValue(), true, 1);
789 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
790 GetBackgroundValue(), true, 1);
793 std::vector<MaskSlicePointType> points2D;
794 std::vector<MaskSlicePointType> centroids1;
795 std::vector<MaskSlicePointType> centroids2;
796 std::vector<MaskSlicePointType> centroids3;
797 std::vector<MaskSlicePointType> centroids4;
798 std::vector<MaskSlicePointType> centroids5;
799 std::vector<MaskSlicePointType> centroids6;
800 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
801 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
802 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
803 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
804 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
805 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
807 // BrachioCephalicVein -> when it is separated into two CCL, we
808 // only consider the most at Right one
809 if (centroids6.size() > 2) {
810 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
811 else centroids6.erase(centroids6.begin()+1);
814 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
815 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
816 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
820 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
821 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
822 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
823 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
824 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
825 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
827 // Sort by angle according to trachea centroid and vertical line,
828 // in polar coordinates :
829 // http://en.wikipedia.org/wiki/Polar_coordinate_system
830 std::vector<MaskSlicePointType> centroids_trachea;
831 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
832 typedef std::pair<MaskSlicePointType, double> PointAngleType;
833 std::vector<PointAngleType> angles;
834 for(unsigned int j=0; j<points2D.size(); j++) {
835 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
836 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
837 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
839 if (x>0) angle = atan(y/x);
840 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
841 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
843 if (y>0) angle = M_PI/2.0;
844 if (y<0) angle = -M_PI/2.0;
847 angle = clitk::rad2deg(angle);
848 // Angle is [-180;180] wrt the X axis. We change the X axis to
849 // be the vertical line, because we want to sort from Right to
850 // Left from Post to Ant.
855 if (angle<0) angle = 360-angle;
858 a.first = points2D[j];
863 // Do nothing if less than 2 points --> n
864 if (points2D.size() < 3) { //continue;
868 // Sort points2D according to polar angles
869 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
870 for(unsigned int j=0; j<angles.size(); j++) {
871 points2D[j] = angles[j].first;
873 // DDV(points2D, points2D.size());
875 /* When vessels are far away, we try to replace the line segment
876 with a curved line that join the two vessels but stay
877 approximately at the same distance from the trachea centroids
881 - let a and c be two successive vessels centroids
882 - id distance(a,c) < threshold, next point
886 - compute mid position between two successive points -
887 compute dist to trachea centroid for the 3 pts - if middle too
890 std::vector<MaskSlicePointType> toadd;
891 unsigned int index = 0;
893 while (index<points2D.size()-1) {
894 MaskSlicePointType a = points2D[index];
895 MaskSlicePointType c = points2D[index+1];
897 double dac = a.EuclideanDistanceTo(c);
900 MaskSlicePointType b;
901 b[0] = a[0]+(c[0]-a[0])/2.0;
902 b[1] = a[1]+(c[1]-a[1])/2.0;
904 // Compute distance to trachea centroid
905 MaskSlicePointType m = centroids_trachea[1];
906 double da = m.EuclideanDistanceTo(a);
907 double db = m.EuclideanDistanceTo(b);
908 //double dc = m.EuclideanDistanceTo(c);
910 // Mean distance, find point on the line from trachea centroid
912 double alpha = (da+db)/2.0;
913 MaskSlicePointType v;
914 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
915 v[0] = (b[0]-m[0])/n;
916 v[1] = (b[1]-m[1])/n;
917 MaskSlicePointType r;
918 r[0] = m[0]+alpha*v[0];
919 r[1] = m[1]+alpha*v[1];
920 points2D.insert(points2D.begin()+index+1, r);
926 // DDV(points2D, points2D.size());
928 // Add some points to close the contour
929 // - H line towards Right
930 MaskSlicePointType p = points2D[0];
932 points2D.insert(points2D.begin(), p);
933 // - corner Right/Post
935 points2D.insert(points2D.begin(), p);
936 // - H line towards Left
939 points2D.push_back(p);
940 // - corner Right/Post
942 points2D.push_back(p);
943 // Close contour with the first point
944 points2D.push_back(points2D[0]);
945 // DDV(points2D, points2D.size());
947 // Build 3D points from the 2D points
948 std::vector<ImagePointType> points3D;
949 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
950 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
952 // Build the mesh from the contour's points
953 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
954 append->AddInput(mesh);
957 // DEBUG: write points3D
958 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
960 // Build the final 3D mesh form the list 2D mesh
962 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
964 // Debug, write the mesh
966 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
968 w->SetFileName("bidon.vtk");
972 // Compute a single binary 3D image from the list of contours
973 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
974 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
975 filter->SetMesh(mesh);
976 filter->SetLikeImage(support);
978 binarizedContour = filter->GetOutput();
981 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
984 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
987 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
988 false, GetBackgroundValue());
990 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
991 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
993 return binarizedContour;
996 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
997 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
998 p[1] += 50; // 50 mm Post from this point must be kept
999 ImageIndexType index;
1000 binarizedContour->TransformPhysicalPointToIndex(p, index);
1001 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1003 // remove from support
1004 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1005 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1006 boolFilter->InPlaceOn();
1007 boolFilter->SetInput1(m_Working_Support);
1008 boolFilter->SetInput2(binarizedContour);
1009 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1010 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1012 boolFilter->SetOperationType(BoolFilterType::And);
1014 boolFilter->SetOperationType(BoolFilterType::AndNot);
1015 boolFilter->Update();
1016 m_Working_Support = boolFilter->GetOutput();
1020 //StopCurrentStep<MaskImageType>(m_Working_Support);
1021 //m_ListOfStations["2R"] = m_Working_Support;
1022 //m_ListOfStations["2L"] = m_Working_Support;
1024 //--------------------------------------------------------------------
1042 //--------------------------------------------------------------------
1043 template <class ImageType>
1044 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1045 clitk::ExtractLymphStationsFilter<ImageType>::
1046 FindAntPostVessels2()
1048 // -----------------------------------------------------
1049 /* Rod says: "The anterior border, as with the Atlas – UM, is
1050 posterior to the vessels (right subclavian vein, left
1051 brachiocephalic vein, right brachiocephalic vein, left subclavian
1052 artery, left common carotid artery and brachiocephalic trunk).
1053 These vessels are not included in the nodal station. The anterior
1054 border is drawn to the midpoint of the vessel and an imaginary
1055 line joins the middle of these vessels. Between the vessels,
1056 station 2 is in contact with station 3a." */
1058 // Check if no already done
1060 MaskImagePointer binarizedContour;
1062 DD("FindAntPostVessels try to get");
1063 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1065 catch(clitk::ExceptionObject e) {
1070 return binarizedContour;
1073 /* Here, we consider the vessels form a kind of anterior barrier. We
1074 link all vessels centroids and remove what is post to it.
1075 - select the list of structure
1076 vessel1 = BrachioCephalicArtery
1077 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1078 vessel3 = CommonCarotidArtery
1079 vessel4 = SubclavianArtery
1082 - crop images as needed
1083 - slice by slice, choose the CCL and extract centroids
1084 - slice by slice, sort according to polar angle wrt Trachea centroid.
1085 - slice by slice, link centroids and close contour
1086 - remove outside this contour
1087 - merge with support
1091 std::map<std::string, MaskImagePointer> MapOfStructures;
1092 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1093 MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1094 MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1095 MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1096 MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1097 MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1098 MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1099 MapOfStructures["Thyroid"] = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1100 MapOfStructures["Aorta"] = GetAFDB()->template GetImage<MaskImageType>("Aorta");
1101 MapOfStructures["Trachea"] = GetAFDB()->template GetImage<MaskImageType>("Trachea");
1103 std::vector<std::string> ListOfStructuresNames;
1105 // Create a temporay support
1106 // From first slice of BrachioCephalicVein to end of 3A
1107 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1108 MaskImagePointType p;
1109 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1110 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1111 GetBackgroundValue(), 2, true, p);
1113 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1114 GetBackgroundValue(), 2, false, p);
1115 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1116 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1117 false, GetBackgroundValue());
1119 // Resize all structures like support
1120 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1121 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1125 typedef std::vector<MaskSlicePointer> SlicesType;
1126 std::map<std::string, SlicesType> MapOfSlices;
1127 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1129 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1130 MapOfSlices[it->first] = s;
1133 unsigned int n= MapOfSlices["Trachea"].size();
1135 // Get the boundaries of one slice
1136 std::vector<MaskSlicePointType> bounds;
1137 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1140 // For all slices, for all structures, find the centroid and build the contour
1141 // List of 3D points (for debug)
1142 std::vector<MaskImagePointType> p3D;
1143 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1144 for(unsigned int i=0; i<n; i++) {
1146 // Labelize the slices
1147 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1148 MaskSlicePointer & s = MapOfSlices[it->first][i];
1149 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1153 std::vector<MaskSlicePointType> points2D;
1154 typedef std::vector<MaskSlicePointType> CentroidsType;
1155 std::map<std::string, CentroidsType> MapOfCentroids;
1156 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1157 std::string structure = it->first;
1158 MaskSlicePointer & s = MapOfSlices[structure][i];
1160 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1161 MapOfCentroids[structure] = c;
1165 // BrachioCephalicVein -> when it is separated into two CCL, we
1166 // only consider the most at Right one
1167 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1169 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1170 else c.erase(c.begin()+1);
1174 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1175 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1176 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1177 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1178 MapOfCentroids["BrachioCephalicVein"].clear();
1182 // Add all 2D points
1183 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1184 std::string structure = it->first;
1185 if (structure != "Trachea") { // do not add centroid of trachea
1186 CentroidsType & c = MapOfCentroids[structure];
1187 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1191 // Sort by angle according to trachea centroid and vertical line,
1192 // in polar coordinates :
1193 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1194 // std::vector<MaskSlicePointType> centroids_trachea;
1195 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1196 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1197 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1198 std::vector<PointAngleType> angles;
1199 for(unsigned int j=0; j<points2D.size(); j++) {
1200 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1201 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1202 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1204 if (x>0) angle = atan(y/x);
1205 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1206 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1208 if (y>0) angle = M_PI/2.0;
1209 if (y<0) angle = -M_PI/2.0;
1210 if (y==0) angle = 0;
1212 angle = clitk::rad2deg(angle);
1213 // Angle is [-180;180] wrt the X axis. We change the X axis to
1214 // be the vertical line, because we want to sort from Right to
1215 // Left from Post to Ant.
1217 angle = (270-angle);
1220 if (angle<0) angle = 360-angle;
1223 a.first = points2D[j];
1225 angles.push_back(a);
1228 // Do nothing if less than 2 points --> n
1229 if (points2D.size() < 3) { //continue;
1233 // Sort points2D according to polar angles
1234 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1235 for(unsigned int j=0; j<angles.size(); j++) {
1236 points2D[j] = angles[j].first;
1238 // DDV(points2D, points2D.size());
1240 /* When vessels are far away, we try to replace the line segment
1241 with a curved line that join the two vessels but stay
1242 approximately at the same distance from the trachea centroids
1246 - let a and c be two successive vessels centroids
1247 - id distance(a,c) < threshold, next point
1251 - compute mid position between two successive points -
1252 compute dist to trachea centroid for the 3 pts - if middle too
1255 std::vector<MaskSlicePointType> toadd;
1256 unsigned int index = 0;
1258 while (index<points2D.size()-1) {
1259 MaskSlicePointType a = points2D[index];
1260 MaskSlicePointType c = points2D[index+1];
1262 double dac = a.EuclideanDistanceTo(c);
1265 MaskSlicePointType b;
1266 b[0] = a[0]+(c[0]-a[0])/2.0;
1267 b[1] = a[1]+(c[1]-a[1])/2.0;
1269 // Compute distance to trachea centroid
1270 MaskSlicePointType m = centroids_trachea[1];
1271 double da = m.EuclideanDistanceTo(a);
1272 double db = m.EuclideanDistanceTo(b);
1273 //double dc = m.EuclideanDistanceTo(c);
1275 // Mean distance, find point on the line from trachea centroid
1277 double alpha = (da+db)/2.0;
1278 MaskSlicePointType v;
1279 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1280 v[0] = (b[0]-m[0])/n;
1281 v[1] = (b[1]-m[1])/n;
1282 MaskSlicePointType r;
1283 r[0] = m[0]+alpha*v[0];
1284 r[1] = m[1]+alpha*v[1];
1285 points2D.insert(points2D.begin()+index+1, r);
1291 // DDV(points2D, points2D.size());
1293 // Add some points to close the contour
1294 // - H line towards Right
1295 MaskSlicePointType p = points2D[0];
1296 p[0] = bounds[0][0];
1297 points2D.insert(points2D.begin(), p);
1298 // - corner Right/Post
1300 points2D.insert(points2D.begin(), p);
1301 // - H line towards Left
1302 p = points2D.back();
1303 p[0] = bounds[2][0];
1304 points2D.push_back(p);
1305 // - corner Right/Post
1307 points2D.push_back(p);
1308 // Close contour with the first point
1309 points2D.push_back(points2D[0]);
1310 // DDV(points2D, points2D.size());
1312 // Build 3D points from the 2D points
1313 std::vector<ImagePointType> points3D;
1314 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1315 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1317 // Build the mesh from the contour's points
1318 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1319 append->AddInput(mesh);
1320 // if (i ==n-1) { // last slice
1321 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1322 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1323 // append->AddInput(mesh);
1327 // DEBUG: write points3D
1328 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1330 // Build the final 3D mesh form the list 2D mesh
1332 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1334 // Debug, write the mesh
1336 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1338 w->SetFileName("bidon.vtk");
1342 // Compute a single binary 3D image from the list of contours
1343 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1344 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1345 filter->SetMesh(mesh);
1346 filter->SetLikeImage(support);
1348 binarizedContour = filter->GetOutput();
1351 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1352 GetForegroundValue(), 2, true, p);
1354 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1355 GetForegroundValue(), 2, false, p);
1356 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1357 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1358 false, GetBackgroundValue());
1361 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1362 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1364 return binarizedContour;
1367 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1368 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1369 p[1] += 50; // 50 mm Post from this point must be kept
1370 ImageIndexType index;
1371 binarizedContour->TransformPhysicalPointToIndex(p, index);
1372 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1374 // remove from support
1375 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1376 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1377 boolFilter->InPlaceOn();
1378 boolFilter->SetInput1(m_Working_Support);
1379 boolFilter->SetInput2(binarizedContour);
1380 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1381 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1383 boolFilter->SetOperationType(BoolFilterType::And);
1385 boolFilter->SetOperationType(BoolFilterType::AndNot);
1386 boolFilter->Update();
1387 m_Working_Support = boolFilter->GetOutput();
1391 //StopCurrentStep<MaskImageType>(m_Working_Support);
1392 //m_ListOfStations["2R"] = m_Working_Support;
1393 //m_ListOfStations["2L"] = m_Working_Support;
1395 //--------------------------------------------------------------------
1399 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX