X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkRelativePositionAnalyzerFilter.txx;h=df12ea33437c42975282662ad74015745d625041;hb=3c2eea5109c2d80123fbacc41670bc77b8043aaf;hp=8b69c6875f981f38fde1adac82d3d35e9034f285;hpb=1023d9d8de9f99a29043d54042ac45082c9a830f;p=clitk.git diff --git a/itk/clitkRelativePositionAnalyzerFilter.txx b/itk/clitkRelativePositionAnalyzerFilter.txx index 8b69c68..df12ea3 100644 --- a/itk/clitkRelativePositionAnalyzerFilter.txx +++ b/itk/clitkRelativePositionAnalyzerFilter.txx @@ -20,14 +20,17 @@ 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); + SetAreaLossTolerance(0.01); + SetSupportSize(0); + SetTargetSize(0); + SetSizeWithThreshold(0); + SetSizeWithReverseThreshold(0); } //-------------------------------------------------------------------- @@ -92,44 +95,20 @@ GenerateOutputInformation() //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -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() { - 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; - // 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()); - // Define filter to compute statics on mask image typedef itk::LabelStatisticsImageFilter StatFilterType; typename StatFilterType::Pointer statFilter = StatFilterType::New(); @@ -138,86 +117,74 @@ GenerateData() statFilter->SetInput(m_Support); statFilter->SetLabelInput(m_Support); statFilter->Update(); - int m_SupportSize = statFilter->GetCount(GetForegroundValue()); - // DD(m_SupportSize); + 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(); - int m_TargetSize = statFilter->GetCount(GetForegroundValue()); - // DD(m_TargetSize); - - // 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)); + 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()); } - - // Loop on all orientations - for(int i=0; i(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 - // this->SetNthOutput(0, working_image); - // this->GraftOutput(working_image); + 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(); } //-------------------------------------------------------------------- @@ -226,24 +193,27 @@ 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(); + // 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); @@ -292,7 +262,8 @@ ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, dou for(int j=0; jGetFrequency(j)/(double)count); if ((!found) && (sum > tolerance)) { - threshold = h->GetBinMin(0,j); + if (j==0) threshold = h->GetBinMin(0,j); + else threshold = h->GetBinMin(0,j-1); // the last before reaching the threshold found = true; } } @@ -304,7 +275,8 @@ ComputeOptimalThresholds(FloatImageType * map, ImageType * target, int bins, dou for(int j=bins-1; j>=0; j--) { sum += ((double)h->GetFrequency(j)/(double)count); if ((!found) && (sum > tolerance)) { - reverseThreshold = h->GetBinMax(0,j); + if (j==bins-1) reverseThreshold = h->GetBinMax(0,j); + else reverseThreshold = h->GetBinMax(0,j+1); found = true; } }