X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=sidebyside;f=itk%2FclitkRelativePositionAnalyzerFilter.txx;h=63af1191987386dfac0b632f0af23559ef967e10;hb=afede09ad631dcf7297e3189aeb1d2288fb25902;hp=c8d126a96d92d1ccf444dd59e08c83282ec77d75;hpb=5bff4d89c8ca42336267c36974e9ceb6fe2ac843;p=clitk.git diff --git a/itk/clitkRelativePositionAnalyzerFilter.txx b/itk/clitkRelativePositionAnalyzerFilter.txx index c8d126a..63af119 100644 --- a/itk/clitkRelativePositionAnalyzerFilter.txx +++ b/itk/clitkRelativePositionAnalyzerFilter.txx @@ -20,18 +20,13 @@ template clitk::RelativePositionAnalyzerFilter:: RelativePositionAnalyzerFilter(): - // clitk::FilterBase(), - clitk::FilterWithAnatomicalFeatureDatabaseManagement(), itk::ImageToImageFilter() { - this->SetNumberOfRequiredInputs(3); // support, object, target - VerboseFlagOff(); + this->SetNumberOfRequiredInputs(3); // Input : support, object, target SetBackgroundValue(0); SetForegroundValue(1); SetNumberOfBins(100); - SetNumberOfAngles(4); SetAreaLossTolerance(0.01); - m_ListOfAngles.clear(); SetSupportSize(0); SetTargetSize(0); SetSizeWithThreshold(0); @@ -106,20 +101,19 @@ void clitk::RelativePositionAnalyzerFilter:: GenerateData() { - this->LoadAFDB(); - - // 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)); - // Remove object from support + // Remove object from support (keep initial image) + m_Support = clitk::Clone(temp); clitk::AndNot(m_Support, m_Object, GetBackgroundValue()); - // Resize object like target (to enable substraction later) - ImagePointer objectLikeTarget = clitk::ResizeImageLike(m_Object, m_Target, GetBackgroundValue()); - + // 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(); @@ -139,82 +133,73 @@ GenerateData() SetTargetSize(statFilter->GetCount(GetForegroundValue())); // DD(GetTargetSize()); - // Build the list of tested orientations - m_ListOfAngles.clear(); - for(uint i=0; i180) a = 180-a; - m_ListOfAngles.push_back(clitk::deg2rad(a)); - RelativePositionOrientationType r; - r.angle1 = clitk::deg2rad(a); - r.angle2 = 0; - r.notFlag = false; - m_ListOfOrientation.push_back(r); - r.notFlag = true; - m_ListOfOrientation.push_back(r); - } - - // Loop on all orientations + // int bins = GetNumberOfBins(); double tolerance = GetAreaLossTolerance(); - for(int i=0; i(map, "fuzzy_"+toString(i)+".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(); - // DD(mThreshold); - // DD(mReverseThreshold); - 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()); - } - - int s2 = GetSupportSize(); - if (mReverseThreshold < 1.0) { - // DD(m_ListOfAngles[1]); - 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()); - } - // Set results values - RelativePositionInformationType r; - r.threshold = mThreshold; - r.sizeAfterThreshold = s1; // DD(s1); - r.sizeBeforeThreshold = GetSupportSize(); - r.sizeReference = GetTargetSize(); - m_ListOfInformation.push_back(r); - - r.threshold = mReverseThreshold; - r.sizeAfterThreshold = s2; // DD(s2); - m_ListOfInformation.push_back(r); - // Print(); - } // end loop on orientations + // 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(); } //-------------------------------------------------------------------- @@ -223,29 +208,36 @@ GenerateData() template typename clitk::RelativePositionAnalyzerFilter::FloatImageType::Pointer clitk::RelativePositionAnalyzerFilter:: -ComputeFuzzyMap(ImageType * object, ImageType * target, double angle) +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(target); + sliceRelPosFilter->SetInput(support); sliceRelPosFilter->SetInputObject(object); sliceRelPosFilter->SetDirection(2); sliceRelPosFilter->SetIntermediateSpacingFlag(false); //sliceRelPosFilter->AddOrientationTypeString(orientation); - sliceRelPosFilter->AddAngles(angle, 0.0); + 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, object, GetForegroundValue(), 0.0, true); + 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"); // end return map; @@ -266,9 +258,10 @@ ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, dou f->SetInput(map); f->SetLabelInput(target); f->UseHistogramsOn(); - f->SetHistogramParameters(bins, 0.0, 1.1); + 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 @@ -280,18 +273,34 @@ ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, dou << "\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)) { - if (j==0) threshold = h->GetBinMin(0,j); - else threshold = h->GetBinMin(0,j-1); // the last before reaching the threshold + // 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; } } @@ -299,28 +308,23 @@ ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, dou sum = 0.0; found = false; reverseThreshold = 1.0; - for(int j=bins-1; j>=0; j--) { + 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); + 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; } } -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -clitk::RelativePositionAnalyzerFilter:: -Print(std::string s, std::ostream & os) -{ - for(int i=0; i