X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkRelativePositionAnalyzerFilter.txx;h=df12ea33437c42975282662ad74015745d625041;hb=39562adceed3b608834959f2ed079b945187bb22;hp=a732aab61e363fda81e3c8244f4509b9507a43ba;hpb=aeb947ddd800ab06cf4916c9371bca3832056b4f;p=clitk.git diff --git a/itk/clitkRelativePositionAnalyzerFilter.txx b/itk/clitkRelativePositionAnalyzerFilter.txx index a732aab..df12ea3 100644 --- a/itk/clitkRelativePositionAnalyzerFilter.txx +++ b/itk/clitkRelativePositionAnalyzerFilter.txx @@ -20,10 +20,17 @@ template clitk::RelativePositionAnalyzerFilter:: RelativePositionAnalyzerFilter(): -itk::ImageToImageFilter() + itk::ImageToImageFilter() { - this->SetNumberOfRequiredInputs(3); // support, object, target - VerboseFlagOff(); + this->SetNumberOfRequiredInputs(3); // Input : support, object, target + SetBackgroundValue(0); + SetForegroundValue(1); + SetNumberOfBins(100); + SetAreaLossTolerance(0.01); + SetSupportSize(0); + SetTargetSize(0); + SetSizeWithThreshold(0); + SetSizeWithReverseThreshold(0); } //-------------------------------------------------------------------- @@ -59,7 +66,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,61 +82,204 @@ 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:: GenerateData() { - DD("Update"); - // Get input pointer - m_Support = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer temp = 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 (keep initial image) + m_Support = clitk::Clone(temp); + clitk::AndNot(m_Support, m_Object, GetBackgroundValue()); - //for(int i=0; i analyze floating point values inside the target area (histo ?) - - // Set threshold and compute new support - - // Compute ratio of area before/after - // http://www.itk.org/Doxygen/html/classitk_1_1LabelStatisticsImageFilter.html#details - /* - typedef itk::LabelStatisticsImageFilter 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); - */ - - // Print results - - //} + // Define filter to compute statics on mask image + typedef itk::LabelStatisticsImageFilter StatFilterType; + typename StatFilterType::Pointer statFilter = StatFilterType::New(); + + // Compute the initial support size + statFilter->SetInput(m_Support); + statFilter->SetLabelInput(m_Support); + statFilter->Update(); + SetSupportSize(statFilter->GetCount(GetForegroundValue())); + // DD(GetSupportSize()); + + // Compute the initial target size + ImagePointer s = clitk::ResizeImageLike(m_Support, m_Target, GetBackgroundValue()); + statFilter->SetInput(s); + statFilter->SetLabelInput(m_Target); + statFilter->Update(); + SetTargetSize(statFilter->GetCount(GetForegroundValue())); + // DD(GetTargetSize()); + + // + int bins = GetNumberOfBins(); + double tolerance = GetAreaLossTolerance(); + + // Compute Fuzzy map + double angle = GetDirection().angle1; + typename FloatImageType::Pointer map = ComputeFuzzyMap(m_Object, m_Target, m_Support, angle); + writeImage(map, "fuzzy_"+toString(clitk::rad2deg(angle))+".mha"); + + // Compute the optimal thresholds (direct and inverse) + double mThreshold=0.0; + double mReverseThreshold=1.0; + ComputeOptimalThresholds(map, m_Target, bins, tolerance, mThreshold, mReverseThreshold); + + // Use the threshold to compute new support + int s1 = GetSupportSize(); + if (mThreshold > 0.0) { + ImagePointer support1 = + clitk::SliceBySliceRelativePosition(m_Support, m_Object, 2, + mThreshold, + angle,false, // inverseFlag + false, // uniqueConnectedComponent + -1, true, + false);//singleObjectCCL + // Compute the new support size + statFilter->SetInput(support1); + statFilter->SetLabelInput(support1); + statFilter->Update(); + s1 = statFilter->GetCount(GetForegroundValue()); + } + + int s2 = GetSupportSize(); + if (mReverseThreshold < 1.0) { + ImagePointer support2 = + clitk::SliceBySliceRelativePosition(m_Support, m_Object, 2, + mReverseThreshold, + angle,true,// inverseFlag + false, // uniqueConnectedComponent + -1, true, + false); //singleObjectCCL + // Compute the new support size + statFilter = StatFilterType::New(); + statFilter->SetInput(support2); + statFilter->SetLabelInput(support2); + statFilter->Update(); + s2 = statFilter->GetCount(GetForegroundValue()); + } + + // Set results values + m_Info.threshold = mThreshold; + m_Info.sizeAfterThreshold = s1; + m_Info.sizeBeforeThreshold = GetSupportSize(); + m_Info.sizeReference = GetTargetSize(); + m_InfoReverse.threshold = mReverseThreshold; + m_InfoReverse.sizeAfterThreshold = s2; + m_InfoReverse.sizeBeforeThreshold = GetSupportSize(); + m_InfoReverse.sizeReference = GetTargetSize(); +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +typename clitk::RelativePositionAnalyzerFilter::FloatImageType::Pointer +clitk::RelativePositionAnalyzerFilter:: +ComputeFuzzyMap(ImageType * object, ImageType * target, ImageType * support, double angle) +{ + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typedef typename SliceRelPosFilterType::FloatImageType FloatImageType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(support); + sliceRelPosFilter->SetInputObject(object); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetIntermediateSpacingFlag(false); + //sliceRelPosFilter->AddOrientationTypeString(orientation); + sliceRelPosFilter->AddAnglesInRad(angle, 0.0); + sliceRelPosFilter->FuzzyMapOnlyFlagOn(); // do not threshold, only compute the fuzzy map + // sliceRelPosFilter->PrintOptions(); + sliceRelPosFilter->Update(); + typename FloatImageType::Pointer map = sliceRelPosFilter->GetFuzzyMap(); + + // Resize map like object to allow SetBackground + map = clitk::ResizeImageLike(map, object, GetBackgroundValue()); + + // 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 - // Final Step -> set output TODO - // this->SetNthOutput(0, working_image); - // this->GraftOutput(working_image); + // 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)) { + if (j==0) threshold = h->GetBinMin(0,j); + else threshold = h->GetBinMin(0,j-1); // the last before reaching the threshold + 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)) { + if (j==bins-1) reverseThreshold = h->GetBinMax(0,j); + else reverseThreshold = h->GetBinMax(0,j+1); + found = true; + } + } } //--------------------------------------------------------------------