2 //--------------------------------------------------------------------
3 template <class TImageType>
5 clitk::ExtractLymphStationsFilter<TImageType>::
6 ExtractStation_7_SI_Limits()
9 MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
10 double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
12 double m_OriginOfRightMiddleLobeBronchusZ = GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2);
13 DD(m_OriginOfRightMiddleLobeBronchusZ);
16 Superior limit = carina
17 Inferior limit = origin right middle lobe bronchus */
18 StartNewStep("[Station7] Inf/Sup mediastinum limits with carina/bronchus");
20 clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
21 m_OriginOfRightMiddleLobeBronchusZ,
23 GetBackgroundValue());
25 Superior limit = carina
26 Inferior limit = origin right middle lobe bronchus*/
28 clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2,
29 m_OriginOfRightMiddleLobeBronchusZ,
31 GetBackgroundValue());
33 StopCurrentStep<MaskImageType>(m_Working_Support);
34 m_ListOfStations["7"] = m_Working_Support;
36 //--------------------------------------------------------------------
39 //--------------------------------------------------------------------
40 template <class TImageType>
42 clitk::ExtractLymphStationsFilter<TImageType>::
43 ExtractStation_7_RL_Limits()
45 // ----------------------------------------------------------------
46 // Separate trachea in two CCL
47 StartNewStep("[Station7] Separate trachea under carina");
49 // Labelize and consider two main labels
50 m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
52 // Carina position must at the first slice that separate the two
53 // main bronchus (not superiorly) Check that upper slice is composed
54 // of at least two labels
55 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
56 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
57 iter.SetFirstDirection(0);
58 iter.SetSecondDirection(1);
59 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
61 while (!iter.IsAtReverseEndOfSlice()) {
62 while (!iter.IsAtReverseEndOfLine()) {
63 if (iter.Get() > maxLabel) maxLabel = iter.Get();
69 clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
72 // Compute 3D centroids of both parts to identify the left from the
74 std::vector<ImagePointType> c;
75 clitk::ComputeCentroids<MaskImageType>(m_working_trachea, GetBackgroundValue(), c);
76 ImagePointType C1 = c[1];
77 ImagePointType C2 = c[2];
79 ImagePixelType leftLabel;
80 ImagePixelType rightLabel;
81 if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
82 else { leftLabel = 2; rightLabel = 1; }
84 StopCurrentStep<MaskImageType>(m_working_trachea);
86 //-----------------------------------------------------
87 // Select LeftLabel (set one label to Backgroundvalue)
89 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
90 rightLabel, GetBackgroundValue(), false);
92 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
93 leftLabel, GetBackgroundValue(), false);
95 StartNewStep("[Station7] Limits with bronchus (slice by slice) : RightTo left bronchus");
97 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
99 GetFuzzyThreshold(), "RightTo",
102 StartNewStep("[Station7] Limits with bronchus (slice by slice) : LeftTo right bronchus");
104 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
106 2, GetFuzzyThreshold(), "LeftTo",
109 StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo left bronchus");
111 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
113 2, GetFuzzyThreshold(), "AntTo",
114 true, 4, true); // NOT
116 StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo right bronchus");
118 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
120 2, GetFuzzyThreshold(), "AntTo",
123 StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo left bronchus");
125 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
127 2, GetFuzzyThreshold(), "PostTo",
130 StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo right bronchus");
132 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
134 2, GetFuzzyThreshold(), "PostTo",
136 m_Station7 = m_Working_Support;
137 StopCurrentStep<MaskImageType>(m_Station7);
139 //--------------------------------------------------------------------
142 //--------------------------------------------------------------------
143 template <class TImageType>
145 clitk::ExtractLymphStationsFilter<TImageType>::
146 ExtractStation_7_Posterior_Limits()
148 StartNewStep("[Station7] Posterior limits -> must be AntTo post wall of the bronchi");
150 // Search for points that are the most left/post/ant and most
151 // right/post/ant of the left and right bronchus
153 // extract, loop slices, label/keep, find extrema x 3
154 FindExtremaPointsInBronchus(m_LeftBronchus, 0, 15,
155 m_RightMostInLeftBronchus,
156 m_AntMostInLeftBronchus,
157 m_PostMostInLeftBronchus);
158 FindExtremaPointsInBronchus(m_RightBronchus, 1, 15,
159 m_LeftMostInRightBronchus,
160 m_AntMostInRightBronchus,
161 m_PostMostInRightBronchus);
163 std::ofstream osrl; openFileForWriting(osrl, "osrl.txt"); osrl << "LANDMARKS1" << std::endl;
164 std::ofstream osal; openFileForWriting(osal, "osal.txt"); osal << "LANDMARKS1" << std::endl;
165 std::ofstream ospl; openFileForWriting(ospl, "ospl.txt"); ospl << "LANDMARKS1" << std::endl;
166 std::ofstream osrr; openFileForWriting(osrr, "osrr.txt"); osrr << "LANDMARKS1" << std::endl;
167 std::ofstream osar; openFileForWriting(osar, "osar.txt"); osar << "LANDMARKS1" << std::endl;
168 std::ofstream ospr; openFileForWriting(ospr, "ospr.txt"); ospr << "LANDMARKS1" << std::endl;
170 for(uint i=0; i<m_RightMostInLeftBronchus.size(); i++) {
171 osrl << i << " " << m_RightMostInLeftBronchus[i][0] << " " << m_RightMostInLeftBronchus[i][1]
172 << " " << m_RightMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
173 osal << i << " " << m_AntMostInLeftBronchus[i][0] << " " << m_AntMostInLeftBronchus[i][1]
174 << " " << m_AntMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
175 ospl << i << " " << m_PostMostInLeftBronchus[i][0] << " " << m_PostMostInLeftBronchus[i][1]
176 << " " << m_PostMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
178 osrr << i << " " << m_LeftMostInRightBronchus[i][0] << " " << m_LeftMostInRightBronchus[i][1]
179 << " " << m_LeftMostInRightBronchus[i][2] << " 0 0 " << std::endl;
180 osar << i << " " << m_AntMostInRightBronchus[i][0] << " " << m_AntMostInRightBronchus[i][1]
181 << " " << m_AntMostInRightBronchus[i][2] << " 0 0 " << std::endl;
182 ospr << i << " " << m_PostMostInRightBronchus[i][0] << " " << m_PostMostInRightBronchus[i][1]
183 << " " << m_PostMostInRightBronchus[i][2] << " 0 0 " << std::endl;
189 // Now uses these points to limit, slice by slice
190 // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
192 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
193 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
194 This will equal zero if the point C is on the line formed by
195 points A and B, and will have a different sign depending on the
196 side. Which side this is depends on the orientation of your (x,y)
197 coordinates, but you can plug test values for A,B and C into this
198 formula to determine whether negative values are to the left or to
200 => to accelerate, start with formula, when change sign -> stop and fill
202 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
203 SliceIteratorType iter = SliceIteratorType(m_Working_Support,
204 m_Working_Support->GetLargestPossibleRegion());
205 iter.SetFirstDirection(0);
206 iter.SetSecondDirection(1);
209 MaskImageType::PointType A;
210 MaskImageType::PointType B;
211 MaskImageType::PointType C;
212 while (!iter.IsAtEnd()) {
213 A = m_PostMostInLeftBronchus[i];
214 B = m_PostMostInRightBronchus[i];
216 C[1] -= 10; // I know I must keep this point
217 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
218 bool isPositive = s<0;
219 while (!iter.IsAtEndOfSlice()) {
220 while (!iter.IsAtEndOfLine()) {
221 // Very slow, I know ... but image should be very small
222 m_Working_Support->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
223 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
224 if (s == 0) iter.Set(2);
226 if (s > 0) iter.Set(GetBackgroundValue());
229 if (s < 0) iter.Set(GetBackgroundValue());
239 //-----------------------------------------------------
240 // StartNewStep("[Station7] Anterior limits");
245 // Station 4R, Station 4L, the right pulmonary artery, and/or the
246 // left superior pulmonary vein
249 //-----------------------------------------------------
250 //-----------------------------------------------------
251 // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
255 m_Station7 = m_Working_Support;
257 //--------------------------------------------------------------------