@matt.mccormick I wrote a function that takes two inputs (float 3D ImageType is the intensity input image and a unsigned char 3D image (the prediction output of my algorithm). Then, I calculated the means of predicted labels using itk::LabelStatisticsImageFilter. I set the parameters numberOfIterations=1000, smothingFactor=2.5, and SetErrorTolerance(1e-7), and the weights for 3D images I set according to this link. However, the output image is a blank image, which I do not know what is the reason causing this issue.
CharImageType::Pointer ApplyMRF(ImageType::Pointer image,CharImageType::Pointer prediction)
{
/*
* This function applies MRF for the correction of
* Since the Markov Random Field algorithm is defined in general for images whose pixels have multiple
components, that is, images of vector type, we must adapt our scalar image in order to satisfy the
interface expected by the MRFImageFilter. We do this by using the itk::ComposeImageFilter.
With this filter we will present our scalar image as a vector image whose vector pixels contain a
single component.
*/
LabelStatisticsImageFilterType::Pointer labelStatisticsImageFilter=LabelStatisticsImageFilterType::New();
labelStatisticsImageFilter->SetLabelInput(prediction);
labelStatisticsImageFilter->SetInput(image);
labelStatisticsImageFilter->Update();
std::cout << "Number of labels: " << labelStatisticsImageFilter->GetNumberOfLabels() << std::endl;
//
// getting the label Statistics
// Link: https://itk.org/Doxygen/html/classitk_1_1LabelStatisticsImageFilter.html
//
//using LabelPixelType=LabelStatisticsImageFilterType::LabelPixelType ; //-
//LabelStatisticsImageFilterType::LabelPixelType LabelPixelType; //-
float labelMeanValue [labelStatisticsImageFilter->GetNumberOfLabels()]; //an array for saving the means of labels
for (auto vIt = labelStatisticsImageFilter->GetValidLabelValues().begin();
vIt != labelStatisticsImageFilter->GetValidLabelValues().end();
++vIt) {
if (labelStatisticsImageFilter->HasLabel(*vIt)) {
//LabelPixelType labelValue = *vIt; //-
labelMeanValue[(int)*vIt]=labelStatisticsImageFilter->GetMean(*vIt); //+
std::cout << "Label: " << (int)*vIt << std::endl;
//std::cout << "\tmean: " <<labelStatisticsImageFilter->GetMean(labelValue)<< std::endl; //-
std::cout << "\tmean: " <<labelMeanValue[(int)*vIt]<< std::endl; //+
}
}
/* printing the sorted means of labels
* int len=sizeof(labelMeanValue)/sizeof(labelMeanValue[0]);
for(int i=0;i<len;++i){ //0 1 2 3 4
std::cout<<labelMeanValue[i]<<"\t"; //0.278904 0.115523 0.496039 0.121493 0.43653
}*/
ScalarToArrayFilterType::Pointer scalerToArrayFilter=ScalarToArrayFilterType::New();
scalerToArrayFilter->SetInput(image);
scalerToArrayFilter->Update();
/*
* With the input image type ImageType and labeled image type CharImageType we instantiate the
type of the itk::MRFImageFilter that will apply the Markov Random Field algorithm in order to
refine the pixel classification
*/
MRFFilterType::Pointer mrfFilter=MRFFilterType::New();
mrfFilter->SetInput(scalerToArrayFilter->GetOutput());
/*
* We set now some of the parameters for the MRF filter. In particular, the number of classes to be
used during the classification, the maximum number of iterations to be run in this filter and the error
tolerance that will be used as a criterion for convergence.
*/
int numberOfIterations=1000;
mrfFilter->SetNumberOfClasses(numberOfClasses);
mrfFilter->SetMaximumNumberOfIterations(numberOfIterations);
mrfFilter->SetErrorTolerance(1e-7);
// The smoothing factor represents the tradeoff between fidelity to the observed image and the smoothness
// of the segmented image. Typical smoothing factors have values between 1-5
float smoothingFactor=2.5;
mrfFilter->SetSmoothingFactor(smoothingFactor);
SupervisedClassifierType::Pointer classifier=SupervisedClassifierType::New();
/*
* The classifier need a decision rule to be set by the user. Note that we must use GetPointer() in the
call of the SetDecisionRule() method because we are passing a SmartPointer, and smart pointer
cannot perform polymorphism, we must then extract the raw pointer that is associated to the smart
pointer. This extraction is done with the GetPointer() method.
*/
typedef itk::Statistics::MinimumDecisionRule DecisionRuleType;
DecisionRuleType::Pointer classificationDecisionRule=DecisionRuleType::New();
classifier->SetDecisionRule(classificationDecisionRule);
//classifier->SetClassifiedImage(prediction); //+
// Instantiating membership function
typedef itk::Statistics::DistanceToCentroidMembershipFunction<ArrayPixelType> MembershipFunctionType;
typedef MembershipFunctionType::Pointer MembershipFunctionPointer;
double meanDistance=0;
MembershipFunctionType::CentroidType centroid(1);
for (unsigned int i=0;i<numberOfClasses;i++) {
MembershipFunctionPointer membershipFunction = MembershipFunctionType::New();
centroid[0] = labelMeanValue[i];
membershipFunction->SetCentroid(centroid);
classifier->AddMembershipFunction(membershipFunction);
meanDistance += static_cast<double>(centroid[0]);
}
if (numberOfClasses>0)
{
meanDistance/=numberOfClasses;
}
else
{
std::cerr <<"Error: numberOfClasses is 0"<<std::endl;
}
//
// we set the neighborhoodradius that will define the size of the clique
//
mrfFilter->SetNeighborhoodRadius(1);
//Start with below slice, setting it as the example mentions
/*
* As far as I understand it, it goes along the rows from the top or bottom slice upwards/downwards…
* The weights have to be the same size as the neighbourhood. For a 3d image with a radius of 1 (3x3x3) It needs 9x3 weights
* Link: https://discourse.itk.org/t/solved-copying-headers-after-performing-operation/1111/3
*/
std::vector<double>weights; //The array is arranged and then passed to the filter by using the method SetMRFNeighborhoodWeight().
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.5); //center
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.7);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.5);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
weights.push_back(1.3);
/*
* We now scale weights so that the smoothing function and the image fidelity functions have comparable value.
* This is necessary since the label image and the input image can have different dynamic ranges.
* The fidelity function is usually computed using a distance function, such as the
itk::DistanceToCentroidMembershipFunction or one of the other membership functions. They
tend to have values in the order of the means specified
*/
double totalWeight=0;
for(std::vector<double>::const_iterator wcIt=weights.begin();
wcIt!=weights.end(); ++wcIt)
{
totalWeight+=*wcIt;
}
for(std::vector<double>::iterator wIt=weights.begin();
wIt!=weights.end();++wIt)
{
*wIt=static_cast<double>((*wIt)*meanDistance/(2*totalWeight));
}
mrfFilter->SetMRFNeighborhoodWeight(weights);
mrfFilter->SetClassifier(classifier);
mrfFilter->Update();
//
// The output image produced by the \doxygen{MRFImageFilter} has the same pixel
// type as the labeled input image. In the following lines we use the
// \code{OutputImageType} in order to instantiate the type of a
// \doxygen{ImageFileWriter}. Then create one, and connect it to the output of
// the classification filter after passing it through an intensity rescaler
// to rescale it to an 8 bit dynamic range
//
using OutputImageType=MRFFilterType::OutputImageType ;
using RescaledOutputImageType=itk::Image<unsigned char, Dimension>; //Similar to the CharImageType
typedef itk::RescaleIntensityImageFilter<RescaledOutputImageType ,OutputImageType > RescalerType;
RescalerType::Pointer intensityRescaler=RescalerType::New();
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
intensityRescaler->SetInput(mrfFilter->GetOutput());
// Software Guide : BeginCodeSnippet
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetInput(intensityRescaler->GetOutput());
writer->SetFileName("MRFOutput.nii.gz");
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We are now ready for triggering the execution of the pipeline. This is done
// by simply invoking the \code{Update()} method in the writer. This call will
// propagate the update request to the reader and then to the MRF filter.
//
// Software Guide : EndLatex
writer->Update();
std::cout << "Number of Iterations : ";
std::cout << mrfFilter->GetNumberOfIterations() << std::endl;
std::cout << "Stop condition: " << std::endl;
std::cout << " (1) Maximum number of iterations " << std::endl;
std::cout << " (2) Error tolerance: " << std::endl;
std::cout << mrfFilter->GetStopCondition() << std::endl;
return intensityRescaler->GetOutput();
}// end of ApplyMRF function