I am coding to extract the deformation vector field stored in a DICOM REG file. This file is generated and exported from ImSim QA.
The size of the grid is 512512141. The DVF value is basically (0,0,9) on each point as I generated it. These values are also shown in ImSim QA perfectly. However, the value on each point shown in Slicer is around (0,0,3). I also tried to extract the DVF in a DIY python script, the result is similar to those shown in Slicer.
Thanks for sharing the example. It would be great if you could share the script you are using to look at the data. I tried loading it in pydicom but the vector data couldn’t be loaded as numpy. If you think the data is valid this could be good to report on the pydicom issue tracker.
d = pydicom.read_file("/tmp/z 9.dcm")
g = d.DeformableRegistrationSequence[0].DeformableRegistrationGridSequence[0]
v = g.VectorGridData
import numpy
a = numpy.array(v)
print(a.min(), a.max())
Traceback (most recent call last):
File "<string>", line 7, in <module>
File "/opt/sr/python-install/lib/python3.9/site-packages/numpy/core/_methods.py", line 45, in _amin
return umr_minimum(a, axis, None, out, keepdims, initial, where)
numpy.core._exceptions._UFuncNoLoopError: ufunc 'minimum' did not contain a loop with signature matching types (dtype('S443547648'), dtype('S443547648')) -> None
However if I look with dcmdump, it appears that the vector data is not all 0,0,9 but actually some other vectors.
A DICOM dump shows that displacement values in the file are approximately all (0, 0, - 3). Therefore, to me it looks like Slicer loaded the transform correctly. We also tested this quite extensively when we developed the importer a number of years ago.
Maybe ImSimQA was confused by somehow wanting to factor in the grid spacing of 3 into the calculations? I would be surprised (and it would be very embarrassing for them), but it is possible. I would recommend to tell Imsim QA developers know about your findings. Report back here if they reply.
Until they get back to you, you can use Slicer to generate displacement fields. You can compute transforms from automatic intensity based image registration (e.g., using Elastix or ANTs) or landmarks (e.g., using Fiducial Registration Wizard module) and then export it to DICOM SRO. You can also do visualization and analysis, apply transforms to structure sets, compute statistics, DVH, DVH metrics, etc. using the Slicer GUI or automate or extend these features using Python scripts.
Thank you so much for the discussion. I think the problem might come from unpacking the data. The data format is little endian.
#get the displacement field
dvf_vals_raw = ds.DeformableRegistrationSequence[0].DeformableRegistrationGridSequence[0].VectorGridData
#unpack it in little endian format
dvf_vals = struct.unpack(f"<{len(dvf_vals_raw)//4}f", dvf_vals_raw)
This is the right explanation. As ImSim QA explained, the DVF stored in VectorGridData is in voxels. ‘‘To convert the DVF values to mm, one must multiply by the voxel size.’’
To me it seems ImSim QA developers made mistake: they think that the vectors are in voxels, while according to the DICOM standard it is in millimeters.
According to the DICOM standard, the displacement vectors are in the “Registered RCS”:
Since the vectors are added, it means X_Start and DeltaX_IJK are in the same coordinate system. The standard defines X_Start as “The standard coordinate, in the Registered RCS, of the deformation grid as specified in Image Position (Patient) (0020, 0032)”.
ImSim QA is a company working in radiation therapy quality assurance, so making such mistake would be an enormous embarrassment for them (who would ever trust their products if they make such elementary mistake?), so I’m slightly hesitant to make a final verdict.
@David_Clunie Could you confirm that Vector Grid Data (0064,0009) vectors are in Registered RCS physical coordinate system (unit is millimeter and not voxel)?
Yes Andras, the encoded displacements should be in mm not voxels, as I understand the standard.
I agree that Equation C.20.3-1 in the standard that describes how to find the vector and then apply it seems pretty clear that the “deformation specified at index (i,j,k) in the deformation grid” is added to the coordinate without there being any multiplication by the voxel size (as distinct from finding the right entry by using the deformation grid resolution).
As HAL would say, “I don’t think there is any question about it.”
For confirmation, I suggest you ask the DICOM WG 7 (RT) and IHE-RO mailing lists, who may want to add a CP to emphasize this in the standard. You might also want to contact the vendor to see if they need to initiate a safety recall.
@David_Clunie thank you very much for the quick answer. I also saw that you already contacted WG7 members, thanks for this amazing responsiveness!
@bqi Please report this to your ImSim QA contacts. Instead of telling their customers to “multiply by voxel size”, they must fix their software and may want to initiate a safety recall.
Not sure if you still need to be on the WG-07 mailing list register to post since they changed the server, but at the very least you need to be signed up for the NEMA discourse site (https://forums.nema.org/).
So I sent a message to them about this just now, since they are meeting F2F in Stockholm this week and can discuss it. I will relay any useful response.
