]> Creatis software - clitk.git/blob - itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
correction orientation, add 'keep object' option
[clitk.git] / itk / clitkAddRelativePositionConstraintToLabelImageFilter.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 // clitk
20 #include "clitkCommon.h"
21 #include "clitkBooleanOperatorLabelImageFilter.h"
22 #include "clitkAutoCropFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
24 #include "clitkBooleanOperatorLabelImageFilter.h"
25
26 // itk
27 #include <deque>
28 #include <itkStatisticsLabelMapFilter.h>
29 #include <itkLabelImageToStatisticsLabelMapFilter.h>
30 #include <itkRegionOfInterestImageFilter.h>
31 #include <itkBinaryThresholdImageFilter.h>
32 #include <itkBinaryErodeImageFilter.h>
33 #include <itkBinaryBallStructuringElement.h>
34 #include <itkAddImageFilter.h>
35 #include <itkDivideByConstantImageFilter.h>
36
37 // itk [Bloch et al] 
38 #include "RelativePositionPropImageFilter.h"
39
40 //--------------------------------------------------------------------
41 template <class ImageType>
42 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
43 AddRelativePositionConstraintToLabelImageFilter():
44   clitk::FilterBase(),
45   itk::ImageToImageFilter<ImageType, ImageType>()
46 {
47   this->SetNumberOfRequiredInputs(2);
48   SetFuzzyThreshold(0.6);
49   SetBackgroundValue(0);
50   SetObjectBackgroundValue(0);
51   ClearOrientationType();
52   ResampleBeforeRelativePositionFilterOn();
53   SetIntermediateSpacing(10);
54   AutoCropFlagOn();
55   InverseOrientationFlagOff();
56   RemoveObjectFlagOn();
57   CombineWithOrFlagOff();
58 }
59 //--------------------------------------------------------------------
60
61
62 //--------------------------------------------------------------------
63 template <class ImageType>
64 void 
65 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
66 SetInput(const ImageType * image) 
67 {
68   // Process object is not const-correct so the const casting is required.
69   this->SetNthInput(0, const_cast<ImageType *>(image));
70 }
71 //--------------------------------------------------------------------
72   
73
74 //--------------------------------------------------------------------
75 template <class ImageType>
76 void 
77 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
78 SetInputObject(const ImageType * image) 
79 {
80   // Process object is not const-correct so the const casting is required.
81   this->SetNthInput(1, const_cast<ImageType *>(image));
82 }
83 //--------------------------------------------------------------------
84   
85
86 //--------------------------------------------------------------------
87 template <class ImageType>
88 void 
89 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
90 ClearOrientationType() 
91 {
92   m_OrientationTypeString.clear();
93   m_OrientationType.clear();
94   m_Angle1.clear();
95   m_Angle2.clear();
96 }
97 //--------------------------------------------------------------------
98
99
100 //--------------------------------------------------------------------
101 template <class ImageType>
102 int
103 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
104 GetNumberOfAngles()
105 {
106   return m_OrientationType.size();
107 }
108 //--------------------------------------------------------------------
109
110
111 //--------------------------------------------------------------------
112 template <class ImageType>
113 void 
114 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
115 AddOrientationTypeString(std::string t) 
116 {
117   m_OrientationTypeString.push_back(t);
118   switch (t[0]) {
119   case 'L' : AddOrientationType(LeftTo); break;
120   case 'R' : AddOrientationType(RightTo);break;
121   case 'A' : AddOrientationType(AntTo);break;
122   case 'P' : AddOrientationType(PostTo);break;
123   case 'S' : AddOrientationType(SupTo);break;
124   case 'I' : AddOrientationType(InfTo);break;
125   default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I");
126   }
127 }
128 //--------------------------------------------------------------------
129   
130
131 //--------------------------------------------------------------------
132 template <class ImageType>
133 void 
134 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
135 GenerateOutputInformation() 
136
137   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
138   ImagePointer outputImage = this->GetOutput(0);
139   outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
140 }
141 //--------------------------------------------------------------------
142
143
144 //--------------------------------------------------------------------
145 template <class ImageType>
146 void 
147 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
148 GenerateInputRequestedRegion() 
149 {
150   // Call default
151   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
152   // Get input pointers and set requested region to common region
153   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
154   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
155   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
156   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
157 }
158 //--------------------------------------------------------------------
159
160   
161 //--------------------------------------------------------------------
162 template <class ImageType>
163 void 
164 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
165 AddAngles(double a, double b) 
166 {
167   AddOrientationTypeString("Angle");
168   m_Angle1.push_back(a);
169   m_Angle2.push_back(b);
170 }
171 //--------------------------------------------------------------------
172
173
174 //--------------------------------------------------------------------
175 template <class ImageType>
176 void 
177 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
178 AddOrientationType(OrientationTypeEnumeration orientation) 
179 {
180   m_OrientationType.push_back(orientation);
181   switch (orientation) {
182   case LeftTo:   
183     m_Angle1.push_back(clitk::deg2rad(0));   
184     m_Angle2.push_back(clitk::deg2rad(0));
185     break;
186   case RightTo:  
187     m_Angle1.push_back(clitk::deg2rad(180)); 
188     m_Angle2.push_back(clitk::deg2rad(0));
189     break;
190   case AntTo:
191     m_Angle1.push_back(clitk::deg2rad(90));
192     m_Angle2.push_back(clitk::deg2rad(0));
193     break;
194   case PostTo:
195     m_Angle1.push_back(clitk::deg2rad(-90)); 
196     m_Angle2.push_back(clitk::deg2rad(0));
197     break;
198   case InfTo:    
199     m_Angle1.push_back(clitk::deg2rad(0));   
200     m_Angle2.push_back(clitk::deg2rad(90));
201     break;
202   case SupTo:    
203     m_Angle1.push_back(clitk::deg2rad(0));   
204     m_Angle2.push_back(clitk::deg2rad(-90));
205     break;
206   case Angle:  break;
207   }
208   /*         A1   A2
209              Left      0    0
210              Right   180    0
211              Ant      90    0
212              Post    -90    0
213              Inf       0   90
214              Sup       0  -90
215   */
216 }
217 //--------------------------------------------------------------------
218
219
220 //--------------------------------------------------------------------
221 template <class ImageType>
222 void 
223 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
224 GenerateData() 
225 {
226   if (GetNumberOfAngles() <1) {
227     clitkExceptionMacro("Add at least one orientation type");
228   }  
229
230   // Get input pointer
231   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
232   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
233
234   //--------------------------------------------------------------------
235   //--------------------------------------------------------------------
236   static const unsigned int dim = ImageType::ImageDimension;
237   StartNewStep("Initial resample");  
238   // Step 1 : resample
239   if (m_ResampleBeforeRelativePositionFilter) {
240     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
241     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
242     resampleFilter->SetInput(object);
243     resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
244     resampleFilter->SetGaussianFilteringEnabled(false);
245     //    resampleFilter->SetVerboseOptions(true);
246     resampleFilter->Update();
247     working_image = resampleFilter->GetOutput();
248   }
249   else {
250     working_image = object;
251   }
252   StopCurrentStep<ImageType>(working_image);
253
254   // Step 2: object pad to input image -> we want to compute the
255   // relative position for each point belonging to the input image
256   // domain, so we have to extend (pad) the object image to fit the
257   // domain size
258   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
259     StartNewStep("Pad object to image size");  
260     typename ImageType::Pointer output = ImageType::New();
261     SizeType size;
262     for(unsigned int i=0; i<dim; i++) {
263       size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
264                        input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
265     }
266
267     // The index of the input is not necessarily zero, so we have to
268     // take it into account (not done)
269     RegionType region;
270     IndexType index = input->GetLargestPossibleRegion().GetIndex();
271     region.SetSize(size);
272     for(unsigned int i=0; i<dim; i++) {
273       if (index[i] != 0) {
274         std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
275         exit(0);
276       }
277     }
278     // output->SetLargestPossibleRegion(region);
279     output->SetRegions(region);
280     output->SetSpacing(working_image->GetSpacing());    
281     PointType origin = input->GetOrigin();
282     for(unsigned int i=0; i<dim; i++) {
283       origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
284     }
285     output->SetOrigin(origin);
286     //    output->SetOrigin(input->GetOrigin());
287
288     output->Allocate();
289     output->FillBuffer(m_BackgroundValue);
290     typename PadFilterType::Pointer padFilter = PadFilterType::New();
291     // typename PadFilterType::InputImageIndexType index;
292     for(unsigned int i=0; i<dim; i++) {
293       index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
294         + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
295     }
296     padFilter->SetSourceImage(working_image);
297     padFilter->SetDestinationImage(output);
298     padFilter->SetDestinationIndex(index);
299     padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
300     padFilter->Update();
301     working_image = padFilter->GetOutput();
302     StopCurrentStep<ImageType>(working_image);
303   }
304   else {
305     // DD("[debug] RelPos : same size and spacing : no padding");
306   }
307   // Keep object image (with resampline and pad)
308   object_resampled = working_image;
309   //  StopCurrentStep<ImageType>(working_image);
310
311   // Step 3: compute rel pos in object
312   StartNewStep("Relative Position Map");  
313   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
314   typename RelPosFilterType::Pointer relPosFilter;
315
316   typename FloatImageType::Pointer m_FuzzyMap;
317   for(int i=0; i<GetNumberOfAngles(); i++) {
318     // Compute fuzzy map
319     relPosFilter = RelPosFilterType::New();
320     relPosFilter->SetInput(working_image);
321     relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
322     relPosFilter->SetAlpha2(m_Angle2[i]);
323     relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
324     relPosFilter->SetFast(true);
325     relPosFilter->SetRadius(1); // seems sufficient in this case
326     // relPosFilter->SetVerboseProgress(true);
327     relPosFilter->Update();
328     relPos = relPosFilter->GetOutput();
329
330     if (GetNumberOfAngles() != 1) {
331       // Creation of the first m_FuzzyMap
332       if (i==0) {
333         m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
334         m_FuzzyMap->FillBuffer(0.0);
335       }
336       
337       // Add to current fuzzy map
338       typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
339       typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
340       addFilter->SetInput1(m_FuzzyMap);
341       addFilter->SetInput2(relPos);
342       addFilter->Update();
343       m_FuzzyMap = addFilter->GetOutput();
344     }
345     else m_FuzzyMap = relPos;
346   }
347
348   // Divide by the number of relpos
349   if (GetNumberOfAngles() != 1) {
350     typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
351     typename DivideFilter::Pointer divideFilter = DivideFilter::New();
352     divideFilter->SetInput(m_FuzzyMap);
353     divideFilter->SetConstant(GetNumberOfAngles());
354     divideFilter->Update();
355     m_FuzzyMap = divideFilter->GetOutput();
356   }
357
358   relPos = m_FuzzyMap;
359   StopCurrentStep<FloatImageType>(relPos);
360                
361   //--------------------------------------------------------------------
362   //--------------------------------------------------------------------
363   StartNewStep("Map Threshold");
364   // Step 1: threshold
365   typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
366   typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
367   thresholdFilter->SetInput(relPos);
368   thresholdFilter->SetOutsideValue(m_BackgroundValue);
369   thresholdFilter->SetInsideValue(m_BackgroundValue+1);
370   thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
371   thresholdFilter->Update();
372   working_image = thresholdFilter->GetOutput();
373   StopCurrentStep<ImageType>(working_image);
374
375   //--------------------------------------------------------------------
376   //--------------------------------------------------------------------
377   StartNewStep("Post Processing: erosion with initial mask");
378   // Step 2 : erosion with initial mask to exclude pixels that were
379   // inside the resampled version and outside the original mask
380   typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType; 
381   StructuringElementType kernel;
382   kernel.SetRadius(1);
383   kernel.CreateStructuringElement();
384   typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
385   typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
386   erodeFilter->SetInput(working_image);
387   erodeFilter->SetKernel(kernel);
388   erodeFilter->SetBackgroundValue(m_BackgroundValue);
389   erodeFilter->SetErodeValue(m_BackgroundValue+1);
390   erodeFilter->Update();
391   working_image = erodeFilter->GetOutput();
392   StopCurrentStep<ImageType>(working_image);
393
394   //--------------------------------------------------------------------
395   //--------------------------------------------------------------------
396   // Step 5: resample to initial spacing
397   if (m_ResampleBeforeRelativePositionFilter) {
398     StartNewStep("Resample to get the same sampling than input");
399     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
400     typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
401     resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
402     resampleFilter->SetInput(working_image);
403     resampleFilter->SetOutputSpacing(input->GetSpacing());
404     resampleFilter->SetGaussianFilteringEnabled(false);
405     // resampleFilter->SetVerboseOptions(true);
406     resampleFilter->Update();
407     working_image = resampleFilter->GetOutput();
408     StopCurrentStep<ImageType>(working_image);
409   }
410
411   //--------------------------------------------------------------------
412   //--------------------------------------------------------------------
413   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
414   //if (!HaveSameSizeAndSpacing(working_image, input)) {
415   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
416     StartNewStep("Pad to get the same size than input");
417     typename ImageType::Pointer temp = ImageType::New();
418     temp->CopyInformation(input);
419     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
420     temp->Allocate();
421     temp->FillBuffer(m_BackgroundValue); 
422     typename PadFilterType::Pointer padFilter2 = PadFilterType::New();
423     padFilter2->SetSourceImage(working_image);
424     padFilter2->SetDestinationImage(temp);
425     // DD(input->GetLargestPossibleRegion().GetIndex());
426     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
427     padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
428     padFilter2->Update();
429     working_image = padFilter2->GetOutput();
430     StopCurrentStep<ImageType>(working_image);
431   }
432   else {
433     //DD("[debug] Rel Pos : no padding after");
434   }
435
436   //--------------------------------------------------------------------
437   //--------------------------------------------------------------------
438   // Step 6: combine input+thresholded relpos
439   StartNewStep("Combine with initial input (boolean And)");
440   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
441   typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
442   combineFilter->SetBackgroundValue(m_BackgroundValue);
443   combineFilter->SetBackgroundValue1(m_BackgroundValue);
444   combineFilter->SetBackgroundValue2(m_BackgroundValue);
445   combineFilter->SetForegroundValue(m_BackgroundValue+1);
446   combineFilter->SetInput1(input);
447   combineFilter->SetInput2(working_image);
448   if (GetInverseOrientationFlag())
449     combineFilter->SetOperationType(BoolFilterType::AndNot);
450   else {
451     if (GetCombineWithOrFlag())
452       combineFilter->SetOperationType(BoolFilterType::Or);
453     else
454       combineFilter->SetOperationType(BoolFilterType::And);
455   }
456   combineFilter->InPlaceOff(); // Do not modify initial input (!)
457   combineFilter->Update(); 
458   working_image = combineFilter->GetOutput();
459
460   // Remove (if needed the object from the support)
461   if (GetRemoveObjectFlag()) {
462     combineFilter = BoolFilterType::New();
463     combineFilter->SetInput1(working_image);
464     combineFilter->SetInput2(object);
465     combineFilter->SetOperationType(BoolFilterType::AndNot);
466     combineFilter->InPlaceOn();
467     combineFilter->Update(); 
468     working_image = combineFilter->GetOutput();
469   }
470   // In the other case, we must *add* the initial object to keep it
471   // but not more than the initial support
472   else { 
473     combineFilter = BoolFilterType::New();
474     combineFilter->SetInput1(working_image);
475     combineFilter->SetInput2(object);
476     combineFilter->SetOperationType(BoolFilterType::Or);
477     combineFilter->InPlaceOn();
478     combineFilter->Update(); 
479     working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
480     combineFilter = BoolFilterType::New();
481     combineFilter->SetInput1(working_image);
482     combineFilter->SetInput2(input);
483     combineFilter->SetOperationType(BoolFilterType::And);
484     combineFilter->InPlaceOn();
485     combineFilter->Update(); 
486     working_image = combineFilter->GetOutput();
487   }
488
489   StopCurrentStep<ImageType>(working_image);
490
491   //--------------------------------------------------------------------
492   //--------------------------------------------------------------------
493   // Step 7: autocrop
494   if (GetAutoCropFlag()) {
495     StartNewStep("Final AutoCrop");
496     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
497     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
498     cropFilter->SetInput(working_image);
499     cropFilter->ReleaseDataFlagOff();
500     cropFilter->Update();   
501     working_image = cropFilter->GetOutput();
502     StopCurrentStep<ImageType>(working_image);
503   }
504
505   //--------------------------------------------------------------------
506   //--------------------------------------------------------------------
507   
508   // Final Step -> set output
509   this->SetNthOutput(0, working_image);
510   //  this->GraftOutput(working_image);
511 }
512 //--------------------------------------------------------------------
513