3D Volume rendering
How to read a dicom dataset¶
import modules¶
In [2]:
import matplotlib.pyplot as plt;
import numpy as np;
from pydicom import dcmread;
from tqdm import tqdm;
example dicom dataset¶
In [3]:
ds = dcmread('/Users/taekjunkim/Documents/UW_PasupathyLab/Experiments/4.3P_Imaging/analysis_code/vol/test.dcm');
Keys in dicom dict dataset is hexadecimal. Therefore, individual item could be accessed as below¶
In [4]:
# to access Dataset.file_meta
print(ds.file_meta[(0x0002,0x0000)])
# to access other elements
print(ds[(0x0008,0x0000)])
(0002,0000) File Meta Information Group Length UL: 152 (0008,0000) Group Length UL: 940
Pixel data¶
Many DICOM SOP classes contain bulk pixel data, which typically represents medical imagery or 2D slices of a 3D volume. This data is most commonly found in the Pixel Data element, however it may be in Float Pixel Data or Double Float Pixel Data instead, depending on the SOP class. The table below lists these possible pixel data containing elements, although it’s important to note that only one may be present in any given dataset.
| Tag | Description | Keyword | VR |
|---|---|---|---|
| (7FE0,0008) | Float Pixel Data | FloatPixelData | OF |
| (7FE0,0009) | Double Pixel Data | DoubleFloatPixelData | OD |
| (7FE0,0010) | Pixel Data | PixelData | OB or OW |
In [5]:
ds.group_dataset(0x7FE0)
Out[5]:
(7FE0,0000) Group Length UL: 13107212 (7FE0,0010) Pixel Data OB: Array of 13107200 elements
A dataset with pixel data should always contain group 0x0028 Image Pixel module elements, which are needed to properly interpret the encoded pixel data byte stream:
the below indicates that this dataset contains 100 grayscale image frames with dimensions 256 x 256
In [6]:
ds.group_dataset(0x0028)
Out[6]:
(0028,0000) Group Length UL: 186 (0028,0002) Samples per Pixel US: 1 (0028,0004) Photometric Interpretation CS: 'MONOCHROME2' (0028,0008) Number of Frames IS: '100' (0028,0010) Rows US: 256 (0028,0011) Columns US: 256 (0028,0100) Bits Allocated US: 16 (0028,0101) Bits Stored US: 12 (0028,0102) High Bit US: 11 (0028,0103) Pixel Representation US: 0 (0028,0301) Burned In Annotation CS: 'NO' (0028,2110) Lossy Image Compression CS: '00' (0028,3003) LUT Explanation LO: 'Philips Real World Value Mapping' (0028,9001) Data Point Rows UL: 1 (0028,9002) Data Point Columns UL: 0
pixel spacing, spacing between slices
Converting the pixel data to a numpy ndarray¶
In [7]:
from pydicom.pixels import pixel_array
arr = pixel_array(ds);
print(np.shape(arr));
print(arr.dtype);
(100, 256, 256) uint16
Create a datacube¶
- Both pixel spacing and slice spacing are taken into account.
- As slice spacing is larger than pixel spacing, interpolation was applied to the data points.
In [8]:
from scipy.interpolate import interpn
# found these information from the meta file
pixel_spacing = [0.5, 0.5];
slice_spacing = 1;
# check x,y,z dimension
Nx, Ny, Nz = np.shape(arr);
Lx = slice_spacing * Nx;
Ly = pixel_spacing[0] * Ny;
Lz = pixel_spacing[1] * Nz;
# x,y,z points
x = np.linspace(-Lx/2, Lx/2, Nx);
y = np.linspace(-Ly/2, Ly/2, Ny);
z = np.linspace(-Lz/2, Lz/2, Nz);
points = (x, y, z);
N = int(np.max([Lx, Ly, Lz]));
c = np.linspace(-N/2, N/2, N);
qx, qy, qz = np.meshgrid(c,c,c);
qi = np.array([qx.ravel(), qy.ravel(), qz.ravel()]).T
datacube = interpn(points, arr, qi, method='linear', bounds_error=False, fill_value=0).reshape((N,N,N));
# normalize datacube (min: 0, max: 1)
datacube = datacube - np.min(datacube);
datacube = datacube / np.max(datacube);
print(np.shape(datacube));
(128, 128, 128)
Display the datacube¶
In [8]:
import matplotlib.animation as animation
from IPython.display import HTML
fig = plt.figure(figsize=(8,3));
plt.subplot(1,3,1);
image1 = plt.imshow(np.squeeze(datacube[0,:,:]),vmin=0,vmax=1,animated=True,origin='lower',cmap='binary');
plt.title('Top view');
plt.axis('off');
plt.subplot(1,3,2);
image2 = plt.imshow(np.squeeze(datacube[:,0,:]),vmin=0,vmax=1,animated=True,origin='lower',cmap='binary');
plt.title('Front view');
plt.axis('off');
plt.subplot(1,3,3);
image3 = plt.imshow(np.squeeze(datacube[:,:,0]),vmin=0,vmax=1,animated=True,origin='lower',cmap='binary');
plt.title('Side view');
plt.axis('off');
plt.close()
def animate(i):
image1.set_data(np.squeeze(datacube[i,:,:]));
image2.set_data(np.squeeze(datacube[:,i,:]));
image3.set_data(np.squeeze(datacube[:,:,i]));
return image1, image2, image3,
anim = animation.FuncAnimation(fig, animate, frames=np.shape(datacube)[0], interval=100, blit=True);
HTML(anim.to_jshtml())
Out[8]:
Histogram of pixel intensity¶
- Plot a histogram of your image's intensity values. This will show you the distribution of pixel intensities.
- Look for Peaks and Valleys: Ideally, you'll see peaks corresponding to different objects or regions in your image, separated by valleys.
- Choose a Valley: A good threshold often lies within a valley between two peaks. This allows you to separate pixels belonging to different groups.
In [9]:
plt.figure(figsize=(6,2))
plt.subplot(1,2,1)
plt.hist(datacube.flatten(),np.arange(0,1.05,0.05))
plt.gca().spines[['right', 'top']].set_visible(False)
plt.gca().tick_params(axis='both',labelsize=7);
plt.title('pixel values: raw signal');
plt.xlabel('0-1 intensity');
plt.ylabel('Number of pixels');
plt.subplot(1,2,2)
plt.hist(np.log(254*datacube+0.001).flatten(), np.arange(-7,6.5,0.5));
plt.gca().spines[['right', 'top']].set_visible(False)
plt.gca().tick_params(axis='both',labelsize=7);
plt.title('pixel values: log signal');
plt.xlabel('log(0-255 intensity)');
plt.ylabel('Number of pixels');
plt.tight_layout()
View the datacube in 3D using matplotlib¶
In [10]:
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np
from skimage import measure
levels = [0.2, 0.4];
fig = plt.figure(figsize=(4, 4))
ax = fig.add_subplot(projection='3d')
for ci, level in enumerate(levels):
p = datacube.transpose(2,1,0);
p = p[::2,::2,::2];
verts, faces, normals, values = measure.marching_cubes(p, level, method='lewiner');
mesh = Poly3DCollection(verts[faces], alpha=0.05)
if ci==0:
face_color = [1, 0.5, 0.5]
elif ci==1:
face_color = [0.5, 1, 0.5]
elif ci==2:
face_color = [0.5, 0.5, 1]
#face_color = [0.5, 0.5, 1]
mesh.set_facecolor(face_color)
ax.add_collection3d(mesh)
ax.set_xlim(0, p.shape[0])
ax.set_ylim(0, p.shape[1])
ax.set_zlim(0, p.shape[2])
ax.tick_params(axis='both',labelsize=7);
plt.close()
def animate_rot(frame):
ax.view_init(elev=10., azim=frame, roll=None);
return fig
# Create the animation
anim = animation.FuncAnimation(fig, animate_rot, frames=np.linspace(0, 360, 18), interval=300)
HTML(anim.to_jshtml())
Out[10]:
View the datacube in 3D using plotly¶
in the
go.Volumefunction
- Values below isomin and above isomax are not rendered
- The range between isomin and isomax is divided into multiple isosurfaces for volume rendering
- The number of surfaces rendered is controlled by the surface_count parameter
the output figure was not included due to its large size.
In [ ]:
import plotly.graph_objects as go
import numpy as np
import plotly.offline as plo;
N = datacube.shape[0]
c = np.linspace(-N/2, N/2, N)
x1, y1, z1 = np.meshgrid(c,c,c)
fig = go.Figure(data=go.Volume(
x=x1[::2,::2,::2].flatten(),
y=y1[::2,::2,::2].flatten(),
z=z1[::2,::2,::2].flatten(),
value=datacube[::2,::2,::2].transpose((1,2,0)).flatten(),
isomin=0.2,
isomax=0.5,
opacity=0.1, # needs to be small to see through all surfaces
surface_count=20,
colorscale='viridis',
))
fig.update_layout(
width = 600,
height = 400,
scene = dict(
xaxis = dict(visible=True),
yaxis = dict(visible=True),
zaxis = dict(visible=True)
),
scene_camera = dict(
eye=dict(x=0.1, y=-1.25, z=0.0)
)
)
fig.update_traces(showscale=True)
fig.show(renderer="notebook")
Volume rendering with a transfer function¶
Transfer functions in volume rendering serve as a crucial mapping mechanism, translating raw voxel data values into visual properties like color and opacity.
This transformation allows for the selective emphasis of specific features within the volumetric dataset.
Got help from this medium post
https://medium.com/swlh/create-your-own-volume-rendering-with-python-655ca839b097
In [12]:
# log intensity centered at 2, 3, 4 as key points
def transferFunction(x):
r = 1.0*np.exp( -(x - 4.5)**2/0.5 ) + 0.2*np.exp( -(x - 3.0)**2/0.5 ) + 0.2*np.exp( -(x - 1.5)**2/0.5 )
g = 1.0*np.exp( -(x - 4.5)**2/0.5 ) + 1.0*np.exp( -(x - 3.0)**2/0.5 ) + 0.2*np.exp( -(x - 1.5)**2/0.5 )
b = 0.2*np.exp( -(x - 4.5)**2/0.5 ) + 0.2*np.exp( -(x - 3.0)**2/0.5 ) + 1.0*np.exp( -(x - 1.5)**2/0.5 )
a = 0.05*np.exp( -(x - 4.5)**2/0.5 ) + 0.03*np.exp( -(x - 3.0)**2/0.5 ) + 0.01*np.exp( -(x - 1.5)**2/0.5 )
return r,g,b,a
using matplotlib¶
In [13]:
N = datacube.shape[0]
c = np.linspace(-N/2, N/2, N)
x1, y1, z1 = np.meshgrid(c,c,c)
points = (c, c, c);
imgs = [];
Nangles = 18
for i in tqdm(np.arange(Nangles)):
# Rotate datacube: slightly looking down
x2 = x1 * np.cos(-2*np.pi/18) - y1 * np.sin(-2*np.pi/18);
y2 = x1 * np.sin(-2*np.pi/18) + y1 * np.cos(-2*np.pi/18);
z2 = z1;
# Camera angle
angle = 2*np.pi * i / Nangles
qxR = x2
qyR = y2 * np.cos(angle) - z2 * np.sin(angle)
qzR = y2 * np.sin(angle) + z2 * np.cos(angle)
qi = np.array([qxR.ravel(), qyR.ravel(), qzR.ravel()]).T
# Interpolate onto Camera Grid
camera_grid = interpn(points, datacube, qi, method='linear', bounds_error=False, fill_value=0).reshape((N,N,N))
# Do Volume Rendering
image = np.ones((camera_grid.shape[1],camera_grid.shape[2],3))
for dataslice in camera_grid:
r,g,b,a = transferFunction(np.log(dataslice*254+0.001))
image[:,:,0] = a*r + (1-a)*image[:,:,0]
image[:,:,1] = a*g + (1-a)*image[:,:,1]
image[:,:,2] = a*b + (1-a)*image[:,:,2]
image = np.clip(image,0.0,1.0)
imgs.append(image);
fig, ax = plt.subplots(figsize=(4,4))
im = ax.imshow(imgs[0],vmin=0,vmax=1,animated=True,origin='lower');
ax.set_axis_off()
plt.close()
def animate(i):
im.set_data(imgs[i])
return im,
anim = animation.FuncAnimation(fig, animate, frames=Nangles, interval=300, blit=True);
HTML(anim.to_jshtml())
0%| | 0/18 [00:00<?, ?it/s]
100%|██████████| 18/18 [00:25<00:00, 1.44s/it]
Out[13]: