X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkRelativePositionAnalyzerFilter.txx;h=63af1191987386dfac0b632f0af23559ef967e10;hb=595624eb4297e747630105d45017de69733caef2;hp=a732aab61e363fda81e3c8244f4509b9507a43ba;hpb=938c4491f8ef70565b124a7eeb4f9e69b1b9baba;p=clitk.git diff --git a/itk/clitkRelativePositionAnalyzerFilter.txx b/itk/clitkRelativePositionAnalyzerFilter.txx index a732aab..63af119 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,249 @@ 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; + ImagePointer temp2 = dynamic_cast(itk::ProcessObject::GetInput(2)); - /* - 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 - - //} + // Remove object from target. Important because sometimes, there is + // overlap between target and object. + m_Target = clitk::Clone(temp2); + clitk::AndNot(m_Target, m_Object, GetBackgroundValue()); + + // 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); + + // DD(mThreshold); + // DD(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()); + } + + // Check threshold, if we gain nothing, we force to max/min thresholds + // DD(GetSupportSize()); + // DD(s1); + // DD(s2); + if (s1 >= GetSupportSize()) mThreshold = 0.0; + if (s2 >= GetSupportSize()) mReverseThreshold = 1.0; + + // 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(); + writeImage(map, "fuzzy_0_"+toString(clitk::rad2deg(angle))+".mha"); + + // Resize object like map to allow SetBackground + ImagePointer temp = clitk::ResizeImageLike(object, map, GetBackgroundValue()); + // writeImage(map, "fuzzy_1_"+toString(clitk::rad2deg(angle))+".mha"); + + // Remove initial object from the fuzzy map + map = clitk::SetBackground(map, temp, GetForegroundValue(), 0.0, true); + writeImage(map, "fuzzy_2_"+toString(clitk::rad2deg(angle))+".mha"); + + // 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 + writeImage(map, "fuzzy_3_"+toString(clitk::rad2deg(angle))+".mha"); - // 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.0/bins), 1.0+(1.0/bins)); + f->Update(); + int count = f->GetCount(GetForegroundValue()); + // DD(count); + 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(); + std::ofstream histogramFile2(std::string("fuzzy_histo_R_"+toString(i)+".txt").c_str()); + for(int j=bins-1; j>=0; j--) { + histogramFile2 << h->GetMeasurement(j,0) + << "\t" << h->GetFrequency(j) + << "\t" << (double)h->GetFrequency(j)/(double)count << std::endl; + } + histogramFile2.close(); + i++; + + // Analyze the histogram (direct) + double sum = 0.0; + bool found = false; + threshold = 0.0; + for(int j=0; jGetFrequency(j)/(double)count); + // DD(j); + // DD(sum); + // DD(threshold); + // DD(h->GetBinMin(0,j)); + // DD(h->GetBinMax(0,j)); + if ((!found) && (sum > tolerance)) { + // We consider as threshold the laste before current, because + if (j==0) + threshold = h->GetBinMin(0,j); + else threshold = h->GetBinMin(0,j-1); // FIXME ? the last before reaching the threshold + // DD(threshold); + found = true; + j = bins; + } + } + + // 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); + // DD(j); + // DD(sum); + // DD(reverseThreshold); + // DD(h->GetBinMin(0,j)); + // DD(h->GetBinMax(0,j)); + if ((!found) && (sum > tolerance)) { + if (j==bins-1) + reverseThreshold = h->GetBinMax(0,j); + else reverseThreshold = h->GetBinMax(0,j-1);// FIXME ? the last before reaching the threshold + // DD(reverseThreshold); + found = true; + j = -1; + } + } + } //--------------------------------------------------------------------