Dose Grid Summation

219 views
Skip to first unread message

Dan Cutright

unread,
Mar 2, 2020, 6:37:34 PM3/2/20
to pydicom
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/qkU2CtYzgLg
I 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.


#!/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))

Darcy Mason

unread,
Mar 2, 2020, 9:37:26 PM3/2/20
to pydicom
I also have an interest in summing dose grids so had a look at this.

I'm by no means completely understanding it, but one thing that worries me, is the xyz_to_ijk(selfxyz) modifies the passed xyz in place as well as returning the altered values, i.e. modifies `self` when used to map the `other` grid.  Trying to follow possible side-effects of that, however, I don't see that there are any.

Second thought is: have you tried a minimal example with some very simple arrays not involving dose grid scaling and so on?  Then add on one extra complexity at a time?


On Monday, March 2, 2020 at 6:37:34 PM UTC-5, Dan Cutright wrote:
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/qkU2CtYzgLg
I 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>

Dan Cutright

unread,
Mar 2, 2020, 10:12:44 PM3/2/20
to pydicom
Hi Darcy, thanks for your response and interest.

I also don't think modifying the xyz_to_ijk is an issue, but I could be wrong... and to your point, I guess I don't need to explicitly return anything.  My test case involves a near whole body scan with a 2mm uniform dose grid (total marrow irradiation case), so I run into memory allocation issues on Windows if I'm not careful (no issues on macOS though).

Below is a simple example working with meshgrid and vstack to get a list points points cycling through x first, then y, then z.  For some reason, I feel like the answer has something do with the points design.

And finally, at the end is a working dose summation that I'm pretty sure is calculated correctly (only validated on one case), but the code isn't as organized and it is super slow.  I think it'd be orders of magnitude faster if I generate points in same order as the dose grid array so I don't have to look up the ijk coordinate for every point and also be able to add all the dose values in one shot.


>>> 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):

print("PlanSum: Using direct summation")
dose_1.direct_sum(dose_2)

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

Dan Cutright

unread,
Mar 3, 2020, 11:05:08 AM3/3/20
to pydicom
So I've had some better luck with scipy.interpolate.RegularGridInterpolator and I think the code posted in the previous message is slow not really due to the interpolation, but due to a combination of the point generation and adding one point at a time (and not in order).... and possibly the frequency of checking for being inside or outside of the grid.  I had done that to avoid memory errors, but I'll have to work on a more efficient way to avoid them.

Unlike map_coordinates, RegularGridInterpolator can directly accept xyz input (which makes it a little easier to understand).  I use bounds_error=False and fill_value=0 to assume 0 dose rather than extrapolating.  The output of calling the interpolator appears to be in the correct order of my dose grid (note that my dose grid = np.swapaxes(ds.pixel_array * ds.DoseGridScaling, 0, 2)

One thing I'm unsure of, what is the consequence of not updating ds.DoseGridScaling?  It looks to me this value is either 1/prescription of 1/max_dose? I'm guessing not updating would only result in rounding errors if the max_dose of the dose_sum is significantly higher than in the original dose grid?

Here is a summary of the code, see attached file for complete code.  Any takers out there for an independent validation?


The actual interpolation and summation:
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

where self and other refer to my DoseGrid class, other values defined below:

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)])

Finally:
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())


dicom_dose_sum.py

Dan Cutright

unread,
Mar 3, 2020, 2:29:29 PM3/3/20
to pydicom
I've figured out that:
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])

is equivalent to:
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()


input:
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())


output:
List method (s): 0.918
Numpy method (s): 0.021
NumPy method was 43.715 times faster
True


Still tracking down an error when I update in chunks instead of calculating all points at once (code not posted here).

Dan Cutright

unread,
Mar 3, 2020, 5:19:22 PM3/3/20
to pydicom
Success!!

I had some red herrings throwing me off with my database software.  But I've validated this code by importing the calculated dose grid into Eclipse to compare with its composite directly.  The results are suspiciously good, so hopefully someone else will be interested in trying this for themselves.  A summation two dose grids (228x112x380 and 228x111x380) took about 10sec on my fairly old MSW7 work computer (8GB RAM, 3.2GHz i5).

The attached file can be used like so:

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)


One could easily modify this code to dose differences too, but I personally have no interest in that... and I didn't feel like dealing trying to store negative doses in DICOM (is that possible?)

I also coded in a MemoryError catch that will do the interpolation in chunks of 50,000 points at a time (the user can edit this).

I suspect scipy's map_coordinates would work just as well too if this is adapted to deal with ijk space again... plus map_coordinates supports b-splines.

And one last note, any insight about the consequences of reusing the same DoseGridScaling?  Is there a convention/standard for how it's calculated?
dicom_dose_sum.py

Steve Gregory

unread,
Mar 4, 2020, 7:57:43 AM3/4/20
to pydicom
Well done! I haven't been through all this, but I have hit an issue with DoseGridScaling before. I think it gets used to keep the stored values within 8-bit or 16-bit etc. but using all the levels to maximise dose resolution. I've seen divide by the max dose, so that it becomes the highest value stored. If the max dose is e.g. 100Gy and using 16-bit then max integer is ~65000,  so 0.2 cGy resolution. But for 8-bit then 40cGy resolution. If you don't recalculate it, and add 2 plans together, then I think the integers might overflow (if you have e.g. 200% in a voxel), unless numpy or pydicom does something internally? Hope that helps.

Steve G, Leeds UK

Dan Cutright

unread,
Mar 4, 2020, 10:36:56 AM3/4/20
to pydicom
Thanks Steve, I think I'll edit my code to set the DoseGridScaling to 1/max_dose.

Not sure if it's appropriate or not, but I'm storing a 32-bit pixel_array with np.uint32().  Eclipse 11, at least, has no issues reading that in.

Dan Cutright

unread,
Mar 4, 2020, 11:28:16 PM3/4/20
to pydicom
One last follow-up. It turns out scipy.nd_image.map_coordinates is significantly faster than scipy.interpolate.RegularGridInterpolator.  For anyone interested in the entire code and any updates, it can be found on the DVH Analytics GitHub page (dvha/tools/dicom_dose_sum.py).

map_coordinates behaves in a similar manner, but it works in ijk space.  The bonus is you can use splines if you specify the order, but it's not clear it much matters... there's way more error/variation in DVH calculations than in dose grid interpolations.

Based on some googling, I think there's potential for another order of magnitude in calculation time reduction, and still on the CPU. But I'm totally happy with 3 sec compared to the 5 min calc time I started with last week!

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()


My test code run on a 2019 MacBook Pro 2.3GHz i9:
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))

Output:
7.105427357601002e-15
RegularGridInterp time (s): 4.324321
map_coordinates time (s): 2.977713
map_coordinates is 45.2% times faster

Dan Cutright

unread,
Mar 4, 2020, 11:59:25 PM3/4/20
to pydicom
I lied...  the above code uses a b-spline (default order=3) for map_coordinates.  To compare apples to apples, I reran with a trilinear interp (order=1) for map_coordinates, the output was:

0.0
RegularGridInterp time (s): 4.31992
map_coordinates time
(s): 1.175025
map_coordinates
is 267.6% times faster

Steve Gregory

unread,
Mar 5, 2020, 4:09:34 AM3/5/20
to pydicom
Just to clarify, DoseGridScaling = (2**32 -1) / max dose
2**32 -1 is the max integer value, also obtained from np.iinfo(np.uint32).max
Steve G, Leeds UK

Dan Cutright

unread,
Mar 5, 2020, 8:35:15 AM3/5/20
to pyd...@googlegroups.com
Ahh ok, that makes much more sense.

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.

Dan Cutright

unread,
Mar 5, 2020, 10:21:29 AM3/5/20
to pydicom
I'm not sure why, but that ends up giving me 0 dose when I calculate DVHs.  Perhaps something other than DoseGridScaling needs to be updated?

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()

Steve Gregory

unread,
Mar 6, 2020, 7:31:45 AM3/6/20
to pydicom
Sorry I wrote it out upside down. Scaling should be max dose / max int.

self.ds.DoseGridScaling = np.max(self.dose_grid) / np.iinfo(np.uint32).max

Such that the max pixel * scaling = max dose. Apologies!

Dan Cutright

unread,
Mar 6, 2020, 9:28:41 AM3/6/20
to pydicom
No worries, I should have noticed that too.  Math is hard :)

Made the swap, all looks good now!
Reply all
Reply to author
Forward
0 new messages