Resampling along the z-axis of binary segmentation stack

caniko · July 23, 2021, 7:09am

I have a binary image stack or space of binarized voxels (either black or white voxels).

Since my data is from electron microscopy the spacing of the pixels are (1, 1, Z) where (x, y, z), and Z is an arbitrary value defined by the limitation of preparatory steps before performing the imaging with the electron microscope.

My current issue is resampling. With this code:

def resample(binary_stack: np.ndarray, z_spacing: float) -> np.ndarray:
    binary_stack_itk = sitk.GetImageFromArray(binary_stack.astype(np.float32))

    new_size = list(binary_stack_itk.GetSize())
    new_size[2] = round(z_spacing * new_size[2])

    return sitk.GetArrayFromImage(sitk.Resample(
        binary_stack_itk,
        new_size,
        sitk.Transform(),
        sitk.sitkLinear,
        binary_stack_itk.GetOrigin(),
        (1.0, 1.0, 1.0 / z_spacing),
        binary_stack_itk.GetDirection(),
        0,
        binary_stack_itk.GetPixelID(),
    )).astype(bool)

If I rotate the mesh 90 degrees (np.rot90(result)) we can inspect the array along the zy-plane. The following is an arbitrary position of somewhere in the middle of the zy-plane in one of my resampled binary stacks:

This is not smooth at all! I tried changing methods from sitk.sitkLinear sitk.sitkBSplineResamplerOrder3, and get this when rotated:

Before resample, segmantation masks stacked on top of each other (same location as the ones above), also rotated (np.rot90(stack)):

The resample is way too simple with respect to what I am looking for. I was hoping to find a resample that would be smoother, natively. I would like the resample to be connecting the edges of the segmentations along zy plane so that they are connected along their edges, not something like the one we see in sitkLinear, where the previous edge is traced to the next edge:

I would like something like this instead:

Where everything above the white line is True!

Any resample method that I missed in the Simple ITK liberary, or is what I am looking for non-existant? I have to avoid functions that would require variable tuning of any sort as these introduce bias to the data.

dzenanz · July 23, 2021, 3:50pm

Perhaps take a look at LabelImageGaussianInterpolateImageFunction and the corresponding example. Does this help you?

zivy · July 23, 2021, 4:22pm

Hello @caniko,

I suspect the issue you are facing is much more complicated than what you think, it isn’t really a resampling problem. Think a cylinder whose main axis is at an angle to the z direction (fiducial that is often used in microscopy). As a result you have a 2D circle in each slice and there is no overlap between their x-y coordinates. You will have to incorporate prior knowledge to correctly reconstruct this object, no resampling algorithm will be able to do it correctly without additional information. The standard resampling assumes that the sampling rate along the axis in which you are performing resampling is sufficient otherwise there will be aliasing (Nyquist theorem).

Some observations with respect to the code, I believe there is a bug in computing the new_size, should be:

new_size[2] = round(z_spacing * (new_size[2]-1) *  binary_stack_itk.GetSpacing()[2])

Also, for visualization and debugging purposes would highly recommend using ITK-SNAP for visualization. This will allow you to look at the binary volumes using x-y, x-z, y-z planes at the same time. See code below:

image_viewer = sitk.ImageViewer()
image_viewer.SetApplication(app='/Applications/ITK-SNAP.app/Contents/MacOS/ITK-SNAP')
image_viewer.Execute(binary_stack_itk)

caniko · July 23, 2021, 8:00pm

That function is similar to sitkLinear when there is only one label, but I might be wrong. Have already tried it by the way, not good enough

caniko · July 23, 2021, 8:04pm

Thank you!

The new_size bug you found is indeed a bug. Moreover, it wasn’t an issue as binary_stack_itk.GetSpacing() -> (1, 1, 1)

I understand my problem more clearly with respect to the Nyquist theorem. Now that I understand it, I am sad… Would you just stick to sitkLinear in that case? Or is there a better interpolation method for me? Right now I am, as you pointed to, performing aliasing using Gaussian filtering followed by threshold (binarization):

def resample(binary_stack: np.ndarray, z_spacing: float) -> np.ndarray:
    binary_stack_itk = sitk.GetImageFromArray(binary_stack.astype(np.float64))

    new_size = list(binary_stack_itk.GetSize())
    new_size[2] = round(z_spacing * (new_size[2] - 1) * binary_stack_itk.GetSpacing()[2])

    resampled = sitk.Resample(
        binary_stack_itk,
        new_size,
        sitk.Transform(),
        sitk.sitkLinear,
        binary_stack_itk.GetOrigin(),
        (1.0, 1.0, 1.0 / z_spacing),
        binary_stack_itk.GetDirection(),
        0,
        binary_stack_itk.GetPixelID(),
    )

    smoother = sitk.SmoothingRecursiveGaussian(
        resampled,
        2.5,    # sigma (x, y, z)
        True    # normalizeAcrossScale
    )

    gauss_array = sitk.GetArrayFromImage(smoother).astype(np.float64)
    binarized = sitk.BinaryThreshold(
        smoother, upperThreshold=1.0,
        lowerThreshold=np.median(gauss_array[gauss_array > 0.05])
    )

    return sitk.GetArrayFromImage(binarized).astype(bool)

zivy · July 23, 2021, 9:06pm

Hello @caniko,

As the interpolation is along the z direction and this is likely not the appropriate direction for interpolating the imaged object you might consider a more sophisticated approach which involves registration and using that information during interpolation.

See this paper for additional details: D. Frakes et al. " A New Method for Registration-Based Medical Image Interpolation", IEEE TMI (available on researchgate).

lassoan · July 25, 2021, 1:54am

If you don’t need to recover high-frequency information along the K axis (it is enough to smoothly interpolate along the slices) then you have a couple of simple options.

A. Anti-aliasing filter on the image. You can use a windowed sinc kernel, with strong smoothing along the K axis and no smoothing along I and J axes.

B. Compute signed distance map and threshold that. This is implemented in the resample module in BRAINSTools.

C. Convert the labelmap to closed surface, apply low-pass filter on the surface (for example, using vtkWindowedSincPolyDataFilter), then convert back to labelmap.

D. Insert empty slices between the segmented slices and run ITKMorphologicalContourInterpolation filter to smoothly interpolate between them.

You can try all these options using Python scripting using ITK and VTK classes, or by a couple of clicks using the GUI of 3D Slicer.

caniko · January 18, 2023, 11:37am

I picked ITKMorphologicalContourInterpolation, thank you!

Do you know any methods that uses GPU? CUDA/OpenCL?

lassoan · January 18, 2023, 3:15pm

ITKMorphologicalContourInterpolation is practically real-time on the CPU.

OpenCL is dead and CUDA-capable GPU is only available for a tiny fraction of users. Therefore, if you want your software to run on most user’s computer then it can only rely on the CPU or commonly available GPU APIs (such as GLSL shaders or WebGPU).