Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update GeodesicSegmentation CLI #1496

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
+23 −10
Open

Update GeodesicSegmentation CLI #1496

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7602df6
Update GeodesicSegmentation.cxx
Getkablue May 28, 2022
470ff95
remove needless comment
Getkablue May 28, 2022
53e4e14
Merge branch 'master' into fix-geodesic
AlexanderGetka-cbica May 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 23 additions & 10 deletions src/applications/GeodesicSegmentation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,15 +309,28 @@ int main(int argc, char **argv)
unsigned short threshold = 0;
if (parser.isPresent("t"))
{
std::string temp;
int temp;
parser.getParameterValue("t", temp);
threshold = std::atoi(temp.c_str());
threshold = temp;
}

std::cout << "Reading inputs.\n";
using ImageType = itk::Image< float, 3 >;
using ImageTypeGeodesic = itk::Image < short, 3 >;
auto inputImage = cbica::ReadImage< ImageType >(inputFile);
auto mask = cbica::ReadImage< ImageType >(drawingFile);
using CastFilterType = itk::CastImageFilter<ImageType, ImageTypeGeodesic>;
auto castFilter = CastFilterType::New();
castFilter->SetInput(inputImage);
castFilter->Update();
auto convertedInputImage = castFilter->GetOutput();
using RescaleFilterType = itk::RescaleIntensityImageFilter<ImageTypeGeodesic, ImageTypeGeodesic>;
auto rescaleFilter = RescaleFilterType::New();
rescaleFilter->SetInput(convertedInputImage);
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);
rescaleFilter->Update();
convertedInputImage = rescaleFilter->GetOutput();

VectorVectorDouble Indices;
typedef itk::ImageRegionIteratorWithIndex <ImageType> IteratorType;
Expand All @@ -338,15 +351,15 @@ int main(int argc, char **argv)
}

std::cout << "Performing Geodesic map estimation.\n";
GeodesicSegmentation< ImageType > segmentationClass;
auto geosOutput = segmentationClass.Run/*< ImageType >*/(inputImage, Indices);
GeodesicSegmentation< ImageTypeGeodesic > segmentationClass;
auto geosOutput = segmentationClass.Run/*< ImageType >*/(convertedInputImage, Indices);
ImageType::PixelType maxVal;

if (normalize)
{
std::cout << "Performing normalization.\n";
maxVal = 255;
auto rescaleFilter = itk::RescaleIntensityImageFilter< ImageType >::New();
auto rescaleFilter = itk::RescaleIntensityImageFilter< ImageTypeGeodesic >::New();
rescaleFilter->SetInput(geosOutput);
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(maxVal);
Expand All @@ -355,7 +368,7 @@ int main(int argc, char **argv)
}
else
{
auto calculator = itk::MinimumMaximumImageCalculator< ImageType >::New();
auto calculator = itk::MinimumMaximumImageCalculator< ImageTypeGeodesic >::New();
calculator->SetImage(geosOutput);
calculator->ComputeMaximum();
maxVal = calculator->GetMaximum();
Expand All @@ -364,18 +377,18 @@ int main(int argc, char **argv)
if (threshold != 0)
{
std::cout << "Performing threshold with threshold = " << threshold << ".\n";
auto thresholder = itk::BinaryThresholdImageFilter< ImageType, ImageType >::New();
auto thresholder = itk::BinaryThresholdImageFilter< ImageTypeGeodesic, ImageTypeGeodesic >::New();
thresholder->SetInput(geosOutput);
thresholder->SetOutsideValue(0);
thresholder->SetInsideValue(1);
thresholder->SetLowerThreshold(threshold);
thresholder->SetUpperThreshold(maxVal);
thresholder->SetLowerThreshold(0);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
thresholder->SetLowerThreshold(0);
thresholder->SetLowerThreshold(threshold);

This should use threshold instead of a hard-coded 0.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your suggestion is what it was before. The geodesic distance map that is created before the thresholding step is represented as increasing values away from the original segmentation/seeds. So setting it as you suggest will yield a 1 label in everything outside the predicted region of interest.

I tested this and compared to the GUI version with my changes, and this makes them consistent.

thresholder->SetUpperThreshold(threshold);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
thresholder->SetUpperThreshold(threshold);
thresholder->SetUpperThreshold(maxVal);

thresholder->Update();
geosOutput = thresholder->GetOutput();
}

std::cout << "Writing Output File:\n" + outputFile + "\n";
cbica::WriteImage< ImageType >(geosOutput, outputFile);
cbica::WriteImage< ImageTypeGeodesic >(geosOutput, outputFile);
}

std::cout << "Finished Successfully.\n";
Expand Down