2 #include <itkBinaryDilateImageFilter.h>
3 #include <itkMirrorPadImageFilter.h>
5 //--------------------------------------------------------------------
6 template <class ImageType>
8 clitk::ExtractLymphStationsFilter<ImageType>::
9 ExtractStationSupports()
11 // Get initial Mediastinum
12 m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
14 // Remove some computed structures
15 this->GetAFDB()->RemoveTag("CarinaZ");
16 this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
18 // Superior and inferior limits.
19 Support_SI_Limit("inferior", "Sup_to_Carina", "inferior", "Carina", 0);
20 Support_SI_Limit("superior", "Inf_to_Carina", "inferior", "Carina", m_Working_Support->GetSpacing()[2]);
22 // Initialise all others supports
23 // m_ListOfSupports["S1R"] = m_ListOfSupports["Sup_to_Carina"];
24 // m_ListOfSupports["S1L"] = m_ListOfSupports["Sup_to_Carina"];
25 // m_ListOfSupports["S2R"] = m_ListOfSupports["Sup_to_Carina"];
26 // m_ListOfSupports["S2L"] = m_ListOfSupports["Sup_to_Carina"];
27 // m_ListOfSupports["S3A"] = m_ListOfSupports["Sup_to_Carina"];
28 // m_ListOfSupports["S3P"] = m_ListOfSupports["Sup_to_Carina"];
29 // m_ListOfSupports["S4R"] = m_ListOfSupports["Sup_to_Carina"];
30 // m_ListOfSupports["S4L"] = m_ListOfSupports["Sup_to_Carina"];
31 // m_ListOfSupports["S5"] = m_Mediastinum; // Not above Carina
32 // m_ListOfSupports["S6"] = m_Mediastinum; // Not above Carina
34 // Read all support limits in a file and apply them
35 ReadSupportLimits(GetSupportLimitsFilename());
36 for(unsigned int i=0; i<m_ListOfSupportLimits.size(); i++) {
37 SupportLimitsType s = m_ListOfSupportLimits[i];
38 Support_SI_Limit(s.station_limit, s.station, s.structure_limit,
39 s.structure, s.offset*m_Working_Support->GetSpacing()[2]);
43 Support_LeftRight_S1R_S1L();
46 Support_LeftRight_S2R_S2L();
49 Support_LeftRight_S4R_S4L();
51 // Post limits of S1,S2,S4
52 Support_Post_S1S2S4();
55 StartNewStep("[Support] Ant limits of S3P with trachea");
56 m_ListOfSupports["S3P"] = LimitsWithTrachea(m_ListOfSupports["S3P"], 1, 0, 10);
59 StartNewStep("[Support] Ant limits of S3A with trachea");
60 m_ListOfSupports["S3A"] = LimitsWithTrachea(m_ListOfSupports["S3A"], 1, 0, -10);
63 // Below Carina S7,8,9,10
64 m_ListOfSupports["S7"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
65 m_ListOfSupports["S8"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
66 m_ListOfSupports["S9"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
67 m_ListOfSupports["S10"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
68 m_ListOfSupports["S11"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
70 // Store image filenames into AFDB
71 WriteImageSupport("S1R"); WriteImageSupport("S1L");
72 WriteImageSupport("S2R"); WriteImageSupport("S2L");
73 WriteImageSupport("S3A"); WriteImageSupport("S3P");
74 WriteImageSupport("S4R"); WriteImageSupport("S4L");
75 WriteImageSupport("S5");
76 WriteImageSupport("S6");
77 WriteImageSupport("S7");
78 WriteImageSupport("S8");
79 WriteImageSupport("S9");
80 WriteImageSupport("S10");
81 WriteImageSupport("S11");
84 //--------------------------------------------------------------------
87 //--------------------------------------------------------------------
88 template <class ImageType>
90 clitk::ExtractLymphStationsFilter<ImageType>::
91 Support_SI_Limit(const std::string station_limit, const std::string station,
92 const std::string structure_limit, const std::string structure,
95 if (!GetCheckSupportFlag())
96 StartNewStep("[Support] "+station_limit+" limit of "+station+" is "+structure_limit+" limit of "+structure);
99 if ((station_limit != "superior") && (station_limit != "inferior")) {
100 clitkExceptionMacro("Error station_limit must be 'inferior' or 'superior', not '"<< station_limit);
102 if ((structure_limit != "superior") && (structure_limit != "inferior")) {
103 clitkExceptionMacro("Error structure_limit must be 'inferior' or 'superior', not '"<< structure_limit);
106 // Get current support
107 if (m_ListOfSupports.find(station) == m_ListOfSupports.end()) {
108 // std::cerr << "Warning: support " << station << " not initialized" << std::endl;
109 m_ListOfSupports[station] = m_Mediastinum;
111 m_Working_Support = m_ListOfSupports[station];
113 // Get structure or structureZ
118 // Try to load structure and compute extrema point
119 if (this->GetAFDB()->TagExist(structure)) {
120 MaskImagePointer Structure = this->GetAFDB()->template GetImage <MaskImageType>(structure);
121 file = this->GetAFDB()->GetTagValue(structure);
122 MaskImagePointType p;
123 p[0] = p[1] = p[2] = 0.0; // to avoid warning
124 if (structure_limit == "superior")
125 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, false, p);
127 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, true, p);
132 // Try to load structureZ
133 if ((found==0) && (this->GetAFDB()->TagExist(structure+"Z"))) {
134 z = this->GetAFDB()->GetDouble(structure+"Z");
138 // Try to load structurePoint
139 if ((found==0) && (this->GetAFDB()->TagExist(structure+"Point"))) {
140 MaskImagePointType p;
141 this->GetAFDB()->GetPoint3D(structure+"Point", p);
146 // Try to see if it is an already computed support
148 if (m_ListOfSupports.find(structure) != m_ListOfSupports.end()) {
149 MaskImagePointer Structure = m_ListOfSupports[structure];
150 MaskImagePointType p;
151 if (structure_limit == "superior")
152 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, false, p);
154 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, true, p);
160 // Try special case : "FindApexOfTheChest"
161 if (structure == "FindApexOfTheChest") {
162 z = FindApexOfTheChest();
165 if (structure == "FindInferiorBorderOfAorticArch") {
166 z = FindInferiorBorderOfAorticArch();
169 if (structure == "FindSuperiorBorderOfAorticArch") {
170 z = FindSuperiorBorderOfAorticArch();
174 // If we find anything
176 std::cerr << "ERROR : I could not find " << structure << " nor " << structure << "Z nor "
177 << structure << "Point" << std::endl;
184 // Remove Lower or greater
185 if (station_limit == "inferior") {
187 clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2, z, true, GetBackgroundValue());
191 clitk::CropImageRemoveGreaterThan<MaskImageType>(m_Working_Support, 2, z, true, GetBackgroundValue());
194 // Check: if reference station is given, display information
195 if (GetCheckSupportFlag()) {
197 MaskImagePointer Ref = this->GetAFDB()->template GetImage <MaskImageType>(station+"_Ref");
198 MaskImagePointType p_support;
199 MaskImagePointType p_ref;
200 if (station_limit == "superior") {
201 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Ref, GetBackgroundValue(), 2, false, p_ref);
202 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(m_Working_Support, GetBackgroundValue(), 2, false, p_support);
205 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Ref, GetBackgroundValue(), 2, true, p_ref);
206 clitk::FindExtremaPointInAGivenDirection<MaskImageType>(m_Working_Support, GetBackgroundValue(), 2, true, p_support);
208 std::ostringstream os;
209 os << "[Support] \t" << station << "\t" << station_limit << " "
210 << "Z = " << z << std::setprecision(2) << std::fixed
211 << "\tSupport = " << p_support[2]
212 << "\tRef = " << p_ref[2]
213 << "\tdiff = " << p_support[2]-p_ref[2] << "\t"
214 << structure << " " << structure_limit;
215 if (found==1) os << " (S "+file+")";
216 if (found==2) os << " (Z)";
217 if (found==3) os << " (P)";
218 if (found==4) os << " (p)";
219 if (found==5) os << " (Apex)";
220 if (found==6) os << " (AorticArch)";
221 StartNewStep(os.str());
222 } catch(clitk::ExceptionObject e) { }
226 m_ListOfSupports[station] = m_Working_Support;
227 StopCurrentStep<MaskImageType>(m_Working_Support);
229 //--------------------------------------------------------------------
232 //--------------------------------------------------------------------
233 template <class ImageType>
235 clitk::ExtractLymphStationsFilter<ImageType>::
236 Support_LeftRight_S1R_S1L()
238 // Step S1RL : Left-Right
239 StartNewStep("[Support] Left-Right S1R S1L");
240 std::vector<ImagePointType> A;
241 std::vector<ImagePointType> B;
242 // Search for centroid positions of trachea
243 MaskImagePointer Trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");
244 MaskImagePointer S1RL = m_ListOfSupports["S1R"];
245 Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, S1RL, GetBackgroundValue());
246 std::vector<MaskSlicePointer> slices;
247 clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices);
248 for(uint i=0; i<slices.size(); i++) {
249 slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
250 std::vector<typename MaskSliceType::PointType> c;
251 clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
253 clitk::PointsUtils<MaskImageType>::Convert2DTo3D(c[1], Trachea, i, a);
259 clitk::WriteListOfLandmarks<MaskImageType>(A, "S1LR_A.txt");
260 clitk::WriteListOfLandmarks<MaskImageType>(B, "S1LR_B.txt");
263 MaskImagePointer S1R = clitk::Clone<MaskImageType>(S1RL);
264 MaskImagePointer S1L = clitk::Clone<MaskImageType>(S1RL);
267 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1R, A, B,
268 GetBackgroundValue(), 0, 10);
269 S1R = clitk::AutoCrop<MaskImageType>(S1R, GetBackgroundValue());
270 m_ListOfSupports["S1R"] = S1R;
273 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1L, A, B,
274 GetBackgroundValue(), 0, -10);
275 S1L = clitk::AutoCrop<MaskImageType>(S1L, GetBackgroundValue());
276 m_ListOfSupports["S1L"] = S1L;
277 StopCurrentStep<MaskImageType>(m_ListOfSupports["S1L"]);
279 //--------------------------------------------------------------------
282 //--------------------------------------------------------------------
283 template <class ImageType>
285 clitk::ExtractLymphStationsFilter<ImageType>::
286 Support_LeftRight_S2R_S2L()
288 // ---------------------------------------------------------------------------
289 /* Step : S2RL LeftRight
290 As for lymph node station 4R, 2R includes nodes extending to the
291 left lateral border of the trachea
292 Rod says: "For station 2 there is a shift, dividing 2R from 2L, from midline
293 to the left paratracheal border."
295 StartNewStep("[Support] Separate 2R/2L according to Trachea");
296 MaskImagePointer S2R = m_ListOfSupports["S2R"];
297 MaskImagePointer S2L = m_ListOfSupports["S2L"];
298 S2R = LimitsWithTrachea(S2R, 0, 1, -10);
299 S2L = LimitsWithTrachea(S2L, 0, 1, 10);
300 S2R = clitk::AutoCrop<MaskImageType>(S2R, GetBackgroundValue());
301 S2L = clitk::AutoCrop<MaskImageType>(S2L, GetBackgroundValue());
302 m_ListOfSupports["S2R"] = S2R;
303 m_ListOfSupports["S2L"] = S2L;
304 this->GetAFDB()->template ReleaseImage<MaskImageType>("Trachea");
306 //--------------------------------------------------------------------
309 //--------------------------------------------------------------------
310 template <class ImageType>
312 clitk::ExtractLymphStationsFilter<ImageType>::
313 Support_LeftRight_S4R_S4L()
315 // ---------------------------------------------------------------------------
316 /* Step : S4RL LeftRight
318 - 4R: includes right paratracheal nodes, and pretracheal nodes
319 extending to the left lateral border of trachea
321 - 4L: includes nodes to the left of the left lateral border of
322 the trachea, medial to the ligamentum arteriosum
326 StartNewStep("[Support] Left Right separation of 4R/4L");
328 MaskImagePointer S4R = m_ListOfSupports["S4R"];
329 MaskImagePointer S4L = m_ListOfSupports["S4L"];
330 S4R = LimitsWithTrachea(S4R, 0, 1, -10);
331 S4L = LimitsWithTrachea(S4L, 0, 1, 10);
332 m_ListOfSupports["S4R"] = S4R;
333 m_ListOfSupports["S4L"] = S4L;
335 //--------------------------------------------------------------------
338 //--------------------------------------------------------------------
339 template <class ImageType>
340 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
341 clitk::ExtractLymphStationsFilter<ImageType>::
342 LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection,
345 MaskImagePointType min, max;
346 GetMinMaxBoundary<MaskImageType>(input, min, max);
347 return LimitsWithTrachea(input, extremaDirection, lineDirection, offset, max[2]);
349 template <class ImageType>
350 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
351 clitk::ExtractLymphStationsFilter<ImageType>::
352 LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection,
353 double offset, double maxSupPosition)
356 Take the input mask, consider the trachea and limit according to
357 Left border of the trachea. Keep at Left or at Right according to
361 MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
363 // Find extrema post positions
364 std::vector<MaskImagePointType> tracheaLeftPositionsA;
365 std::vector<MaskImagePointType> tracheaLeftPositionsB;
367 // Crop Trachea only on the Sup-Inf axes, without autocrop
368 // Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, input, GetBackgroundValue());
369 MaskImagePointType min, max;
370 GetMinMaxBoundary<MaskImageType>(input, min, max);
371 Trachea = clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2, min[2], max[2],
372 false, GetBackgroundValue());
374 // Select the main CCL (because of bronchus)
375 Trachea = SliceBySliceKeepMainCCL<MaskImageType>(Trachea, GetBackgroundValue(), GetForegroundValue());
377 // Slice by slice, build the separation line
378 clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
379 GetBackgroundValue(), 2,
380 extremaDirection, false, // Left
381 lineDirection, // Vertical line
383 tracheaLeftPositionsA,
384 tracheaLeftPositionsB);
385 // Do not consider trachea above the limit
386 int indexMax=tracheaLeftPositionsA.size();
387 for(uint i=0; i<tracheaLeftPositionsA.size(); i++) {
388 if (tracheaLeftPositionsA[i][2] > maxSupPosition) {
390 i = tracheaLeftPositionsA.size(); // stop loop
393 tracheaLeftPositionsA.erase(tracheaLeftPositionsA.begin()+indexMax, tracheaLeftPositionsA.end());
394 tracheaLeftPositionsB.erase(tracheaLeftPositionsB.begin()+indexMax, tracheaLeftPositionsB.end());
396 // Cut post to this line
397 clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(input,
398 tracheaLeftPositionsA,
399 tracheaLeftPositionsB,
400 GetBackgroundValue(),
401 extremaDirection, offset);
402 MaskImagePointer output = clitk::AutoCrop<MaskImageType>(input, GetBackgroundValue());
405 //--------------------------------------------------------------------
408 //--------------------------------------------------------------------
409 template <class ImageType>
411 clitk::ExtractLymphStationsFilter<ImageType>::
412 Support_Post_S1S2S4()
414 StartNewStep("[Support] Post limits of S1RL, S2RL, S4RL");
416 double m_ApexOfTheChest = FindApexOfTheChest();
418 // Post limits with Trachea
419 MaskImagePointer S1R = m_ListOfSupports["S1R"];
420 MaskImagePointer S1L = m_ListOfSupports["S1L"];
421 MaskImagePointer S2R = m_ListOfSupports["S2R"];
422 MaskImagePointer S2L = m_ListOfSupports["S2L"];
423 MaskImagePointer S4R = m_ListOfSupports["S4R"];
424 MaskImagePointer S4L = m_ListOfSupports["S4L"];
425 m_ListOfSupports["S1R"] = LimitsWithTrachea(S1L, 1, 0, -10, m_ApexOfTheChest);
426 m_ListOfSupports["S1L"] = LimitsWithTrachea(S1R, 1, 0, -10, m_ApexOfTheChest);
427 m_ListOfSupports["S2R"] = LimitsWithTrachea(S2R, 1, 0, -10, m_ApexOfTheChest);
428 m_ListOfSupports["S2L"] = LimitsWithTrachea(S2L, 1, 0, -10, m_ApexOfTheChest);
429 m_ListOfSupports["S4R"] = LimitsWithTrachea(S4R, 1, 0, -10, m_ApexOfTheChest);
430 m_ListOfSupports["S4L"] = LimitsWithTrachea(S4L, 1, 0, -10, m_ApexOfTheChest);
432 //--------------------------------------------------------------------