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");
131 // Extract Station2RL
132 StartNewStep("Station 2RL");
134 ExtractStation_2RL();
138 StartNewStep("Station 7");
143 if (0) { // temporary suppress
144 // Extract Station4RL
145 StartNewStep("Station 4RL");
147 //ExtractStation_4RL();
151 //--------------------------------------------------------------------
154 //--------------------------------------------------------------------
155 template <class TImageType>
157 clitk::ExtractLymphStationsFilter<TImageType>::
158 GenerateInputRequestedRegion() {
159 //DD("GenerateInputRequestedRegion (nothing?)");
161 //--------------------------------------------------------------------
164 //--------------------------------------------------------------------
165 template <class TImageType>
167 clitk::ExtractLymphStationsFilter<TImageType>::
169 // Final Step -> graft output (if SetNthOutput => redo)
170 this->GraftOutput(m_ListOfStations["8"]);
172 //--------------------------------------------------------------------
175 //--------------------------------------------------------------------
176 template <class TImageType>
178 clitk::ExtractLymphStationsFilter<TImageType>::
179 CheckForStation(std::string station)
181 // Compute Station name
182 std::string s = "Station"+station;
185 // Check if station already exist in DB
187 if (GetAFDB()->TagExist(s)) {
188 m_ListOfStations[station] = GetAFDB()->template GetImage<MaskImageType>(s);
192 // Define the starting support
193 if (found && GetComputeStation(station)) {
194 std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
196 if (!found || GetComputeStation(station)) {
197 m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
201 std::cout << "Station " << station << " found. I used it" << std::endl;
205 //--------------------------------------------------------------------
208 //--------------------------------------------------------------------
209 template <class TImageType>
211 clitk::ExtractLymphStationsFilter<TImageType>::
212 GetComputeStation(std::string station)
214 return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
216 //--------------------------------------------------------------------
219 //--------------------------------------------------------------------
220 template <class TImageType>
222 clitk::ExtractLymphStationsFilter<TImageType>::
223 AddComputeStation(std::string station)
225 m_ComputeStationMap[station] = true;
227 //--------------------------------------------------------------------
231 //--------------------------------------------------------------------
232 template<class PointType>
233 class comparePointsX {
235 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
237 //--------------------------------------------------------------------
240 //--------------------------------------------------------------------
241 template<class PairType>
242 class comparePointsWithAngle {
244 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
246 //--------------------------------------------------------------------
249 //--------------------------------------------------------------------
251 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
252 std::vector<itk::Point<double, Dim-1> > previous;
253 HypercubeCorners<Dim-1>(previous);
254 out.resize(previous.size()*2);
255 for(unsigned int i=0; i<out.size(); i++) {
256 itk::Point<double, Dim> p;
257 if (i<previous.size()) p[0] = 0;
259 for(int j=0; j<Dim-1; j++)
261 p[j+1] = previous[i%previous.size()][j];
266 //--------------------------------------------------------------------
269 //--------------------------------------------------------------------
271 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
276 //--------------------------------------------------------------------
279 //--------------------------------------------------------------------
280 template<class ImageType>
281 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
282 std::vector<typename ImageType::PointType> & bounds)
284 // Get image max/min coordinates
285 const unsigned int dim=ImageType::ImageDimension;
286 typedef typename ImageType::PointType PointType;
287 typedef typename ImageType::IndexType IndexType;
288 PointType min_c, max_c;
289 IndexType min_i, max_i;
290 min_i = image->GetLargestPossibleRegion().GetIndex();
291 for(unsigned int i=0; i<dim; i++)
292 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
293 image->TransformIndexToPhysicalPoint(min_i, min_c);
294 image->TransformIndexToPhysicalPoint(max_i, max_c);
296 // Get corners coordinates
297 HypercubeCorners<ImageType::ImageDimension>(bounds);
298 for(unsigned int i=0; i<bounds.size(); i++) {
299 for(unsigned int j=0; j<dim; j++) {
300 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
301 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
305 //--------------------------------------------------------------------
308 //--------------------------------------------------------------------
309 template <class TImageType>
311 clitk::ExtractLymphStationsFilter<TImageType>::
312 ExtractStation_4RL() {
316 WARNING ONLY 4R FIRST !!! (not same inf limits)
318 ExtractStation_4RL_SI_Limits();
319 ExtractStation_4RL_LR_Limits();
320 ExtractStation_4RL_AP_Limits();
322 //--------------------------------------------------------------------
325 //--------------------------------------------------------------------
326 template <class ImageType>
328 clitk::ExtractLymphStationsFilter<ImageType>::
329 Remove_Structures(std::string station, std::string s)
332 StartNewStep("[Station"+station+"] Remove "+s);
333 MaskImagePointer Structure = GetAFDB()->template GetImage<MaskImageType>(s);
334 clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
336 catch(clitk::ExceptionObject e) {
337 std::cout << s << " not found, skip." << std::endl;
340 //--------------------------------------------------------------------
343 //--------------------------------------------------------------------
344 template <class TImageType>
346 clitk::ExtractLymphStationsFilter<TImageType>::
347 SetFuzzyThreshold(std::string station, std::string tag, double value)
349 m_FuzzyThreshold[station][tag] = value;
351 //--------------------------------------------------------------------
354 //--------------------------------------------------------------------
355 template <class TImageType>
357 clitk::ExtractLymphStationsFilter<TImageType>::
358 SetThreshold(std::string station, std::string tag, double value)
360 m_Threshold[station][tag] = value;
362 //--------------------------------------------------------------------
365 //--------------------------------------------------------------------
366 template <class TImageType>
368 clitk::ExtractLymphStationsFilter<TImageType>::
369 GetFuzzyThreshold(std::string station, std::string tag)
371 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
372 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
376 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
377 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
381 return m_FuzzyThreshold[station][tag];
383 //--------------------------------------------------------------------
386 //--------------------------------------------------------------------
387 template <class TImageType>
389 clitk::ExtractLymphStationsFilter<TImageType>::
390 GetThreshold(std::string station, std::string tag)
392 if (m_Threshold.find(station) == m_Threshold.end()) {
393 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
397 if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
398 clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
402 return m_Threshold[station][tag];
404 //--------------------------------------------------------------------
407 //--------------------------------------------------------------------
408 template <class TImageType>
410 clitk::ExtractLymphStationsFilter<TImageType>::
411 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
413 // Create line from A to B with
414 // A = upper border of LLL at left
415 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
418 GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
419 GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
421 catch(clitk::ExceptionObject & o) {
423 DD("FindLineForS7S8Separation");
424 // Load LeftLowerLobeBronchus and get centroid point
425 MaskImagePointer LeftLowerLobeBronchus =
426 GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
427 std::vector<MaskImagePointType> c;
428 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
431 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
432 MaskImagePointer RightMiddleLobeBronchus =
433 GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
434 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
435 GetBackgroundValue(),
438 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
441 // Insert into the DB
442 GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
443 GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
446 //--------------------------------------------------------------------
449 //--------------------------------------------------------------------
450 template <class TImageType>
452 clitk::ExtractLymphStationsFilter<TImageType>::
457 z = GetAFDB()->GetDouble("CarinaZ");
459 catch(clitk::ExceptionObject e) {
460 DD("FindCarinaSlicePosition");
462 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
464 // Get Centroid and Z value
465 std::vector<MaskImagePointType> centroids;
466 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
468 // We dont need Carina structure from now
469 GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
471 // Put inside the AFDB
472 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
473 GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
479 //--------------------------------------------------------------------
482 //--------------------------------------------------------------------
483 template <class TImageType>
485 clitk::ExtractLymphStationsFilter<TImageType>::
490 z = GetAFDB()->GetDouble("ApexOfTheChestZ");
492 catch(clitk::ExceptionObject e) {
493 DD("FindApexOfTheChestPosition");
494 MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
495 MaskImagePointType p;
496 p[0] = p[1] = p[2] = 0.0; // to avoid warning
497 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
499 // We dont need Lungs structure from now
500 GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
502 // Put inside the AFDB
503 GetAFDB()->SetPoint3D("ApexOfTheChest", p);
504 p[2] -= 5; // We consider 5 mm lower
505 GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
511 //--------------------------------------------------------------------
514 //--------------------------------------------------------------------
515 template <class TImageType>
517 clitk::ExtractLymphStationsFilter<TImageType>::
518 FindLeftAndRightBronchi()
521 m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
522 m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
524 catch(clitk::ExceptionObject & o) {
526 DD("FindLeftAndRightBronchi");
527 // The goal is to separate the trachea inferiorly to the carina into
528 // a Left and Right bronchus.
531 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
533 // Get the Carina position
534 double m_CarinaZ = FindCarina();
536 // Consider only inferiorly to the Carina
537 MaskImagePointer m_Working_Trachea =
538 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
539 GetBackgroundValue());
541 // Labelize the trachea
542 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
544 // Carina position must at the first slice that separate the two
545 // main bronchus (not superiorly). We thus first check that the
546 // upper slice is composed of at least two labels
547 MaskImagePointer RightBronchus;
548 MaskImagePointer LeftBronchus;
549 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
550 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
551 iter.SetFirstDirection(0);
552 iter.SetSecondDirection(1);
553 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
555 while (!iter.IsAtReverseEndOfSlice()) {
556 while (!iter.IsAtReverseEndOfLine()) {
557 if (iter.Get() > maxLabel) maxLabel = iter.Get();
563 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
566 // Compute 3D centroids of both parts to identify the left from the
568 std::vector<ImagePointType> c;
569 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
570 ImagePointType C1 = c[1];
571 ImagePointType C2 = c[2];
573 ImagePixelType rightLabel;
574 ImagePixelType leftLabel;
575 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
576 else { rightLabel = 2; leftLabel = 1; }
578 // Select LeftLabel (set one label to Backgroundvalue)
580 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
581 GetBackgroundValue(), GetForegroundValue());
583 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
584 leftLabel, GetBackgroundValue(), false);
586 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
587 GetBackgroundValue(), GetForegroundValue());
589 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
590 rightLabel, GetBackgroundValue(), false);
594 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
595 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
597 // Insert int AFDB if need after
598 GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
599 RightBronchus, true);
600 GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
604 //--------------------------------------------------------------------
607 //--------------------------------------------------------------------
608 template <class TImageType>
610 clitk::ExtractLymphStationsFilter<TImageType>::
611 FindSuperiorBorderOfAorticArch()
615 z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
617 catch(clitk::ExceptionObject e) {
618 DD("FindSuperiorBorderOfAorticArch");
619 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
620 MaskImagePointType p;
621 p[0] = p[1] = p[2] = 0.0; // to avoid warning
622 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
623 p[2] += Aorta->GetSpacing()[2]; // the slice above
625 // We dont need Lungs structure from now
626 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
628 // Put inside the AFDB
629 GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
630 GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
636 //--------------------------------------------------------------------
639 //--------------------------------------------------------------------
640 template <class TImageType>
642 clitk::ExtractLymphStationsFilter<TImageType>::
643 FindInferiorBorderOfAorticArch()
647 z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
649 catch(clitk::ExceptionObject e) {
650 DD("FindInferiorBorderOfAorticArch");
651 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
652 std::vector<MaskSlicePointer> slices;
653 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
655 uint i = slices.size()-1;
657 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
658 std::vector<typename MaskSliceType::PointType> c;
659 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
667 MaskImageIndexType index;
668 index[0] = index[1] = 0.0;
670 MaskImagePointType lower;
671 Aorta->TransformIndexToPhysicalPoint(index, lower);
673 // We dont need Lungs structure from now
674 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
676 // Put inside the AFDB
677 GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
678 GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
684 //--------------------------------------------------------------------
687 //--------------------------------------------------------------------
688 template <class ImageType>
689 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
690 clitk::ExtractLymphStationsFilter<ImageType>::
693 // -----------------------------------------------------
694 /* Rod says: "The anterior border, as with the Atlas – UM, is
695 posterior to the vessels (right subclavian vein, left
696 brachiocephalic vein, right brachiocephalic vein, left subclavian
697 artery, left common carotid artery and brachiocephalic trunk).
698 These vessels are not included in the nodal station. The anterior
699 border is drawn to the midpoint of the vessel and an imaginary
700 line joins the middle of these vessels. Between the vessels,
701 station 2 is in contact with station 3a." */
703 // Check if no already done
705 MaskImagePointer binarizedContour;
707 DD("FindAntPostVessels try to get");
708 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
710 catch(clitk::ExceptionObject e) {
715 return binarizedContour;
718 /* Here, we consider the vessels form a kind of anterior barrier. We
719 link all vessels centroids and remove what is post to it.
720 - select the list of structure
721 vessel1 = BrachioCephalicArtery
722 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
723 vessel3 = CommonCarotidArtery
724 vessel4 = SubclavianArtery
727 - crop images as needed
728 - slice by slice, choose the CCL and extract centroids
729 - slice by slice, sort according to polar angle wrt Trachea centroid.
730 - slice by slice, link centroids and close contour
731 - remove outside this contour
736 MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
737 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
738 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
739 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
740 MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
741 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
742 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
744 // Create a temporay support
745 // From first slice of BrachioCephalicVein to end of 3A
746 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
747 MaskImagePointType p;
748 p[0] = p[1] = p[2] = 0.0; // to avoid warning
749 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
751 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
752 GetBackgroundValue(), 2, false, p);
754 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
755 false, GetBackgroundValue());
757 // Resize all structures like support
758 BrachioCephalicArtery =
759 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
760 CommonCarotidArtery =
761 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
763 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
765 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
767 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
768 BrachioCephalicVein =
769 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
771 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
774 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
775 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
776 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
777 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
778 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
779 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
780 std::vector<MaskSlicePointer> slices_SubclavianArtery;
781 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
782 std::vector<MaskSlicePointer> slices_Thyroid;
783 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
784 std::vector<MaskSlicePointer> slices_Aorta;
785 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
786 std::vector<MaskSlicePointer> slices_Trachea;
787 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
788 unsigned int n= slices_BrachioCephalicArtery.size();
790 // Get the boundaries of one slice
791 std::vector<MaskSlicePointType> bounds;
792 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
794 // For all slices, for all structures, find the centroid and build the contour
795 // List of 3D points (for debug)
796 std::vector<MaskImagePointType> p3D;
798 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
799 for(unsigned int i=0; i<n; i++) {
800 // Labelize the slices
801 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
802 GetBackgroundValue(), true, 1);
803 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
804 GetBackgroundValue(), true, 1);
805 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
806 GetBackgroundValue(), true, 1);
807 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
808 GetBackgroundValue(), true, 1);
809 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
810 GetBackgroundValue(), true, 1);
811 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
812 GetBackgroundValue(), true, 1);
815 std::vector<MaskSlicePointType> points2D;
816 std::vector<MaskSlicePointType> centroids1;
817 std::vector<MaskSlicePointType> centroids2;
818 std::vector<MaskSlicePointType> centroids3;
819 std::vector<MaskSlicePointType> centroids4;
820 std::vector<MaskSlicePointType> centroids5;
821 std::vector<MaskSlicePointType> centroids6;
822 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
823 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
824 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
825 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
826 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
827 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
829 // BrachioCephalicVein -> when it is separated into two CCL, we
830 // only consider the most at Right one
831 if (centroids6.size() > 2) {
832 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
833 else centroids6.erase(centroids6.begin()+1);
836 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
837 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
838 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
842 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
843 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
844 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
845 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
846 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
847 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
849 // Sort by angle according to trachea centroid and vertical line,
850 // in polar coordinates :
851 // http://en.wikipedia.org/wiki/Polar_coordinate_system
852 std::vector<MaskSlicePointType> centroids_trachea;
853 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
854 typedef std::pair<MaskSlicePointType, double> PointAngleType;
855 std::vector<PointAngleType> angles;
856 for(unsigned int j=0; j<points2D.size(); j++) {
857 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
858 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
859 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
861 if (x>0) angle = atan(y/x);
862 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
863 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
865 if (y>0) angle = M_PI/2.0;
866 if (y<0) angle = -M_PI/2.0;
869 angle = clitk::rad2deg(angle);
870 // Angle is [-180;180] wrt the X axis. We change the X axis to
871 // be the vertical line, because we want to sort from Right to
872 // Left from Post to Ant.
877 if (angle<0) angle = 360-angle;
880 a.first = points2D[j];
885 // Do nothing if less than 2 points --> n
886 if (points2D.size() < 3) { //continue;
890 // Sort points2D according to polar angles
891 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
892 for(unsigned int j=0; j<angles.size(); j++) {
893 points2D[j] = angles[j].first;
895 // DDV(points2D, points2D.size());
897 /* When vessels are far away, we try to replace the line segment
898 with a curved line that join the two vessels but stay
899 approximately at the same distance from the trachea centroids
903 - let a and c be two successive vessels centroids
904 - id distance(a,c) < threshold, next point
908 - compute mid position between two successive points -
909 compute dist to trachea centroid for the 3 pts - if middle too
912 std::vector<MaskSlicePointType> toadd;
913 unsigned int index = 0;
915 while (index<points2D.size()-1) {
916 MaskSlicePointType a = points2D[index];
917 MaskSlicePointType c = points2D[index+1];
919 double dac = a.EuclideanDistanceTo(c);
922 MaskSlicePointType b;
923 b[0] = a[0]+(c[0]-a[0])/2.0;
924 b[1] = a[1]+(c[1]-a[1])/2.0;
926 // Compute distance to trachea centroid
927 MaskSlicePointType m = centroids_trachea[1];
928 double da = m.EuclideanDistanceTo(a);
929 double db = m.EuclideanDistanceTo(b);
930 //double dc = m.EuclideanDistanceTo(c);
932 // Mean distance, find point on the line from trachea centroid
934 double alpha = (da+db)/2.0;
935 MaskSlicePointType v;
936 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
937 v[0] = (b[0]-m[0])/n;
938 v[1] = (b[1]-m[1])/n;
939 MaskSlicePointType r;
940 r[0] = m[0]+alpha*v[0];
941 r[1] = m[1]+alpha*v[1];
942 points2D.insert(points2D.begin()+index+1, r);
948 // DDV(points2D, points2D.size());
950 // Add some points to close the contour
951 // - H line towards Right
952 MaskSlicePointType p = points2D[0];
954 points2D.insert(points2D.begin(), p);
955 // - corner Right/Post
957 points2D.insert(points2D.begin(), p);
958 // - H line towards Left
961 points2D.push_back(p);
962 // - corner Right/Post
964 points2D.push_back(p);
965 // Close contour with the first point
966 points2D.push_back(points2D[0]);
967 // DDV(points2D, points2D.size());
969 // Build 3D points from the 2D points
970 std::vector<ImagePointType> points3D;
971 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
972 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
974 // Build the mesh from the contour's points
975 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
976 append->AddInput(mesh);
979 // DEBUG: write points3D
980 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
982 // Build the final 3D mesh form the list 2D mesh
984 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
986 // Debug, write the mesh
988 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
990 w->SetFileName("bidon.vtk");
994 // Compute a single binary 3D image from the list of contours
995 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
996 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
997 filter->SetMesh(mesh);
998 filter->SetLikeImage(support);
1000 binarizedContour = filter->GetOutput();
1003 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1006 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1009 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1010 false, GetBackgroundValue());
1012 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1013 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1015 return binarizedContour;
1018 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1019 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1020 p[1] += 50; // 50 mm Post from this point must be kept
1021 ImageIndexType index;
1022 binarizedContour->TransformPhysicalPointToIndex(p, index);
1023 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1025 // remove from support
1026 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1027 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1028 boolFilter->InPlaceOn();
1029 boolFilter->SetInput1(m_Working_Support);
1030 boolFilter->SetInput2(binarizedContour);
1031 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1032 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1034 boolFilter->SetOperationType(BoolFilterType::And);
1036 boolFilter->SetOperationType(BoolFilterType::AndNot);
1037 boolFilter->Update();
1038 m_Working_Support = boolFilter->GetOutput();
1042 //StopCurrentStep<MaskImageType>(m_Working_Support);
1043 //m_ListOfStations["2R"] = m_Working_Support;
1044 //m_ListOfStations["2L"] = m_Working_Support;
1046 //--------------------------------------------------------------------
1064 //--------------------------------------------------------------------
1065 template <class ImageType>
1066 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1067 clitk::ExtractLymphStationsFilter<ImageType>::
1068 FindAntPostVessels2()
1070 // -----------------------------------------------------
1071 /* Rod says: "The anterior border, as with the Atlas – UM, is
1072 posterior to the vessels (right subclavian vein, left
1073 brachiocephalic vein, right brachiocephalic vein, left subclavian
1074 artery, left common carotid artery and brachiocephalic trunk).
1075 These vessels are not included in the nodal station. The anterior
1076 border is drawn to the midpoint of the vessel and an imaginary
1077 line joins the middle of these vessels. Between the vessels,
1078 station 2 is in contact with station 3a." */
1080 // Check if no already done
1082 MaskImagePointer binarizedContour;
1084 DD("FindAntPostVessels try to get");
1085 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1087 catch(clitk::ExceptionObject e) {
1092 return binarizedContour;
1095 /* Here, we consider the vessels form a kind of anterior barrier. We
1096 link all vessels centroids and remove what is post to it.
1097 - select the list of structure
1098 vessel1 = BrachioCephalicArtery
1099 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1100 vessel3 = CommonCarotidArtery
1101 vessel4 = SubclavianArtery
1104 - crop images as needed
1105 - slice by slice, choose the CCL and extract centroids
1106 - slice by slice, sort according to polar angle wrt Trachea centroid.
1107 - slice by slice, link centroids and close contour
1108 - remove outside this contour
1109 - merge with support
1113 std::map<std::string, MaskImagePointer> MapOfStructures;
1114 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1115 MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1116 MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1117 MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1118 MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1119 MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1120 MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1121 MapOfStructures["Thyroid"] = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1122 MapOfStructures["Aorta"] = GetAFDB()->template GetImage<MaskImageType>("Aorta");
1123 MapOfStructures["Trachea"] = GetAFDB()->template GetImage<MaskImageType>("Trachea");
1125 std::vector<std::string> ListOfStructuresNames;
1127 // Create a temporay support
1128 // From first slice of BrachioCephalicVein to end of 3A
1129 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1130 MaskImagePointType p;
1131 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1132 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1133 GetBackgroundValue(), 2, true, p);
1135 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1136 GetBackgroundValue(), 2, false, p);
1137 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1138 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1139 false, GetBackgroundValue());
1141 // Resize all structures like support
1142 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1143 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1147 typedef std::vector<MaskSlicePointer> SlicesType;
1148 std::map<std::string, SlicesType> MapOfSlices;
1149 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1151 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1152 MapOfSlices[it->first] = s;
1155 unsigned int n= MapOfSlices["Trachea"].size();
1157 // Get the boundaries of one slice
1158 std::vector<MaskSlicePointType> bounds;
1159 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1162 // For all slices, for all structures, find the centroid and build the contour
1163 // List of 3D points (for debug)
1164 std::vector<MaskImagePointType> p3D;
1165 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1166 for(unsigned int i=0; i<n; i++) {
1168 // Labelize the slices
1169 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1170 MaskSlicePointer & s = MapOfSlices[it->first][i];
1171 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1175 std::vector<MaskSlicePointType> points2D;
1176 typedef std::vector<MaskSlicePointType> CentroidsType;
1177 std::map<std::string, CentroidsType> MapOfCentroids;
1178 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1179 std::string structure = it->first;
1180 MaskSlicePointer & s = MapOfSlices[structure][i];
1182 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1183 MapOfCentroids[structure] = c;
1187 // BrachioCephalicVein -> when it is separated into two CCL, we
1188 // only consider the most at Right one
1189 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1191 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1192 else c.erase(c.begin()+1);
1196 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1197 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1198 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1199 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1200 MapOfCentroids["BrachioCephalicVein"].clear();
1204 // Add all 2D points
1205 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1206 std::string structure = it->first;
1207 if (structure != "Trachea") { // do not add centroid of trachea
1208 CentroidsType & c = MapOfCentroids[structure];
1209 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1213 // Sort by angle according to trachea centroid and vertical line,
1214 // in polar coordinates :
1215 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1216 // std::vector<MaskSlicePointType> centroids_trachea;
1217 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1218 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1219 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1220 std::vector<PointAngleType> angles;
1221 for(unsigned int j=0; j<points2D.size(); j++) {
1222 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1223 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1224 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1226 if (x>0) angle = atan(y/x);
1227 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1228 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1230 if (y>0) angle = M_PI/2.0;
1231 if (y<0) angle = -M_PI/2.0;
1232 if (y==0) angle = 0;
1234 angle = clitk::rad2deg(angle);
1235 // Angle is [-180;180] wrt the X axis. We change the X axis to
1236 // be the vertical line, because we want to sort from Right to
1237 // Left from Post to Ant.
1239 angle = (270-angle);
1242 if (angle<0) angle = 360-angle;
1245 a.first = points2D[j];
1247 angles.push_back(a);
1250 // Do nothing if less than 2 points --> n
1251 if (points2D.size() < 3) { //continue;
1255 // Sort points2D according to polar angles
1256 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1257 for(unsigned int j=0; j<angles.size(); j++) {
1258 points2D[j] = angles[j].first;
1260 // DDV(points2D, points2D.size());
1262 /* When vessels are far away, we try to replace the line segment
1263 with a curved line that join the two vessels but stay
1264 approximately at the same distance from the trachea centroids
1268 - let a and c be two successive vessels centroids
1269 - id distance(a,c) < threshold, next point
1273 - compute mid position between two successive points -
1274 compute dist to trachea centroid for the 3 pts - if middle too
1277 std::vector<MaskSlicePointType> toadd;
1278 unsigned int index = 0;
1280 while (index<points2D.size()-1) {
1281 MaskSlicePointType a = points2D[index];
1282 MaskSlicePointType c = points2D[index+1];
1284 double dac = a.EuclideanDistanceTo(c);
1287 MaskSlicePointType b;
1288 b[0] = a[0]+(c[0]-a[0])/2.0;
1289 b[1] = a[1]+(c[1]-a[1])/2.0;
1291 // Compute distance to trachea centroid
1292 MaskSlicePointType m = centroids_trachea[1];
1293 double da = m.EuclideanDistanceTo(a);
1294 double db = m.EuclideanDistanceTo(b);
1295 //double dc = m.EuclideanDistanceTo(c);
1297 // Mean distance, find point on the line from trachea centroid
1299 double alpha = (da+db)/2.0;
1300 MaskSlicePointType v;
1301 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1302 v[0] = (b[0]-m[0])/n;
1303 v[1] = (b[1]-m[1])/n;
1304 MaskSlicePointType r;
1305 r[0] = m[0]+alpha*v[0];
1306 r[1] = m[1]+alpha*v[1];
1307 points2D.insert(points2D.begin()+index+1, r);
1313 // DDV(points2D, points2D.size());
1315 // Add some points to close the contour
1316 // - H line towards Right
1317 MaskSlicePointType p = points2D[0];
1318 p[0] = bounds[0][0];
1319 points2D.insert(points2D.begin(), p);
1320 // - corner Right/Post
1322 points2D.insert(points2D.begin(), p);
1323 // - H line towards Left
1324 p = points2D.back();
1325 p[0] = bounds[2][0];
1326 points2D.push_back(p);
1327 // - corner Right/Post
1329 points2D.push_back(p);
1330 // Close contour with the first point
1331 points2D.push_back(points2D[0]);
1332 // DDV(points2D, points2D.size());
1334 // Build 3D points from the 2D points
1335 std::vector<ImagePointType> points3D;
1336 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1337 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1339 // Build the mesh from the contour's points
1340 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1341 append->AddInput(mesh);
1342 // if (i ==n-1) { // last slice
1343 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1344 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1345 // append->AddInput(mesh);
1349 // DEBUG: write points3D
1350 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1352 // Build the final 3D mesh form the list 2D mesh
1354 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1356 // Debug, write the mesh
1358 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1360 w->SetFileName("bidon.vtk");
1364 // Compute a single binary 3D image from the list of contours
1365 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1366 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1367 filter->SetMesh(mesh);
1368 filter->SetLikeImage(support);
1370 binarizedContour = filter->GetOutput();
1373 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1374 GetForegroundValue(), 2, true, p);
1376 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1377 GetForegroundValue(), 2, false, p);
1378 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1379 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1380 false, GetBackgroundValue());
1383 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1384 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1386 return binarizedContour;
1389 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1390 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1391 p[1] += 50; // 50 mm Post from this point must be kept
1392 ImageIndexType index;
1393 binarizedContour->TransformPhysicalPointToIndex(p, index);
1394 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1396 // remove from support
1397 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1398 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1399 boolFilter->InPlaceOn();
1400 boolFilter->SetInput1(m_Working_Support);
1401 boolFilter->SetInput2(binarizedContour);
1402 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1403 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1405 boolFilter->SetOperationType(BoolFilterType::And);
1407 boolFilter->SetOperationType(BoolFilterType::AndNot);
1408 boolFilter->Update();
1409 m_Working_Support = boolFilter->GetOutput();
1413 //StopCurrentStep<MaskImageType>(m_Working_Support);
1414 //m_ListOfStations["2R"] = m_Working_Support;
1415 //m_ListOfStations["2L"] = m_Working_Support;
1417 //--------------------------------------------------------------------
1421 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX