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(unsigned int 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 unsigned int 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(unsigned int i=0; i<dim; i++)
75 max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
76 image->TransformIndexToPhysicalPoint(min_i, min_c);
77 image->TransformIndexToPhysicalPoint(max_i, max_c);
79 // Get corners coordinates
80 HypercubeCorners<ImageType::ImageDimension>(bounds);
81 for(unsigned int i=0; i<bounds.size(); i++) {
82 for(unsigned int j=0; j<dim; j++) {
83 if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
84 if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
88 //--------------------------------------------------------------------
91 //--------------------------------------------------------------------
92 template <class ImageType>
94 clitk::ExtractLymphStationsFilter<ImageType>::
95 ExtractStation_2RL_SetDefaultValues()
97 SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7);
98 SetFuzzyThreshold("2RL", "BrachioCephalicArtery", 0.7);
99 SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3);
100 SetFuzzyThreshold("2RL", "Aorta", 0.7);
101 SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5);
102 SetFuzzyThreshold("2RL", "SubclavianArteryLeft", 0.8);
104 //--------------------------------------------------------------------
107 //--------------------------------------------------------------------
108 template <class ImageType>
110 clitk::ExtractLymphStationsFilter<ImageType>::
111 ExtractStation_2RL_SI_Limits()
113 // Apex of the chest or Sternum & Carina.
114 StartNewStep("[Station 2RL] Inf/Sup limits with Sternum and TopOfAorticArch/CaudalMarginOfLeftBrachiocephalicVein");
116 /* Rod says: "For the inferior border, unlike in Atlas – UM, there
117 is now a difference between 2R and 2L. 2R stops at the
118 intersection of the caudal margin of the innominate vein with the
119 trachea. 2L extends less inferiorly to the superior border of the
122 /* Get TopOfAorticArch and CaudalMarginOfLeftBrachiocephalicVein
123 - TopOfAorticArch -> can be obtain from Aorta -> most sup part.
124 - CaudalMarginOfLeftBrachiocephalicVein -> must inf part of BrachicephalicVein
126 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
127 MaskImagePointType p = Aorta->GetOrigin(); // initialise to avoid warning
128 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
129 double TopOfAorticArchZ=p[2];
130 GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ);
132 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
133 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
134 double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2];
135 GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ);
137 // First, cut on the most inferior part. Add one slice because this
138 // inf slice should not be included.
139 double inf = std::min(CaudalMarginOfLeftBrachiocephalicVeinZ, TopOfAorticArchZ) + m_Working_Support->GetSpacing()[2];
141 // Get Sternum and search for the upper position
142 MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
144 // Search most sup point
145 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, p);
146 double m_SternumZ = p[2];
147 // +Sternum->GetSpacing()[2]; // One more slice, because it is below this point // NOT HERE ?? WHY DIFFERENT FROM 3A ??
151 clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
152 inf, m_SternumZ, true,
153 GetBackgroundValue());
155 StopCurrentStep<MaskImageType>(m_Working_Support);
156 m_ListOfStations["2R"] = m_Working_Support;
157 m_ListOfStations["2L"] = m_Working_Support;
159 //--------------------------------------------------------------------
162 //--------------------------------------------------------------------
163 template <class ImageType>
165 clitk::ExtractLymphStationsFilter<ImageType>::
166 ExtractStation_2RL_Ant_Limits()
168 // -----------------------------------------------------
169 /* Rod says: "The anterior border, as with the Atlas – UM, is
170 posterior to the vessels (right subclavian vein, left
171 brachiocephalic vein, right brachiocephalic vein, left subclavian
172 artery, left common carotid artery and brachiocephalic trunk).
173 These vessels are not included in the nodal station. The anterior
174 border is drawn to the midpoint of the vessel and an imaginary
175 line joins the middle of these vessels. Between the vessels,
176 station 2 is in contact with station 3a." */
178 // -----------------------------------------------------
179 StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery");
181 // Get CommonCarotidArtery
182 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
184 // Remove Ant to CommonCarotidArtery
185 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
186 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
187 sliceRelPosFilter->VerboseStepFlagOff();
188 sliceRelPosFilter->WriteStepFlagOff();
189 sliceRelPosFilter->SetInput(m_Working_Support);
190 sliceRelPosFilter->SetInputObject(CommonCarotidArtery);
191 sliceRelPosFilter->SetDirection(2);
192 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "CommonCarotidArtery"));
193 sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
194 sliceRelPosFilter->IntermediateSpacingFlagOn();
195 sliceRelPosFilter->SetIntermediateSpacing(2);
196 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
197 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
198 sliceRelPosFilter->AutoCropFlagOn();
199 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
200 sliceRelPosFilter->RemoveObjectFlagOff();
201 sliceRelPosFilter->Update();
202 m_Working_Support = sliceRelPosFilter->GetOutput();
205 StopCurrentStep<MaskImageType>(m_Working_Support);
206 m_ListOfStations["2R"] = m_Working_Support;
207 m_ListOfStations["2L"] = m_Working_Support;
209 // -----------------------------------------------------
210 // Remove Ant to H line from the Ant most part of the
211 // CommonCarotidArtery until we reach the first slice of
212 // BrachioCephalicArtery
213 StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line");
215 // First, find the first slice of BrachioCephalicArtery
216 MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
217 MaskImagePointType p = BrachioCephalicArtery->GetOrigin(); // initialise to avoid warning
218 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicArtery, GetBackgroundValue(), 2, false, p);
219 double TopOfBrachioCephalicArteryZ=p[2] + BrachioCephalicArtery->GetSpacing()[2]; // Add one slice
221 // Remove CommonCarotidArtery below this point
222 CommonCarotidArtery = clitk::CropImageRemoveLowerThan<MaskImageType>(CommonCarotidArtery, 2, TopOfBrachioCephalicArteryZ, true, GetBackgroundValue());
224 // Find most Ant points
225 std::vector<MaskImagePointType> ccaAntPositionA;
226 std::vector<MaskImagePointType> ccaAntPositionB;
227 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(CommonCarotidArtery,
228 GetBackgroundValue(), 2,
230 0, // Horizontal line
234 // Cut ant to this line
235 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
238 GetBackgroundValue(), 1, 10);
241 StopCurrentStep<MaskImageType>(m_Working_Support);
242 m_ListOfStations["2R"] = m_Working_Support;
243 m_ListOfStations["2L"] = m_Working_Support;
245 // -----------------------------------------------------
246 // Ant limit with the BrachioCephalicArtery
247 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery line");
249 // Remove Ant to BrachioCephalicArtery
251 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicArtery, 2,
252 GetFuzzyThreshold("2RL", "BrachioCephalicArtery"), "NotAntTo", false, 2, true, false);
254 StopCurrentStep<MaskImageType>(m_Working_Support);
255 m_ListOfStations["2R"] = m_Working_Support;
256 m_ListOfStations["2L"] = m_Working_Support;
258 // -----------------------------------------------------
259 // Ant limit with the BrachioCephalicArtery H line
260 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery, Horizontal line");
262 // Find most Ant points
263 std::vector<MaskImagePointType> bctAntPositionA;
264 std::vector<MaskImagePointType> bctAntPositionB;
265 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(BrachioCephalicArtery,
266 GetBackgroundValue(), 2,
268 0, // Horizontal line
272 // Cut ant to this line
273 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
276 GetBackgroundValue(), 1, 10);
278 StopCurrentStep<MaskImageType>(m_Working_Support);
279 m_ListOfStations["2R"] = m_Working_Support;
280 m_ListOfStations["2L"] = m_Working_Support;
282 // -----------------------------------------------------
283 // Ant limit with the BrachioCephalicVein
284 StartNewStep("[Station 2RL] Ant limits with BrachioCephalicVein");
286 // Remove Ant to BrachioCephalicVein
287 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
289 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicVein, 2,
290 GetFuzzyThreshold("2RL", "BrachioCephalicVein"), "NotAntTo", false, 2, true, false);
292 StopCurrentStep<MaskImageType>(m_Working_Support);
293 m_ListOfStations["2R"] = m_Working_Support;
294 m_ListOfStations["2L"] = m_Working_Support;
296 //--------------------------------------------------------------------
299 //--------------------------------------------------------------------
300 template <class ImageType>
302 clitk::ExtractLymphStationsFilter<ImageType>::
303 ExtractStation_2RL_Ant_Limits2()
305 // -----------------------------------------------------
306 /* Rod says: "The anterior border, as with the Atlas – UM, is
307 posterior to the vessels (right subclavian vein, left
308 brachiocephalic vein, right brachiocephalic vein, left subclavian
309 artery, left common carotid artery and brachiocephalic trunk).
310 These vessels are not included in the nodal station. The anterior
311 border is drawn to the midpoint of the vessel and an imaginary
312 line joins the middle of these vessels. Between the vessels,
313 station 2 is in contact with station 3a." */
315 // -----------------------------------------------------
316 StartNewStep("[Station 2RL] Ant limits with vessels centroids");
318 /* Here, we consider the vessels form a kind of anterior barrier. We
319 link all vessels centroids and remove what is post to it.
320 - select the list of structure
321 vessel1 = BrachioCephalicArtery
322 vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
323 vessel3 = CommonCarotidArtery
324 vessel4 = SubclavianArtery
327 - crop images as needed
328 - slice by slice, choose the CCL and extract centroids
329 - slice by slice, sort according to polar angle wrt Trachea centroid.
330 - slice by slice, link centroids and close contour
331 - remove outside this contour
336 MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
337 MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
338 MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
339 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
340 MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
341 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
342 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
344 // Resize all structures like support
345 BrachioCephalicArtery =
346 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue());
347 CommonCarotidArtery =
348 clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
350 clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
352 clitk::ResizeImageLike<MaskImageType>(Thyroid, m_Working_Support, GetBackgroundValue());
354 clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
355 BrachioCephalicVein =
356 clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, m_Working_Support, GetBackgroundValue());
358 clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
361 std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
362 clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
363 std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
364 clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
365 std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
366 clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
367 std::vector<MaskSlicePointer> slices_SubclavianArtery;
368 clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
369 std::vector<MaskSlicePointer> slices_Thyroid;
370 clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
371 std::vector<MaskSlicePointer> slices_Aorta;
372 clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
373 std::vector<MaskSlicePointer> slices_Trachea;
374 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
375 unsigned int n= slices_BrachioCephalicArtery.size();
377 // Get the boundaries of one slice
378 std::vector<MaskSlicePointType> bounds;
379 ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
381 // For all slices, for all structures, find the centroid and build the contour
382 // List of 3D points (for debug)
383 std::vector<MaskImagePointType> p3D;
385 vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
386 for(unsigned int i=0; i<n; i++) {
387 // Labelize the slices
388 slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
389 GetBackgroundValue(), true, 1);
390 slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
391 GetBackgroundValue(), true, 1);
392 slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
393 GetBackgroundValue(), true, 1);
394 slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
395 GetBackgroundValue(), true, 1);
396 slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
397 GetBackgroundValue(), true, 1);
398 slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
399 GetBackgroundValue(), true, 1);
402 std::vector<MaskSlicePointType> points2D;
403 std::vector<MaskSlicePointType> centroids1;
404 std::vector<MaskSlicePointType> centroids2;
405 std::vector<MaskSlicePointType> centroids3;
406 std::vector<MaskSlicePointType> centroids4;
407 std::vector<MaskSlicePointType> centroids5;
408 std::vector<MaskSlicePointType> centroids6;
409 ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
410 ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
411 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
412 ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
413 ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
414 ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
416 // BrachioCephalicVein -> when it is separated into two CCL, we
417 // only consider the most at Right one
418 if (centroids6.size() > 2) {
419 if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
420 else centroids6.erase(centroids6.begin()+1);
423 // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
424 // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
425 if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
429 for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
430 for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
431 for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
432 for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
433 for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
434 for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
436 // Sort by angle according to trachea centroid and vertical line,
437 // in polar coordinates :
438 // http://en.wikipedia.org/wiki/Polar_coordinate_system
439 std::vector<MaskSlicePointType> centroids_trachea;
440 ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
441 typedef std::pair<MaskSlicePointType, double> PointAngleType;
442 std::vector<PointAngleType> angles;
443 for(unsigned int j=0; j<points2D.size(); j++) {
444 //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
445 double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
446 double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
448 if (x>0) angle = atan(y/x);
449 if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
450 if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
452 if (y>0) angle = M_PI/2.0;
453 if (y<0) angle = -M_PI/2.0;
456 angle = clitk::rad2deg(angle);
457 // Angle is [-180;180] wrt the X axis. We change the X axis to
458 // be the vertical line, because we want to sort from Right to
459 // Left from Post to Ant.
464 if (angle<0) angle = 360-angle;
467 a.first = points2D[j];
472 // Do nothing if less than 2 points --> n
473 if (points2D.size() < 3) { //continue;
477 // Sort points2D according to polar angles
478 std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
479 for(unsigned int j=0; j<angles.size(); j++) {
480 points2D[j] = angles[j].first;
482 // DDV(points2D, points2D.size());
484 /* When vessels are far away, we try to replace the line segment
485 with a curved line that join the two vessels but stay
486 approximately at the same distance from the trachea centroids
490 - let a and c be two successive vessels centroids
491 - id distance(a,c) < threshold, next point
495 - compute mid position between two successive points -
496 compute dist to trachea centroid for the 3 pts - if middle too
499 std::vector<MaskSlicePointType> toadd;
500 unsigned int index = 0;
502 while (index<points2D.size()-1) {
503 MaskSlicePointType a = points2D[index];
504 MaskSlicePointType c = points2D[index+1];
506 double dac = a.EuclideanDistanceTo(c);
509 MaskSlicePointType b;
510 b[0] = a[0]+(c[0]-a[0])/2.0;
511 b[1] = a[1]+(c[1]-a[1])/2.0;
513 // Compute distance to trachea centroid
514 MaskSlicePointType m = centroids_trachea[1];
515 double da = m.EuclideanDistanceTo(a);
516 double db = m.EuclideanDistanceTo(b);
517 //double dc = m.EuclideanDistanceTo(c);
519 // Mean distance, find point on the line from trachea centroid
521 double alpha = (da+db)/2.0;
522 MaskSlicePointType v;
523 double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
524 v[0] = (b[0]-m[0])/n;
525 v[1] = (b[1]-m[1])/n;
526 MaskSlicePointType r;
527 r[0] = m[0]+alpha*v[0];
528 r[1] = m[1]+alpha*v[1];
529 points2D.insert(points2D.begin()+index+1, r);
535 // DDV(points2D, points2D.size());
537 // Add some points to close the contour
538 // - H line towards Right
539 MaskSlicePointType p = points2D[0];
541 points2D.insert(points2D.begin(), p);
542 // - corner Right/Post
544 points2D.insert(points2D.begin(), p);
545 // - H line towards Left
548 points2D.push_back(p);
549 // - corner Right/Post
551 points2D.push_back(p);
552 // Close contour with the first point
553 points2D.push_back(points2D[0]);
554 // DDV(points2D, points2D.size());
556 // Build 3D points from the 2D points
557 std::vector<ImagePointType> points3D;
558 clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D);
559 for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
561 // Build the mesh from the contour's points
562 vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
563 append->AddInput(mesh);
566 // DEBUG: write points3D
567 clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
569 // Build the final 3D mesh form the list 2D mesh
571 vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
573 // Debug, write the mesh
575 vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
577 w->SetFileName("bidon.vtk");
581 // Compute a single binary 3D image from the list of contours
582 clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
583 clitk::MeshToBinaryImageFilter<MaskImageType>::New();
584 filter->SetMesh(mesh);
585 filter->SetLikeImage(m_Working_Support);
587 MaskImagePointer binarizedContour = filter->GetOutput();
589 // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
590 ImagePointType p = p3D[2]; // This is the first centroid of the first slice
591 p[1] += 50; // 50 mm Post from this point must be kept
592 ImageIndexType index;
593 binarizedContour->TransformPhysicalPointToIndex(p, index);
594 bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
596 // remove from support
597 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
598 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
599 boolFilter->InPlaceOn();
600 boolFilter->SetInput1(m_Working_Support);
601 boolFilter->SetInput2(binarizedContour);
602 boolFilter->SetBackgroundValue1(GetBackgroundValue());
603 boolFilter->SetBackgroundValue2(GetBackgroundValue());
605 boolFilter->SetOperationType(BoolFilterType::And);
607 boolFilter->SetOperationType(BoolFilterType::AndNot);
608 boolFilter->Update();
609 m_Working_Support = boolFilter->GetOutput();
612 StopCurrentStep<MaskImageType>(m_Working_Support);
613 m_ListOfStations["2R"] = m_Working_Support;
614 m_ListOfStations["2L"] = m_Working_Support;
616 //--------------------------------------------------------------------
619 //--------------------------------------------------------------------
620 template <class ImageType>
622 clitk::ExtractLymphStationsFilter<ImageType>::
623 ExtractStation_2RL_Post_Limits()
625 StartNewStep("[Station 2RL] Post limits with post wall of Trachea");
628 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
630 // Resize like the current support (to have the same number of slices)
631 Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
633 // Find extrema post positions
634 std::vector<MaskImagePointType> tracheaPostPositionsA;
635 std::vector<MaskImagePointType> tracheaPostPositionsB;
636 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
637 GetBackgroundValue(), 2,
639 0, // Horizontal line
641 tracheaPostPositionsA,
642 tracheaPostPositionsB);
643 // Cut post to this line
644 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
645 tracheaPostPositionsA,
646 tracheaPostPositionsB,
647 GetBackgroundValue(), 1, -10);
649 StopCurrentStep<MaskImageType>(m_Working_Support);
650 m_ListOfStations["2R"] = m_Working_Support;
651 m_ListOfStations["2L"] = m_Working_Support;
653 //--------------------------------------------------------------------
656 //--------------------------------------------------------------------
657 // Build a vtk mesh from a list of slice number/closed-contours
658 template <class ImageType>
659 vtkSmartPointer<vtkPolyData>
660 clitk::ExtractLymphStationsFilter<ImageType>::
661 Build3DMeshFrom2DContour(const std::vector<ImagePointType> & points)
663 // create a contour, polydata.
664 vtkSmartPointer<vtkPolyData> mesh = vtkSmartPointer<vtkPolyData>::New();
665 mesh->Allocate(); //for cell structures
666 mesh->SetPoints(vtkPoints::New());
668 int point_number = points.size();
669 for (unsigned int i=0; i<points.size(); i++) {
670 mesh->GetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]);
672 ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
673 mesh->GetLines()->InsertNextCell(2,ids);
678 //--------------------------------------------------------------------
681 //--------------------------------------------------------------------
682 template <class ImageType>
684 clitk::ExtractLymphStationsFilter<ImageType>::
685 ExtractStation_2RL_LR_Limits()
687 // ---------------------------------------------------------------------------
688 StartNewStep("[Station 2RL] Left/Right limits with Aorta");
689 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
690 // DD(GetFuzzyThreshold("2RL", "BrachioCephalicVein"));
692 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Aorta, 2,
693 GetFuzzyThreshold("2RL", "Aorta"),
694 "RightTo", false, 2, true, false);
696 StopCurrentStep<MaskImageType>(m_Working_Support);
697 m_ListOfStations["2R"] = m_Working_Support;
698 m_ListOfStations["2L"] = m_Working_Support;
700 // ---------------------------------------------------------------------------
701 StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Right)");
703 // SliceBySliceRelativePosition + select CCL most at Right
704 MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
705 typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
706 typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
707 sliceRelPosFilter->VerboseStepFlagOff();
708 sliceRelPosFilter->WriteStepFlagOff();
709 sliceRelPosFilter->SetInput(m_Working_Support);
710 sliceRelPosFilter->SetInputObject(SubclavianArtery);
711 sliceRelPosFilter->SetDirection(2);
712 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryRight"));
713 sliceRelPosFilter->AddOrientationTypeString("NotRightTo");
714 sliceRelPosFilter->IntermediateSpacingFlagOn();
715 sliceRelPosFilter->SetIntermediateSpacing(2);
716 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
717 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
719 sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
720 sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
721 sliceRelPosFilter->SetCCLSelectionDirection(-1); // select most at Right
722 sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOn(); // ignore if only one CCL
724 sliceRelPosFilter->AutoCropFlagOn();
725 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
726 sliceRelPosFilter->RemoveObjectFlagOff();
727 sliceRelPosFilter->Update();
728 m_Working_Support = sliceRelPosFilter->GetOutput();
731 StopCurrentStep<MaskImageType>(m_Working_Support);
732 m_ListOfStations["2R"] = m_Working_Support;
733 m_ListOfStations["2L"] = m_Working_Support;
736 // ---------------------------------------------------------------------------
737 StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Left)");
739 // SliceBySliceRelativePosition + select CCL most at Right
740 sliceRelPosFilter = SliceRelPosFilterType::New();
741 sliceRelPosFilter->VerboseStepFlagOff();
742 sliceRelPosFilter->WriteStepFlagOff();
743 sliceRelPosFilter->SetInput(m_Working_Support);
744 sliceRelPosFilter->SetInputObject(SubclavianArtery);
745 sliceRelPosFilter->SetDirection(2);
746 sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryLeft"));
747 sliceRelPosFilter->AddOrientationTypeString("NotLeftTo");
748 sliceRelPosFilter->IntermediateSpacingFlagOn();
749 sliceRelPosFilter->SetIntermediateSpacing(2);
750 sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
751 sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
753 sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
754 sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
755 sliceRelPosFilter->SetCCLSelectionDirection(+1); // select most at Left
756 sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOff(); // do not ignore if only one CCL
758 sliceRelPosFilter->AutoCropFlagOn();
759 sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
760 sliceRelPosFilter->RemoveObjectFlagOff();
761 sliceRelPosFilter->Update();
762 m_Working_Support = sliceRelPosFilter->GetOutput();
765 StopCurrentStep<MaskImageType>(m_Working_Support);
766 m_ListOfStations["2R"] = m_Working_Support;
767 m_ListOfStations["2L"] = m_Working_Support;
769 //--------------------------------------------------------------------
771 //--------------------------------------------------------------------
772 template <class ImageType>
774 clitk::ExtractLymphStationsFilter<ImageType>::
775 ExtractStation_2RL_Remove_Structures()
777 Remove_Structures("2RL", "BrachioCephalicVein");
778 Remove_Structures("2RL", "CommonCarotidArtery");
779 Remove_Structures("2RL", "SubclavianArtery");
780 Remove_Structures("2RL", "Thyroid");
781 Remove_Structures("2RL", "Aorta");
784 StopCurrentStep<MaskImageType>(m_Working_Support);
785 m_ListOfStations["2R"] = m_Working_Support;
786 m_ListOfStations["2L"] = m_Working_Support;
788 //--------------------------------------------------------------------
791 //--------------------------------------------------------------------
792 template <class ImageType>
794 clitk::ExtractLymphStationsFilter<ImageType>::
795 ExtractStation_2RL_SeparateRL()
797 // ---------------------------------------------------------------------------
798 StartNewStep("[Station 2RL] Separate 2R/2L according to Trachea");
802 "For station 2 there is a shift, dividing 2R from 2L, from midline
803 to the left paratracheal border."
806 - Consider Trachea SliceBySlice
807 - find extrema at Left
808 - add margins towards Right
809 - remove what is at Left of this line
813 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
815 // Resize like the current support (to have the same number of slices)
816 Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
818 // Find extrema post positions
819 std::vector<MaskImagePointType> tracheaLeftPositionsA;
820 std::vector<MaskImagePointType> tracheaLeftPositionsB;
821 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
822 GetBackgroundValue(), 2,
826 tracheaLeftPositionsA,
827 tracheaLeftPositionsB);
828 // Copy support for R and L
829 typedef itk::ImageDuplicator<MaskImageType> DuplicatorType;
830 DuplicatorType::Pointer duplicator = DuplicatorType::New();
831 duplicator->SetInputImage(m_Working_Support);
832 duplicator->Update();
833 MaskImageType::Pointer m_Working_Support2 = duplicator->GetOutput();
835 // Cut post to this line for Right part
836 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
837 tracheaLeftPositionsA,
838 tracheaLeftPositionsB,
839 GetBackgroundValue(), 0, -10);
840 writeImage<MaskImageType>(m_Working_Support, "R.mhd");
842 // Cut post to this line for Left part
843 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support2,
844 tracheaLeftPositionsA,
845 tracheaLeftPositionsB,
846 GetBackgroundValue(), 0, +10);
847 writeImage<MaskImageType>(m_Working_Support2, "L.mhd");
850 StopCurrentStep<MaskImageType>(m_Working_Support);
851 m_ListOfStations["2R"] = m_Working_Support;
852 m_ListOfStations["2L"] = m_Working_Support2;
854 //--------------------------------------------------------------------