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);
- SetNumberOfAngles(4);
SetAreaLossTolerance(0.01);
- m_ListOfAngles.clear();
SetSupportSize(0);
SetTargetSize(0);
SetSizeWithThreshold(0);
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();
SetTargetSize(statFilter->GetCount(GetForegroundValue()));
// DD(GetTargetSize());
- // Build the list of tested orientations
- m_ListOfAngles.clear();
- for(uint i=0; i<GetNumberOfAngles(); i++) {
- double a = i*360.0/GetNumberOfAngles();
- if (a>180) 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<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=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<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());
- }
-
- int s2 = GetSupportSize();
- if (mReverseThreshold < 1.0) {
- // DD(m_ListOfAngles[1]);
- 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());
- }
- // 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<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());
+ }
+
+ // 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)) {
- 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;
}
}
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 <class ImageType>
-void
-clitk::RelativePositionAnalyzerFilter<ImageType>::
-Print(std::string s, std::ostream & os)
-{
- for(int i=0; i<m_ListOfOrientation.size(); i++) {
- os << s << " ";
- m_ListOfOrientation[i].Print(os);
- m_ListOfInformation[i].Println(os);
- }
}
//--------------------------------------------------------------------
+