+ typename RelPosFilterType::Pointer relPosFilter;
+
+ for(int i=0; i<GetNumberOfAngles(); i++) {
+ // Compute fuzzy map
+ relPosFilter = RelPosFilterType::New();
+ relPosFilter->SetFast(GetFastFlag());
+ relPosFilter->SetRadius(GetRadius());
+ relPosFilter->SetInput(working_image);
+ relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
+ relPosFilter->SetAlpha2(m_Angle2[i]);
+ relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
+
+ // relPosFilter->SetFast(true);
+ // relPosFilter->SetRadius(1); // seems sufficient in this case
+
+ // relPosFilter->SetVerboseProgress(true);
+ relPosFilter->Update();
+ relPos = relPosFilter->GetOutput();
+
+ if (GetNumberOfAngles() != 1) {
+ // Creation of the first m_FuzzyMap
+ if (i==0) {
+ m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
+ m_FuzzyMap->FillBuffer(0.0);
+ }
+
+ // Add to current fuzzy map
+ typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
+ typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
+ addFilter->SetInput1(m_FuzzyMap);
+ addFilter->SetInput2(relPos);
+ addFilter->Update();
+ m_FuzzyMap = addFilter->GetOutput();
+ }
+ else m_FuzzyMap = relPos;
+ }
+
+ // Divide by the number of relpos
+ if (GetNumberOfAngles() != 1) {
+#if ITK_VERSION_MAJOR >= 4
+ typedef itk::DivideImageFilter<FloatImageType, FloatImageType, FloatImageType> DivideFilter;
+ typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+ divideFilter->SetConstant2(GetNumberOfAngles());
+#else
+ typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
+ typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+ divideFilter->SetConstant(GetNumberOfAngles());
+#endif
+ divideFilter->SetInput(m_FuzzyMap);
+ divideFilter->Update();
+ m_FuzzyMap = divideFilter->GetOutput();
+ }
+
+ relPos = m_FuzzyMap;
+ StopCurrentStep<FloatImageType>(relPos);
+ if (GetFuzzyMapOnlyFlag()) {
+ // If the spacing is used, recompute fuzzy map with the input size/spacing
+ if (m_IntermediateSpacingFlag) {
+ StartNewStep("Resample FuzzyMap to come back to the same sampling than input");
+ typedef clitk::ResampleImageWithOptionsFilter<FloatImageType> ResampleFilterType;
+ typename ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
+ resampleFilter->SetDefaultPixelValue(0.0); //Default fuzzymap value FIXME
+ resampleFilter->SetInput(relPos);
+ resampleFilter->SetOutputSpacing(input->GetSpacing());
+ resampleFilter->SetGaussianFilteringEnabled(false);
+ resampleFilter->Update();
+ relPos = m_FuzzyMap = resampleFilter->GetOutput();
+
+ // Need to put exactly the same size
+ if (relPos->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
+ StartNewStep("Pad to get the same size than input");
+ typename FloatImageType::Pointer temp = FloatImageType::New();
+ temp->CopyInformation(input);
+ temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
+ temp->Allocate();
+ temp->FillBuffer(0.0); // Default fuzzymap value FIXME
+ typename PasteFloatFilterType::Pointer padFilter2 = PasteFloatFilterType::New();
+ padFilter2->SetSourceImage(relPos);
+ padFilter2->SetDestinationImage(temp);
+ padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
+ padFilter2->SetSourceRegion(relPos->GetLargestPossibleRegion());
+ padFilter2->Update();
+ relPos = padFilter2->GetOutput();
+ StopCurrentStep<FloatImageType>(relPos);
+ m_FuzzyMap = relPos;
+ }
+ }
+ else {
+ StopCurrentStep<FloatImageType>(relPos);
+ }
+ return;
+ }