+
+//--------------------------------------------------------------------
+template <class ImageType>
+typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
+clitk::RelativePositionAnalyzerFilter<ImageType>::
+ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
+{
+ typedef clitk::SliceBySliceRelativePositionFilter<ImageType> SliceRelPosFilterType;
+ typedef typename SliceRelPosFilterType::FloatImageType FloatImageType;
+ typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
+ sliceRelPosFilter->VerboseStepFlagOff();
+ sliceRelPosFilter->WriteStepFlagOff();
+ sliceRelPosFilter->SetInput(target);
+ sliceRelPosFilter->SetInputObject(object);
+ sliceRelPosFilter->SetDirection(2);
+ sliceRelPosFilter->SetIntermediateSpacingFlag(false);
+ //sliceRelPosFilter->AddOrientationTypeString(orientation);
+ sliceRelPosFilter->AddAngles(angle, 0.0);
+ sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map
+ // sliceRelPosFilter->PrintOptions();
+ sliceRelPosFilter->Update();
+ typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap();
+
+ // Remove initial object from the fuzzy map
+ map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
+
+ // Resize the fuzzy map like the target, put 2.0 when outside
+ map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0); // Put 2.0 when out of initial map
+
+ // end
+ return map;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::RelativePositionAnalyzerFilter<ImageType>::
+ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance,
+ double & threshold, double & reverseThreshold)
+{
+ // Get the histogram of fuzzy values inside the target image
+ typedef itk::LabelStatisticsImageFilter<FloatImageType, ImageType> FloatStatFilterType;
+ typename FloatStatFilterType::Pointer f = FloatStatFilterType::New();
+ f->SetInput(map);
+ f->SetLabelInput(target);
+ f->UseHistogramsOn();
+ f->SetHistogramParameters(bins, 0.0, 1.1);
+ f->Update();
+ int count = f->GetCount(GetForegroundValue());
+ typename FloatStatFilterType::HistogramPointer h = f->GetHistogram(GetForegroundValue());
+
+ // Debug : dump histogram
+ static int i=0;
+ std::ofstream histogramFile(std::string("fuzzy_histo_"+toString(i)+".txt").c_str());
+ for(int j=0; j<bins; j++) {
+ histogramFile << h->GetMeasurement(j,0)
+ << "\t" << h->GetFrequency(j)
+ << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl;
+ }
+ histogramFile.close();
+ i++;
+
+ // Analyze the histogram (direct)
+ double sum = 0.0;
+ bool found = false;
+ threshold = 0.0;
+ for(int j=0; j<bins; j++) {
+ sum += ((double)h->GetFrequency(j)/(double)count);
+ if ((!found) && (sum > tolerance)) {
+ threshold = h->GetBinMin(0,j);
+ found = true;
+ }
+ }
+
+ // Analyze the histogram (reverse)
+ sum = 0.0;
+ found = false;
+ reverseThreshold = 1.0;
+ for(int j=bins-1; j>=0; j--) {
+ sum += ((double)h->GetFrequency(j)/(double)count);
+ if ((!found) && (sum > tolerance)) {
+ reverseThreshold = h->GetBinMax(0,j);
+ found = true;
+ }
+ }
+}
+//--------------------------------------------------------------------
+