Calculate summed DVH of Two Treatment Plans using Plansum
238 views
Skip to first unread message
Laura Stenzel
Jun 22, 2019, 6:08:41 PM6/22/19
to dicompyler
Hello,
I am a Medical Physics MSc student currently undertaking a research project where I would like to calculate a DVH which combines the dose received by a region of interest from a previous treatment plan and current treatment plan i.e. recurrent cancer cases.
This means I would have to calculate the dose received from both the old and new treatment plans and sum them together in order to create the summed DVH. The dose grid co-ordinate system parameters, ImagePositionPatient, GridFrameOffsetVector, Pixel_Array of the old and new plans vary and the Plansum plugin looks like it considers these discrepancies.
From my understanding, Plansum creates a new dose grid based on the co-ordinate system parameters of both the old and new treatment plans and using this new grid system, a tri-linear interpolation method is used to calculate the summed dose while correcting the discrepancies between the old and new treatment plan dose grids. Please correct me if I am mistaken.
My question is, how does Plansum relate structure information to the summed dose? To me, Plansum looks like it calculates the summed dose of the whole dose grid whereas I would be interested in calculating the summed dose for a region of interest such as the lung or heart.
I would also like to apply The EQD2 parameter in my calculations of the summed DVH in order to investigate radiobiological effects.
Any information or advice on how to tackle this research would be greatly appreciated.
Kind regards,
Laura
Aditya Panchal
Jun 24, 2019, 11:59:07 PM6/24/19
to dicom...@googlegroups.com
Hi Laura,
I will reply to your points by section.
This means I would have to calculate the dose received from both the old and new treatment plans and sum them together in order to create the summed DVH. The dose grid co-ordinate system parameters, ImagePositionPatient, GridFrameOffsetVector, Pixel_Array of the old and new plans vary and the Plansum plugin looks like it considers these discrepancies.
You are correct, it does account for these discrepancies. However the plugin assumes these plans were calculated in the same patient space (essentially same CT scan).
From my understanding, Plansum creates a new dose grid based on the co-ordinate system parameters of both the old and new treatment plans and using this new grid system, a tri-linear interpolation method is used to calculate the summed dose while correcting the discrepancies between the old and new treatment plan dose grids. Please correct me if I am mistaken.
Correct again.
My question is, how does Plansum relate structure information to the summed dose? To me, Plansum looks like it calculates the summed dose of the whole dose grid whereas I would be interested in calculating the summed dose for a region of interest such as the lung or heart.
The plansum plugin actually won't do anything with the structure information. The main thing to understand with DICOM RT is that the Structure Set and the RT Dose have to be on the same frame of reference. As stated above, as long as the treatment plans were calculated on the same CT (essentially same frame of reference) then any structure set from these treatment plans can be used with the summed RT dose as well.
I would also like to apply The EQD2 parameter in my calculations of the summed DVH in order to investigate radiobiological effects.
Essentially, you would just use plansum to sum your dose grids and either use dicompyler or python code to recalculate the DVH. You could then apply EQD2 to the dose grid prior to calculating the DVH.
Hope that is a good starting point.
Adit
--
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/d67d7e83-ab76-410c-8fcf-2fcfff4be6c5%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Laura Stenzel
Jun 25, 2019, 3:43:58 PM6/25/19
to dicompyler
Hi Aditya,
Thank you very much for your reply. I will investigate your advice during my research.
I have a follow up question with regards to your statement:
"As stated above, as long as the treatment plans were calculated on the same CT (essentially same frame of reference) then any structure set from these treatment plans can be used with the summed RT dose as well."
When both treatment plans are calculated on the same CT/patient space, will the parameters, ImagePositionPatient, GridFrameOffsetVector, Pixel_Array and Pixel Spacing be equal for dose files associated with both treatments?
In Plansum, I am aware that there is an if statement for when these parameters are not equal in both plans. I am curious when this would be the case when both plans are calculated on the same CT set? It was my understanding that these parameters would be equal when both plans are calculated in the same CT set.
On a sidenote, thank you for your work done on Dicompyler, it is a very impressive piece of software.
Kind regards,
Laura
On Tuesday, June 25, 2019 at 4:59:07 AM UTC+1, Aditya Panchal wrote:
Hi Laura,I will reply to your points by section.This means I would have to calculate the dose received from both the old and new treatment plans and sum them together in order to create the summed DVH. The dose grid co-ordinate system parameters, ImagePositionPatient, GridFrameOffsetVector, Pixel_Array of the old and new plans vary and the Plansum plugin looks like it considers these discrepancies.You are correct, it does account for these discrepancies. However the plugin assumes these plans were calculated in the same patient space (essentially same CT scan).From my understanding, Plansum creates a new dose grid based on the co-ordinate system parameters of both the old and new treatment plans and using this new grid system, a tri-linear interpolation method is used to calculate the summed dose while correcting the discrepancies between the old and new treatment plan dose grids. Please correct me if I am mistaken.Correct again.My question is, how does Plansum relate structure information to the summed dose? To me, Plansum looks like it calculates the summed dose of the whole dose grid whereas I would be interested in calculating the summed dose for a region of interest such as the lung or heart.The plansum plugin actually won't do anything with the structure information. The main thing to understand with DICOM RT is that the Structure Set and the RT Dose have to be on the same frame of reference. As stated above, as long as the treatment plans were calculated on the same CT (essentially same frame of reference) then any structure set from these treatment plans can be used with the summed RT dose as well.I would also like to apply The EQD2 parameter in my calculations of the summed DVH in order to investigate radiobiological effects.Essentially, you would just use plansum to sum your dose grids and either use dicompyler or python code to recalculate the DVH. You could then apply EQD2 to the dose grid prior to calculating the DVH.Hope that is a good starting point.Adit
On Sat, Jun 22, 2019 at 5:08 PM Laura Stenzel <laura...@gmail.com> wrote:
Hello,--I am a Medical Physics MSc student currently undertaking a research project where I would like to calculate a DVH which combines the dose received by a region of interest from a previous treatment plan and current treatment plan i.e. recurrent cancer cases.This means I would have to calculate the dose received from both the old and new treatment plans and sum them together in order to create the summed DVH. The dose grid co-ordinate system parameters, ImagePositionPatient, GridFrameOffsetVector, Pixel_Array of the old and new plans vary and the Plansum plugin looks like it considers these discrepancies.From my understanding, Plansum creates a new dose grid based on the co-ordinate system parameters of both the old and new treatment plans and using this new grid system, a tri-linear interpolation method is used to calculate the summed dose while correcting the discrepancies between the old and new treatment plan dose grids. Please correct me if I am mistaken.My question is, how does Plansum relate structure information to the summed dose? To me, Plansum looks like it calculates the summed dose of the whole dose grid whereas I would be interested in calculating the summed dose for a region of interest such as the lung or heart.I would also like to apply The EQD2 parameter in my calculations of the summed DVH in order to investigate radiobiological effects.Any information or advice on how to tackle this research would be greatly appreciated.Kind regards,Laura
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicom...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicom...@googlegroups.com.
Stephen Terry
Jun 25, 2019, 3:55:16 PM6/25/19
to dicom...@googlegroups.com
Laura,
Sorry for not getting back to you sooner. I've been on vacation.
The thing to remember is that the resolution of the dose grid does not necessarily match up with the resolution of the CT data set. So, for example, you could have a CT series with slices every 1 mm, but calculate a dose grid with points every 10 mm. So two different data sets can have different dicom parameters in regards to spatial positioning. Even if spacing is the same, the range of the dose grid could be different. In your TPS you can set the volume in which to calculate dose, and those would be different for an initial prostate plan with nodes than for a boost plan treating only the prostate.
I'll have to go back through the code to be sure, but I think the plansum plugin takes the overlap between the two dose regions and calculates the sum only in that overlap region. That prevents having to extrapolate dose outside the calculation region.
Steve
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/6f2c22d8-4f39-4712-a69d-2fc3b7f21cc5%40googlegroups.com.
Dan Cutright
Jun 25, 2019, 8:16:19 PM6/25/19
to dicompyler
Hi Steve,
Error: The length of the pixel data in the dataset doesn't match the expected amount (23674208 vs. 23878296 bytes). The dataset may be corrupted or there may be an issue with the pixel data handler.
I verified that the dimensions of old, new (initially), and new (final) objects are the same (e.g., new.pixel_array.shape).
While troubleshooting, I got a two-site plan from Eclipse (on the same CT and Structure set) to have the same dose grid position and dimensions. But for some reason the ImagePositionPatient was still off (by ~0.2mm). So I forced SumPlan to do a direct sum anyway (by editing the if condition), which worked as expected, and I verified the DVH against the composite in Eclipse. However, these two dose grids generated the error above if summed with interpolation.
Any suggestions?
Laura,
If I understand your project correctly, I think you also need a method to align the dose grids based on anatomical information prior to a dose sum (e.g., a static or deformable CT fusion)?
Thanks everyone,
Dan
Dan Cutright
Jun 25, 2019, 8:55:20 PM6/25/19
to dicompyler
Here's some more context:
the function sum_two_dose_grids is basically SumPlan without the threading or callback function. I can supply that code if you like.
import pydicom
from tools.dicom_dose_sum import sum_two_dose_grids
dose1 = pydicom.read_file('dose1.dcm')
dose2 = pydicom.read_file('dose2.dcm')
dose3 = sum_two_dose_grids(dose1, dose2)
from dicompylercore.dvhcalc import get_dvh
a = get_dvh('structure2.dcm', dose3, 8)
Traceback (most recent call last):
File "<input>", line 1, in <module>
File "/usr/local/lib/python3.7/site-packages/dicompylercore/dvhcalc.py", line 75, in get_dvh
callback)
File "/usr/local/lib/python3.7/site-packages/dicompylercore/dvhcalc.py", line 125, in calculate_dvh
dd = dose.GetDoseData()
File "/usr/local/lib/python3.7/site-packages/dicompylercore/dicomparser.py", line 768, in GetDoseData
data['dosemax'] = float(self.ds.pixel_array.max())
File "/usr/local/lib/python3.7/site-packages/pydicom/dataset.py", line 949, in pixel_array
self.convert_pixel_data()
File "/usr/local/lib/python3.7/site-packages/pydicom/dataset.py", line 895, in convert_pixel_data
raise last_exception
File "/usr/local/lib/python3.7/site-packages/pydicom/dataset.py", line 863, in convert_pixel_data
arr = handler.get_pixeldata(self)
File "/usr/local/lib/python3.7/site-packages/pydicom/pixel_data_handlers/numpy_handler.py", line 307, in get_pixeldata
.format(actual_length, expected_len + expected_len % 2)
ValueError: The length of the pixel data in the dataset doesn't match the expected amount (9702464 vs. 9766296 bytes). The dataset may be corrupted or there may be an issue with the pixel data handler.
Laura Stenzel
Jun 26, 2019, 4:11:19 AM6/26/19
to dicompyler
Hi Steve,
Thank you for your reply.
I was not aware of the phenomena that the resolution of dose grids of both plans may not match up with the resolution of the CT set even when both plans are calculated on the same CT set (if my understanding of your reply is correct).
I will have to investigate further if it is possible to implement Plansum during my project.
Thank you again for your insight.
Kind regards,
Laura
Laura Stenzel
Jun 26, 2019, 4:14:31 AM6/26/19
to dicompyler
Hi Dan,
Thank you for your advice.
A deformable CT fusion software known as Mirada will be installed into my department soon so I hope to investigate the use of this during my project.
I'm unable to help you with your query as I am still at a learning stage with Python but I hope you find a solution!
Thank you again!
Kind regards,
Laura
Dan Cutright
Jul 2, 2019, 1:26:45 PM7/2/19
to dicompyler
Hi Steve, I was just curious if you've had a chance to take a look? I was hoping to use this code in my database software (https://github.com/cutright/DVH-Analytics).
Stephen Terry
Jul 5, 2019, 9:34:12 AM7/5/19
to dicom...@googlegroups.com
I haven't, sorry. It's probably something simple regarding setting the size of the new array in pydicom. Next week should settle down a little bit, and I'll take a look. If I can't reproduce it, would it be possible to get an anonymized copy of your data? Thanks.
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/20ba9d07-71d8-4d2d-86bd-0975836a0615%40googlegroups.com.
Dan Cutright
Jul 5, 2019, 1:35:03 PM7/5/19
to dicompyler
Sure thing, just let me know if you need it.
On Friday, July 5, 2019 at 8:34:12 AM UTC-5, Stephen Terry wrote:
I haven't, sorry. It's probably something simple regarding setting the size of the new array in pydicom. Next week should settle down a little bit, and I'll take a look. If I can't reproduce it, would it be possible to get an anonymized copy of your data? Thanks.
Dan Cutright
Dec 8, 2019, 1:33:26 PM12/8/19
to dicompyler
Hi Steve, I don't suppose you were able to take a look?
Stephen Terry
Dec 10, 2019, 10:08:52 AM12/10/19
to dicom...@googlegroups.com
Dan,
I took a look back through the code, and I can't seem to remember where I left things. Would you mind sending me an anonymized version of your data so I can take a look? Sorry for dropping the ball earlier.
Steve
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/62375397-b29d-489b-880a-513383ebb401%40googlegroups.com.
Dan Cutright
Dec 10, 2019, 3:23:03 PM12/10/19
to dicompyler
No worries at all, thanks Steve... I'll send you a private message with a link when it's ready.
I just tried to reproduce the issue, but now I'm seeing a different error. I'll poke around a little before sending you on a wild goose chase.
Dan Cutright
Jan 14, 2020, 4:43:49 PM1/14/20
to dicompyler
Hi Steve,
Sorry for the delay...
It looks like the dose sum works when providing only two DICOM dose objects
Are you able to add the output of your dose sum to a third DICOM dose? I get the following error:
Traceback (most recent call last):
File "<input>", line 1, in <module>
File "C:\Users\dcutright\PycharmProjects\DVH-Analytics\dvha\tools\dicom_dose_sum.py", line 108, in sum_two_dose_grids
dose_sum = interpolate_image(np.swapaxes(old.pixel_array, 0, 2), scale_old, old.ImagePositionPatient, sum_xyz_coords) * old.DoseGridScaling + \
File "C:\Users\dcutright\AppData\Local\Programs\Python\Python37-32\lib\site-packages\pydicom\dataset.py", line 949, in pixel_array
self.convert_pixel_data()
File "C:\Users\dcutright\AppData\Local\Programs\Python\Python37-32\lib\site-packages\pydicom\dataset.py", line 895, in convert_pixel_data
raise last_exception
File "C:\Users\dcutright\AppData\Local\Programs\Python\Python37-32\lib\site-packages\pydicom\dataset.py", line 863, in convert_pixel_data
arr = handler.get_pixeldata(self)
File "C:\Users\dcutright\AppData\Local\Programs\Python\Python37-32\lib\site-packages\pydicom\pixel_data_handlers\numpy_handler.py", line 307, in get_pixeldata
.format(actual_length, expected_len + expected_len % 2)
ValueError: The length of the pixel data in the dataset doesn't match the expected amount (29098640 vs. 29303560 bytes). The dataset may be corrupted or there may be an issue with the pixel data handler.
On Tuesday, December 10, 2019 at 9:08:52 AM UTC-6, Stephen Terry wrote:
Dan,I took a look back through the code, and I can't seem to remember where I left things. Would you mind sending me an anonymized version of your data so I can take a look? Sorry for dropping the ball earlier.Steve
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/62375397-b29d-489b-880a-513383ebb401%40googlegroups.com.
Dan Cutright
Jan 14, 2020, 4:48:38 PM1/14/20
to dicompyler
False alarm, it looks like a regular sum causes an issue. I'll send you a PM with the two dose grids anonymized.
Dan Cutright
Feb 10, 2020, 12:45:31 PM2/10/20
to dicompyler
Problem solved!
The interpolate_image function was returning a different sized array than was assumed for the new Rows and Columns attributes of the pydicom dataset.
Currently at https://github.com/dicompyler/dicompyler-plugins/blob/master/plugins/plansum/plansum.py
Lines 186 to 210:
x_vals = np.arange(x0,x1,scale_sum[0])
y_vals = np.arange(y0,y1,scale_sum[1])
z_vals = np.arange(z0,z1,scale_sum[2])
#Create a 3 x i x j x k array of xyz coordinates for the interpolation.
sum_xyz_coords = np.array([i*scale_sum[0] + sum_ip[0],
j*scale_sum[1] + sum_ip[1],
k*scale_sum[2] + sum_ip[2]])
#Dicom pixel_array objects seem to have the z axis in the first index
#(zyx). The x and z axes are swapped before interpolation to coincide
#with the xyz ordering of ImagePositionPatient
sum = interpolate_image(np.swapaxes(old.pixel_array,0,2), scale_old,
old.ImagePositionPatient, sum_xyz_coords,
progressFunc)*old.DoseGridScaling + \
interpolate_image(np.swapaxes(new.pixel_array,0,2), scale_new,
new.ImagePositionPatient, sum_xyz_coords,
progressFunc)*new.DoseGridScaling
#Swap the x and z axes back
sum = np.swapaxes(sum, 0, 2)
sum_dcm.ImagePositionPatient = list(sum_ip)
sum_dcm.Rows = len(y_vals)
sum_dcm.Columns = len(x_vals)
sum_dcm.NumberOfFrames = len(z_vals)
If I change Lines 208 to 210 to the following, issues resolved:
sum_dcm.Rows = sum.shape[2]
sum_dcm.Columns = sum.shape[1]
sum_dcm.NumberOfFrames = sum.shape[0]
Stephen Terry
Feb 10, 2020, 1:41:29 PM2/10/20
to dicom...@googlegroups.com
Great! Thanks for tracking that down. I've been fighting with Anaconda and Pydev, trying to get them to work together. When I get them working, I'll verify this on my end then upload the fix.
--
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/95ead741-e656-4c45-b07f-26daca3e02f2%40googlegroups.com.
Dan Cutright
Feb 10, 2020, 6:22:22 PM2/10/20
to dicompyler
Sure thing. This plugin is very helpful... no idea why Eclipse can't export a composite dose grid!
Dan Cutright
Feb 13, 2020, 5:25:26 PM2/13/20
to dicompyler
I came across another indexing error, this time in the trilinear_interp function. This only showed up when I was adding two very large dose grids (total marrow irradiation plans), and each grid had a different resolution (a 2mm and a 3mm). If this doesn't make sense, I can anonymize the DICOM dose files that do this for me.
You had capped the x1, y1, and z1 values, but I was getting values too large in z0.
The other changes you see in the commit link below are related to optimizing memory allocation.
Dan Cutright
Feb 27, 2020, 11:43:59 AM2/27/20
to dicompyler
Hi Stephen,
As an aside, I had tried writing my own summation using scipy's RegularGridInterpolator, and it was insanely slow. Took about 10 minutes for one summation!
Below is the code I'm currently using in DVH Analytics. Some of the arithmetic looks a little funny, but it prevents MemoryErrors. Also, I wonder if it would be more accurate (and faster) to interpolate one dose grid into the same space as the other, so you only interpolate one grid?
# Slightly modified from https://github.com/dicompyler/dicompyler-plugins/blob/master/plugins/plansum/plansum.py
def sum_two_dose_grids(old, new):
""" Given two Dicom RTDose objects, returns a summed RTDose object"""
"""The summed RTDose object will consist of pixels inside the region of
overlap between the two pixel_arrays. The pixel spacing will be the
coarser of the two objects in each direction. The new DoseGridScaling
tag will be the sum of the tags of the two objects.
interp_method: A string that is one of ['scipy','weave','python'].
This forces SumPlan to use a particular interpolation method, even if
the dose objects could be directly summed. Used for unit testing."""
# Recycle the new Dicom object to store the summed dose values
sum_dcm = new
# Test if dose grids are coincident. If so, we can directly sum the
# pixel arrays.
if (old.ImagePositionPatient == new.ImagePositionPatient and
old.pixel_array.shape == new.pixel_array.shape and
old.PixelSpacing == new.PixelSpacing and
old.GridFrameOffsetVector == new.GridFrameOffsetVector):
print("PlanSum: Using direct summation")
dose_sum = old.pixel_array * old.DoseGridScaling + new.pixel_array * new.DoseGridScaling
else:
print("PlanSum: Using interpolation")
# Compute mapping from xyz (physical) space to ijk (index) space
scale_old = np.array([old.PixelSpacing[0], old.PixelSpacing[1],
old.GridFrameOffsetVector[1] - old.GridFrameOffsetVector[0]])
scale_new = np.array([new.PixelSpacing[0], new.PixelSpacing[1],
new.GridFrameOffsetVector[1] - new.GridFrameOffsetVector[0]])
scale_sum = np.maximum(scale_old, scale_new)
# Find region of overlap
xmin = np.array([old.ImagePositionPatient[0], new.ImagePositionPatient[0]])
ymin = np.array([old.ImagePositionPatient[1], new.ImagePositionPatient[1]])
zmin = np.array([old.ImagePositionPatient[2], new.ImagePositionPatient[2]])
xmax = np.array([old.ImagePositionPatient[0] + old.PixelSpacing[0] * old.Columns,
new.ImagePositionPatient[0] + new.PixelSpacing[0] * new.Columns])
ymax = np.array([old.ImagePositionPatient[1] + old.PixelSpacing[1] * old.Rows,
new.ImagePositionPatient[1] + new.PixelSpacing[1] * new.Rows])
zmax = np.array([old.ImagePositionPatient[2] + scale_old[2] * len(old.GridFrameOffsetVector),
new.ImagePositionPatient[2] + scale_new[2] * len(new.GridFrameOffsetVector)])
x0 = xmin[np.argmin(abs(xmin))]
x1 = xmax[np.argmin(abs(xmax))]
y0 = ymin[np.argmin(abs(ymin))]
y1 = ymax[np.argmin(abs(ymax))]
z0 = zmin[np.argmin(abs(zmin))]
z1 = zmax[np.argmin(abs(zmax))]
sum_ip = np.array([x0, y0, z0])
# Create index grid for the sum array
i, j, k = np.mgrid[0:int((x1 - x0) / scale_sum[0]),
0:int((y1 - y0) / scale_sum[1]),
0:int((z1 - z0) / scale_sum[2])]
# x_vals = np.arange(x0, x1, scale_sum[0])
# y_vals = np.arange(y0, y1, scale_sum[1])
z_vals = np.arange(z0, z1, scale_sum[2])
# Create a 3 x i x j x k array of xyz coordinates for the interpolation.
sum_xyz_coords = np.array([i * scale_sum[0] + sum_ip[0],
j * scale_sum[1] + sum_ip[1],
k * scale_sum[2] + sum_ip[2]])
# Dicom pixel_array objects seem to have the z axis in the first index
# (zyx). The x and z axes are swapped before interpolation to coincide
# with the xyz ordering of ImagePositionPatient
dose_sum = interpolate_image(np.swapaxes(old.pixel_array, 0, 2), scale_old,
old.ImagePositionPatient, sum_xyz_coords) * old.DoseGridScaling
dose_interp_2 = interpolate_image(np.swapaxes(new.pixel_array, 0, 2), scale_new,
new.ImagePositionPatient, sum_xyz_coords) * new.DoseGridScaling
dose_sum += dose_interp_2
# Swap the x and z axes back
dose_sum = np.swapaxes(dose_sum, 0, 2)
sum_dcm.ImagePositionPatient = list(sum_ip)
sum_dcm.Rows = dose_sum.shape[2]
sum_dcm.Columns = dose_sum.shape[1]
sum_dcm.NumberOfFrames = dose_sum.shape[0]
sum_dcm.PixelSpacing = [scale_sum[0], scale_sum[1]]
sum_dcm.GridFrameOffsetVector = list(z_vals - sum_ip[2])
sum_scaling = old.DoseGridScaling + new.DoseGridScaling
dose_sum = dose_sum / sum_scaling
dose_sum = np.uint32(dose_sum)
# sum_dcm.pixel_array = sum
sum_dcm.BitsAllocated = 32
sum_dcm.BitsStored = 32
sum_dcm.HighBit = 31
sum_dcm.PixelData = dose_sum.tostring()
sum_dcm.DoseGridScaling = sum_scaling
return sum_dcm
def interpolate_image(input_array, scale, offset, xyz_coords):
"""Interpolates an array at the xyz coordinates given"""
"""Parameters:
input_array: a 3D numpy array
scale: a list of 3 floats which give the pixel spacing in xyz
offset: the xyz coordinates of the origin of the input_array
xyz_coordinates: the coordinates at which the input is evaluated.
The purpose of this function is to convert the xyz coordinates to
index space, which trilinear_interp
requires. Following the scipy convention, the xyz_coordinates array is
an array of three i x j x k element arrays. The first element contains
the x axis value for each of the i x j x k elements, the second contains
the y axis value, etc."""
indices = np.empty(xyz_coords.shape)
indices[0] = (xyz_coords[0] - offset[0]) / scale[0]
indices[1] = (xyz_coords[1] - offset[1]) / scale[1]
indices[2] = (xyz_coords[2] - offset[2]) / scale[2]
return trilinear_interp(input_array, indices)
def trilinear_interp(input_array, indices):
"""Evaluate the input_array data at the indices given"""
# output = np.empty(indices[0].shape)
x_indices = indices[0]
y_indices = indices[1]
z_indices = indices[2]
x0 = x_indices.astype(np.integer)
y0 = y_indices.astype(np.integer)
z0 = z_indices.astype(np.integer)
# Check if xyz0 is beyond array boundary:
x0[np.where(x0 >= input_array.shape[0])] = np.shape(input_array)[0] - 1
y0[np.where(y0 >= input_array.shape[1])] = np.shape(input_array)[1] - 1
z0[np.where(z0 >= input_array.shape[2])] = np.shape(input_array)[2] - 1
x1 = x0 + 1
y1 = y0 + 1
z1 = z0 + 1
# Check if xyz1 is beyond array boundary:
x1[np.where(x1 >= input_array.shape[0])] = x0.max()
y1[np.where(y1 >= input_array.shape[1])] = y0.max()
z1[np.where(z1 >= input_array.shape[2])] = z0.max()
x = x_indices - x0
y = y_indices - y0
z = z_indices - z0
output = input_array[x0, y0, z0] * tri_linear_factor(x, y, z, 0, 0, 0)
for zi, zz in enumerate([z0, z1]):
for yi, yy in enumerate([y0, y1]):
for xi, xx in enumerate([x0, x1]):
if xi + yi + zi: # ignore 0 0 0 since it was created before loops
output += input_array[xx, yy, zz] * tri_linear_factor(x, y, z, xi, yi, zi)
return output
def tri_linear_factor(x, y, z, xi, yi, zi):
ans = x if xi else (1 - x)
ans = ans * y if yi else ans * (1 - y)
ans = ans * z if zi else ans * (1 - z)
return ans
Dan Cutright
Feb 27, 2020, 7:15:21 PM2/27/20
to dicompyler
The code in the previous post gives some nonsensical results, including very much incorrect structure volumes... so I'm leaning towards the overlapping grid is calculated incorrectly. I'm guessing my edits to avoid indexing and memory errors have caused the issue, but I can't test that unfortunately.
Below is an interpolator I wrote which is generating correct DVHs, but it's insanely slow. I'll look into using something like dolo for interpolation instead of scipy... unless there are other suggestions out there?
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
Mar 3, 2020, 5:52:31 PM3/3/20
to dicompyler
The good news is I have an interpolator that works for me. I've posted the code on the pydicom user group here: https://groups.google.com/d/msg/pydicom/JRViIu5PbFM/0k72RQR3AgAJ
The bad news is that it uses scipy which I don't believe is a part of the dicompyler GUI? But it allows me to calculate the interpolations in chunks to avoid memory allocation errors. I suspect Stephen's method could be adapted for that, but I don't understand the method well enough to suggest how.
Thanks to Stephen for sharing his code, it helped me out immensely.
Stephen Terry
Mar 3, 2020, 6:05:36 PM3/3/20
to dicom...@googlegroups.com
Sure. I apologize for not being very helpful. I'm swamped at the moment and trying to work with Pydev is killing me at the moment.
--
-- You received this message because you are subscribed to the Google Groups dicompyler group. To post to this group, send email to dicom...@googlegroups.com. To unsubscribe from this group, send email to dicompyler+...@googlegroups.com. For more options, visit this group at https://groups.google.com/d/forum/dicompyler?hl=en
---
You received this message because you are subscribed to the Google Groups "dicompyler" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dicompyler+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dicompyler/781833b2-209e-4fe8-bd94-7ab3f67cdcf4%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages