]> Creatis software - clitk.git/blob - segmentation/clitkExtractLymphStationsFilter.old.txx
Add one slice for the S1RL support
[clitk.git] / segmentation / clitkExtractLymphStationsFilter.old.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
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
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17   ======================================================================-====*/
18
19 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
21
22 // clitk
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"
30
31 // itk
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>
39
40 // itk ENST
41 #include "RelativePositionPropImageFilter.h"
42
43 //--------------------------------------------------------------------
44 template <class TImageType>
45 clitk::ExtractLymphStationsFilter<TImageType>::
46 ExtractLymphStationsFilter():
47   clitk::FilterBase(),
48   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
49   itk::ImageToImageFilter<TImageType, MaskImageType>()
50 {
51   this->SetNumberOfRequiredInputs(1);
52   SetBackgroundValue(0);
53   SetForegroundValue(1);
54   SetFuzzyThreshold(0.5);
55 }
56 //--------------------------------------------------------------------
57
58
59 //--------------------------------------------------------------------
60 template <class TImageType>
61 void 
62 clitk::ExtractLymphStationsFilter<TImageType>::
63 GenerateOutputInformation() { 
64   //  Superclass::GenerateOutputInformation();
65   
66   // Get input
67   LoadAFDB();
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);
74   ImagePointType mlb;
75   GetAFDB()->GetPoint3D("middleLobeBronchus", mlb);
76   m_MiddleLobeBronchusZPositionInMM = mlb[2];
77   DD(m_MiddleLobeBronchusZPositionInMM);
78
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;
91
92   // ----------------------------------------------------------------
93   // ----------------------------------------------------------------
94   // Superior limit = carina
95   // Inferior limit = origin right middle lobe bronchus
96   StartNewStep("Inf/Sup trachea limits with carina/bronchus");
97   m_working_trachea = 
98     clitk::CropImageAlongOneAxis<MaskImageType>(m_trachea, 2, 
99                                                 m_MiddleLobeBronchusZPositionInMM, 
100                                                 m_CarinaZPositionInMM, true,
101                                                 GetBackgroundValue());
102   StopCurrentStep<MaskImageType>(m_working_trachea);
103
104   // ----------------------------------------------------------------
105   // ----------------------------------------------------------------
106   // Separate trachea in two CCL
107   StartNewStep("Separate trachea under carina");
108
109   // Labelize and consider two main labels
110   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
111
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)
119   int maxLabel=0;
120   while (!iter.IsAtReverseEndOfSlice()) {
121     while (!iter.IsAtReverseEndOfLine()) {    
122       //  DD(iter.GetIndex());
123       if (iter.Get() > maxLabel) maxLabel = iter.Get();
124       --iter;
125     }
126     iter.PreviousLine();
127   }
128   if (maxLabel < 2) {
129     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
130   }
131
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();
146
147   ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
148   ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
149
150   ImagePixelType leftLabel;
151   ImagePixelType rightLabel;  
152   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
153   else { leftLabel = 2; rightLabel = 1; }
154   DD(leftLabel);
155   DD(rightLabel);
156
157   StopCurrentStep<MaskImageType>(m_working_trachea);
158
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");
170
171   m_working_mediastinum = 
172     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
173                                                        leftBronchus, 2, 
174                                                        GetFuzzyThreshold(), "RightTo", 
175                                                        true, 4);
176   m_working_mediastinum = 
177     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
178                                                        rightBronchus, 
179                                                        2, GetFuzzyThreshold(), "LeftTo", 
180                                                        true, 4);
181   m_working_mediastinum = 
182     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
183                                                        leftBronchus, 
184                                                        2, GetFuzzyThreshold(), "AntTo", 
185                                                        true, 4, true); // NOT
186   m_working_mediastinum = 
187     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
188                                                        rightBronchus, 
189                                                        2, GetFuzzyThreshold(), "AntTo", 
190                                                        true, 4, true);
191   m_working_mediastinum = 
192     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
193                                                        leftBronchus, 
194                                                        2, GetFuzzyThreshold(), "PostTo", 
195                                                        true, 4, true);
196   m_working_mediastinum = 
197     clitk::SliceBySliceRelativePosition<MaskImageType>(m_working_mediastinum, 
198                                                        rightBronchus, 
199                                                        2, GetFuzzyThreshold(), "PostTo", 
200                                                        true, 4, true);
201   m_output = m_working_mediastinum;
202   StopCurrentStep<MaskImageType>(m_output);
203
204   //-----------------------------------------------------
205   //-----------------------------------------------------
206   StartNewStep("Posterior limits");  
207
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);
218   
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);
226   
227   assert(leftBronchusSlices.size() == rightBronchusSlices.size());
228   
229   std::vector<MaskImageType::PointType> leftPoints;
230   std::vector<MaskImageType::PointType> rightPoints;
231   for(uint i=0; i<leftBronchusSlices.size(); i++) {
232     // Keep main CCL
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;
242
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(), 
249                                                           0, false, a, 0);
250     // find post most point in left, not far away from rightMost
251     SliceType::PointType p = 
252       clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i], 
253                                                           GetBackgroundValue(), 
254                                                           1, false, rightMost, 
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);
260
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(), 
266                                                           0, true, a, 0);
267     // find post most point in left, not far away from leftMost
268     p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i], 
269                                                             GetBackgroundValue(), 
270                                                             1, false, leftMost, 
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);
275   }
276
277   // DEBUG
278   std::ofstream osl;
279   openFileForWriting(osl, "left.txt");
280   osl << "LANDMARKS1" << std::endl;
281   std::ofstream osr;
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;
289   }
290   osl.close();
291   osr.close();
292
293   // Now uses these points to limit, slice by slice 
294   // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
295   /*
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.
299
300     => to accelerate, start with formula, when change sign -> stop and fill
301   */
302   //  typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
303   iter = SliceIteratorType(m_working_mediastinum, m_working_mediastinum->GetLargestPossibleRegion());
304   iter.SetFirstDirection(0);
305   iter.SetSecondDirection(1);
306   iter.GoToBegin();
307   int i=0;
308   MaskImageType::PointType A;
309   MaskImageType::PointType B;
310   MaskImageType::PointType C;
311   while (!iter.IsAtEnd()) {
312     A = leftPoints[i];
313     B = rightPoints[i];
314     C = A;
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);
324         if (isPositive) {
325           if (s > 0) iter.Set(GetBackgroundValue());
326         }
327         else {
328           if (s < 0) iter.Set(GetBackgroundValue());
329         }
330         ++iter;
331       }
332       iter.NextLine();
333     }
334     iter.NextSlice();
335     ++i;
336   }
337
338   //-----------------------------------------------------
339   //-----------------------------------------------------
340   // StartNewStep("Anterior limits");  
341  
342
343   // MISSING FROM NOW 
344   
345   // Station 4R, Station 4L, the right pulmonary artery, and/or the
346   // left superior pulmonary vein
347
348
349   //-----------------------------------------------------
350   //-----------------------------------------------------
351   // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
352
353   
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());
359
360 }
361 //--------------------------------------------------------------------
362
363
364 //--------------------------------------------------------------------
365 template <class TImageType>
366 void 
367 clitk::ExtractLymphStationsFilter<TImageType>::
368 GenerateInputRequestedRegion() {
369   DD("GenerateInputRequestedRegion");
370   // Call default
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));
375   DD("set reg");
376   m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion());
377   m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion());
378   DD("end");
379 }
380 //--------------------------------------------------------------------
381
382
383 //--------------------------------------------------------------------
384 template <class TImageType>
385 void 
386 clitk::ExtractLymphStationsFilter<TImageType>::
387 GenerateData() {
388   DD("GenerateData");
389   // Final Step -> graft output (if SetNthOutput => redo)
390   this->GraftOutput(m_output);
391 }
392 //--------------------------------------------------------------------
393   
394
395 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX