marcomusy / vedo

A python module for scientific analysis of 3D data based on VTK and Numpy
https://vedo.embl.es
MIT License
1.98k stars 257 forks source link

Slice rendering without smoothing #1109

Closed sudmat closed 1 month ago

sudmat commented 1 month ago

Is there a way to render a slice without smoothing? The first picture is a slice rendered using vedo, second is the slice shown in ITK-SNAP. I would like something like the second. Thanks for the help. image

image

The code in vedo is below (I have already used the nearest interpolation), image

marcomusy commented 1 month ago

You may try with

slice.map_points_to_cells()
print(slice)
slice.cmap("binary_r", on="cells")
sudmat commented 1 month ago

You may try with

slice.map_points_to_cells()
print(slice)
slice.cmap("binary_r", on="cells")

Thank you, this works as expected. I have a related question, when I try to overlay a binary mask on an image, the mask color will “leak”. I assumed it was because the rendering was doing interpolation around the voxels, so the voxels between value=2 and value=0 will be something close to 1, therefore be colored the same as value=1. But even if I use cmap on cells, it still has the problem.Any ideas of this? Appreciate your help.

Current vedo rendering image

Expected (ITK-snap) image

LuT Building for the binary mask image

Rendering mask using binary_r cmap in vedo. I am sure that the mask only have values of 0, 1, and 2, no values in between. image

marcomusy commented 1 month ago

Hi, this happens because there is a different convention in the vtk for the Volume (wrapping vtkImageData). In this dataset the values are associated to the vertices (points) and not to the cells or voxels. For most practical application the difference is irrelevant but It can indeed lead to mismatches like the one you show.. Have a look at this toy example which should illustrate the problem and how to solve it:

import numpy as np
from vedo import Volume, CellCenters, show

# make up some fake data
X, Y, Z = np.mgrid[:4, :4, :2]
scalar_field = (X - 2) ** 2 + (Y - 2) ** 2 + (Z - 2) ** 2

vol = Volume(scalar_field.astype(int))
spacing = vol.spacing() # get the voxel size
# print("spacing", spacing)
# print('numpy array from Volume:', vol.pointdata)
# print('input_scalars', vol.pointdata['input_scalars'])

# extract a z-slice at index k=1
zslice = vol.zslice(k=1)
zslice.cmap("hot_r").lw(1).alpha(0.9).add_scalarbar3d()
# print("input_scalars", zslice.pointdata["input_scalars"])

# create a set of points at the cell centers and associate the scalar value
# corresponding to the bottom left corner of each voxel
cc = CellCenters(zslice).shift([-spacing[0] / 2, -spacing[1] / 2, 0])
cc.resample_data_from(zslice)

zslice2 = zslice.clone()
zslice2.celldata["pixel_value"] = cc.pointdata["input_scalars"]
print(zslice2.celldata["pixel_value"])

lego = vol.legosurface(vmin=0, vmax=10).wireframe()

show(
    [
        (vol, lego, CellCenters(vol)),
        (
            lego,
            cc,
            zslice,
            zslice.labels("id"),
            zslice.labels("cellid"),
        ),
        zslice2,
    ],
    N=3,
    axes=1,
)

Screenshot from 2024-05-06 22-04-40

note that this solution is only meaningful if you slice your volume along one of the cartesian axes.

sudmat commented 1 month ago

Excellent illustration. This solves my problem.