1 //--------------------------------------------------------------------
2 template <class TImageType>
4 clitk::ExtractLymphStationsFilter<TImageType>::
5 ExtractStation_7_SI_Limits()
8 double m_CarinaZPositionInMM;
9 double m_MiddleLobeBronchusZPositionInMM;
12 m_Trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
13 m_CarinaZPositionInMM = GetAFDB()->GetPoint3D("carina", 2);
14 DD(m_CarinaZPositionInMM);
15 m_MiddleLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightMiddleLobeBronchus", 2);
16 DD(m_MiddleLobeBronchusZPositionInMM);
19 Superior limit = carina
20 Inferior limit = origin right middle lobe bronchus */
21 StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
23 clitk::CropImageAlongOneAxis<MaskImageType>(m_Support, 2,
24 m_MiddleLobeBronchusZPositionInMM,
25 m_CarinaZPositionInMM, true,
26 GetBackgroundValue());
27 StopCurrentStep<MaskImageType>(m_Working_Support);
30 Superior limit = carina
31 Inferior limit = origin right middle lobe bronchus*/
32 StartNewStep("Inf/Sup trachea limits with carina/bronchus");
34 clitk::CropImageAlongOneAxis<MaskImageType>(m_Trachea, 2,
35 m_MiddleLobeBronchusZPositionInMM,
36 m_CarinaZPositionInMM, true,
37 GetBackgroundValue());
38 StopCurrentStep<MaskImageType>(m_working_trachea);
40 m_Station7 = m_Working_Support;
42 //--------------------------------------------------------------------
44 //--------------------------------------------------------------------
45 template <class TImageType>
47 clitk::ExtractLymphStationsFilter<TImageType>::
48 ExtractStation_7_RL_Limits()
50 // ----------------------------------------------------------------
51 // Separate trachea in two CCL
52 StartNewStep("Separate trachea under carina");
54 // Labelize and consider two main labels
55 m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
57 // Carina position must at the first slice that separate the two main bronchus (not superiorly)
58 // Check that upper slice is composed of at least two labels
59 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
60 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
61 iter.SetFirstDirection(0);
62 iter.SetSecondDirection(1);
63 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
65 while (!iter.IsAtReverseEndOfSlice()) {
66 while (!iter.IsAtReverseEndOfLine()) {
67 // DD(iter.GetIndex());
68 if (iter.Get() > maxLabel) maxLabel = iter.Get();
74 clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
77 // Compute centroid of both parts to identify the left from the right bronchus
78 typedef long LabelType;
79 static const unsigned int Dim = ImageType::ImageDimension;
80 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
81 typedef itk::LabelMap< LabelObjectType > LabelMapType;
82 typedef itk::LabelImageToLabelMapFilter<MaskImageType, LabelMapType> ImageToMapFilterType;
83 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
84 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskImageType> ShapeFilterType;
85 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
86 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
87 imageToLabelFilter->SetInput(m_working_trachea);
88 statFilter->SetInput(imageToLabelFilter->GetOutput());
90 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
92 ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
93 ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
95 ImagePixelType leftLabel;
96 ImagePixelType rightLabel;
97 if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
98 else { leftLabel = 2; rightLabel = 1; }
102 StopCurrentStep<MaskImageType>(m_working_trachea);
104 //-----------------------------------------------------
105 StartNewStep("Left limits with bronchus (slice by slice)");
106 // Select LeftLabel (set one label to Backgroundvalue)
108 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
109 rightLabel, GetBackgroundValue());
111 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
112 leftLabel, GetBackgroundValue());
113 writeImage<MaskImageType>(m_LeftBronchus, "left.mhd");
114 writeImage<MaskImageType>(m_RightBronchus, "right.mhd");
117 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
119 GetFuzzyThreshold(), "RightTo",
122 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
124 2, GetFuzzyThreshold(), "LeftTo",
127 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
129 2, GetFuzzyThreshold(), "AntTo",
130 true, 4, true); // NOT
132 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
134 2, GetFuzzyThreshold(), "AntTo",
137 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
139 2, GetFuzzyThreshold(), "PostTo",
142 clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
144 2, GetFuzzyThreshold(), "PostTo",
146 m_Station7 = m_Working_Support;
147 StopCurrentStep<MaskImageType>(m_Station7);
149 //--------------------------------------------------------------------
152 //--------------------------------------------------------------------
153 template <class TImageType>
155 clitk::ExtractLymphStationsFilter<TImageType>::
156 ExtractStation_7_Posterior_Limits()
158 StartNewStep("Posterior limits");
160 // Left Bronchus slices
161 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
162 typedef typename ExtractSliceFilterType::SliceType SliceType;
163 typename ExtractSliceFilterType::Pointer
164 extractSliceFilter = ExtractSliceFilterType::New();
165 extractSliceFilter->SetInput(m_LeftBronchus);
166 extractSliceFilter->SetDirection(2);
167 extractSliceFilter->Update();
168 std::vector<typename SliceType::Pointer> leftBronchusSlices;
169 extractSliceFilter->GetOutputSlices(leftBronchusSlices);
171 // Right Bronchus slices
172 extractSliceFilter = ExtractSliceFilterType::New();
173 extractSliceFilter->SetInput(m_RightBronchus);
174 extractSliceFilter->SetDirection(2);
175 extractSliceFilter->Update();
176 std::vector<typename SliceType::Pointer> rightBronchusSlices ;
177 extractSliceFilter->GetOutputSlices(rightBronchusSlices);
179 assert(leftBronchusSlices.size() == rightBronchusSlices.size());
181 std::vector<MaskImageType::PointType> leftPoints;
182 std::vector<MaskImageType::PointType> rightPoints;
183 for(uint i=0; i<leftBronchusSlices.size(); i++) {
185 leftBronchusSlices[i] = Labelize<SliceType>(leftBronchusSlices[i], 0, true, 10);
186 leftBronchusSlices[i] = KeepLabels<SliceType>(leftBronchusSlices[i],
187 GetBackgroundValue(),
188 GetForegroundValue(), 1, 1, true);
189 rightBronchusSlices[i] = Labelize<SliceType>(rightBronchusSlices[i], 0, true, 10);
190 rightBronchusSlices[i] = KeepLabels<SliceType>(rightBronchusSlices[i],
191 GetBackgroundValue(),
192 GetForegroundValue(), 1, 1, true);
193 double distance_max_from_center_point = 15;
195 // ------- Find point in left Bronchus -------
196 // find right most point in left = rightMost
197 SliceType::PointType a;
198 SliceType::PointType rightMost =
199 clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
200 GetBackgroundValue(),
202 // find post most point in left, not far away from rightMost
203 SliceType::PointType p =
204 clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
205 GetBackgroundValue(),
207 distance_max_from_center_point);
208 MaskImageType::PointType pp;
209 pp[0] = p[0]; pp[1] = p[1];
210 pp[2] = i*m_LeftBronchus->GetSpacing()[2] + m_LeftBronchus->GetOrigin()[2];
211 leftPoints.push_back(pp);
213 // ------- Find point in right Bronchus -------
214 // find left most point in right = leftMost
215 SliceType::PointType leftMost =
216 clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
217 GetBackgroundValue(),
219 // find post most point in left, not far away from leftMost
220 p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
221 GetBackgroundValue(),
223 distance_max_from_center_point);
224 pp[0] = p[0]; pp[1] = p[1];
225 pp[2] = i*m_RightBronchus->GetSpacing()[2] + m_RightBronchus->GetOrigin()[2];
226 rightPoints.push_back(pp);
231 openFileForWriting(osl, "left.txt");
232 osl << "LANDMARKS1" << std::endl;
234 openFileForWriting(osr, "right.txt");
235 osr << "LANDMARKS1" << std::endl;
236 for(uint i=0; i<leftBronchusSlices.size(); i++) {
237 osl << i << " " << leftPoints[i][0] << " " << leftPoints[i][1]
238 << " " << leftPoints[i][2] << " 0 0 " << std::endl;
239 osr << i << " " << rightPoints[i][0] << " " << rightPoints[i][1]
240 << " " << rightPoints[i][2] << " 0 0 " << std::endl;
245 // Now uses these points to limit, slice by slice
246 // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
248 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
249 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
250 This will equal zero if the point C is on the line formed by
251 points A and B, and will have a different sign depending on the
252 side. Which side this is depends on the orientation of your (x,y)
253 coordinates, but you can plug test values for A,B and C into this
254 formula to determine whether negative values are to the left or to
256 => to accelerate, start with formula, when change sign -> stop and fill
258 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
259 SliceIteratorType iter = SliceIteratorType(m_Working_Support,
260 m_Working_Support->GetLargestPossibleRegion());
261 iter.SetFirstDirection(0);
262 iter.SetSecondDirection(1);
265 MaskImageType::PointType A;
266 MaskImageType::PointType B;
267 MaskImageType::PointType C;
268 while (!iter.IsAtEnd()) {
272 C[1] -= 10; // I know I must keep this point
273 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
274 bool isPositive = s<0;
275 while (!iter.IsAtEndOfSlice()) {
276 while (!iter.IsAtEndOfLine()) {
277 // Very slow, I know ... but image should be very small
278 m_Working_Support->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
279 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
280 if (s == 0) iter.Set(2);
282 if (s > 0) iter.Set(GetBackgroundValue());
285 if (s < 0) iter.Set(GetBackgroundValue());
295 //-----------------------------------------------------
296 // StartNewStep("Anterior limits");
301 // Station 4R, Station 4L, the right pulmonary artery, and/or the
302 // left superior pulmonary vein
305 //-----------------------------------------------------
306 //-----------------------------------------------------
307 // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
311 m_Station7 = m_Working_Support;
313 //--------------------------------------------------------------------