#!/usr/bin/env python
# -*- coding: utf-8 -*-
# tools.dicom_dose_sum.py
"""
Class for summing dose grids
"""
# This file is part of DVH Analytics, released under a BSD license.
# See the file LICENSE included with this distribution, also
# available at https://github.com/cutright/DVH-Analytics
import numpy as np
from os.path import isfile
import pydicom
from scipy.ndimage import map_coordinates
class DoseGrid:
"""
Class to easily access commonly used attributes of a DICOM dose grid and perform summations
Example: Add two dose grids
grid_1 = DoseGrid(dose_file_1)
grid_2 = DoseGrid(dose_file_2)
grid_1.add(grid_2)
grid_1.save_dcm(some_file_path)
"""
def __init__(self, rt_dose):
"""
:param rt_dose: an RT Dose dicom dataset or file_path
:type rt_dose: pydicom.FileDataset
"""
self.ds = self.__validate_input(rt_dose)
if self.ds:
self.x_axis = np.arange(self.ds.Columns) * self.ds.PixelSpacing[0] + self.ds.ImagePositionPatient[0]
self.y_axis = np.arange(self.ds.Rows) * self.ds.PixelSpacing[1] + self.ds.ImagePositionPatient[1]
self.z_axis = np.array(self.ds.GridFrameOffsetVector) + self.ds.ImagePositionPatient[2]
# x and z are swapped in the pixel_array
self.dose_grid = np.swapaxes(self.ds.pixel_array * self.ds.DoseGridScaling, 0, 2)
@staticmethod
def __validate_input(rt_dose):
"""Ensure provided input is either an RT Dose pydicom.FileDataset or a file_path to one"""
if type(rt_dose) is pydicom.FileDataset:
if rt_dose.Modality.lower() == 'rtdose':
return rt_dose
print("The provided pydicom.FileDataset is not RTDOSE")
return
elif isfile(rt_dose):
try:
rt_dose_ds = pydicom.read_file(rt_dose)
if rt_dose_ds.Modality.lower() == 'rtdose':
return rt_dose_ds
print('The provided file_path points to a DICOM file, but it is not an RT Dose file.')
except Exception as e:
print(e)
print('The provided input is neither a pydicom.FileDataset nor could it be read by pydicom.')
return
####################################################
# Basic properties
####################################################
@property
def shape(self):
"""Get the x, y, z dimensions of the dose grid"""
return tuple([self.ds.Columns, self.ds.Rows, len(self.ds.GridFrameOffsetVector)])
@property
def axes(self):
"""Get the x, y, z axes of the dose grid (in mm)"""
return [self.x_axis, self.y_axis, self.z_axis]
@property
def scale(self):
"""Get the dose grid resolution (xyz)"""
return np.array([self.ds.PixelSpacing[0],
self.ds.PixelSpacing[1],
self.ds.GridFrameOffsetVector[1] - self.ds.GridFrameOffsetVector[0]])
@property
def offset(self):
"""Get the coordinates of the dose grid origin (mm)"""
return np.array(self.ds.ImagePositionPatient, dtype='float')
@property
def points(self):
"""Get all of the points in the dose grid"""
y, z, x = np.meshgrid(self.y_axis, self.z_axis, self.x_axis) # iterate through x, then y, then z
points = np.vstack((x.ravel(), y.ravel(), z.ravel()))
return points
####################################################
# Tools
####################################################
def is_coincident(self, other):
"""Check dose grid coincidence, if True a direct summation is appropriate"""
return self.ds.ImagePositionPatient == other.ds.ImagePositionPatient and \
self.ds.pixel_array.shape == other.ds.pixel_array.shape and \
self.ds.PixelSpacing == other.ds.PixelSpacing and \
self.ds.GridFrameOffsetVector == other.ds.GridFrameOffsetVector
def set_pixel_data(self, pixel_data):
"""
Update the PixelData in the pydicom.FileDataset
:param pixel_data: a 3D numpy array with the same dimensions as self.ds.pixel_array
:type pixel_data: np.array
"""
self.ds.BitsAllocated = 32
self.ds.BitsStored = 32
self.ds.HighBit = 31
self.ds.PixelData = np.uint32(pixel_data / self.ds.DoseGridScaling).tostring()
def save_dcm(self, file_path):
"""Save the pydicom.FileDataset to file"""
self.ds.save_as(file_path)
def xyz_to_ijk(self, xyz):
"""Convert xyz points (e.g., self.points) into ijk space"""
for dim in range(3):
xyz[dim] = (xyz[dim] - self.offset[dim]) / self.scale[dim]
return xyz
####################################################
# Dose Summation
####################################################
def add(self, other, interp_order=1):
"""
Add another 3D dose grid to this 3D dose grid, with interpolation if needed
:param other: another DoseGrid
:type other: DoseGrid
:param interp_order: the order to be passed to scipy.ndimage.map_coordinates, default is trilinear interpolation
:type interp_order: int
"""
if self.is_coincident(other):
self.direct_sum(other)
else:
self.interp_sum(other, order=interp_order)
def direct_sum(self, other):
"""Directly sum two dose grids (only works if both are coincident)"""
dose_sum = self.ds.pixel_array * self.ds.DoseGridScaling + other.ds.pixel_array * other.ds.DoseGridScaling
self.set_pixel_data(dose_sum)
def interp_sum(self, other, order=1):
"""
Interpolate the other dose grid to this dose grid's axes, then directly sum
:param other: another DoseGrid
:type other: DoseGrid
:param order: interpolation order: default 0: nearest point
1: linear
2 to 5: spline
:type order: int
"""
ijk_points = other.xyz_to_ijk(self.points)
other_grid = map_coordinates(input=other.dose_grid,
coordinates=ijk_points,
order=order).reshape(self.shape)
self.dose_grid += other_grid
self.set_pixel_data(np.swapaxes(self.dose_grid, 0, 2))
Does anyone know of any open-source code for dose grid summation with interpolation?I've tried the dicompyler-plugin, but it crashes for my dcm files... we're still working through it here: https://groups.google.com/forum/#!topic/dicompyler/qkU2CtYzgLgI posted some different code in that thread, which works correctly (using scipy.interpolate.RegularGridInterpolator), but it's about ~400x slower than the code below. Even though the code below is incorrect, it seems to be doing the same number of calculations as it should... so I'm hopeful the corrected code would be just as fast.Is it possible the code below is correct and the scipy.ndimage.map_coordinates interpolation is just way off? scipy.interpolate.RegularGridInterpolator came very close to Eclipse.Note that I swapped the order of the meshgrid arguments in DoseGrid.points to get the vstack to iterate though x values, then y, then z. Using np.meshgrid(self.x_axis, self.y_axis, self.z_axis) didn't work either.
<code>
>>> import numpy as np
>>> x, y, z = [0, 1, 2], [-3, -2, -1], [10, 11, 12]
>>> Y, Z, X = np.meshgrid(y, z, x)
>>> points = np.vstack((X.ravel(), Y.ravel(), Z.ravel()))
>>> points.transpose()
array([[ 0, -3, 10],
[ 1, -3, 10],
[ 2, -3, 10],
[ 0, -2, 10],
[ 1, -2, 10],
[ 2, -2, 10],
[ 0, -1, 10],
[ 1, -1, 10],
[ 2, -1, 10],
[ 0, -3, 11],
[ 1, -3, 11],
[ 2, -3, 11],
[ 0, -2, 11],
[ 1, -2, 11],
[ 2, -2, 11],
[ 0, -1, 11],
[ 1, -1, 11],
[ 2, -1, 11],
[ 0, -3, 12],
[ 1, -3, 12],
[ 2, -3, 12],
[ 0, -2, 12],
[ 1, -2, 12],
[ 2, -2, 12],
[ 0, -1, 12],
[ 1, -1, 12],
[ 2, -1, 12]])
import numpy as np
import pydicom
from scipy.interpolate import RegularGridInterpolator
class DoseGrid:
def __init__(self, ds):
self.ds = ds
self.__set_axes()
# x and z are swapped in the pixel_array
self.dose_grid = np.swapaxes(self.ds.pixel_array * self.ds.DoseGridScaling, 0, 2)
def __set_axes(self):
self.x_axis = np.arange(self.ds.Columns) * self.ds.PixelSpacing[0] + self.ds.ImagePositionPatient[0]
self.y_axis = np.arange(self.ds.Rows) * self.ds.PixelSpacing[1] + self.ds.ImagePositionPatient[1]
self.z_axis = np.array(self.ds.GridFrameOffsetVector) + self.ds.ImagePositionPatient[2]
def is_point_inside_grid(self, xyz):
for a, axis in enumerate(self.axes):
if not (axis[0] <= xyz[a] <= axis[-1]):
return False
return True
def is_coincident(self, other):
return self.ds.ImagePositionPatient == other.ds.ImagePositionPatient and \
self.ds.pixel_array.shape == other.ds.pixel_array.shape and \
self.ds.PixelSpacing == other.ds.PixelSpacing and \
self.ds.GridFrameOffsetVector == other.ds.GridFrameOffsetVector
def direct_sum(self, other):
dose_sum = self.ds.pixel_array * self.ds.DoseGridScaling + other.ds.pixel_array * other.ds.DoseGridScaling
self.set_pixel_array(dose_sum)
def set_pixel_array(self, pixel_data):
self.ds.BitsAllocated = 32
self.ds.BitsStored = 32
self.ds.HighBit = 31
self.ds.PixelData = np.uint32(pixel_data / self.ds.DoseGridScaling).tostring()
def update_pixel_array(self):
self.set_pixel_array(np.swapaxes(self.dose_grid, 0, 2))
@property
def grid_size(self):
return [self.ds.Columns, self.ds.Rows, len(self.ds.GridFrameOffsetVector)]
@property
def axes(self):
return [self.x_axis, self.y_axis, self.z_axis]
def get_xyz(self, ijk):
"""Convert an ijk coordinate into xyz space"""
i, j, k = ijk[0], ijk[1], ijk[2]
return np.array([self.x_axis[i], self.y_axis[j], self.z_axis[k]])
def get_ijk(self, xyz):
"""Convert an xyz coordinate into ijk space"""
return [int(np.interp(xyz[a], axis, np.arange(self.grid_size[a]))) for a, axis in enumerate(self.axes)]
def add_dose(self, xyz, dose):
if self.is_point_inside_grid(xyz):
ijk = self.get_ijk(xyz)
self.dose_grid[ijk[0], ijk[1], ijk[2]] += dose
class DoseInterpolator(DoseGrid):
"""Inherit DoseGrid, separate class so a RegularGridInterpolator isn't created for every DoseGrid"""
def __init__(self, ds):
DoseGrid.__init__(self, ds)
self.interpolator = RegularGridInterpolator(points=self.axes, values=self.dose_grid, bounds_error=False)
def get_dose(self, xyz):
return self.interpolator([xyz])[0]
def get_doses(self, list_of_xyz):
list_of_xyz = [xyz for xyz in list_of_xyz if self.is_point_inside_grid(xyz)]
return self.interpolator(list_of_xyz), list_of_xyz
def add_dose_grids(ds1, ds2):
if type(ds1) is not pydicom.FileDataset:
ds1 = pydicom.read_file(ds1)
if type(ds2) is not pydicom.FileDataset:
ds2 = pydicom.read_file(ds2)
dose_1 = DoseGrid(ds1)
dose_2 = DoseInterpolator(ds2)
if dose_1.is_coincident(dose_2):dose_1.direct_sum(dose_2)
print("PlanSum: Using direct summation")
else:
slice_count = float(len(dose_1.z_axis))
for slice_counter, z in enumerate(dose_1.z_axis):
print("%0.1f%%" % (100 * slice_counter / slice_count))
points = []
for x in dose_1.x_axis:
for y in dose_1.y_axis:
points.append([x, y, z])
doses, points = dose_2.get_doses(points)
for p, point in enumerate(points):
dose_1.add_dose(point, doses[p])
dose_1.update_pixel_array()
return dose_1.ds
interpolator = RegularGridInterpolator(points=other.axes, values=other.dose_grid, bounds_error=False, fill_value=0)
other_grid = interpolator(self.points).reshape(self.shape)
self.dose_grid += other_grid
self.x_axis = np.arange(self.ds.Columns) * self.ds.PixelSpacing[0] + self.ds.ImagePositionPatient[0]
self.y_axis = np.arange(self.ds.Rows) * self.ds.PixelSpacing[1] + self.ds.ImagePositionPatient[1]
self.z_axis = np.array(self.ds.GridFrameOffsetVector) + self.ds.ImagePositionPatient[2]
self.axes =[self.x_axis, self.y_axis, self.z_axis]
self.points = []for x in self.x_axis: for y in self.y_axis: for z in self.z_axis: self.points.append([x, y, z])
self.shape = tuple([self.ds.Columns, self.ds.Rows, len(self.ds.GridFrameOffsetVector)])
self.ds.BitsAllocated = 32
self.ds.BitsStored = 32
self.ds.HighBit = 31
self.ds.PixelData = np.uint32(np.swapaxes(self.dose_grid/ self.ds.DoseGridScaling, 0, 2).tostring())
self.points = []
for x in self.x_axis:
for y in self.y_axis:
for z in self.z_axis:
self.points.append([x, y, z])
y, x, z = np.meshgrid(self.y_axis, self.x_axis, self.z_axis)
stack = np.vstack((x.ravel(), y.ravel(), z.ravel()))
self.points = stack.transpose()
import numpy as np
from datetime import datetime
X = np.linspace(0, 100, 101)
Y = np.linspace(100, 201, 101)
Z = np.linspace(200, 301, 101)
start_list = datetime.now()
points = []
for x in X:
for y in Y:
for z in Z:
points.append([x, y, z])
end_list = datetime.now()
start_np = datetime.now()
y, x, z = np.meshgrid(Y, X, Z)
stack = np.vstack((x.ravel(), y.ravel(), z.ravel()))
points_np = stack.transpose()
end_np = datetime.now()
delta = {'list': end_list - start_list,
'numpy': end_np - start_np}
factor = delta['list'].total_seconds() / delta['numpy'].total_seconds()
print('List method (s): %0.3f' % delta['list'].total_seconds())
print('Numpy method (s): %0.3f' % delta['numpy'].total_seconds())
print('NumPy method was %0.3f times faster' % factor)
print(points == points_np.tolist())
List method (s): 0.918
Numpy method (s): 0.021
NumPy method was 43.715 times faster
True
from dicom_dose_sum import DoseGrid
grid_1 = DoseGrid(dose_file_1)
grid_2 = DoseGrid(dose_file_2)
grid_sum = grid_1 + grid_2 # or grid_1.add(grid_1) to reuse grid_1
grid_sum.save_dcm(some_file_path)
def interp_sum_new(self, other):
axes = [(np.array(axis) - other.offset[a]) / other.scale[a] for a, axis in enumerate(self.axes)]
j, i, k = np.meshgrid(axes[1], axes[0], axes[2])
points = np.vstack((i.ravel(), j.ravel(), k.ravel()))
other_grid = map_coordinates(input=other.dose_grid, coordinates=points).reshape(self.shape)
self.dose_grid += other_grid
self.set_pixel_data()
from dvha.tools.dicom_dose_sum import DoseGrid
from datetime import datetime
import numpy as np
# RegularGridInterpolator method
dose_1 = DoseGrid('ANON22831/dose1.dcm')
dose_2 = DoseGrid('ANON22831/dose2.dcm')
start_reg_grid_interp = datetime.now()
dose_1.interp_sum(dose_2)
end_reg_grid_interp = datetime.now()
# map_coordinates method
dose_1_new = DoseGrid('ANON22831/dose1.dcm')
dose_2_new = DoseGrid('ANON22831/dose2.dcm')
start_map_coord = datetime.now()
dose_1_new.interp_sum_new(dose_2_new)
end_map_coord = datetime.now()
diff = dose_1.dose_grid - dose_1_new.dose_grid
print(np.max(diff))
rgi_time = (end_reg_grid_interp - start_reg_grid_interp).total_seconds()
mc_time = (end_map_coord - start_map_coord).total_seconds()
print('RegularGridInterp time (s):', rgi_time)
print('map_coordinates time (s):', mc_time)
print('map_coordinates is %0.1f%% times faster' % (100 * (rgi_time - mc_time) / mc_time))
7.105427357601002e-15RegularGridInterp time (s): 4.324321map_coordinates time (s): 2.977713map_coordinates is 45.2% times faster
0.0
RegularGridInterp time (s): 4.31992
map_coordinates time (s): 1.175025
map_coordinates is 267.6% times faster
On Mar 5, 2020, at 3:09 AM, Steve Gregory <stephen....@gmail.com> wrote:
--
You received this message because you are subscribed to the Google Groups "pydicom" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pydicom+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pydicom/91359b2e-c891-472c-ae94-383baf845cf5%40googlegroups.com.
self.ds.BitsAllocated = 32
self.ds.BitsStored = 32
self.ds.HighBit = 31
self.ds.DoseGridScaling = np.iinfo(np.uint32).max / np.max(self.dose_grid)
pixel_data = np.swapaxes(self.dose_grid, 0, 2) / self.ds.DoseGridScaling
self.ds.PixelData = np.uint32(pixel_data).tostring()
self.ds.DoseGridScaling = np.max(self.dose_grid) / np.iinfo(np.uint32).max
Such that the max pixel * scaling = max dose. Apologies!