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