# Resampling along the z-axis of binary segmentation stack
**URL:** https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286
**Category:** Beginner Questions
**Tags:** python
**Created:** [July 23, 2021, 7:09am UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286 "2021-07-23T07:09:26Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@caniko](https://discourse.itk.org/u/caniko)
#### Post date: [July 23, 2021, 7:09am UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/1 "2021-07-23T07:09:26Z")
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:
```python
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.**
---
### Author:  [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [July 23, 2021, 3:50pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/2 "2021-07-23T15:50:51Z")
Perhaps take a look at [LabelImageGaussianInterpolateImageFunction](https://itk.org/Doxygen/html/classitk_1_1LabelImageGaussianInterpolateImageFunction.html) and the corresponding [example](https://itk.org/ITKExamples/src/Core/ImageFunction/ResampleSegmentedImage/Documentation.html). Does this help you?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [July 23, 2021, 4:22pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/3 "2021-07-23T16:22:25Z")
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:
```auto
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](http://www.itksnap.org/pmwiki/pmwiki.php) 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:
```auto
image_viewer = sitk.ImageViewer()
image_viewer.SetApplication(app='/Applications/ITK-SNAP.app/Contents/MacOS/ITK-SNAP')
image_viewer.Execute(binary_stack_itk)
```
---
### Author:  [@caniko](https://discourse.itk.org/u/caniko)
#### Post date: [July 23, 2021, 8:00pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/4 "2021-07-23T20:00:13Z")
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 🙂
---
### Author:  [@caniko](https://discourse.itk.org/u/caniko)
#### Post date: [July 23, 2021, 8:04pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/5 "2021-07-23T20:04:24Z")
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):
```python
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)
```
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [July 23, 2021, 9:06pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/6 "2021-07-23T21:06:13Z")
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](https://www.researchgate.net/publication/5518683_A_New_Method_for_Registration-Based_Medical_Image_Interpolation)).
---
### Author:  [@lassoan](https://discourse.itk.org/u/lassoan)
#### Post date: [July 25, 2021, 1:54am UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/7 "2021-07-25T01:54:29Z")
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](https://github.com/BRAINSia/BRAINSTools/tree/master/BRAINSResample).
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](https://github.com/KitwareMedical/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.
---
### Author:  [@caniko](https://discourse.itk.org/u/caniko)
#### Post date: [January 18, 2023, 11:37am UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/8 "2023-01-18T11:37:21Z")
I picked [ITKMorphologicalContourInterpolation](https://github.com/KitwareMedical/ITKMorphologicalContourInterpolation), thank you!
Do you know any methods that uses GPU? CUDA/OpenCL?
---
### Author:  [@lassoan](https://discourse.itk.org/u/lassoan)
#### Post date: [January 18, 2023, 3:15pm UTC](https://discourse.itk.org/t/resampling-along-the-z-axis-of-binary-segmentation-stack/4286/9 "2023-01-18T15:15:20Z")
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).