template <class ImageType>
clitk::RelativePositionAnalyzerFilter<ImageType>::
RelativePositionAnalyzerFilter():
- clitk::FilterBase(),
- clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
itk::ImageToImageFilter<ImageType, ImageType>()
{
- 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);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-//--------------------------------------------------------------------
-template <class ImageType>
-void
-clitk::RelativePositionAnalyzerFilter<ImageType>::
-GenerateInputRequestedRegion()
-{
- // Call default
- itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
- // Get input pointers and set requested region to common region
- ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
- ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
- ImagePointer input3 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
- input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
- input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
- input3->SetRequestedRegion(input3->GetLargestPossibleRegion());
-}
-//--------------------------------------------------------------------
-
//--------------------------------------------------------------------
template <class ImageType>
void
clitk::RelativePositionAnalyzerFilter<ImageType>::
GenerateData()
{
- this->LoadAFDB();
-
- // Get input pointer
- m_Support = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
+ ImagePointer temp = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
m_Object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
- m_Target = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
- static const unsigned int dim = ImageType::ImageDimension;
+ ImagePointer temp2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
- // Remove object from support
+ // Remove object from support (keep initial image)
+ m_Support = clitk::Clone<ImageType>(temp);
clitk::AndNot<ImageType>(m_Support, m_Object, GetBackgroundValue());
- // Resize object like target (to enable substraction later)
- ImagePointer objectLikeTarget = clitk::ResizeImageLike<ImageType>(m_Object, m_Target, GetBackgroundValue());
-
+ // Remove object from target. Important because sometimes, there is
+ // overlap between target and object.
+ m_Target = clitk::Clone<ImageType>(temp2);
+ clitk::AndNot<ImageType>(m_Target, m_Object, GetBackgroundValue());
+
// Define filter to compute statics on mask image
typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
typename StatFilterType::Pointer statFilter = StatFilterType::New();
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<ImageType>(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<double> m_ListOfAngles;
- int m_NumberOfAngles = this->GetAFDB()->GetDouble("NumberOfAngles");
- for(uint i=0; i<m_NumberOfAngles; i++) {
- double a = i*360.0/m_NumberOfAngles;
- if (a>180) 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<FloatImageType>(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<ImageType>(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<ImageType>(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());
}
-
- // Loop on all orientations
- for(int i=0; i<m_ListOfAngles.size(); i++) {
-
- // Compute Fuzzy map
- typename FloatImageType::Pointer map = ComputeFuzzyMap(objectLikeTarget, m_Target, m_ListOfAngles[i]);
- writeImage<FloatImageType>(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<ImageType>(m_Support, m_Object, 2,
- mThreshold,
- m_ListOfAngles[i],false,
- false, -1, true, false);
- writeImage<ImageType>(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<ImageType>(m_Support, m_Object, 2,
- mReverseThreshold,
- m_ListOfAngles[i],true,
- false, -1, true, false);
- writeImage<ImageType>(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);
+ // 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 <class ImageType>
typename clitk::RelativePositionAnalyzerFilter<ImageType>::FloatImageType::Pointer
clitk::RelativePositionAnalyzerFilter<ImageType>::
-ComputeFuzzyMap(ImageType * object, ImageType * target, double angle)
+ComputeFuzzyMap(ImageType * object, ImageType * target, ImageType * support, double angle)
{
typedef clitk::SliceBySliceRelativePositionFilter<ImageType> 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<FloatImageType>(map, "fuzzy_0_"+toString(clitk::rad2deg(angle))+".mha");
+ // Resize object like map to allow SetBackground
+ ImagePointer temp = clitk::ResizeImageLike<ImageType>(object, map, GetBackgroundValue());
+ // writeImage<FloatImageType>(map, "fuzzy_1_"+toString(clitk::rad2deg(angle))+".mha");
+
// Remove initial object from the fuzzy map
- map = clitk::SetBackground<FloatImageType, ImageType>(map, object, GetForegroundValue(), 0.0, true);
+ map = clitk::SetBackground<FloatImageType, ImageType>(map, temp, GetForegroundValue(), 0.0, true);
+ writeImage<FloatImageType>(map, "fuzzy_2_"+toString(clitk::rad2deg(angle))+".mha");
// Resize the fuzzy map like the target, put 2.0 when outside
map = clitk::ResizeImageLike<FloatImageType>(map, target, 2.0); // Put 2.0 when out of initial map
+ writeImage<FloatImageType>(map, "fuzzy_3_"+toString(clitk::rad2deg(angle))+".mha");
// end
return map;
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
<< "\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; j<bins; j++) {
+ for(int j=0; j<bins-1; j++) {
sum += ((double)h->GetFrequency(j)/(double)count);
+ // DD(j);
+ // DD(sum);
+ // DD(threshold);
+ // DD(h->GetBinMin(0,j));
+ // DD(h->GetBinMax(0,j));
if ((!found) && (sum > tolerance)) {
- threshold = h->GetBinMin(0,j);
+ // 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;
}
}
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);// FIXME ? the last before reaching the threshold
+ // DD(reverseThreshold);
found = true;
+ j = -1;
}
}
+
}
//--------------------------------------------------------------------