Working with large field medical images

Beginner Questions
1

Hello,

I’m working on a preprocessing and a segmentation pipeline for vessel extraction from micro-scan medical images. The image volumes are around 2048x2048x1400 and 2560x2560x1400, with 5,5 Gb and 8,5 Gb in 8 bits. In my pipeline i’m using some Frangi , k-mean and morphological filters and i’m struggling with processing such volumes.
I’ve made some research and tried to use stream divisions in the filters to solve the problem but nothing changed ( maybe i’m using it wrong). Apparently, loading and processing the whole image is the problem, so i tried to split in different stacks but filters such as Frangi has different results working like that.

So, i was wondering what options could be more interesting to deal with such volumes of medical images. Could you guys help me on this?
I’m using Python 3.7.6 on windows 10.

Thank you for the attention!

Lucca

2

Hello Lucca,

We are actively working on improving support for this use case. Can you share a sample dataset and a simple, self-contained set of code for your particular pipeline?

Thanks,
Matt

3

Hello Matt,

Unfortunately, the dataset that i’m working on is confidential, so i don’t have permission to send it. But i can send you a part of my code. The image is a .tiff file from a micro-ct scan and i’m using a Dell Precision Tower with 16Gb RAM.

Thanks!

import itk

inputFileName = 'filename'
PixelType = itk.UC
ImageType = itk.Image[PixelType, 3]
reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFileName)
reader.SetUseStreaming(True)        #Streaming to chunk image

sigmoidFilter = itk.SigmoidImageFilter[ImageType, ImageType].New()
sigmoidFilter.SetInput(reader.GetOutput())
sigmoidFilter.SetOutputMinimum(0)
sigmoidFilter.SetOutputMaximum(255)
sigmoidFilter.SetAlpha(7)
sigmoidFilter.SetBeta(63)
# sigmoidFilter.Update()        

'''Frangi Multi-scale vesselness filtering'''
FloatImageType = itk.Image[itk.F, 3]
castImageFilter = itk.CastImageFilter[ImageType, FloatImageType].New()    
castImageFilter.SetInput(reader.GetOutput())

HessianPixelType = itk.SymmetricSecondRankTensor[itk.D, 3]       
HessianImageType = itk.Image[HessianPixelType, 3]

objectness_filter = itk.HessianToObjectnessMeasureImageFilter[HessianImageType, FloatImageType].New()
objectness_filter.SetBrightObject(False)                        
objectness_filter.SetScaleObjectnessMeasure(True)
objectness_filter.SetObjectDimension(1)
objectness_filter.SetAlpha(1.0)                               
objectness_filter.SetBeta(1.0)
objectness_filter.SetGamma(5.0)                                

multi_scale_filter = itk.MultiScaleHessianBasedMeasureImageFilter[FloatImageType, HessianImageType, FloatImageType].New()
multi_scale_filter.SetInput(castImageFilter.GetOutput())
multi_scale_filter.SetHessianToMeasureFilter(objectness_filter)
multi_scale_filter.SetSigmaStepMethodToLogarithmic()
multi_scale_filter.SetSigmaMinimum(0.2)               
multi_scale_filter.SetSigmaMaximum(3.0)
multi_scale_filter.SetNumberOfSigmaSteps(5)

rescale_filter = itk.RescaleIntensityImageFilter[FloatImageType, ImageType].New()       
rescale_filter.SetInput(multi_scale_filter.GetOutput())
# rescale_filter.Update()


sigmoidFilter2 = itk.SigmoidImageFilter[ImageType, ImageType].New()
sigmoidFilter2.SetInput(rescale_filter.GetOutput())
sigmoidFilter2.SetOutputMinimum(0)
sigmoidFilter2.SetOutputMaximum(255)
sigmoidFilter2.SetAlpha(36)
sigmoidFilter2.SetBeta(130)

ImageType = itk.Image[itk.UC, 3]
'''K-mean pipeline for segmentation'''    
kmeansFilter = itk.ScalarImageKmeansImageFilter[ImageType, ImageType].New()
kmeansFilter.SetInput(sigmoidFilter2.GetOutput())
kmeansFilter.SetUseNonContiguousLabels(True)
userMeans = [10, 100]
for k in range(2):
    kmeansFilter.AddClassWithInitialMean(userMeans[k])
# kmeansFilter.Update()

'''Thresholding to make binary image'''
thresholdFilter = itk.ThresholdImageFilter[ImageType].New()
thresholdFilter.SetInput(kmeansFilter.GetOutput())
thresholdFilter.ThresholdAbove(50)
thresholdFilter.SetOutsideValue(255) 
# thresholdFilter.Update()

'''Closing to fill holes'''
StructuringElementType = itk.FlatStructuringElement[3]
structuringElement = StructuringElementType.Ball(2)
closingFilter = itk.BinaryMorphologicalClosingImageFilter[ImageType,ImageType,StructuringElementType].New()
closingFilter.SetInput(thresholdFilter.GetOutput())
closingFilter.SetKernel(structuringElement)
# closingFilter.Update()

'''Opening to erase unconnected structures'''
openfilter = itk.BinaryShapeOpeningImageFilter[ImageType].New()
openfilter.SetInput(closingFilter.GetOutput())
openfilter.SetFullyConnected(True)
openfilter.SetBackgroundValue(0)
openfilter.SetAttribute("NumberOfPixels")
openfilter.SetLambda(100)

'''Writing image on file'''
imwrite = itk.ImageFileWriter[ImageType].New()
imwrite.SetInput(openfilter.GetOutput())
imwrite.SetNumberOfStreamDivisions(150)         #Streaming to chunk image
imwrite.SetFileName('filename')
imwrite.Update()
1 Like
4

Hi Lucca,

The pipeline looks great!

A few quick fixes to limit memory consumption:

  • Use a file format for the input and output that supports streaming, e.g. MetaImage, .mha extension.
  • Separate filters, that do not support streaming, e.g. the morphological filters, into separate pipelines calls. If the filters that do not support streaming are at the end of the pipeline, then all memory will be required for the preceeding steps of the pipeline.
3 Likes
5

Hello,

The ITK filter to generate the Hessian is the HessianRecursiveGaussianImageFilter, this filter is not streamable and needs the whole image to generate its output. With a 8-bit scalar as input, the 64-bit double Hessian matrix would be 48x times the size of the input image.

In the SimpleITKFilters remote module there is a filter HessianImageFilter which directly computes the Hessian using discrete derivatives. This filter can enable the computation of Hessian based objecteness measures in a streamed fashion.

4 Likes
6

Thanks very much for the help Matt, i’ll do those changes to improve the code. Just one more question, tiff files and the k-mean filter they doesn’t suport streaming too?

Lucca

7

Hello,

Thank you very much for the help! I’ll try to implement this change to improve the code.

Lucca

8

@blowekamp did a lot of work on the itk.TIFFImageIO and is the expert on this – @blowekamp is their streaming support for TIFF’s?

For the k-means filter, SetImageRegion can be called to limit the input region required.

1 Like
9

Hrmm… It looks the the itk::TIFFImageIO does not support streaming after reviewing the code. I recall stream whole slices at a time. Perhaps I used separate tools to split multi-splice images into separate files, then used the ImageSeriesReader.

I believe the SCIFIOImageIO can read many TIFF files and has better support for streaming.

3 Likes