From 1023d9d8de9f99a29043d54042ac45082c9a830f Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 24 Oct 2011 08:27:19 +0200 Subject: [PATCH] First version of RelativePositionAnalyzer --- itk/clitkRelativePositionAnalyzerFilter.h | 38 +-- itk/clitkRelativePositionAnalyzerFilter.txx | 242 +++++++++++++++++--- 2 files changed, 235 insertions(+), 45 deletions(-) diff --git a/itk/clitkRelativePositionAnalyzerFilter.h b/itk/clitkRelativePositionAnalyzerFilter.h index 2933958..7d70fe7 100644 --- a/itk/clitkRelativePositionAnalyzerFilter.h +++ b/itk/clitkRelativePositionAnalyzerFilter.h @@ -20,28 +20,35 @@ #define CLITKRELATIVEPOSITIONANALYZERFILTER_H // clitk -#include "clitkCommon.h" +#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkFilterBase.h" +#include "clitkSliceBySliceRelativePositionFilter.h" // itk #include +#include namespace clitk { //-------------------------------------------------------------------- /* - TODO + Analyze the relative position of a Target (mask image) according + to some Object, in a given Support. Indicate the main important + position of this Target according the Object. */ //-------------------------------------------------------------------- template - class ITK_EXPORT RelativePositionAnalyzerFilter: - public itk::ImageToImageFilter + class RelativePositionAnalyzerFilter: + public virtual FilterBase, + public clitk::FilterWithAnatomicalFeatureDatabaseManagement, + public itk::ImageToImageFilter { public: /** Standard class typedefs. */ typedef itk::ImageToImageFilter Superclass; - typedef RelativePositionAnalyzerFilter Self; + typedef RelativePositionAnalyzerFilter Self; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; @@ -63,6 +70,7 @@ namespace clitk { /** ImageDimension constants */ itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + FILTERBASE_INIT; typedef itk::Image FloatImageType; /** Input : initial image and object */ @@ -71,10 +79,6 @@ namespace clitk { void SetInputTarget(const ImageType * image); // Options - itkSetMacro(VerboseFlag, bool); - itkGetConstMacro(VerboseFlag, bool); - itkBooleanMacro(VerboseFlag); - itkGetConstMacro(BackgroundValue, PixelType); itkSetMacro(BackgroundValue, PixelType); @@ -84,19 +88,29 @@ namespace clitk { // For debug void PrintOptions(); - protected: + // I dont want to verify inputs information + virtual void VerifyInputInformation() { } + + protected: RelativePositionAnalyzerFilter(); virtual ~RelativePositionAnalyzerFilter() {} - bool m_VerboseFlag; PixelType m_BackgroundValue; PixelType m_ForegroundValue; ImagePointer m_Support; ImagePointer m_Object; ImagePointer m_Target; + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); virtual void GenerateData(); + typename FloatImageType::Pointer + ComputeFuzzyMap(ImageType * object, ImageType * target, double angle); + + void + ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, double tolerance, + double & threshold, double & reverseThreshold); private: RelativePositionAnalyzerFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented @@ -107,8 +121,6 @@ namespace clitk { } // end namespace clitk //-------------------------------------------------------------------- -#ifndef ITK_MANUAL_INSTANTIATION #include "clitkRelativePositionAnalyzerFilter.txx" -#endif #endif diff --git a/itk/clitkRelativePositionAnalyzerFilter.txx b/itk/clitkRelativePositionAnalyzerFilter.txx index a732aab..8b69c68 100644 --- a/itk/clitkRelativePositionAnalyzerFilter.txx +++ b/itk/clitkRelativePositionAnalyzerFilter.txx @@ -20,10 +20,14 @@ template clitk::RelativePositionAnalyzerFilter:: RelativePositionAnalyzerFilter(): -itk::ImageToImageFilter() + clitk::FilterBase(), + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + itk::ImageToImageFilter() { this->SetNumberOfRequiredInputs(3); // support, object, target VerboseFlagOff(); + SetBackgroundValue(0); + SetForegroundValue(1); } //-------------------------------------------------------------------- @@ -59,7 +63,7 @@ clitk::RelativePositionAnalyzerFilter:: SetInputTarget(const ImageType * image) { // Process object is not const-correct so the const casting is required. - this->SetNthInput(1, const_cast(image)); + this->SetNthInput(2, const_cast(image)); } //-------------------------------------------------------------------- @@ -75,56 +79,140 @@ PrintOptions() //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionAnalyzerFilter:: +GenerateOutputInformation() +{ + ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer outputImage = this->GetOutput(0); + outputImage->SetRegions(outputImage->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionAnalyzerFilter:: +GenerateInputRequestedRegion() +{ + // Call default + itk::ImageToImageFilter::GenerateInputRequestedRegion(); + // Get input pointers and set requested region to common region + ImagePointer input1 = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer input2 = dynamic_cast(itk::ProcessObject::GetInput(1)); + ImagePointer input3 = dynamic_cast(itk::ProcessObject::GetInput(2)); + input1->SetRequestedRegion(input1->GetLargestPossibleRegion()); + input2->SetRequestedRegion(input2->GetLargestPossibleRegion()); + input3->SetRequestedRegion(input3->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- + //-------------------------------------------------------------------- template void clitk::RelativePositionAnalyzerFilter:: GenerateData() { - DD("Update"); + this->LoadAFDB(); + // Get input pointer m_Support = dynamic_cast(itk::ProcessObject::GetInput(0)); m_Object = dynamic_cast(itk::ProcessObject::GetInput(1)); m_Target = dynamic_cast(itk::ProcessObject::GetInput(2)); static const unsigned int dim = ImageType::ImageDimension; - /* - Prerequisite: - - target fully inside support ? - */ + // Remove object from support + clitk::AndNot(m_Support, m_Object, GetBackgroundValue()); - //for(int i=0; i(m_Object, m_Target, GetBackgroundValue()); + + // Define filter to compute statics on mask image + typedef itk::LabelStatisticsImageFilter StatFilterType; + typename StatFilterType::Pointer statFilter = StatFilterType::New(); - // Get minimal value in the target area + // Compute the initial support size + statFilter->SetInput(m_Support); + statFilter->SetLabelInput(m_Support); + statFilter->Update(); + int m_SupportSize = statFilter->GetCount(GetForegroundValue()); + // DD(m_SupportSize); + + // Compute the initial target size + ImagePointer s = clitk::ResizeImageLike(m_Support, m_Target, GetBackgroundValue()); + statFilter->SetInput(s); + statFilter->SetLabelInput(m_Target); + statFilter->Update(); + int m_TargetSize = statFilter->GetCount(GetForegroundValue()); + // DD(m_TargetSize); - // Or -> analyze floating point values inside the target area (histo ?) + // Build the list of tested orientations + std::vector m_ListOfAngles; + int m_NumberOfAngles = this->GetAFDB()->GetDouble("NumberOfAngles"); + for(uint i=0; i180) a = 180-a; + m_ListOfAngles.push_back(clitk::deg2rad(a)); + } - // Set threshold and compute new support + // Loop on all orientations + for(int i=0; i StatisticsImageFilterType; - typename StatisticsImageFilterType::Pointer statisticsFilter = StatisticsImageFilterType::New(); - statisticsFilter->SetInput(m_Input); - statisticsFilter->SetLabelInput(m_Input); - statisticsFilter->Update(); - int n = labelStatisticsImageFilter->GetCount(GetForegroundValue()); - DD(n); - statisticsFilter = StatisticsImageFilterType::New(); - statisticsFilter->SetInput(m_Output); - statisticsFilter->SetLabelInput(m_Output); - statisticsFilter->Update(); - int m = labelStatisticsImageFilter->GetCount(GetForegroundValue()); - DD(m); - */ + // Compute Fuzzy map + typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]); + writeImage(map, "fuzzy_"+toString(i)+".mha"); + + // Compute the optimal thresholds (direct and inverse) + double mThreshold; + double mReverseThreshold; + int bins = this->GetAFDB()->GetDouble("bins"); + double tolerance = this->GetAFDB()->GetDouble("TargetAreaLossTolerance"); + ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold); + + // Use the threshold to compute new support + int s1; + if (mThreshold > 0.0) { + ImagePointer support1 = + clitk::SliceBySliceRelativePosition(m_Support, m_Object, 2, + mThreshold, + m_ListOfAngles[i],false, + false, -1, true, false); + writeImage(support1, "sup_"+toString(i)+".mha"); + // Compute the new support size + statFilter->SetInput(support1); + statFilter->SetLabelInput(support1); + statFilter->Update(); + s1 = statFilter->GetCount(GetForegroundValue()); + } + else s1 = m_SupportSize; + + int s2; + if (mReverseThreshold < 1.0) { + ImagePointer support2 = + clitk::SliceBySliceRelativePosition(m_Support, m_Object, 2, + mReverseThreshold, + m_ListOfAngles[i],true, + false, -1, true, false); + writeImage(support2, "sup_rev_"+toString(i)+".mha"); + // Compute the new support size + statFilter = StatFilterType::New(); + statFilter->SetInput(support2); + statFilter->SetLabelInput(support2); + statFilter->Update(); + s2 = statFilter->GetCount(GetForegroundValue()); + } + else s2 =m_SupportSize; // Print results + std::cout << i << " " << clitk::rad2deg(m_ListOfAngles[i]) << "\t" + << m_SupportSize << " " << m_TargetSize << "\t" + << s1/(double)m_SupportSize << " " << s2/(double)m_SupportSize << "\t" + << mThreshold << " " << mReverseThreshold << std::endl; - //} + } // end loop on orientations // Final Step -> set output TODO @@ -133,3 +221,93 @@ GenerateData() } //-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +typename clitk::RelativePositionAnalyzerFilter::FloatImageType::Pointer +clitk::RelativePositionAnalyzerFilter:: +ComputeFuzzyMap(ImageType * object, ImageType * target, double angle) +{ + typedef clitk::SliceBySliceRelativePositionFilter 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(map, object, GetForegroundValue(), 0.0, true); + + // Resize the fuzzy map like the target, put 2.0 when outside + map = clitk::ResizeImageLike(map, target, 2.0); // Put 2.0 when out of initial map + + // end + return map; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionAnalyzerFilter:: +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 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; jGetMeasurement(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; jGetFrequency(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; + } + } +} +//-------------------------------------------------------------------- + -- 2.47.1