1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
23 #include "clitkCommon.h"
24 #include "clitkExtractLymphStationsFilter.h"
25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
26 #include "clitkSegmentationUtils.h"
27 #include "clitkAutoCropFilter.h"
28 #include "clitkSegmentationUtils.h"
29 #include "clitkSliceBySliceRelativePositionFilter.h"
32 #include <itkStatisticsLabelMapFilter.h>
33 #include <itkLabelImageToStatisticsLabelMapFilter.h>
34 #include <itkRegionOfInterestImageFilter.h>
35 #include <itkBinaryThresholdImageFilter.h>
36 #include <itkImageSliceConstIteratorWithIndex.h>
37 #include <itkImageSliceIteratorWithIndex.h>
38 #include <itkBinaryThinningImageFilter.h>
41 #include "RelativePositionPropImageFilter.h"
43 //--------------------------------------------------------------------
44 template <class TImageType>
45 clitk::ExtractLymphStationsFilter<TImageType>::
46 ExtractLymphStationsFilter():
48 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
49 itk::ImageToImageFilter<TImageType, MaskImageType>()
51 this->SetNumberOfRequiredInputs(1);
52 SetBackgroundValue(0);
53 SetForegroundValue(1);
54 SetFuzzyThreshold(0.5);
56 //--------------------------------------------------------------------
59 //--------------------------------------------------------------------
60 template <class TImageType>
62 clitk::ExtractLymphStationsFilter<TImageType>::
63 GenerateOutputInformation() {
64 // Superclass::GenerateOutputInformation();
68 m_mediastinum = GetAFDB()->template GetImage <MaskImageType>("mediastinum");
69 m_trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
70 ImagePointType carina;
71 GetAFDB()->GetPoint3D("carina", carina);
72 m_CarinaZPositionInMM = carina[2];
73 DD(m_CarinaZPositionInMM);
75 GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
76 m_MiddleLobeBronchusZPositionInMM = mlb[2];
77 DD(m_MiddleLobeBronchusZPositionInMM);
79 // ----------------------------------------------------------------
80 // ----------------------------------------------------------------
81 // Superior limit = carina
82 // Inferior limit = origin right middle lobe bronchus
83 StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
84 m_working_mediastinum =
85 clitk::CropImageAlongOneAxis<MaskImageType>(m_mediastinum, 2,
86 m_MiddleLobeBronchusZPositionInMM,
87 m_CarinaZPositionInMM, true,
88 GetBackgroundValue());
89 StopCurrentStep<MaskImageType>(m_working_mediastinum);
90 m_output = m_working_mediastinum;
92 // ----------------------------------------------------------------
93 // ----------------------------------------------------------------
94 // Superior limit = carina
95 // Inferior limit = origin right middle lobe bronchus
96 StartNewStep("Inf/Sup trachea limits with carina/bronchus");
98 clitk::CropImageAlongOneAxis<MaskImageType>(m_trachea, 2,
99 m_MiddleLobeBronchusZPositionInMM,
100 m_CarinaZPositionInMM, true,
101 GetBackgroundValue());
102 StopCurrentStep<MaskImageType>(m_working_trachea);
104 // ----------------------------------------------------------------
105 // ----------------------------------------------------------------
106 // Separate trachea in two CCL
107 StartNewStep("Separate trachea under carina");
109 // Labelize and consider two main labels
110 m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
112 // Carina position must at the first slice that separate the two main bronchus (not superiorly)
113 // Check that upper slice is composed of at least two labels
114 typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
115 SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
116 iter.SetFirstDirection(0);
117 iter.SetSecondDirection(1);
118 iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
120 while (!iter.IsAtReverseEndOfSlice()) {
121 while (!iter.IsAtReverseEndOfLine()) {
122 // DD(iter.GetIndex());
123 if (iter.Get() > maxLabel) maxLabel = iter.Get();
129 clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
132 // Compute centroid of both parts to identify the left from the right bronchus
133 typedef long LabelType;
134 static const unsigned int Dim = ImageType::ImageDimension;
135 typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
136 typedef itk::LabelMap< LabelObjectType > LabelMapType;
137 typedef itk::LabelImageToLabelMapFilter<MaskImageType, LabelMapType> ImageToMapFilterType;
138 typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
139 typedef itk::ShapeLabelMapFilter<LabelMapType, MaskImageType> ShapeFilterType;
140 typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
141 imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
142 imageToLabelFilter->SetInput(m_working_trachea);
143 statFilter->SetInput(imageToLabelFilter->GetOutput());
144 statFilter->Update();
145 typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
147 ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
148 ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
150 ImagePixelType leftLabel;
151 ImagePixelType rightLabel;
152 if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
153 else { leftLabel = 2; rightLabel = 1; }
157 StopCurrentStep<MaskImageType>(m_working_trachea);
159 //-----------------------------------------------------
160 StartNewStep("Left limits with bronchus (slice by slice)");
161 // Select LeftLabel (set one label to Backgroundvalue)
162 MaskImagePointer leftBronchus =
163 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
164 rightLabel, GetBackgroundValue());
165 MaskImagePointer rightBronchus =
166 SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
167 leftLabel, GetBackgroundValue());
168 writeImage<MaskImageType>(leftBronchus, "left.mhd");
169 writeImage<MaskImageType>(rightBronchus, "right.mhd");
171 m_working_mediastinum =
172 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
174 GetFuzzyThreshold(), "RightTo",
176 m_working_mediastinum =
177 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
179 2, GetFuzzyThreshold(), "LeftTo",
181 m_working_mediastinum =
182 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
184 2, GetFuzzyThreshold(), "AntTo",
185 true, 4, true); // NOT
186 m_working_mediastinum =
187 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
189 2, GetFuzzyThreshold(), "AntTo",
191 m_working_mediastinum =
192 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
194 2, GetFuzzyThreshold(), "PostTo",
196 m_working_mediastinum =
197 clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum,
199 2, GetFuzzyThreshold(), "PostTo",
201 m_output = m_working_mediastinum;
202 StopCurrentStep<MaskImageType>(m_output);
204 //-----------------------------------------------------
205 //-----------------------------------------------------
206 StartNewStep("Posterior limits");
208 // Left Bronchus slices
209 typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
210 typedef typename ExtractSliceFilterType::SliceType SliceType;
211 typename ExtractSliceFilterType::Pointer
212 extractSliceFilter = ExtractSliceFilterType::New();
213 extractSliceFilter->SetInput(leftBronchus);
214 extractSliceFilter->SetDirection(2);
215 extractSliceFilter->Update();
216 std::vector<typename SliceType::Pointer> leftBronchusSlices;
217 extractSliceFilter->GetOutputSlices(leftBronchusSlices);
219 // Right Bronchus slices
220 extractSliceFilter = ExtractSliceFilterType::New();
221 extractSliceFilter->SetInput(rightBronchus);
222 extractSliceFilter->SetDirection(2);
223 extractSliceFilter->Update();
224 std::vector<typename SliceType::Pointer> rightBronchusSlices ;
225 extractSliceFilter->GetOutputSlices(rightBronchusSlices);
227 assert(leftBronchusSlices.size() == rightBronchusSlices.size());
229 std::vector<MaskImageType::PointType> leftPoints;
230 std::vector<MaskImageType::PointType> rightPoints;
231 for(uint i=0; i<leftBronchusSlices.size(); i++) {
233 leftBronchusSlices[i] = Labelize<SliceType>(leftBronchusSlices[i], 0, true, 10);
234 leftBronchusSlices[i] = KeepLabels<SliceType>(leftBronchusSlices[i],
235 GetBackgroundValue(),
236 GetForegroundValue(), 1, 1, true);
237 rightBronchusSlices[i] = Labelize<SliceType>(rightBronchusSlices[i], 0, true, 10);
238 rightBronchusSlices[i] = KeepLabels<SliceType>(rightBronchusSlices[i],
239 GetBackgroundValue(),
240 GetForegroundValue(), 1, 1, true);
241 double distance_max_from_center_point = 15;
243 // ------- Find point in left Bronchus -------
244 // find right most point in left = rightMost
245 SliceType::PointType a;
246 SliceType::PointType rightMost =
247 clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
248 GetBackgroundValue(),
250 // find post most point in left, not far away from rightMost
251 SliceType::PointType p =
252 clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
253 GetBackgroundValue(),
255 distance_max_from_center_point);
256 MaskImageType::PointType pp;
257 pp[0] = p[0]; pp[1] = p[1];
258 pp[2] = i*leftBronchus->GetSpacing()[2] + leftBronchus->GetOrigin()[2];
259 leftPoints.push_back(pp);
261 // ------- Find point in right Bronchus -------
262 // find left most point in right = leftMost
263 SliceType::PointType leftMost =
264 clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
265 GetBackgroundValue(),
267 // find post most point in left, not far away from leftMost
268 p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
269 GetBackgroundValue(),
271 distance_max_from_center_point);
272 pp[0] = p[0]; pp[1] = p[1];
273 pp[2] = i*rightBronchus->GetSpacing()[2] + rightBronchus->GetOrigin()[2];
274 rightPoints.push_back(pp);
279 openFileForWriting(osl, "left.txt");
280 osl << "LANDMARKS1" << std::endl;
282 openFileForWriting(osr, "right.txt");
283 osr << "LANDMARKS1" << std::endl;
284 for(uint i=0; i<leftBronchusSlices.size(); i++) {
285 osl << i << " " << leftPoints[i][0] << " " << leftPoints[i][1]
286 << " " << leftPoints[i][2] << " 0 0 " << std::endl;
287 osr << i << " " << rightPoints[i][0] << " " << rightPoints[i][1]
288 << " " << rightPoints[i][2] << " 0 0 " << std::endl;
293 // Now uses these points to limit, slice by slice
294 // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
296 Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
297 (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
298 This will equal zero if the point C is on the line formed by points A and B, and will have a different sign depending on the side. Which side this is depends on the orientation of your (x,y) coordinates, but you can plug test values for A,B and C into this formula to determine whether negative values are to the left or to the right.
300 => to accelerate, start with formula, when change sign -> stop and fill
302 // typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
303 iter = SliceIteratorType(m_working_mediastinum, m_working_mediastinum->GetLargestPossibleRegion());
304 iter.SetFirstDirection(0);
305 iter.SetSecondDirection(1);
308 MaskImageType::PointType A;
309 MaskImageType::PointType B;
310 MaskImageType::PointType C;
311 while (!iter.IsAtEnd()) {
315 C[1] -= 10; // I know I must keep this point
316 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
317 bool isPositive = s<0;
318 while (!iter.IsAtEndOfSlice()) {
319 while (!iter.IsAtEndOfLine()) {
320 // Very slow, I know ... but image should be very small
321 m_working_mediastinum->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
322 double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
323 if (s == 0) iter.Set(2);
325 if (s > 0) iter.Set(GetBackgroundValue());
328 if (s < 0) iter.Set(GetBackgroundValue());
338 //-----------------------------------------------------
339 //-----------------------------------------------------
340 // StartNewStep("Anterior limits");
345 // Station 4R, Station 4L, the right pulmonary artery, and/or the
346 // left superior pulmonary vein
349 //-----------------------------------------------------
350 //-----------------------------------------------------
351 // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
354 // Set output image information (required)
355 MaskImagePointer outputImage = this->GetOutput(0);
356 outputImage->SetRegions(m_working_mediastinum->GetLargestPossibleRegion());
357 outputImage->SetOrigin(m_working_mediastinum->GetOrigin());
358 outputImage->SetRequestedRegion(m_working_mediastinum->GetLargestPossibleRegion());
361 //--------------------------------------------------------------------
364 //--------------------------------------------------------------------
365 template <class TImageType>
367 clitk::ExtractLymphStationsFilter<TImageType>::
368 GenerateInputRequestedRegion() {
369 DD("GenerateInputRequestedRegion");
371 // Superclass::GenerateInputRequestedRegion();
372 // Following needed because output region can be greater than input (trachea)
373 //ImagePointer mediastinum = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
374 //ImagePointer trachea = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
376 m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
377 m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
380 //--------------------------------------------------------------------
383 //--------------------------------------------------------------------
384 template <class TImageType>
386 clitk::ExtractLymphStationsFilter<TImageType>::
389 // Final Step -> graft output (if SetNthOutput => redo)
390 this->GraftOutput(m_output);
392 //--------------------------------------------------------------------
395 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX