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://www.centreleonberard.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 ===========================================================================**/
20 #include "clitkCommon.h"
21 #include "clitkBooleanOperatorLabelImageFilter.h"
22 #include "clitkAutoCropFilter.h"
23 #include "clitkResampleImageWithOptionsFilter.h"
24 #include "clitkBooleanOperatorLabelImageFilter.h"
25 #include "clitkCropLikeImageFilter.h"
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>
39 #include <itkDivideByConstantImageFilter.h>
43 #include "RelativePositionPropImageFilter.h"
45 //--------------------------------------------------------------------
46 template <class ImageType>
47 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
48 AddRelativePositionConstraintToLabelImageFilter():
50 itk::ImageToImageFilter<ImageType, ImageType>()
52 this->SetNumberOfRequiredInputs(2);
53 SetFuzzyThreshold(0.6);
54 SetBackgroundValue(0);
55 SetObjectBackgroundValue(0);
56 ClearOrientationType();
57 IntermediateSpacingFlagOn();
58 SetIntermediateSpacing(10);
60 InverseOrientationFlagOff();
62 CombineWithOrFlagOff();
65 FuzzyMapOnlyFlagOff();
68 SetK1(vcl_acos(-1.0)/2);
70 //--------------------------------------------------------------------
73 //--------------------------------------------------------------------
74 template <class ImageType>
76 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
77 SetInput(const ImageType * image)
79 // Process object is not const-correct so the const casting is required.
80 this->SetNthInput(0, const_cast<ImageType *>(image));
82 //--------------------------------------------------------------------
85 //--------------------------------------------------------------------
86 template <class ImageType>
88 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
89 SetInputObject(const ImageType * image)
91 // Process object is not const-correct so the const casting is required.
92 this->SetNthInput(1, const_cast<ImageType *>(image));
94 //--------------------------------------------------------------------
97 //--------------------------------------------------------------------
98 template <class ImageType>
100 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
101 ClearOrientationType()
103 m_OrientationTypeString.clear();
104 m_OrientationType.clear();
108 //--------------------------------------------------------------------
111 //--------------------------------------------------------------------
112 template <class ImageType>
114 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
117 return m_OrientationType.size();
119 //--------------------------------------------------------------------
122 //--------------------------------------------------------------------
123 template <class ImageType>
125 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
126 AddOrientationTypeString(std::string t)
128 m_OrientationTypeString.push_back(t);
130 if (t == "LeftTo") { AddOrientationType(LeftTo); return; }
131 if (t == "RightTo") { AddOrientationType(RightTo); return; }
132 if (t == "AntTo") { AddOrientationType(AntTo); return; }
133 if (t == "PostTo") { AddOrientationType(PostTo); return; }
134 if (t == "SupTo") { AddOrientationType(SupTo); return; }
135 if (t == "InfTo") { AddOrientationType(InfTo); return; }
137 if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; }
138 if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; }
139 if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; }
140 if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; }
141 if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; }
142 if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; }
144 if (t == "Angle") return;
146 clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t);
148 //--------------------------------------------------------------------
151 //--------------------------------------------------------------------
152 template <class ImageType>
154 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
155 GenerateOutputInformation()
157 ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
158 ImagePointer outputImage = this->GetOutput(0);
159 outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
161 //--------------------------------------------------------------------
164 //--------------------------------------------------------------------
165 template <class ImageType>
167 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
168 GenerateInputRequestedRegion()
171 itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
172 // Get input pointers and set requested region to common region
173 ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
174 ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
175 input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
176 input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
178 //--------------------------------------------------------------------
181 //--------------------------------------------------------------------
182 template <class ImageType>
184 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
185 AddAnglesInRad(double a, double b)
187 m_OrientationTypeString.push_back("Angle");
188 m_OrientationType.push_back(Angle);
189 m_Angle1.push_back(a);
190 m_Angle2.push_back(b);
192 //--------------------------------------------------------------------
195 //--------------------------------------------------------------------
196 template <class ImageType>
198 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
199 AddAnglesInDeg(double a, double b)
201 AddAnglesInRad(clitk::deg2rad(a), clitk::deg2rad(b));
203 //--------------------------------------------------------------------
206 //--------------------------------------------------------------------
207 template <class ImageType>
209 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
210 AddOrientationType(OrientationTypeEnumeration orientation)
212 m_OrientationType.push_back(orientation);
213 switch (orientation) {
215 m_Angle1.push_back(clitk::deg2rad(0));
216 m_Angle2.push_back(clitk::deg2rad(0));
219 m_Angle1.push_back(clitk::deg2rad(180));
220 m_Angle2.push_back(clitk::deg2rad(0));
223 m_Angle1.push_back(clitk::deg2rad(90));
224 m_Angle2.push_back(clitk::deg2rad(0));
227 m_Angle1.push_back(clitk::deg2rad(-90));
228 m_Angle2.push_back(clitk::deg2rad(0));
231 m_Angle1.push_back(clitk::deg2rad(0));
232 m_Angle2.push_back(clitk::deg2rad(90));
235 m_Angle1.push_back(clitk::deg2rad(0));
236 m_Angle2.push_back(clitk::deg2rad(-90));
249 //--------------------------------------------------------------------
252 //--------------------------------------------------------------------
253 template <class ImageType>
255 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
258 DD((int)this->GetBackgroundValue());
259 DD((int)this->GetObjectBackgroundValue());
260 DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles());
261 DD(this->GetIntermediateSpacingFlag());
262 DD(this->GetIntermediateSpacing());
263 DD(this->GetFuzzyThreshold());
264 DD(this->GetAutoCropFlag());
265 DD(this->GetInverseOrientationFlag());
266 DD(this->GetRemoveObjectFlag());
267 DD(this->GetCombineWithOrFlag());
269 //--------------------------------------------------------------------
272 //--------------------------------------------------------------------
273 template <class ImageType>
275 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
278 if (GetNumberOfAngles() <1) {
279 clitkExceptionMacro("Add at least one orientation type");
282 if (GetVerboseOptionFlag()) {
283 for(int i=0; i<GetNumberOfAngles(); i++) {
284 std::cout << "Orientation \t:" << GetOrientationTypeString(i) << std::endl;
286 std::cout << "Interm Spacing \t:" << GetIntermediateSpacingFlag() << " " << GetIntermediateSpacing() << "mm" << std::endl;
287 std::cout << "Fuzzy threshold\t:" << GetFuzzyThreshold() << std::endl;
288 std::cout << "AutoCrop \t:" << GetAutoCropFlag() << std::endl;
289 std::cout << "InverseOrient \t:" << GetInverseOrientationFlag() << std::endl;
290 std::cout << "RemoveObject \t:" << GetRemoveObjectFlag() << std::endl;
291 std::cout << "CombineWithOr \t:" << GetCombineWithOrFlag() << std::endl;
295 input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
296 object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
297 static const unsigned int dim = ImageType::ImageDimension;
299 // Step 2: object pad to input image -> we want to compute the
300 // relative position for each point belonging to the input image
301 // domain, so we have to extend (pad) the object image to fit the
303 working_image = object;
304 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
305 StartNewStep("Pad (resize) object to input size");
307 if (0) { // OLD VERSION (TO REMOVE)
308 StartNewStep("Pad object to image size");
309 typename ImageType::Pointer output = ImageType::New();
311 for(unsigned int i=0; i<dim; i++) {
312 size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
313 input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
316 // The index of the input is not necessarily zero, so we have to
317 // take it into account (not done)
319 IndexType index = input->GetLargestPossibleRegion().GetIndex();
320 region.SetSize(size);
321 for(unsigned int i=0; i<dim; i++) {
323 std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
327 // output->SetLargestPossibleRegion(region);
328 output->SetRegions(region);
329 output->SetSpacing(working_image->GetSpacing());
330 PointType origin = input->GetOrigin();
331 for(unsigned int i=0; i<dim; i++) {
332 origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
334 output->SetOrigin(origin);
335 // output->SetOrigin(input->GetOrigin());
338 output->FillBuffer(m_BackgroundValue);
339 typename PasteFilterType::Pointer padFilter = PasteFilterType::New();
340 // typename PasteFilterType::InputImageIndexType index;
341 for(unsigned int i=0; i<dim; i++) {
342 index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
343 + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
345 padFilter->SetSourceImage(working_image);
346 padFilter->SetDestinationImage(output);
347 padFilter->SetDestinationIndex(index);
348 padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
350 working_image = padFilter->GetOutput();
353 // Resize object like input
354 working_image = clitk::ResizeImageLike<ImageType>(working_image, input, GetBackgroundValue());
355 StopCurrentStep<ImageType>(working_image);
358 //--------------------------------------------------------------------
359 //--------------------------------------------------------------------
361 if (m_IntermediateSpacingFlag) {
362 StartNewStep("Resample object to intermediate spacing (" + toString(m_IntermediateSpacing) + ")");
363 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
364 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
365 resampleFilter->SetInput(working_image);
366 resampleFilter->SetDefaultPixelValue(0);
367 resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
368 resampleFilter->SetGaussianFilteringEnabled(false);
369 //resampleFilter->SetVerboseOptions(true);
370 resampleFilter->Update();
371 working_image = resampleFilter->GetOutput();
372 StopCurrentStep<ImageType>(working_image);
375 // Keep object image (with resampline and pad)
376 object_resampled = working_image;
378 // Step 3: compute rel pos in object
379 StartNewStep("Relative Position Map");
380 typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
381 typename RelPosFilterType::Pointer relPosFilter;
383 for(int i=0; i<GetNumberOfAngles(); i++) {
385 relPosFilter = RelPosFilterType::New();
386 relPosFilter->SetFast(GetFastFlag());
387 relPosFilter->SetRadius(GetRadius());
388 relPosFilter->SetInput(working_image);
389 relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
390 relPosFilter->SetAlpha2(m_Angle2[i]);
391 relPosFilter->SetK1(GetK1());// M_PI/2.0); // Opening parameter, default = pi/2
392 // relPosFilter->SetVerboseProgress(true);
394 relPosFilter->Update();
395 relPos = relPosFilter->GetOutput();
397 if (GetNumberOfAngles() != 1) {
398 // Creation of the first m_FuzzyMap
400 m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
401 m_FuzzyMap->FillBuffer(0.0);
404 // Add to current fuzzy map
405 typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
406 typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
407 addFilter->SetInput1(m_FuzzyMap);
408 addFilter->SetInput2(relPos);
410 m_FuzzyMap = addFilter->GetOutput();
412 else m_FuzzyMap = relPos;
415 // Divide by the number of relpos
416 if (GetNumberOfAngles() != 1) {
417 #if ITK_VERSION_MAJOR >= 4
418 typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
419 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
420 divideFilter->SetConstant2(GetNumberOfAngles());
422 typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
423 typename DivideFilter::Pointer divideFilter = DivideFilter::New();
424 divideFilter->SetConstant(GetNumberOfAngles());
426 divideFilter->SetInput(m_FuzzyMap);
427 divideFilter->Update();
428 m_FuzzyMap = divideFilter->GetOutput();
432 StopCurrentStep<FloatImageType>(relPos);
433 if (GetFuzzyMapOnlyFlag()) {
434 // If the spacing is used, recompute fuzzy map with the input size/spacing
435 if (m_IntermediateSpacingFlag) {
436 StartNewStep("Resample FuzzyMap to come back to the same sampling than input");
437 typedef clitk::ResampleImageWithOptionsFilter<FloatImageType> ResampleFilterType;
438 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
439 resampleFilter->SetDefaultPixelValue(0.0); //Default fuzzymap value FIXME
440 resampleFilter->SetInput(relPos);
441 resampleFilter->SetOutputSpacing(input->GetSpacing());
442 resampleFilter->SetGaussianFilteringEnabled(false);
443 resampleFilter->Update();
444 relPos = m_FuzzyMap = resampleFilter->GetOutput();
445 StopCurrentStep<FloatImageType>(relPos);
447 // Need to put exactly the same size
448 if (relPos->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
449 StartNewStep("Pad to get the same size than input");
450 typename FloatImageType::Pointer temp = FloatImageType::New();
451 temp->CopyInformation(input);
452 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
454 temp->FillBuffer(0.0); // Default fuzzymap value FIXME
455 typename PasteFloatFilterType::Pointer padFilter2 = PasteFloatFilterType::New();
456 padFilter2->SetSourceImage(relPos);
457 padFilter2->SetDestinationImage(temp);
458 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
459 padFilter2->SetSourceRegion(relPos->GetLargestPossibleRegion());
460 padFilter2->Update();
461 relPos = padFilter2->GetOutput();
462 StopCurrentStep<FloatImageType>(relPos);
467 StopCurrentStep<FloatImageType>(relPos);
472 //--------------------------------------------------------------------
473 //--------------------------------------------------------------------
474 StartNewStep("Map Threshold");
476 typedef itk::BinaryThresholdImageFilter<FloatImageType, ImageType> BinaryThresholdImageFilterType;
477 typename BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();
478 thresholdFilter->SetInput(relPos);
479 thresholdFilter->SetOutsideValue(m_BackgroundValue);
480 thresholdFilter->SetInsideValue(m_BackgroundValue+1);
481 thresholdFilter->SetLowerThreshold(m_FuzzyThreshold);
482 thresholdFilter->Update();
483 working_image = thresholdFilter->GetOutput();
484 StopCurrentStep<ImageType>(working_image);
486 //--------------------------------------------------------------------
487 //--------------------------------------------------------------------
488 StartNewStep("Post Processing: erosion with initial mask");
489 // Step 2 : erosion with initial mask to exclude pixels that were
490 // inside the resampled version and outside the original mask
491 typedef itk::BinaryBallStructuringElement<unsigned int, ImageDimension> StructuringElementType;
492 StructuringElementType kernel;
494 kernel.CreateStructuringElement();
495 typedef itk::BinaryErodeImageFilter<ImageType, ImageType, StructuringElementType> ErodeFilterType;
496 typename ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
497 erodeFilter->SetInput(working_image);
498 erodeFilter->SetKernel(kernel);
499 erodeFilter->SetBackgroundValue(m_BackgroundValue);
500 erodeFilter->SetErodeValue(m_BackgroundValue+1);
501 erodeFilter->Update();
502 working_image = erodeFilter->GetOutput();
503 StopCurrentStep<ImageType>(working_image);
505 //--------------------------------------------------------------------
506 //--------------------------------------------------------------------
507 // Step 5: resample to initial spacing
508 if (m_IntermediateSpacingFlag) {
509 StartNewStep("Resample to come back to the same sampling than input");
510 typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
511 typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
512 resampleFilter->SetDefaultPixelValue(m_BackgroundValue);
513 resampleFilter->SetInput(working_image);
514 resampleFilter->SetOutputSpacing(input->GetSpacing());
515 resampleFilter->SetGaussianFilteringEnabled(false);
516 // resampleFilter->SetVerboseOptions(true);
517 resampleFilter->Update();
518 working_image = resampleFilter->GetOutput();
519 StopCurrentStep<ImageType>(working_image);
522 //--------------------------------------------------------------------
523 //--------------------------------------------------------------------
524 // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
525 //if (!HaveSameSizeAndSpacing(working_image, input)) {
526 if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
527 StartNewStep("Pad to get the same size than input");
528 typename ImageType::Pointer temp = ImageType::New();
529 temp->CopyInformation(input);
530 temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
532 temp->FillBuffer(m_BackgroundValue);
533 typename PasteFilterType::Pointer padFilter2 = PasteFilterType::New();
534 padFilter2->SetSourceImage(working_image);
535 padFilter2->SetDestinationImage(temp);
536 padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
537 padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
538 padFilter2->Update();
539 working_image = padFilter2->GetOutput();
540 StopCurrentStep<ImageType>(working_image);
543 //DD("[debug] Rel Pos : no padding after");
546 //--------------------------------------------------------------------
547 //--------------------------------------------------------------------
548 // Step 6: combine input+thresholded relpos
549 StartNewStep("Combine with initial input (boolean And)");
550 typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
551 typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
552 combineFilter->SetBackgroundValue(m_BackgroundValue);
553 combineFilter->SetBackgroundValue1(m_BackgroundValue);
554 combineFilter->SetBackgroundValue2(m_BackgroundValue);
555 combineFilter->SetForegroundValue(m_BackgroundValue+1);
556 combineFilter->SetInput1(input);
557 combineFilter->SetInput2(working_image);
558 if (GetInverseOrientationFlag())
559 combineFilter->SetOperationType(BoolFilterType::AndNot);
561 if (GetCombineWithOrFlag())
562 combineFilter->SetOperationType(BoolFilterType::Or);
564 combineFilter->SetOperationType(BoolFilterType::And);
566 combineFilter->InPlaceOff(); // Do not modify initial input (!)
567 combineFilter->Update();
568 working_image = combineFilter->GetOutput();
570 // Remove (if needed the object from the support)
571 if (GetRemoveObjectFlag()) {
572 combineFilter = BoolFilterType::New();
573 combineFilter->SetInput1(working_image);
574 combineFilter->SetInput2(object);
575 combineFilter->SetOperationType(BoolFilterType::AndNot);
576 combineFilter->InPlaceOn();
577 combineFilter->Update();
578 working_image = combineFilter->GetOutput();
580 // In the other case, we must *add* the initial object to keep it
581 // but not more than the initial support
583 combineFilter = BoolFilterType::New();
584 combineFilter->SetInput1(working_image);
585 combineFilter->SetInput2(object);
586 combineFilter->SetOperationType(BoolFilterType::Or);
587 combineFilter->InPlaceOn();
588 combineFilter->Update();
589 working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
590 combineFilter = BoolFilterType::New();
591 combineFilter->SetInput1(working_image);
592 combineFilter->SetInput2(input);
593 combineFilter->SetOperationType(BoolFilterType::And);
594 combineFilter->InPlaceOn();
595 combineFilter->Update();
596 working_image = combineFilter->GetOutput();
599 StopCurrentStep<ImageType>(working_image);
601 //--------------------------------------------------------------------
602 //--------------------------------------------------------------------
604 if (GetAutoCropFlag()) {
605 StartNewStep("Final AutoCrop");
606 typedef clitk::AutoCropFilter<ImageType> CropFilterType;
607 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
608 cropFilter->SetInput(working_image);
609 cropFilter->ReleaseDataFlagOff();
610 cropFilter->Update();
611 working_image = cropFilter->GetOutput();
612 StopCurrentStep<ImageType>(working_image);
615 //--------------------------------------------------------------------
616 //--------------------------------------------------------------------
618 // Final Step -> set output
619 this->SetNthOutput(0, working_image);
620 // this->GraftOutput(working_image);
622 //--------------------------------------------------------------------