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 GetFuzzyThreshold(std::string station, std::string tag)
360 if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
361 clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
365 if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
366 clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
370 return m_FuzzyThreshold[station][tag];
372 //--------------------------------------------------------------------
375 //--------------------------------------------------------------------
376 template <class TImageType>
378 clitk::ExtractLymphStationsFilter<TImageType>::
379 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
381 // Create line from A to B with
382 // A = upper border of LLL at left
383 // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
386 GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
387 GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
389 catch(clitk::ExceptionObject & o) {
391 DD("FindLineForS7S8Separation");
392 // Load LeftLowerLobeBronchus and get centroid point
393 MaskImagePointer LeftLowerLobeBronchus =
394 GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
395 std::vector<MaskImagePointType> c;
396 clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
399 // Load RightMiddleLobeBronchus and get superior point (not centroid here)
400 MaskImagePointer RightMiddleLobeBronchus =
401 GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
402 bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
403 GetBackgroundValue(),
406 clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
409 // Insert into the DB
410 GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
411 GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
414 //--------------------------------------------------------------------
417 //--------------------------------------------------------------------
418 template <class TImageType>
420 clitk::ExtractLymphStationsFilter<TImageType>::
425 z = GetAFDB()->GetDouble("CarinaZ");
427 catch(clitk::ExceptionObject e) {
428 DD("FindCarinaSlicePosition");
430 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
432 // Get Centroid and Z value
433 std::vector<MaskImagePointType> centroids;
434 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
436 // We dont need Carina structure from now
437 GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
439 // Put inside the AFDB
440 GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
441 GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
447 //--------------------------------------------------------------------
450 //--------------------------------------------------------------------
451 template <class TImageType>
453 clitk::ExtractLymphStationsFilter<TImageType>::
458 z = GetAFDB()->GetDouble("ApexOfTheChestZ");
460 catch(clitk::ExceptionObject e) {
461 DD("FindApexOfTheChestPosition");
462 MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
463 MaskImagePointType p;
464 p[0] = p[1] = p[2] = 0.0; // to avoid warning
465 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
467 // We dont need Lungs structure from now
468 GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
470 // Put inside the AFDB
471 GetAFDB()->SetPoint3D("ApexOfTheChest", p);
472 p[2] -= 5; // We consider 5 mm lower
473 GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
479 //--------------------------------------------------------------------
482 //--------------------------------------------------------------------
483 template <class TImageType>
485 clitk::ExtractLymphStationsFilter<TImageType>::
486 FindLeftAndRightBronchi()
489 m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
490 m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
492 catch(clitk::ExceptionObject & o) {
494 DD("FindLeftAndRightBronchi");
495 // The goal is to separate the trachea inferiorly to the carina into
496 // a Left and Right bronchus.
499 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
501 // Get the Carina position
502 double m_CarinaZ = FindCarina();
504 // Consider only inferiorly to the Carina
505 MaskImagePointer m_Working_Trachea =
506 clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
507 GetBackgroundValue());
509 // Labelize the trachea
510 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
512 // Carina position must at the first slice that separate the two
513 // main bronchus (not superiorly). We thus first check that the
514 // upper slice is composed of at least two labels
515 MaskImagePointer RightBronchus;
516 MaskImagePointer LeftBronchus;
517 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
518 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
519 iter.SetFirstDirection(0);
520 iter.SetSecondDirection(1);
521 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
523 while (!iter.IsAtReverseEndOfSlice()) {
524 while (!iter.IsAtReverseEndOfLine()) {
525 if (iter.Get() > maxLabel) maxLabel = iter.Get();
531 clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
534 // Compute 3D centroids of both parts to identify the left from the
536 std::vector<ImagePointType> c;
537 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
538 ImagePointType C1 = c[1];
539 ImagePointType C2 = c[2];
541 ImagePixelType rightLabel;
542 ImagePixelType leftLabel;
543 if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
544 else { rightLabel = 2; leftLabel = 1; }
546 // Select LeftLabel (set one label to Backgroundvalue)
548 clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
549 GetBackgroundValue(), GetForegroundValue());
551 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
552 leftLabel, GetBackgroundValue(), false);
554 LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
555 GetBackgroundValue(), GetForegroundValue());
557 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
558 rightLabel, GetBackgroundValue(), false);
562 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
563 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
565 // Insert int AFDB if need after
566 GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
567 RightBronchus, true);
568 GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
572 //--------------------------------------------------------------------
575 //--------------------------------------------------------------------
576 template <class TImageType>
578 clitk::ExtractLymphStationsFilter<TImageType>::
579 FindSuperiorBorderOfAorticArch()
583 z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
585 catch(clitk::ExceptionObject e) {
586 DD("FindSuperiorBorderOfAorticArch");
587 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
588 MaskImagePointType p;
589 p[0] = p[1] = p[2] = 0.0; // to avoid warning
590 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
591 p[2] += Aorta->GetSpacing()[2]; // the slice above
593 // We dont need Lungs structure from now
594 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
596 // Put inside the AFDB
597 GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
598 GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
604 //--------------------------------------------------------------------
607 //--------------------------------------------------------------------
608 template <class TImageType>
610 clitk::ExtractLymphStationsFilter<TImageType>::
611 FindInferiorBorderOfAorticArch()
615 z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
617 catch(clitk::ExceptionObject e) {
618 DD("FindInferiorBorderOfAorticArch");
619 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
620 std::vector<MaskSlicePointer> slices;
621 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
623 uint i = slices.size()-1;
625 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
626 std::vector<typename MaskSliceType::PointType> c;
627 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
635 MaskImageIndexType index;
636 index[0] = index[1] = 0.0;
638 MaskImagePointType lower;
639 Aorta->TransformIndexToPhysicalPoint(index, lower);
641 // We dont need Lungs structure from now
642 GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
644 // Put inside the AFDB
645 GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
646 GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
652 //--------------------------------------------------------------------
655 //--------------------------------------------------------------------
656 template <class ImageType>
657 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
658 clitk::ExtractLymphStationsFilter<ImageType>::
661 // -----------------------------------------------------
662 /* Rod says: "The anterior border, as with the Atlas – UM, is
663 posterior to the vessels (right subclavian vein, left
664 brachiocephalic vein, right brachiocephalic vein, left subclavian
665 artery, left common carotid artery and brachiocephalic trunk).
666 These vessels are not included in the nodal station. The anterior
667 border is drawn to the midpoint of the vessel and an imaginary
668 line joins the middle of these vessels. Between the vessels,
669 station 2 is in contact with station 3a." */
671 // Check if no already done
673 MaskImagePointer binarizedContour;
675 DD("FindAntPostVessels try to get");
676 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
678 catch(clitk::ExceptionObject e) {
683 return binarizedContour;
686 /* Here, we consider the vessels form a kind of anterior barrier. We
687 link all vessels centroids and remove what is post to it.
688 - select the list of structure
689 vessel1 = BrachioCephalicArtery
690 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
691 vessel3 = CommonCarotidArtery
692 vessel4 = SubclavianArtery
695 - crop images as needed
696 - slice by slice, choose the CCL and extract centroids
697 - slice by slice, sort according to polar angle wrt Trachea centroid.
698 - slice by slice, link centroids and close contour
699 - remove outside this contour
704 MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
705 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
706 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
707 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
708 MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
709 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
710 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
712 // Create a temporay support
713 // From first slice of BrachioCephalicVein to end of 3A
714 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
715 MaskImagePointType p;
716 p[0] = p[1] = p[2] = 0.0; // to avoid warning
717 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
719 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
720 GetBackgroundValue(), 2, false, p);
722 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
723 false, GetBackgroundValue());
725 // Resize all structures like support
726 BrachioCephalicArtery =
727 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
728 CommonCarotidArtery =
729 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
731 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
733 clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
735 clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
736 BrachioCephalicVein =
737 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
739 clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
742 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
743 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
744 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
745 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
746 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
747 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
748 std::vector<MaskSlicePointer> slices_SubclavianArtery;
749 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
750 std::vector<MaskSlicePointer> slices_Thyroid;
751 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
752 std::vector<MaskSlicePointer> slices_Aorta;
753 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
754 std::vector<MaskSlicePointer> slices_Trachea;
755 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
756 unsigned int n= slices_BrachioCephalicArtery.size();
758 // Get the boundaries of one slice
759 std::vector<MaskSlicePointType> bounds;
760 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
762 // For all slices, for all structures, find the centroid and build the contour
763 // List of 3D points (for debug)
764 std::vector<MaskImagePointType> p3D;
766 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
767 for(unsigned int i=0; i<n; i++) {
768 // Labelize the slices
769 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
770 GetBackgroundValue(), true, 1);
771 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
772 GetBackgroundValue(), true, 1);
773 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
774 GetBackgroundValue(), true, 1);
775 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
776 GetBackgroundValue(), true, 1);
777 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
778 GetBackgroundValue(), true, 1);
779 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
780 GetBackgroundValue(), true, 1);
783 std::vector<MaskSlicePointType> points2D;
784 std::vector<MaskSlicePointType> centroids1;
785 std::vector<MaskSlicePointType> centroids2;
786 std::vector<MaskSlicePointType> centroids3;
787 std::vector<MaskSlicePointType> centroids4;
788 std::vector<MaskSlicePointType> centroids5;
789 std::vector<MaskSlicePointType> centroids6;
790 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
791 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
792 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
793 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
794 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
795 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
797 // BrachioCephalicVein -> when it is separated into two CCL, we
798 // only consider the most at Right one
799 if (centroids6.size() > 2) {
800 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
801 else centroids6.erase(centroids6.begin()+1);
804 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
805 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
806 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
810 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
811 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
812 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
813 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
814 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
815 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
817 // Sort by angle according to trachea centroid and vertical line,
818 // in polar coordinates :
819 // http://en.wikipedia.org/wiki/Polar_coordinate_system
820 std::vector<MaskSlicePointType> centroids_trachea;
821 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
822 typedef std::pair<MaskSlicePointType, double> PointAngleType;
823 std::vector<PointAngleType> angles;
824 for(unsigned int j=0; j<points2D.size(); j++) {
825 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
826 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
827 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
829 if (x>0) angle = atan(y/x);
830 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
831 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
833 if (y>0) angle = M_PI/2.0;
834 if (y<0) angle = -M_PI/2.0;
837 angle = clitk::rad2deg(angle);
838 // Angle is [-180;180] wrt the X axis. We change the X axis to
839 // be the vertical line, because we want to sort from Right to
840 // Left from Post to Ant.
845 if (angle<0) angle = 360-angle;
848 a.first = points2D[j];
853 // Do nothing if less than 2 points --> n
854 if (points2D.size() < 3) { //continue;
858 // Sort points2D according to polar angles
859 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
860 for(unsigned int j=0; j<angles.size(); j++) {
861 points2D[j] = angles[j].first;
863 // DDV(points2D, points2D.size());
865 /* When vessels are far away, we try to replace the line segment
866 with a curved line that join the two vessels but stay
867 approximately at the same distance from the trachea centroids
871 - let a and c be two successive vessels centroids
872 - id distance(a,c) < threshold, next point
876 - compute mid position between two successive points -
877 compute dist to trachea centroid for the 3 pts - if middle too
880 std::vector<MaskSlicePointType> toadd;
881 unsigned int index = 0;
883 while (index<points2D.size()-1) {
884 MaskSlicePointType a = points2D[index];
885 MaskSlicePointType c = points2D[index+1];
887 double dac = a.EuclideanDistanceTo(c);
890 MaskSlicePointType b;
891 b[0] = a[0]+(c[0]-a[0])/2.0;
892 b[1] = a[1]+(c[1]-a[1])/2.0;
894 // Compute distance to trachea centroid
895 MaskSlicePointType m = centroids_trachea[1];
896 double da = m.EuclideanDistanceTo(a);
897 double db = m.EuclideanDistanceTo(b);
898 //double dc = m.EuclideanDistanceTo(c);
900 // Mean distance, find point on the line from trachea centroid
902 double alpha = (da+db)/2.0;
903 MaskSlicePointType v;
904 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
905 v[0] = (b[0]-m[0])/n;
906 v[1] = (b[1]-m[1])/n;
907 MaskSlicePointType r;
908 r[0] = m[0]+alpha*v[0];
909 r[1] = m[1]+alpha*v[1];
910 points2D.insert(points2D.begin()+index+1, r);
916 // DDV(points2D, points2D.size());
918 // Add some points to close the contour
919 // - H line towards Right
920 MaskSlicePointType p = points2D[0];
922 points2D.insert(points2D.begin(), p);
923 // - corner Right/Post
925 points2D.insert(points2D.begin(), p);
926 // - H line towards Left
929 points2D.push_back(p);
930 // - corner Right/Post
932 points2D.push_back(p);
933 // Close contour with the first point
934 points2D.push_back(points2D[0]);
935 // DDV(points2D, points2D.size());
937 // Build 3D points from the 2D points
938 std::vector<ImagePointType> points3D;
939 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
940 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
942 // Build the mesh from the contour's points
943 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
944 append->AddInput(mesh);
947 // DEBUG: write points3D
948 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
950 // Build the final 3D mesh form the list 2D mesh
952 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
954 // Debug, write the mesh
956 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
958 w->SetFileName("bidon.vtk");
962 // Compute a single binary 3D image from the list of contours
963 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
964 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
965 filter->SetMesh(mesh);
966 filter->SetLikeImage(support);
968 binarizedContour = filter->GetOutput();
971 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
974 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
977 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
978 false, GetBackgroundValue());
980 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
981 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
983 return binarizedContour;
986 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
987 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
988 p[1] += 50; // 50 mm Post from this point must be kept
989 ImageIndexType index;
990 binarizedContour->TransformPhysicalPointToIndex(p, index);
991 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
993 // remove from support
994 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
995 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
996 boolFilter->InPlaceOn();
997 boolFilter->SetInput1(m_Working_Support);
998 boolFilter->SetInput2(binarizedContour);
999 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1000 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1002 boolFilter->SetOperationType(BoolFilterType::And);
1004 boolFilter->SetOperationType(BoolFilterType::AndNot);
1005 boolFilter->Update();
1006 m_Working_Support = boolFilter->GetOutput();
1010 //StopCurrentStep<MaskImageType>(m_Working_Support);
1011 //m_ListOfStations["2R"] = m_Working_Support;
1012 //m_ListOfStations["2L"] = m_Working_Support;
1014 //--------------------------------------------------------------------
1032 //--------------------------------------------------------------------
1033 template <class ImageType>
1034 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1035 clitk::ExtractLymphStationsFilter<ImageType>::
1036 FindAntPostVessels2()
1038 // -----------------------------------------------------
1039 /* Rod says: "The anterior border, as with the Atlas – UM, is
1040 posterior to the vessels (right subclavian vein, left
1041 brachiocephalic vein, right brachiocephalic vein, left subclavian
1042 artery, left common carotid artery and brachiocephalic trunk).
1043 These vessels are not included in the nodal station. The anterior
1044 border is drawn to the midpoint of the vessel and an imaginary
1045 line joins the middle of these vessels. Between the vessels,
1046 station 2 is in contact with station 3a." */
1048 // Check if no already done
1050 MaskImagePointer binarizedContour;
1052 DD("FindAntPostVessels try to get");
1053 binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1055 catch(clitk::ExceptionObject e) {
1060 return binarizedContour;
1063 /* Here, we consider the vessels form a kind of anterior barrier. We
1064 link all vessels centroids and remove what is post to it.
1065 - select the list of structure
1066 vessel1 = BrachioCephalicArtery
1067 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1068 vessel3 = CommonCarotidArtery
1069 vessel4 = SubclavianArtery
1072 - crop images as needed
1073 - slice by slice, choose the CCL and extract centroids
1074 - slice by slice, sort according to polar angle wrt Trachea centroid.
1075 - slice by slice, link centroids and close contour
1076 - remove outside this contour
1077 - merge with support
1081 std::map<std::string, MaskImagePointer> MapOfStructures;
1082 typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1083 MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1084 MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1085 MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1086 MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1087 MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1088 MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1089 MapOfStructures["Thyroid"] = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1090 MapOfStructures["Aorta"] = GetAFDB()->template GetImage<MaskImageType>("Aorta");
1091 MapOfStructures["Trachea"] = GetAFDB()->template GetImage<MaskImageType>("Trachea");
1093 std::vector<std::string> ListOfStructuresNames;
1095 // Create a temporay support
1096 // From first slice of BrachioCephalicVein to end of 3A
1097 MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1098 MaskImagePointType p;
1099 p[0] = p[1] = p[2] = 0.0; // to avoid warning
1100 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1101 GetBackgroundValue(), 2, true, p);
1103 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1104 GetBackgroundValue(), 2, false, p);
1105 double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1106 support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1107 false, GetBackgroundValue());
1109 // Resize all structures like support
1110 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1111 it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1115 typedef std::vector<MaskSlicePointer> SlicesType;
1116 std::map<std::string, SlicesType> MapOfSlices;
1117 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1119 clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1120 MapOfSlices[it->first] = s;
1123 unsigned int n= MapOfSlices["Trachea"].size();
1125 // Get the boundaries of one slice
1126 std::vector<MaskSlicePointType> bounds;
1127 ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1130 // For all slices, for all structures, find the centroid and build the contour
1131 // List of 3D points (for debug)
1132 std::vector<MaskImagePointType> p3D;
1133 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1134 for(unsigned int i=0; i<n; i++) {
1136 // Labelize the slices
1137 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1138 MaskSlicePointer & s = MapOfSlices[it->first][i];
1139 s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1143 std::vector<MaskSlicePointType> points2D;
1144 typedef std::vector<MaskSlicePointType> CentroidsType;
1145 std::map<std::string, CentroidsType> MapOfCentroids;
1146 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1147 std::string structure = it->first;
1148 MaskSlicePointer & s = MapOfSlices[structure][i];
1150 clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1151 MapOfCentroids[structure] = c;
1155 // BrachioCephalicVein -> when it is separated into two CCL, we
1156 // only consider the most at Right one
1157 CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1159 if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1160 else c.erase(c.begin()+1);
1164 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1165 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1166 if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1167 && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1168 MapOfCentroids["BrachioCephalicVein"].clear();
1172 // Add all 2D points
1173 for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1174 std::string structure = it->first;
1175 if (structure != "Trachea") { // do not add centroid of trachea
1176 CentroidsType & c = MapOfCentroids[structure];
1177 for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1181 // Sort by angle according to trachea centroid and vertical line,
1182 // in polar coordinates :
1183 // http://en.wikipedia.org/wiki/Polar_coordinate_system
1184 // std::vector<MaskSlicePointType> centroids_trachea;
1185 //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1186 CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1187 typedef std::pair<MaskSlicePointType, double> PointAngleType;
1188 std::vector<PointAngleType> angles;
1189 for(unsigned int j=0; j<points2D.size(); j++) {
1190 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1191 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1192 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1194 if (x>0) angle = atan(y/x);
1195 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1196 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1198 if (y>0) angle = M_PI/2.0;
1199 if (y<0) angle = -M_PI/2.0;
1200 if (y==0) angle = 0;
1202 angle = clitk::rad2deg(angle);
1203 // Angle is [-180;180] wrt the X axis. We change the X axis to
1204 // be the vertical line, because we want to sort from Right to
1205 // Left from Post to Ant.
1207 angle = (270-angle);
1210 if (angle<0) angle = 360-angle;
1213 a.first = points2D[j];
1215 angles.push_back(a);
1218 // Do nothing if less than 2 points --> n
1219 if (points2D.size() < 3) { //continue;
1223 // Sort points2D according to polar angles
1224 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1225 for(unsigned int j=0; j<angles.size(); j++) {
1226 points2D[j] = angles[j].first;
1228 // DDV(points2D, points2D.size());
1230 /* When vessels are far away, we try to replace the line segment
1231 with a curved line that join the two vessels but stay
1232 approximately at the same distance from the trachea centroids
1236 - let a and c be two successive vessels centroids
1237 - id distance(a,c) < threshold, next point
1241 - compute mid position between two successive points -
1242 compute dist to trachea centroid for the 3 pts - if middle too
1245 std::vector<MaskSlicePointType> toadd;
1246 unsigned int index = 0;
1248 while (index<points2D.size()-1) {
1249 MaskSlicePointType a = points2D[index];
1250 MaskSlicePointType c = points2D[index+1];
1252 double dac = a.EuclideanDistanceTo(c);
1255 MaskSlicePointType b;
1256 b[0] = a[0]+(c[0]-a[0])/2.0;
1257 b[1] = a[1]+(c[1]-a[1])/2.0;
1259 // Compute distance to trachea centroid
1260 MaskSlicePointType m = centroids_trachea[1];
1261 double da = m.EuclideanDistanceTo(a);
1262 double db = m.EuclideanDistanceTo(b);
1263 //double dc = m.EuclideanDistanceTo(c);
1265 // Mean distance, find point on the line from trachea centroid
1267 double alpha = (da+db)/2.0;
1268 MaskSlicePointType v;
1269 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1270 v[0] = (b[0]-m[0])/n;
1271 v[1] = (b[1]-m[1])/n;
1272 MaskSlicePointType r;
1273 r[0] = m[0]+alpha*v[0];
1274 r[1] = m[1]+alpha*v[1];
1275 points2D.insert(points2D.begin()+index+1, r);
1281 // DDV(points2D, points2D.size());
1283 // Add some points to close the contour
1284 // - H line towards Right
1285 MaskSlicePointType p = points2D[0];
1286 p[0] = bounds[0][0];
1287 points2D.insert(points2D.begin(), p);
1288 // - corner Right/Post
1290 points2D.insert(points2D.begin(), p);
1291 // - H line towards Left
1292 p = points2D.back();
1293 p[0] = bounds[2][0];
1294 points2D.push_back(p);
1295 // - corner Right/Post
1297 points2D.push_back(p);
1298 // Close contour with the first point
1299 points2D.push_back(points2D[0]);
1300 // DDV(points2D, points2D.size());
1302 // Build 3D points from the 2D points
1303 std::vector<ImagePointType> points3D;
1304 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1305 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1307 // Build the mesh from the contour's points
1308 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1309 append->AddInput(mesh);
1310 // if (i ==n-1) { // last slice
1311 // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1312 // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1313 // append->AddInput(mesh);
1317 // DEBUG: write points3D
1318 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1320 // Build the final 3D mesh form the list 2D mesh
1322 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1324 // Debug, write the mesh
1326 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1328 w->SetFileName("bidon.vtk");
1332 // Compute a single binary 3D image from the list of contours
1333 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1334 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1335 filter->SetMesh(mesh);
1336 filter->SetLikeImage(support);
1338 binarizedContour = filter->GetOutput();
1341 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1342 GetForegroundValue(), 2, true, p);
1344 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1345 GetForegroundValue(), 2, false, p);
1346 sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1347 binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1348 false, GetBackgroundValue());
1351 writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1352 GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1354 return binarizedContour;
1357 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1358 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1359 p[1] += 50; // 50 mm Post from this point must be kept
1360 ImageIndexType index;
1361 binarizedContour->TransformPhysicalPointToIndex(p, index);
1362 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1364 // remove from support
1365 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1366 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1367 boolFilter->InPlaceOn();
1368 boolFilter->SetInput1(m_Working_Support);
1369 boolFilter->SetInput2(binarizedContour);
1370 boolFilter->SetBackgroundValue1(GetBackgroundValue());
1371 boolFilter->SetBackgroundValue2(GetBackgroundValue());
1373 boolFilter->SetOperationType(BoolFilterType::And);
1375 boolFilter->SetOperationType(BoolFilterType::AndNot);
1376 boolFilter->Update();
1377 m_Working_Support = boolFilter->GetOutput();
1381 //StopCurrentStep<MaskImageType>(m_Working_Support);
1382 //m_ListOfStations["2R"] = m_Working_Support;
1383 //m_ListOfStations["2L"] = m_Working_Support;
1385 //--------------------------------------------------------------------
1389 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX