2 #include <itkBinaryDilateImageFilter.h>
3 #include <itkMirrorPadImageFilter.h>
5 //--------------------------------------------------------------------
6 template <class ImageType>
8 clitk::ExtractLymphStationsFilter<ImageType>::
9 ExtractStation_8_SetDefaultValues()
11 SetDistanceMaxToAnteriorPartOfTheSpine(10);
13 p[0] = 15; p[1] = 2; p[2] = 1;
14 SetEsophagusDiltationForAnt(p);
15 p[0] = 5; p[1] = 10; p[2] = 1;
16 SetEsophagusDiltationForRight(p);
17 SetFuzzyThresholdForS8(0.5);
19 //--------------------------------------------------------------------
21 //--------------------------------------------------------------------
22 template <class ImageType>
24 clitk::ExtractLymphStationsFilter<ImageType>::
25 ExtractStation_8_SI_Limits()
28 Station 8: paraeosphageal nodes
30 "Superiorly, Station 8 begins at the level of the carina and,
31 therefore, abuts Station 3P (Fig. 3C). Inferiorly, the lower limit
32 of Station 8 was not specified in the Mountain and Dresler
33 classification (1). We delineate this volume until it reaches the
34 gastroesphogeal junction. "
36 => new classification IASCL 2009:
38 "Lower border: the diaphragm"
40 "Upper border: the upper border of the lower lobe bronchus on the
41 left; the lower border of the bronchus intermedius on the right"
44 StartNewStep("[Station8] Inf/Sup mediastinum limits with carina/diaphragm junction");
46 // Get Carina Z position
47 MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
49 // CHANGE TO -> OriginOfRightMiddleLobeBronchus ???
51 std::vector<MaskImagePointType> centroids;
52 clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
53 m_CarinaZ = centroids[1][2];
54 // add one slice to include carina ?
55 m_CarinaZ += m_Mediastinum->GetSpacing()[2];
56 // We dont need Carina structure from now
58 GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
60 // Get left lower lobe bronchus (L), take the upper border
61 // Get right bronchus intermedius (RML), take the lower border
64 double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) +
65 m_Mediastinum->GetSpacing()[2]; // add one slice to include carina
66 // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
69 // Found most inferior part of the lung
70 MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
71 // It should be already croped, so I took the origin and add 10mm above
72 m_DiaphragmInferiorLimit = Lungs->GetOrigin()[2]+10;
73 // Lungs->Delete(); // we don't need it, release memory -> it we want to release, also free in AFDB
74 clitk::PrintMemory(GetVerboseMemoryFlag(), "after reading lungs");
75 GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
76 clitk::PrintMemory(GetVerboseMemoryFlag(), "after release lungs");
79 Superior limit = carina
80 Inferior limit = DiaphragmInferiorLimit (old=Gastroesphogeal junction) */
82 clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
83 m_DiaphragmInferiorLimit, //GOjunctionZ,
85 GetBackgroundValue());
87 // Remove some structures that we know are excluded
88 MaskImagePointer VertebralBody =
89 GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
90 MaskImagePointer Aorta =
91 GetAFDB()->template GetImage <MaskImageType>("Aorta");
93 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
94 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
95 boolFilter->InPlaceOn();
96 boolFilter->SetInput1(m_Working_Support);
97 boolFilter->SetInput2(VertebralBody);
98 boolFilter->SetOperationType(BoolFilterType::AndNot);
100 boolFilter->SetInput1(boolFilter->GetOutput());
101 boolFilter->SetInput2(Aorta);
102 boolFilter->SetOperationType(BoolFilterType::AndNot);
103 boolFilter->Update();
104 m_Working_Support = boolFilter->GetOutput();
107 StopCurrentStep<MaskImageType>(m_Working_Support);
108 m_ListOfStations["8"] = m_Working_Support;
110 //--------------------------------------------------------------------
113 //--------------------------------------------------------------------
114 template <class ImageType>
116 clitk::ExtractLymphStationsFilter<ImageType>::
117 ExtractStation_8_Post_Limits()
120 Station 8: paraeosphageal nodes
122 Anteriorly, it is in contact with Station 7 and the
123 left main stem bronchus in its superior aspect (Fig. 3C–H) and
124 with the heart more inferiorly. Posteriorly, Station 8 abuts the
125 descending aorta and the anterior aspect of the vertebral body
126 until an imaginary horizontal line running 1 cm posterior to the
127 anterior border of the vertebral body (Fig. 3C).
129 New classification IASCL 2009 :
131 "Nodes lying adjacent to the wall of the esophagus and to the
132 right or left of the midline, excluding subcarinal nodes (S7)
137 StartNewStep("[Station8] Post limits with vertebral body");
138 MaskImagePointer VertebralBody =
139 GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
141 // Consider vertebral body slice by slice
142 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
143 typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
144 typedef typename ExtractSliceFilterType::SliceType SliceType;
145 std::vector<typename SliceType::Pointer> vertebralSlices;
146 extractSliceFilter->SetInput(VertebralBody);
147 extractSliceFilter->SetDirection(2);
148 extractSliceFilter->Update();
149 extractSliceFilter->GetOutputSlices(vertebralSlices);
151 // For each slice, compute the most anterior point
152 std::map<int, typename SliceType::PointType> vertebralAntPositionBySlice;
153 for(uint i=0; i<vertebralSlices.size(); i++) {
154 typename SliceType::PointType p;
155 bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vertebralSlices[i],
156 GetBackgroundValue(),
159 vertebralAntPositionBySlice[i] = p;
162 // no Foreground in this slice (it appends generally at the
163 // beginning and the end of the volume). Do nothing in this
168 // Convert 2D points in slice into 3D points
169 std::vector<MaskImagePointType> vertebralAntPositions;
170 clitk::PointsUtils<MaskImageType>::Convert2DTo3DList(vertebralAntPositionBySlice,
172 vertebralAntPositions);
174 // DEBUG : write list of points
175 clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
176 "vertebralMostAntPositions.txt");
178 // Cut support posteriorly 1cm the most anterior point of the
179 // VertebralBody. Do nothing below and above the VertebralBody. To
180 // do that compute several region, slice by slice and fill.
181 MaskImageRegionType region;
182 MaskImageSizeType size;
183 MaskImageIndexType start;
184 size[2] = 1; // a single slice
185 start[0] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[0];
186 size[0] = m_Working_Support->GetLargestPossibleRegion().GetSize()[0];
187 for(uint i=0; i<vertebralAntPositions.size(); i++) {
188 typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
190 IteratorType(m_Working_Support, m_Working_Support->GetLargestPossibleRegion());
191 MaskImageIndexType index;
192 // Consider some cm posterior to most anterior positions (usually
194 vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheSpine();
195 // Get index of this point
196 m_Working_Support->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
197 // Compute region (a single slice)
199 start[1] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[1]+index[1];
200 size[1] = m_Working_Support->GetLargestPossibleRegion().GetSize()[1]-start[1];
201 region.SetSize(size);
202 region.SetIndex(start);
204 if (m_Working_Support->GetLargestPossibleRegion().IsInside(start)) {
205 itk::ImageRegionIterator<MaskImageType> it(m_Working_Support, region);
207 while (!it.IsAtEnd()) {
208 it.Set(GetBackgroundValue());
214 StopCurrentStep<MaskImageType>(m_Working_Support);
215 m_ListOfStations["8"] = m_Working_Support;
217 //--------------------------------------------------------------------
220 //--------------------------------------------------------------------
221 template <class ImageType>
223 clitk::ExtractLymphStationsFilter<ImageType>::
224 ExtractStation_8_Ant_Limits()
226 //--------------------------------------------------------------------
227 StartNewStep("[Station8] Ant limits with S7 above Carina");
229 Anteriorly, it is in contact with Station 7 and the
230 left main stem bronchus in its superior aspect (Fig. 3C–H) and
231 with the heart more inferiorly.
233 1) line post wall bronchi (S7), till originRMLB
234 - LUL bronchus : to detect in trachea. But not needed here ??
235 2) heart ! -> not delineated.
236 ==> below S7, crop CT not to far away (ant), then try with threshold
238 1) ==> how to share with S7 ? Prepare both support at the same time !
239 NEED vector of initial support for each station ? No -> map if it exist before, take it !!
243 // Ant limit from carina (start) to end of S7 = originRMLB
245 MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
247 MaskImagePointer m_Working_Trachea =
248 clitk::CropImageAbove<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
249 GetBackgroundValue());
251 // Seprate into two main bronchi
252 MaskImagePointer LeftBronchus;
253 MaskImagePointer RightBronchus;
255 // Labelize and consider the two first (main) labels
256 m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
258 // Carina position must at the first slice that separate the two
259 // main bronchus (not superiorly). We thus first check that the
260 // upper slice is composed of at least two labels
261 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
262 SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
263 iter.SetFirstDirection(0);
264 iter.SetSecondDirection(1);
265 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
267 while (!iter.IsAtReverseEndOfSlice()) {
268 while (!iter.IsAtReverseEndOfLine()) {
269 if (iter.Get() > maxLabel) maxLabel = iter.Get();
275 clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
278 // Compute 3D centroids of both parts to identify the left from the
280 std::vector<ImagePointType> c;
281 clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
282 ImagePointType C1 = c[1];
283 ImagePointType C2 = c[2];
285 ImagePixelType leftLabel;
286 ImagePixelType rightLabel;
287 if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
288 else { leftLabel = 2; rightLabel = 1; }
290 // Select LeftLabel (set one label to Backgroundvalue)
292 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
293 rightLabel, GetBackgroundValue(), false);
295 SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
296 leftLabel, GetBackgroundValue(), false);
298 LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
299 RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
301 // Insert int AFDB if need after
302 GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
304 GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
305 RightBronchus, true);
307 // Now crop below OriginOfRightMiddleLobeBronchusZ
308 // It is not done before to keep entire bronchi.
310 MaskImagePointer OriginOfRightMiddleLobeBronchus =
311 GetAFDB()->template GetImage<MaskImageType>("OriginOfRightMiddleLobeBronchus");
312 std::vector<MaskImagePointType> centroids;
313 clitk::ComputeCentroids<MaskImageType>(OriginOfRightMiddleLobeBronchus, GetBackgroundValue(), centroids);
314 m_OriginOfRightMiddleLobeBronchusZ = centroids[1][2];
315 // add one slice to include carina ?
316 m_OriginOfRightMiddleLobeBronchusZ += LeftBronchus->GetSpacing()[2];
317 // We dont need Carina structure from now
318 OriginOfRightMiddleLobeBronchus->Delete();
321 clitk::CropImageBelow<MaskImageType>(LeftBronchus, 2,
322 m_OriginOfRightMiddleLobeBronchusZ,
324 GetBackgroundValue());
326 clitk::CropImageBelow<MaskImageType>(RightBronchus, 2,
327 m_OriginOfRightMiddleLobeBronchusZ,
329 GetBackgroundValue());
331 // Search for points that are the most left/post/ant and most
332 // right/post/ant of the left and right bronchus
333 // 15 = not 15 mm more distance than the middle point.
334 FindExtremaPointsInBronchus(LeftBronchus, 0, 10, m_RightMostInLeftBronchus,
335 m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
336 FindExtremaPointsInBronchus(RightBronchus, 1, 10, m_LeftMostInRightBronchus,
337 m_AntMostInRightBronchus, m_PostMostInRightBronchus);
339 // DEBUG : write the list of points
341 v.reserve(m_RightMostInLeftBronchus.size()+m_AntMostInLeftBronchus.size()+
342 m_PostMostInLeftBronchus.size());
343 v.insert(v.end(), m_RightMostInLeftBronchus.begin(), m_RightMostInLeftBronchus.end() );
344 v.insert(v.end(), m_AntMostInLeftBronchus.begin(), m_AntMostInLeftBronchus.end() );
345 v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() );
346 clitk::WriteListOfLandmarks<MaskImageType>(v, "LeftBronchusPoints.txt");
349 v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+
350 m_PostMostInRightBronchus.size());
351 v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() );
352 v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() );
353 v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() );
354 clitk::WriteListOfLandmarks<MaskImageType>(v, "RightBronchusPoints.txt");
357 // Now uses these points to limit, slice by slice
358 // line is mainly horizontal, so mainDirection=1
359 writeImage<MaskImageType>(m_Working_Support, "before.mhd");
360 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
361 m_PostMostInLeftBronchus,
362 m_PostMostInRightBronchus,
363 GetBackgroundValue(), 1, 10);
364 writeImage<MaskImageType>(m_Working_Support, "after.mhd");
366 // Keep main 3D CCL :
367 m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
368 m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
369 GetBackgroundValue(),
370 GetForegroundValue(), 1, 1, true);
373 m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
376 StopCurrentStep<MaskImageType>(m_Working_Support);
377 // m_ListOfStations["8"] = m_Working_Support;
381 //--------------------------------------------------------------------
382 template <class ImageType>
384 clitk::ExtractLymphStationsFilter<ImageType>::
385 ExtractStation_8_LR_Limits()
388 //--------------------------------------------------------------------
389 StartNewStep("[Station8] Ant part (not post to Esophagus)");
391 Consider Esophagus, dilate it and remove ant part. It remains part
392 on L & R, than can be partly removed by cutting what remains at
393 right of vertebral body.
397 MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
398 clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus");
400 // In images from the original article, Atlas – UM, the oesophagus
401 //was included in nodal stations 3p and 8. Having said that, in the
402 //description for station 8, it indicates that “the delineation of
403 //station 8 is limited to the soft tissues surrounding the
404 //oesophagus”. In the recent article, The IASLC Lung Cancer Staging
405 //Project (J Thorac Oncol 4:5, 568-77), the images drawn by
406 //Dr. Aletta Frasier exclude this structure. From an oncological
407 //prospective, the oesophagus should be excluded from these nodal
409 typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
410 typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
411 boolFilter->InPlaceOn();
412 boolFilter->SetInput1(m_Working_Support);
413 boolFilter->SetInput2(Esophagus);
414 boolFilter->SetOperationType(BoolFilterType::AndNot);
415 boolFilter->Update();
416 m_Working_Support = boolFilter->GetOutput();
418 // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ
420 clitk::CropImageAbove<MaskImageType>(Esophagus, 2,
421 m_OriginOfRightMiddleLobeBronchusZ,
423 GetBackgroundValue());
425 // Dilate to keep only not Anterior positions
426 MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
427 Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
428 Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
430 GetBackgroundValue(),
431 GetForegroundValue(), true);
432 writeImage<MaskImageType>(Esophagus, "enlarged-eso.mhd");
434 // Remove Anterior part according to this dilatated esophagus. Note:
435 // because we crop Esophagus with ORML, the support will also be
436 // croped in the same way. Here it is a desired feature. If we dont
437 // want, use SetIgnoreEmptySliceObject(true)
439 // In the new IASCL definition, it is not clear if sup limits is
440 // around carina or On the right, it is “the lower border of the
441 // bronchus intermedius”, indicated on the image set as a point
442 // (“lower border of the bronchus intermedius”)
444 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
445 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
446 relPosFilter->VerboseStepFlagOff();
447 relPosFilter->SetInput(m_Working_Support);
448 relPosFilter->SetInputObject(Esophagus);
449 relPosFilter->AddOrientationTypeString("P");
450 relPosFilter->InverseOrientationFlagOff();
451 relPosFilter->SetDirection(2); // Z axis
452 relPosFilter->UniqueConnectedComponentBySliceOff();
453 relPosFilter->SetIntermediateSpacing(3);
454 relPosFilter->ResampleBeforeRelativePositionFilterOn();
455 relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8());
456 relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated
457 relPosFilter->Update();
458 m_Working_Support = relPosFilter->GetOutput();
460 // AutoCrop (OR SbS ?)
461 m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
463 StopCurrentStep<MaskImageType>(m_Working_Support);
464 //--------------------------------------------------------------------
465 StartNewStep("[Station8] Left and Right limits arround esophagus (below Carina)");
467 // Estract slices for current support for slice by slice processing
468 std::vector<typename MaskSliceType::Pointer> slices;
469 clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
471 // Estract slices of Esophagus (resize like support before to have the same set of slices)
472 MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(Esophagus, m_Working_Support, GetBackgroundValue());
473 std::vector<typename MaskSliceType::Pointer> eso_slices;
474 clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
476 // Estract slices of Vertebral (resize like support before to have the same set of slices)
477 MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
478 VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
479 std::vector<typename MaskSliceType::Pointer> vert_slices;
480 clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
482 // Estract slices of Aorta (resize like support before to have the same set of slices)
483 MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
484 Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
485 std::vector<typename MaskSliceType::Pointer> aorta_slices;
486 clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
488 // Extract slices of Mediastinum (resize like support before to have the same set of slices)
489 m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
490 m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
491 std::vector<typename MaskSliceType::Pointer> mediast_slices;
492 clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
495 std::vector<MaskImagePointType> p_RightMostAnt;
496 std::vector<MaskImagePointType> p_RightMostPost;
497 std::vector<MaskImagePointType> p_LeftMostAnt;
498 std::vector<MaskImagePointType> p_LeftMostPost;
499 std::vector<MaskImagePointType> p_AllPoints;
502 In the following, we search for the LeftRight limits. We search
503 for the most Right points in Esophagus and in VertebralBody and
504 consider a line between those to most right points. All points in
505 the support which are most right to this line are discarded. Same
506 for the left part. The underlying assumption is that the support
507 is concave between Eso/VertebralBody. Esophagus is a bit
508 dilatated. On VertebralBody we go right (or left) until we reach
509 the lung (but no more 20 mm).
513 MaskImagePointType p;
514 MaskImagePointType pp;
515 for(uint i=0; i<slices.size() ; i++) {
516 // Declare all needed points (sp = slice point)
517 typename MaskSliceType::PointType sp_maxRight_Eso;
518 typename MaskSliceType::PointType sp_maxRight_Aorta;
519 typename MaskSliceType::PointType sp_maxRight_Vertebra;
520 typename MaskSliceType::PointType sp_maxLeft_Eso;
521 typename MaskSliceType::PointType sp_maxLeft_Aorta;
522 typename MaskSliceType::PointType sp_maxLeft_Vertebra;
524 // Right is at left on screen, coordinate decrease
525 // Left is at right on screen, coordinate increase
527 // Find limit of Vertebral -> only at most Post part of current
528 // slice support. First found most ant point in VertebralBody
529 typedef MaskSliceType SliceType;
530 typename SliceType::PointType p_slice_ant;
531 bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[i], GetBackgroundValue(), 1, true, p_slice_ant);
533 // It should not happen ! But sometimes, a contour is missing or
534 // the VertebralBody is not delineated enough inferiorly ... in
535 // those cases, we consider the first found slice.
536 std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
540 found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[j], GetBackgroundValue(), 1, true, p_slice_ant);
541 //clitkExceptionMacro("No foreground pixels in this VertebralBody slices ??");
545 p_slice_ant[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
547 // The, find most Right and Left points on that AP position
548 typename SliceType::IndexType indexR;
549 typename SliceType::IndexType indexL;
550 vert_slices[i]->TransformPhysicalPointToIndex(p_slice_ant, indexR);
552 // Check that is inside the mask
553 indexR[1] = std::min(indexR[1], (long)vert_slices[i]->GetLargestPossibleRegion().GetSize()[1]-1);
554 indexL[1] = indexR[1];
555 while (vert_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
556 indexR[0] --; // Go to the right
558 while (vert_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
559 indexL[0] ++; // Go to the left
561 vert_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
562 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, VertebralBody, i, p);
563 p_AllPoints.push_back(p);
564 vert_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
565 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, VertebralBody, i, p);
566 p_AllPoints.push_back(p);
568 // Find last point out of the mediastinum on this line, Right :
569 mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxRight_Vertebra, indexR);
570 double distance = 0.0;
571 while (mediast_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
573 distance += mediast_slices[i]->GetSpacing()[0];
575 if (distance < 30) { // Ok in this case, we found limit with lung
576 mediast_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
577 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
579 else { // in that case, we are probably below the diaphragm, so we
580 // add aribtrarly few mm
581 sp_maxRight_Vertebra[0] -= 2; // Leave 2 mm around the VertebralBody
582 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
584 p_RightMostPost.push_back(p);
585 p_AllPoints.push_back(p);
587 // Find last point out of the mediastinum on this line, Left :
588 mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxLeft_Vertebra, indexL);
590 while (mediast_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
592 distance += mediast_slices[i]->GetSpacing()[0];
594 if (distance < 30) { // Ok in this case, we found limit with lung
595 mediast_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
596 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
598 else { // in that case, we are probably below the diaphragm, so we
599 // add aribtrarly few mm
600 sp_maxLeft_Vertebra[0] += 2; // Leave 2 mm around the VertebralBody
601 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
603 p_LeftMostPost.push_back(p);
604 p_AllPoints.push_back(p);
606 // Find right limit of Esophagus and Aorta
607 clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, true, sp_maxRight_Eso);
608 bool f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, true, sp_maxRight_Aorta);
609 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Eso, EsophagusForSlice, i, p);
610 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Aorta, Aorta, i, pp);
611 pp[0] -= 2; // Add a margin of 2 mm to include the Aorta 'wall'
612 p_AllPoints.push_back(p);
614 p_AllPoints.push_back(pp);
615 MaskImagePointType A = p_RightMostPost.back();
616 MaskImagePointType B = p;
617 MaskImagePointType C = pp;
618 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
619 if (s>0) p_RightMostAnt.push_back(p);
620 else p_RightMostAnt.push_back(pp);
623 p_RightMostAnt.push_back(p);
626 // --------------------------------------------------------------------------
627 // Find the limit on the Left: most left point between Eso and
628 // Vertebra. (Left is left on screen, coordinate increase)
630 // Find left limit of Esophagus
631 clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Eso);
632 f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Aorta);
633 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Eso, EsophagusForSlice, i, p);
634 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Aorta, Aorta, i, pp);
635 p_AllPoints.push_back(p);
636 pp[0] += 2; // Add a margin of 2 mm to include the 'wall'
637 if (f) { // not below Aorta
638 p_AllPoints.push_back(pp);
639 MaskImagePointType A = p_LeftMostPost.back();
640 MaskImagePointType B = p;
641 MaskImagePointType C = pp;
642 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
643 if (s<0) p_LeftMostAnt.push_back(p); // Insert point most at Left
644 else p_LeftMostAnt.push_back(pp);
647 p_LeftMostAnt.push_back(p);
649 } // End of slice loop
651 clitk::WriteListOfLandmarks<MaskImageType>(p_AllPoints, "LR-Eso-Vert.txt");
653 // Now uses these points to limit, slice by slice
654 // Line is mainly vertical, so mainDirection=0
655 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
656 p_RightMostAnt, p_RightMostPost,
657 GetBackgroundValue(), 0, 10);
658 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
659 p_LeftMostAnt, p_LeftMostPost,
660 GetBackgroundValue(), 0, -10);
662 StopCurrentStep<MaskImageType>(m_Working_Support);
663 m_ListOfStations["8"] = m_Working_Support;
666 //--------------------------------------------------------------------
669 //--------------------------------------------------------------------
670 template <class ImageType>
671 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
672 clitk::ExtractLymphStationsFilter<ImageType>::
673 EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
675 // Check if Esophagus is delineated at least until GOjunction. Now,
676 // because we use AutoCrop, Origin[2] gives this max inferior
678 // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
680 if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) {
681 std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it."
683 double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit;
685 // Pad with few more slices
686 typedef itk::MirrorPadImageFilter<MaskImageType, MaskImageType> PadFilterType;
687 typename PadFilterType::Pointer padFilter = PadFilterType::New();
688 padFilter->SetInput(Esophagus);
690 b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1;
691 padFilter->SetPadLowerBound(b);
693 Esophagus = padFilter->GetOutput();
697 //--------------------------------------------------------------------
700 //--------------------------------------------------------------------
701 template <class ImageType>
703 clitk::ExtractLymphStationsFilter<ImageType>::
704 ExtractStation_8_LR_Limits_old()
707 Station 8: paraeosphageal nodes
709 Laterally, it is within the pleural envelope and again abuts the
710 descending aorta on the left. Reasonably, the delineation of
711 Station 8 is limited to the soft tissue surrounding the esophagus
715 StartNewStep("[Station8] Right limits (around esophagus)");
717 MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
719 // Autocrop to get first slice with starting Esophagus
720 Esophagus = clitk::AutoCrop<MaskImageType>(Esophagus, GetBackgroundValue());
723 // LR dilatation -> large to keep point inside
724 // AP dilatation -> few mm
725 // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ)
726 MaskImagePointType radiusInMM = GetEsophagusDiltationForRight();
727 Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
728 Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
730 GetBackgroundValue(),
731 GetForegroundValue(), true);
732 writeImage<MaskImageType>(Esophagus, "dilateEso2.mhd");
734 writeImage<MaskImageType>(m_Working_Support, "before-relpos2.mhd");
736 // Remove Right (Left on the screen) part according to this
737 // dilatated esophagus
738 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
739 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
740 relPosFilter->VerboseStepFlagOff();
741 relPosFilter->WriteStepFlagOff();
742 relPosFilter->SetInput(m_Working_Support);
743 relPosFilter->SetInputObject(Esophagus);
744 relPosFilter->AddOrientationTypeString("L");
745 relPosFilter->InverseOrientationFlagOn(); // Not Left to
746 relPosFilter->SetDirection(2); // Z axis
747 relPosFilter->UniqueConnectedComponentBySliceOff(); // important
748 relPosFilter->SetIntermediateSpacing(4);
749 relPosFilter->ResampleBeforeRelativePositionFilterOn();
750 relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
751 relPosFilter->RemoveObjectFlagOff();
752 relPosFilter->Update();
753 m_Working_Support = relPosFilter->GetOutput();
755 // Get a single 3D CCL
756 m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
757 m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
758 GetBackgroundValue(),
759 GetForegroundValue(), 1, 1, true);
763 // Re-Add post to Esophagus -> sometimes previous relpos remove some
764 // correct part below esophagus.
765 MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
766 EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
767 writeImage<MaskImageType>(Esophagus, "e-again.mhd");
768 typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
769 typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
770 relPosFilter->VerboseStepOff();
771 relPosFilter->WriteStepOff();
772 relPosFilter->SetInput(m_Working_Support);
773 relPosFilter->SetInputObject(Esophagus);
774 relPosFilter->SetOrientationTypeString("P");
775 relPosFilter->InverseOrientationFlagOff();
776 relPosFilter->SetDirection(2); // Z axis
777 relPosFilter->UniqueConnectedComponentBySliceOff(); // important
778 relPosFilter->SetIntermediateSpacing(4);
779 relPosFilter->ResampleBeforeRelativePositionFilterOn();
780 relPosFilter->CombineWithOrFlagOn();
781 relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
782 relPosFilter->RemoveObjectFlagOff();
783 relPosFilter->Update();
784 m_Working_Support = relPosFilter->GetOutput();
787 StopCurrentStep<MaskImageType>(m_Working_Support);
788 m_ListOfStations["8"] = m_Working_Support;
790 //--------------------------------------------------------------------