4 #include <vtkAppendPolyData.h>
5 #include <vtkPolyDataWriter.h>
6 #include <vtkCellArray.h>
9 #include "clitkMeshToBinaryImageFilter.h"
12 #include <itkImageDuplicator.h>
14 //--------------------------------------------------------------------
15 template<class PointType>
16 class comparePointsX {
18 bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
20 //--------------------------------------------------------------------
23 //--------------------------------------------------------------------
24 template<class PairType>
25 class comparePointsWithAngle {
27 bool operator() (PairType i, PairType j) { return (i.second < j.second); }
29 //--------------------------------------------------------------------
32 //--------------------------------------------------------------------
34 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
35 std::vector<itk::Point<double, Dim-1> > previous;
36 HypercubeCorners<Dim-1>(previous);
37 out.resize(previous.size()*2);
38 for(uint i=0; i<out.size(); i++) {
39 itk::Point<double, Dim> p;
40 if (i<previous.size()) p[0] = 0;
42 for(int j=0; j<Dim-1; j++)
44 p[j+1] = previous[i%previous.size()][j];
49 //--------------------------------------------------------------------
52 //--------------------------------------------------------------------
54 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
59 //--------------------------------------------------------------------
62 //--------------------------------------------------------------------
63 template<class ImageType>
64 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
65 std::vector<typename ImageType::PointType> & bounds)
67 // Get image max/min coordinates
68 const uint dim=ImageType::ImageDimension;
69 typedef typename ImageType::PointType PointType;
70 typedef typename ImageType::IndexType IndexType;
71 PointType min_c, max_c;
72 IndexType min_i, max_i;
73 min_i = image->GetLargestPossibleRegion().GetIndex();
74 for(uint i=0; i<dim; i++) max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
75 image->TransformIndexToPhysicalPoint(min_i, min_c);
76 image->TransformIndexToPhysicalPoint(max_i, max_c);
78 // Get corners coordinates
79 HypercubeCorners<ImageType::ImageDimension>(bounds);
80 for(uint i=0; i<bounds.size(); i++) {
81 for(uint j=0; j<dim; j++) {
82 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
83 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template <class ImageType>
93 clitk::ExtractLymphStationsFilter<ImageType>::
94 ExtractStation_2RL_SetDefaultValues()
96 SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7);
97 SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", 0.7);
98 SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3);
99 SetFuzzyThreshold("2RL", "Aorta", 0.7);
100 SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5);
101 SetFuzzyThreshold("2RL", "SubclavianArteryLeft", 0.8);
103 //--------------------------------------------------------------------
106 //--------------------------------------------------------------------
107 template <class ImageType>
109 clitk::ExtractLymphStationsFilter<ImageType>::
110 ExtractStation_2RL_SI_Limits()
112 // Apex of the chest or Sternum & Carina.
113 StartNewStep("[Station 2RL] Inf/Sup limits with Sternum and TopOfAorticArch/CaudalMarginOfLeftBrachiocephalicVein");
115 /* Rod says: "For the inferior border, unlike in Atlas – UM, there
116 is now a difference between 2R and 2L. 2R stops at the
117 intersection of the caudal margin of the innominate vein with the
118 trachea. 2L extends less inferiorly to the superior border of the
121 /* Get TopOfAorticArch and CaudalMarginOfLeftBrachiocephalicVein
122 - TopOfAorticArch -> can be obtain from Aorta -> most sup part.
123 - CaudalMarginOfLeftBrachiocephalicVein -> must inf part of BrachicephalicVein
125 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
126 MaskImagePointType p = Aorta->GetOrigin(); // initialise to avoid warning
127 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
128 double TopOfAorticArchZ=p[2];
129 GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ);
131 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
132 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
133 double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2];
134 GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ);
136 // First, cut on the most inferior part. Add one slice because this
137 // inf slice should not be included.
138 double inf = std::min(CaudalMarginOfLeftBrachiocephalicVeinZ, TopOfAorticArchZ) + m_Working_Support->GetSpacing()[2];
140 // Get Sternum and search for the upper position
141 MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
143 // Search most sup point
144 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, p);
145 double m_SternumZ = p[2];
146 // +Sternum->GetSpacing()[2]; // One more slice, because it is below this point // NOT HERE ?? WHY DIFFERENT FROM 3A ??
150 clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
151 inf, m_SternumZ, true,
152 GetBackgroundValue());
154 StopCurrentStep<MaskImageType>(m_Working_Support);
155 m_ListOfStations["2R"] = m_Working_Support;
156 m_ListOfStations["2L"] = m_Working_Support;
158 //--------------------------------------------------------------------
161 //--------------------------------------------------------------------
162 template <class ImageType>
164 clitk::ExtractLymphStationsFilter<ImageType>::
165 ExtractStation_2RL_Ant_Limits()
167 // -----------------------------------------------------
168 /* Rod says: "The anterior border, as with the Atlas – UM, is
169 posterior to the vessels (right subclavian vein, left
170 brachiocephalic vein, right brachiocephalic vein, left subclavian
171 artery, left common carotid artery and brachiocephalic trunk).
172 These vessels are not included in the nodal station. The anterior
173 border is drawn to the midpoint of the vessel and an imaginary
174 line joins the middle of these vessels. Between the vessels,
175 station 2 is in contact with station 3a." */
177 // -----------------------------------------------------
178 StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery");
180 // Get CommonCarotidArtery
181 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
183 // Remove Ant to CommonCarotidArtery
184 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
185 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
186 sliceRelPosFilter->VerboseStepFlagOff();
187 sliceRelPosFilter->WriteStepFlagOff();
188 sliceRelPosFilter->SetInput(m_Working_Support);
189 sliceRelPosFilter->SetInputObject(CommonCarotidArtery);
190 sliceRelPosFilter->SetDirection(2);
191 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "CommonCarotidArtery"));
192 sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
193 sliceRelPosFilter->IntermediateSpacingFlagOn();
194 sliceRelPosFilter->SetIntermediateSpacing(2);
195 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
196 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
197 sliceRelPosFilter->AutoCropFlagOn();
198 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
199 sliceRelPosFilter->RemoveObjectFlagOff();
200 sliceRelPosFilter->Update();
201 m_Working_Support = sliceRelPosFilter->GetOutput();
204 StopCurrentStep<MaskImageType>(m_Working_Support);
205 m_ListOfStations["2R"] = m_Working_Support;
206 m_ListOfStations["2L"] = m_Working_Support;
208 // -----------------------------------------------------
209 // Remove Ant to H line from the Ant most part of the
210 // CommonCarotidArtery until we reach the first slice of
211 // BrachioCephalicTrunk
212 StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line");
214 // First, find the first slice of BrachioCephalicTrunk
215 MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicTrunk");
216 MaskImagePointType p = BrachioCephalicTrunk->GetOrigin(); // initialise to avoid warning
217 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicTrunk, GetBackgroundValue(), 2, false, p);
218 double TopOfBrachioCephalicTrunkZ=p[2] + BrachioCephalicTrunk->GetSpacing()[2]; // Add one slice
220 // Remove CommonCarotidArtery below this point
221 CommonCarotidArtery = clitk::CropImageBelow<MaskImageType>(CommonCarotidArtery, 2, TopOfBrachioCephalicTrunkZ, true, GetBackgroundValue());
223 // Find most Ant points
224 std::vector<MaskImagePointType> ccaAntPositionA;
225 std::vector<MaskImagePointType> ccaAntPositionB;
226 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(CommonCarotidArtery,
227 GetBackgroundValue(), 2,
229 0, // Horizontal line
233 // Cut ant to this line
234 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
237 GetBackgroundValue(), 1, 10);
240 StopCurrentStep<MaskImageType>(m_Working_Support);
241 m_ListOfStations["2R"] = m_Working_Support;
242 m_ListOfStations["2L"] = m_Working_Support;
244 // -----------------------------------------------------
245 // Ant limit with the BrachioCephalicTrunk
246 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk line");
248 // Remove Ant to BrachioCephalicTrunk
250 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicTrunk, 2,
251 GetFuzzyThreshold("2RL", "BrachioCephalicTrunk"), "NotAntTo", false, 2, true, false);
253 StopCurrentStep<MaskImageType>(m_Working_Support);
254 m_ListOfStations["2R"] = m_Working_Support;
255 m_ListOfStations["2L"] = m_Working_Support;
257 // -----------------------------------------------------
258 // Ant limit with the BrachioCephalicTrunk H line
259 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk, Horizontal line");
261 // Find most Ant points
262 std::vector<MaskImagePointType> bctAntPositionA;
263 std::vector<MaskImagePointType> bctAntPositionB;
264 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(BrachioCephalicTrunk,
265 GetBackgroundValue(), 2,
267 0, // Horizontal line
271 // Cut ant to this line
272 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
275 GetBackgroundValue(), 1, 10);
277 StopCurrentStep<MaskImageType>(m_Working_Support);
278 m_ListOfStations["2R"] = m_Working_Support;
279 m_ListOfStations["2L"] = m_Working_Support;
281 // -----------------------------------------------------
282 // Ant limit with the BrachioCephalicVein
283 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicVein");
285 // Remove Ant to BrachioCephalicVein
286 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
288 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicVein, 2,
289 GetFuzzyThreshold("2RL", "BrachioCephalicVein"), "NotAntTo", false, 2, true, false);
291 StopCurrentStep<MaskImageType>(m_Working_Support);
292 m_ListOfStations["2R"] = m_Working_Support;
293 m_ListOfStations["2L"] = m_Working_Support;
295 //--------------------------------------------------------------------
298 //--------------------------------------------------------------------
299 template <class ImageType>
301 clitk::ExtractLymphStationsFilter<ImageType>::
302 ExtractStation_2RL_Ant_Limits2()
304 // -----------------------------------------------------
305 /* Rod says: "The anterior border, as with the Atlas – UM, is
306 posterior to the vessels (right subclavian vein, left
307 brachiocephalic vein, right brachiocephalic vein, left subclavian
308 artery, left common carotid artery and brachiocephalic trunk).
309 These vessels are not included in the nodal station. The anterior
310 border is drawn to the midpoint of the vessel and an imaginary
311 line joins the middle of these vessels. Between the vessels,
312 station 2 is in contact with station 3a." */
314 // -----------------------------------------------------
315 StartNewStep("[Station 2RL] Ant limits with vessels centroids");
317 /* Here, we consider the vessels form a kind of anterior barrier. We
318 link all vessels centroids and remove what is post to it.
319 - select the list of structure
320 vessel1 = BrachioCephalicTrunk
321 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
322 vessel3 = CommonCarotidArtery
323 vessel4 = SubclavianArtery
326 - crop images as needed
327 - slice by slice, choose the CCL and extract centroids
328 - slice by slice, sort according to polar angle wrt Trachea centroid.
329 - slice by slice, link centroids and close contour
330 - remove outside this contour
335 MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicTrunk");
336 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
337 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
338 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
339 MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
340 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
341 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
343 // Resize all structures like support
344 BrachioCephalicTrunk =
345 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicTrunk, m_Working_Support, GetBackgroundValue());
346 CommonCarotidArtery =
347 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
349 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
351 clitk::ResizeImageLike<MaskImageType>(Thyroid, m_Working_Support, GetBackgroundValue());
353 clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
354 BrachioCephalicVein =
355 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, m_Working_Support, GetBackgroundValue());
357 clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
360 std::vector<MaskSlicePointer> slices_BrachioCephalicTrunk;
361 clitk::ExtractSlices<MaskImageType>(BrachioCephalicTrunk, 2, slices_BrachioCephalicTrunk);
362 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
363 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
364 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
365 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
366 std::vector<MaskSlicePointer> slices_SubclavianArtery;
367 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
368 std::vector<MaskSlicePointer> slices_Thyroid;
369 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
370 std::vector<MaskSlicePointer> slices_Aorta;
371 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
372 std::vector<MaskSlicePointer> slices_Trachea;
373 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
374 uint n= slices_BrachioCephalicTrunk.size();
376 // Get the boundaries of one slice
377 std::vector<MaskSlicePointType> bounds;
378 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicTrunk[0], bounds);
380 // For all slices, for all structures, find the centroid and build the contour
381 // List of 3D points (for debug)
382 std::vector<MaskImagePointType> p3D;
384 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
385 for(uint i=0; i<n; i++) {
386 // Labelize the slices
387 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
388 GetBackgroundValue(), true, 1);
389 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
390 GetBackgroundValue(), true, 1);
391 slices_BrachioCephalicTrunk[i] = Labelize<MaskSliceType>(slices_BrachioCephalicTrunk[i],
392 GetBackgroundValue(), true, 1);
393 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
394 GetBackgroundValue(), true, 1);
395 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
396 GetBackgroundValue(), true, 1);
397 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
398 GetBackgroundValue(), true, 1);
401 std::vector<MaskSlicePointType> points2D;
402 std::vector<MaskSlicePointType> centroids1;
403 std::vector<MaskSlicePointType> centroids2;
404 std::vector<MaskSlicePointType> centroids3;
405 std::vector<MaskSlicePointType> centroids4;
406 std::vector<MaskSlicePointType> centroids5;
407 std::vector<MaskSlicePointType> centroids6;
408 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
409 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
410 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicTrunk[i], GetBackgroundValue(), centroids3);
411 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
412 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
413 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
415 // BrachioCephalicVein -> when it is separated into two CCL, we
416 // only consider the most at Right one
417 if (centroids6.size() > 2) {
418 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
419 else centroids6.erase(centroids6.begin()+1);
422 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
423 // (BrachioCephalicTrunk is divided) -> forget BrachioCephalicVein
424 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
428 for(uint j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
429 for(uint j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
430 for(uint j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
431 for(uint j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
432 for(uint j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
433 for(uint j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
435 // Sort by angle according to trachea centroid and vertical line,
436 // in polar coordinates :
437 // http://en.wikipedia.org/wiki/Polar_coordinate_system
438 std::vector<MaskSlicePointType> centroids_trachea;
439 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
440 typedef std::pair<MaskSlicePointType, double> PointAngleType;
441 std::vector<PointAngleType> angles;
442 for(uint j=0; j<points2D.size(); j++) {
443 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
444 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
445 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
447 if (x>0) angle = atan(y/x);
448 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
449 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
451 if (y>0) angle = M_PI/2.0;
452 if (y<0) angle = -M_PI/2.0;
455 angle = clitk::rad2deg(angle);
456 // Angle is [-180;180] wrt the X axis. We change the X axis to
457 // be the vertical line, because we want to sort from Right to
458 // Left from Post to Ant.
463 if (angle<0) angle = 360-angle;
466 a.first = points2D[j];
471 // Do nothing if less than 2 points --> n
472 if (points2D.size() < 3) { //continue;
476 // Sort points2D according to polar angles
477 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
478 for(uint j=0; j<angles.size(); j++) {
479 points2D[j] = angles[j].first;
481 // DDV(points2D, points2D.size());
483 /* When vessels are far away, we try to replace the line segment
484 with a curved line that join the two vessels but stay
485 approximately at the same distance from the trachea centroids
489 - let a and c be two successive vessels centroids
490 - id distance(a,c) < threshold, next point
494 - compute mid position between two successive points -
495 compute dist to trachea centroid for the 3 pts - if middle too
498 std::vector<MaskSlicePointType> toadd;
501 while (index<points2D.size()-1) {
502 MaskSlicePointType a = points2D[index];
503 MaskSlicePointType c = points2D[index+1];
505 double dac = a.EuclideanDistanceTo(c);
508 MaskSlicePointType b;
509 b[0] = a[0]+(c[0]-a[0])/2.0;
510 b[1] = a[1]+(c[1]-a[1])/2.0;
512 // Compute distance to trachea centroid
513 MaskSlicePointType m = centroids_trachea[1];
514 double da = m.EuclideanDistanceTo(a);
515 double db = m.EuclideanDistanceTo(b);
516 //double dc = m.EuclideanDistanceTo(c);
518 // Mean distance, find point on the line from trachea centroid
520 double alpha = (da+db)/2.0;
521 MaskSlicePointType v;
522 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
523 v[0] = (b[0]-m[0])/n;
524 v[1] = (b[1]-m[1])/n;
525 MaskSlicePointType r;
526 r[0] = m[0]+alpha*v[0];
527 r[1] = m[1]+alpha*v[1];
528 points2D.insert(points2D.begin()+index+1, r);
534 // DDV(points2D, points2D.size());
536 // Add some points to close the contour
537 // - H line towards Right
538 MaskSlicePointType p = points2D[0];
540 points2D.insert(points2D.begin(), p);
541 // - corner Right/Post
543 points2D.insert(points2D.begin(), p);
544 // - H line towards Left
547 points2D.push_back(p);
548 // - corner Right/Post
550 points2D.push_back(p);
551 // Close contour with the first point
552 points2D.push_back(points2D[0]);
553 // DDV(points2D, points2D.size());
555 // Build 3D points from the 2D points
556 std::vector<ImagePointType> points3D;
557 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D);
558 for(uint x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
560 // Build the mesh from the contour's points
561 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
562 append->AddInput(mesh);
565 // DEBUG: write points3D
566 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
568 // Build the final 3D mesh form the list 2D mesh
570 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
572 // Debug, write the mesh
574 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
576 w->SetFileName("bidon.vtk");
580 // Compute a single binary 3D image from the list of contours
581 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
582 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
583 filter->SetMesh(mesh);
584 filter->SetLikeImage(m_Working_Support);
586 MaskImagePointer binarizedContour = filter->GetOutput();
588 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
589 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
590 p[1] += 50; // 50 mm Post from this point must be kept
591 ImageIndexType index;
592 binarizedContour->TransformPhysicalPointToIndex(p, index);
593 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
595 // remove from support
596 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
597 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
598 boolFilter->InPlaceOn();
599 boolFilter->SetInput1(m_Working_Support);
600 boolFilter->SetInput2(binarizedContour);
601 boolFilter->SetBackgroundValue1(GetBackgroundValue());
602 boolFilter->SetBackgroundValue2(GetBackgroundValue());
604 boolFilter->SetOperationType(BoolFilterType::And);
606 boolFilter->SetOperationType(BoolFilterType::AndNot);
607 boolFilter->Update();
608 m_Working_Support = boolFilter->GetOutput();
611 StopCurrentStep<MaskImageType>(m_Working_Support);
612 m_ListOfStations["2R"] = m_Working_Support;
613 m_ListOfStations["2L"] = m_Working_Support;
615 //--------------------------------------------------------------------
618 //--------------------------------------------------------------------
619 template <class ImageType>
621 clitk::ExtractLymphStationsFilter<ImageType>::
622 ExtractStation_2RL_Post_Limits()
624 StartNewStep("[Station 2RL] Post limits with post wall of Trachea");
627 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
629 // Resize like the current support (to have the same number of slices)
630 Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
632 // Find extrema post positions
633 std::vector<MaskImagePointType> tracheaPostPositionsA;
634 std::vector<MaskImagePointType> tracheaPostPositionsB;
635 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
636 GetBackgroundValue(), 2,
638 0, // Horizontal line
640 tracheaPostPositionsA,
641 tracheaPostPositionsB);
642 // Cut post to this line
643 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
644 tracheaPostPositionsA,
645 tracheaPostPositionsB,
646 GetBackgroundValue(), 1, -10);
648 StopCurrentStep<MaskImageType>(m_Working_Support);
649 m_ListOfStations["2R"] = m_Working_Support;
650 m_ListOfStations["2L"] = m_Working_Support;
652 //--------------------------------------------------------------------
655 //--------------------------------------------------------------------
656 // Build a vtk mesh from a list of slice number/closed-contours
657 template <class ImageType>
658 vtkSmartPointer<vtkPolyData>
659 clitk::ExtractLymphStationsFilter<ImageType>::
660 Build3DMeshFrom2DContour(const std::vector<ImagePointType> & points)
662 // create a contour, polydata.
663 vtkSmartPointer<vtkPolyData> mesh = vtkSmartPointer<vtkPolyData>::New();
664 mesh->Allocate(); //for cell structures
665 mesh->SetPoints(vtkPoints::New());
667 int point_number = points.size();
668 for (unsigned int i=0; i<points.size(); i++) {
669 mesh->GetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]);
671 ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
672 mesh->GetLines()->InsertNextCell(2,ids);
677 //--------------------------------------------------------------------
680 //--------------------------------------------------------------------
681 template <class ImageType>
683 clitk::ExtractLymphStationsFilter<ImageType>::
684 ExtractStation_2RL_LR_Limits()
686 // ---------------------------------------------------------------------------
687 StartNewStep("[Station 2RL] Left/Right limits with Aorta");
688 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
689 // DD(GetFuzzyThreshold("2RL", "BrachioCephalicVein"));
691 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Aorta, 2,
692 GetFuzzyThreshold("2RL", "Aorta"),
693 "RightTo", false, 2, true, false);
695 StopCurrentStep<MaskImageType>(m_Working_Support);
696 m_ListOfStations["2R"] = m_Working_Support;
697 m_ListOfStations["2L"] = m_Working_Support;
699 // ---------------------------------------------------------------------------
700 StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Right)");
702 // SliceBySliceRelativePosition + select CCL most at Right
703 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
704 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
705 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
706 sliceRelPosFilter->VerboseStepFlagOff();
707 sliceRelPosFilter->WriteStepFlagOff();
708 sliceRelPosFilter->SetInput(m_Working_Support);
709 sliceRelPosFilter->SetInputObject(SubclavianArtery);
710 sliceRelPosFilter->SetDirection(2);
711 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryRight"));
712 sliceRelPosFilter->AddOrientationTypeString("NotRightTo");
713 sliceRelPosFilter->IntermediateSpacingFlagOn();
714 sliceRelPosFilter->SetIntermediateSpacing(2);
715 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
716 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
718 sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
719 sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
720 sliceRelPosFilter->SetCCLSelectionDirection(-1); // select most at Right
721 sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOn(); // ignore if only one CCL
723 sliceRelPosFilter->AutoCropFlagOn();
724 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
725 sliceRelPosFilter->RemoveObjectFlagOff();
726 sliceRelPosFilter->Update();
727 m_Working_Support = sliceRelPosFilter->GetOutput();
730 StopCurrentStep<MaskImageType>(m_Working_Support);
731 m_ListOfStations["2R"] = m_Working_Support;
732 m_ListOfStations["2L"] = m_Working_Support;
735 // ---------------------------------------------------------------------------
736 StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Left)");
738 // SliceBySliceRelativePosition + select CCL most at Right
739 sliceRelPosFilter = SliceRelPosFilterType::New();
740 sliceRelPosFilter->VerboseStepFlagOff();
741 sliceRelPosFilter->WriteStepFlagOff();
742 sliceRelPosFilter->SetInput(m_Working_Support);
743 sliceRelPosFilter->SetInputObject(SubclavianArtery);
744 sliceRelPosFilter->SetDirection(2);
745 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryLeft"));
746 sliceRelPosFilter->AddOrientationTypeString("NotLeftTo");
747 sliceRelPosFilter->IntermediateSpacingFlagOn();
748 sliceRelPosFilter->SetIntermediateSpacing(2);
749 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
750 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
752 sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
753 sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
754 sliceRelPosFilter->SetCCLSelectionDirection(+1); // select most at Left
755 sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOff(); // do not ignore if only one CCL
757 sliceRelPosFilter->AutoCropFlagOn();
758 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
759 sliceRelPosFilter->RemoveObjectFlagOff();
760 sliceRelPosFilter->Update();
761 m_Working_Support = sliceRelPosFilter->GetOutput();
764 StopCurrentStep<MaskImageType>(m_Working_Support);
765 m_ListOfStations["2R"] = m_Working_Support;
766 m_ListOfStations["2L"] = m_Working_Support;
768 //--------------------------------------------------------------------
770 //--------------------------------------------------------------------
771 template <class ImageType>
773 clitk::ExtractLymphStationsFilter<ImageType>::
774 ExtractStation_2RL_Remove_Structures()
776 Remove_Structures("2RL", "BrachioCephalicVein");
777 Remove_Structures("2RL", "CommonCarotidArtery");
778 Remove_Structures("2RL", "SubclavianArtery");
779 Remove_Structures("2RL", "Thyroid");
780 Remove_Structures("2RL", "Aorta");
783 StopCurrentStep<MaskImageType>(m_Working_Support);
784 m_ListOfStations["2R"] = m_Working_Support;
785 m_ListOfStations["2L"] = m_Working_Support;
787 //--------------------------------------------------------------------
790 //--------------------------------------------------------------------
791 template <class ImageType>
793 clitk::ExtractLymphStationsFilter<ImageType>::
794 ExtractStation_2RL_SeparateRL()
796 // ---------------------------------------------------------------------------
797 StartNewStep("[Station 2RL] Separate 2R/2L according to Trachea");
801 "For station 2 there is a shift, dividing 2R from 2L, from midline
802 to the left paratracheal border."
805 - Consider Trachea SliceBySlice
806 - find extrema at Left
807 - add margins towards Right
808 - remove what is at Left of this line
812 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
814 // Resize like the current support (to have the same number of slices)
815 Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
817 // Find extrema post positions
818 std::vector<MaskImagePointType> tracheaLeftPositionsA;
819 std::vector<MaskImagePointType> tracheaLeftPositionsB;
820 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
821 GetBackgroundValue(), 2,
825 tracheaLeftPositionsA,
826 tracheaLeftPositionsB);
827 // Copy support for R and L
828 typedef itk::ImageDuplicator<MaskImageType> DuplicatorType;
829 DuplicatorType::Pointer duplicator = DuplicatorType::New();
830 duplicator->SetInputImage(m_Working_Support);
831 duplicator->Update();
832 MaskImageType::Pointer m_Working_Support2 = duplicator->GetOutput();
834 // Cut post to this line for Right part
835 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
836 tracheaLeftPositionsA,
837 tracheaLeftPositionsB,
838 GetBackgroundValue(), 0, -10);
839 writeImage<MaskImageType>(m_Working_Support, "R.mhd");
841 // Cut post to this line for Left part
842 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support2,
843 tracheaLeftPositionsA,
844 tracheaLeftPositionsB,
845 GetBackgroundValue(), 0, +10);
846 writeImage<MaskImageType>(m_Working_Support2, "L.mhd");
849 StopCurrentStep<MaskImageType>(m_Working_Support);
850 m_ListOfStations["2R"] = m_Working_Support;
851 m_ListOfStations["2L"] = m_Working_Support2;
853 //--------------------------------------------------------------------