Thresholding largest/smallest component in ITK

Hello,

I’m trying to remove the largest island and also the smallest island.
I can get the components correctly, and they are sorted using relabel filter. but I want to threshold the smallest components, and keep the largest, and vice versa, but it seems its not working at all.

constexpr unsigned int Dimension = 3;
    using PixelType = unsigned int;
    using ImageType_cc = itk::Image<PixelType, Dimension>;

    using OutputImageType = itk::Image<unsigned int, Dimension>;

    vtkNew<vtkImageCast> castToUint;
    castToUint->SetInputData(input);
    castToUint->SetOutputScalarTypeToUnsignedInt();
    castToUint->Update();

    using FilterType_vtk = itk::VTKImageToImageFilter<ImageType_cc>;
    auto filter = FilterType_vtk::New();
    filter->SetInput(castToUint->GetOutput());

    filter->Update();


    using ConnectedComponentImageFilterType = itk::ConnectedComponentImageFilter<ImageType_cc, OutputImageType>;
    typedef itk::ConnectedComponentImageFilter<ImageType_cc, ImageType_cc> ConnectedComponentType;
    typename ConnectedComponentType::Pointer ccfilter = ConnectedComponentType::New();
    typedef itk::RelabelComponentImageFilter<ImageType_cc, ImageType_cc> RelabelComponentType;
    typename RelabelComponentType::Pointer relabel = RelabelComponentType::New();

    ccfilter->SetFullyConnected(false);
    ccfilter->SetInput(filter->GetOutput());
    relabel->SetInput(ccfilter->GetOutput());
    relabel->SetMinimumObjectSize(20);
    relabel->Update();

    std::cout << relabel->GetNumberOfObjects() << endl;
    std::cout << relabel->GetOriginalNumberOfObjects() << endl;

    using ImageTypeOutput = itk::Image<unsigned char, 3>;

    using FilterType = itk::BinaryThresholdImageFilter<ImageType_cc, ImageTypeOutput>;
    auto thresholder = FilterType::New();
    thresholder->SetInput(relabel->GetOutput());
    thresholder->SetLowerThreshold(0.0);
    thresholder->SetUpperThreshold(1.0);
    
    thresholder->Update();


    typename ImageTypeOutput::Pointer output = thresholder->GetOutput();

This is the issue. Relabel gives objects integer labels (1, 2, 3, …, n). Thus the output of this thresholder is going to be maybe just the largest object, or probably just an empty image (also known as all zero or all black). You might want to set the lower threshold equal to the number of large objects you want to filter out, and keep the upper threshold at the default (max value).

@dzenanz Can you show me sample code of your idea ? I didn’t get it
Also how to remove largest/smallest islands?

    unsigned nOriginal = relabel->GetOriginalNumberOfObjects();
    unsigned nObjects = relabel->GetNumberOfObjects();
    unsigned nLargeToFilter = 1 + nObjects / 5;
    // ...
    thresholder->SetInput(relabel->GetOutput());
    thresholder->SetLowerThreshold(nLargeToFilter);
    // thresholder->SetUpperThreshold(MAXINT); // do not set

Is that for largest or smallest ?

This is the original number of objects. Some/most of them will be removed by the relabel filter, due to relabel->SetMinimumObjectSize(20);. In my code fragment, this number is completely ignored.

1 Like

You should remove relabel->SetMinimumObjectSize(20);, if you want to be able to get them using this thresholding approach.

1 Like

@dzenanz Thanks so much, but Still don’t know how to keep either largest/ and remove smallest, If I want to make two functions one for Keeping Largest, and one for Removing Smallest.

1 Like

@dzenanz Should I reverse if I want to keep largest, remove smallest, the min, max of thresholding ?
using your code fragments, all the points are removed. :frowning:

// do not invoke SetMinimumObjectSize
if (removeLargest)
  thresholder->SetLowerThreshold(nToFilter);
else // removeSmallest
  thresholder->SetUpperThreshold(nToFilter);
1 Like

@dzenanz Thanks so much! Much apperciated, and what about KEEPing Largest?

@dzenanz I get NumberOfObjects = 1, is that normal ? I have a large volume with many blobs

SetMinimumObjectSize removes a lot of small objects. Remove this call if you want to get a better feel of what is going on.

I removed it, still getting numberofobjects = 1

3D images frequently have strong topological connectivity. You could explore your image using 3D Slicer. It’s Islands operations use ITK’s connected component and relabel filters.

@dzenanz
Does

 typedef itk::RelabelComponentImageFilter<ImageType_cc, ImageType_cc> RelabelComponentType;
    typename RelabelComponentType::Pointer relabel = RelabelComponentType::New();

Work with unsigned int /short ?

It works with whatever you instantiate it with, including short. But since the number of islands can be large, it is recommended to use a larger type, otherwise you might run into overflow.

@dzenanz it seems the thresholding is not reliable approach, can you show me an example of keeping largest, removing smallest ? I mean another way or how to use relable components

@dzenanz
Can you review the full pipeline please ? It doesn’t work

template <typename TPixel, typename TVolume, unsigned int VImageDimension>
bool VolumeProcessingItk<TPixel, TVolume, VImageDimension>::RemoveLargestislands(std::shared_ptr<DNA_Mask> mask)
{
    constexpr unsigned int Dimension = 3;

    using OutputImageType = itk::Image<unsigned int, Dimension>;

    using PixelType = unsigned char;
    using ImageType_cc = itk::Image<TPixel, Dimension>;

    using FilterType_threshold_operation = itk::BinaryThresholdImageFilter<ImageType_cc, ImageType_cc>;
    auto thresholder = FilterType_threshold_operation::New();

    thresholder->SetInput(mImage);
    thresholder->SetLowerThreshold(0);
    thresholder->SetUpperThreshold(500);
    thresholder->Update();


    using ConnectedComponentImageFilterType = itk::ConnectedComponentImageFilter<ImageType_cc, OutputImageType>;
    typedef itk::ConnectedComponentImageFilter<ImageType_cc, OutputImageType> ConnectedComponentType;
    typename ConnectedComponentType::Pointer ccfilter = ConnectedComponentType::New();
    typedef itk::RelabelComponentImageFilter<OutputImageType, ImageType_cc> RelabelComponentType;
    typename RelabelComponentType::Pointer relabel = RelabelComponentType::New();

    ccfilter->SetFullyConnected(false);
    ccfilter->SetInput(thresholder->GetOutput());
    relabel->SetInput(ccfilter->GetOutput());
   // relabel->SetMinimumObjectSize(500);
    relabel->Update();

    std::cout << relabel->GetNumberOfObjects() << endl;
    std::cout << relabel->GetOriginalNumberOfObjects() << endl;

    using FilterType_threshold_operation_final = itk::BinaryThresholdImageFilter<ImageType_cc, ImageType_cc>;
    auto thresholder_final = FilterType_threshold_operation_final::New();
    unsigned nOriginal = relabel->GetOriginalNumberOfObjects();
    unsigned nObjects = relabel->GetNumberOfObjects();
    unsigned nLargeToFilter = 1 + nObjects / 5;

    thresholder_final->SetInput(relabel->GetOutput());
    thresholder->SetLowerThreshold(nLargeToFilter);
   // thresholder_final->SetUpperThreshold(0);
    thresholder_final->Update();


    typename ImageType_cc::Pointer output = thresholder_final->GetOutput();

    DnaVolHelper<BinaryType, DNA_Mask> volInfo;

    typename ImageType_cc::RegionType region_filtered = output->GetLargestPossibleRegion();
    typename ImageType_cc::RegionType::SizeType  size_filtered = region_filtered.GetSize();

    volInfo.SetWidth((unsigned int)size_filtered[0]);
    volInfo.SetHeight((unsigned int)size_filtered[1]);
    volInfo.SetDepth((unsigned int)size_filtered[2]);

    typename ImageType_cc::SpacingType spacing;
    spacing = output->GetSpacing();
    volInfo.SetXSpacing(spacing[0]);
    volInfo.SetYSpacing(spacing[1]);
    volInfo.SetZSpacing(spacing[2]);

    itk::Point<itk::SpacePrecisionType, 3> origin = output->GetOrigin();
    volInfo.SetXOrigin(origin[0]);
    volInfo.SetYOrigin(origin[1]);
    volInfo.SetZOrigin(origin[2]);

    if (!volInfo.Populate(mask))
        return false;

    std::shared_ptr<BinaryType[]> pData = volInfo.GetData();

    using IteratorType = itk::ImageRegionIterator<ImageType_cc>;

    IteratorType  it_in(output, region_filtered); // select and image+region to visit
    auto it_out = pData.get();
    it_in.GoToBegin();  // move the iterator to the starting point
    while (!it_in.IsAtEnd())
    {
        *it_out = it_in.Get();  // get the value
        ++it_in;
        ++it_out;
    }

    return true;
}