dose actor attached to patient CT data (ImageNestedParametrisedVolume)

[djboersma]

When I tried to configure my dose actor such that the dosel geometry would be the same as the patient voxel geometry by setting the same voxelsize, I got a 511x511x55 dose distribution while the patient is actually 512x512x55, and the "elementspacing" in the dose MHD file does not match the voxelsize which I configured, and it is also different from the patient voxelsize; I think that the only quantity that was preserved is the total size of the image volume (I think you call that the "bounding box"). In Z the spacing is a round number, namely 3mm, but in X and Y it's 1.03516mm, so I'm guessing there is some rounding issue going on. If instead of the VoxelSize I set the Resolution to 512x512x55 this problem seems to be solved (or rather: worked around).

I have browsed a bit through the dose and image actor code but it's quite a deep rabbit hole, so I can't really claim yet that I found a bug. I can look into it later (and learn a lot about ITK and imaging). But maybe also someone else can have a look at this who has more experience with imaging code and has a clearer idea where/why this might happen. I can share the details for reproducing the problem, if needed.

Another solution which I tried is to omit both voxel size and resolution. I was hoping that it would take those settings from the patient data, but instead the dose actor regarded the entire patient image as a single voxel. I thought that was odd.

I was going to clarify the DoseActor documentation on the wiki, based on these experiences, but after some more pondering I thought that this is actually something that needs fixing, or at least discussion.

So this issue has two parts:

(1) There seems to be a rounding problem with setVoxelSize in the dose actor. (2) Why doesn't the dose actor by default use the same voxelization as the patient data?

[ThomasDeschler]

Hello, This is a known issue which, I think, does not depend on dose actor. Can you try the branch dicom-dev and tell me if the problem remain ? Thomas

[djboersma]

Hi Tom, thanks for pointing me to dicom-dev, I'll try it out! Probably also very useful for many other things...

[BishopWolf]

Refer to https://github.com/OpenGATE/Gate/pull/29 for discussion. The line that does this is the 155 in the GateVImageActor.cc. So you want to change the default behaviour from

mResolution = G4ThreeVector(1.0, 1.0, 1.0);

But this requires a definition tweak, the Voxel actor can be used in ANY kind of geometry, not only voxelized ones, so what is "the same voxelization as the patient data" in a case where the are no patient at all? I think current behaviour is more general

[BishopWolf]

"If instead of the VoxelSize I set the Resolution to 512x512x55 this problem seems to be solved (or rather: worked around)."

This is related with the lines: if (mVoxelSizeIsSet) { mResolution.setX(std::max(1, (int)floor(mHalfSize.x()*2/mVoxelSize.x()))); mResolution.setY(std::max(1, (int)floor(mHalfSize.y()*2/mVoxelSize.y()))); mResolution.setZ(std::max(1, (int)floor(mHalfSize.z()*2/mVoxelSize.z()))); }

so instead of floor we might need to change these lines to use round. What are your thoughts??

[djboersma]

Yes, I saw those lines. I thought that floor was chosen with good reasons, because when someone deliberately sets a different voxelsize then you want to avoid that any part of a dosel sticks out of the volume described by the patient data. A hackish "solution" would be to inspect mHalfSize.x()*2/mVoxelSize.x() and if it is very close to an integer (NNN.99999) then you use round, otherwise floor. The most hackish part in this "solution" is choosing an arbitrary margin value to define what "very close" means.

But even if this hack would work, I am still surprised that the voxel size that I set in my macro sometimes gets overridden. I configured 1.03516, but it got overridden with 1.03719. This made me think that there must be a bug somewhere. Note that 511*1.03719=530.00409 while 512*1.03516=530.00192, which makes me think that the following happened: (1) code sees that patient size is 530mm and gets 1.03516 for voxelsize (2) since 512*1.03516 is slightly larger than 530, the number of voxels becomes 511 instead, thanks to floor. (3) I would now think that the volume described by the dose actor would simply be 511*1.03516=528.96676mm, but it looks like instead the code decides to use the same total size as for the patient volume, namely 530mm (4) in order to make the total size (530mm), the number of voxels (511) and the voxelsize match, the voxelsize now gets adjust to 530.0/511=1.03718 (I got 1.03719 in the dose actor output MHD file).

The way I understand issue #27, the user should be able to define 2 out of 3 variables (total size, number of voxels, voxelsize) and then the third gets computed from the other two. But if my guess in (3)+(4) is right, then this does not always happen; in particular the voxelsize as provided by the user gets overridden. That would be a bug.

Let me add that I did not explicitly set the total size in the dose actor: I assumed that this would be taken from the patient data to which it is attached. So in the dose actor I only specified the voxelsize. I think that therefore in line 143 of GateVImageActor.cc ConstructBoundingBox gets called to compute the total size. Then in lines 167-171 the resolution gets set, as you wrote above. Right after that a call mImage.SetResolutionAndVoxelSize(mResolution, mVoxelSize); happens, and as far as I can see (in GateVImage.cc) that does the correct thing (use the specified resolution and voxelsize, update the total size). And also in the rest of GateDoseActor.cc code I cannot really find anything that would be responsible for the behavior that I insinuated/deduced above in points (3)+(4), it just uses the geometry as configured/constructed by GateVImage. So I'm confused now.

[BishopWolf]

First of all a Definition Phase:

Actors and Input Geometries are different things, and they shall remain different things. So never expect to have an actor defaulted to input geometry properties.

With that said, your assumption is wrong (" this would be taken from the patient data to which it is attached"), again imagine a case with no patient at all. So you have to entirely define your actor properties.

BUT, If something overrides user definitions in actors then it is a bug! I currently dont see anything that triggers your behaviour if we respect the above.

[ThomasDeschler]

This type of mhd bug should be fixed in the dicom-dev branch. Tell me if it is not.

[dsarrut]

Are you sure Tom ? You think it is related to floating point representation ?

@djboersma @BishopWolf : yes the actors and input volumes are different. But the actor, if the user did not specify the size, assume the size of the bounding box of the volume it is attached to. Note that the volume could be anything (a sphere etc), so some physical space of this bounding box may be outside the real volume space ; and in case of a particle moving there, the actor will not record event (which is expected).

Code is here: https://github.com/OpenGATE/Gate/blob/develop/source/digits_hits/src/GateVImageActor.cc#L133

Let me know if things change with the dicom-dev branch ...

[ThomasDeschler]

No I'm not sure, but I have found that floating point representation can cause a false rounding of the resolution (and size) of the mhd images when they are opened by Gate.

Maybe it is not correlated with this problem but I think, on this case, dicom-dev worth a try.

[djboersma]

OK, I checked out the dicomdev branch and it looks like that for my current case, it performs how I was expecting Gate to behave: when I specify the voxelsize values taken directly from the patient MHD file, but not the total size nor the resolution, then the output from dose actor has the same voxelsize, resolution and offset as the input figure; while, as reported above, dose actor in gate v7.2 changes things a bit (running exactly same script, changing only the Gate executable and the output directory path):

$ diff output/srb_8_4_voxeltest_dicomdev/dose_pat53_srb_8_4_3-Dose.mhd output/srb_8_4_voxeltest_gate72/dose_pat53_srb_8_4_3-Dose.mhd
7c7
< Offset = -265 -148 -637
---
> Offset = -264.999 -147.999 -637
9,10c9,10
< ElementSpacing = 1.03516 1.03516 3
< DimSize = 512 512 55
---
> ElementSpacing = 1.03719 1.03719 3
> DimSize = 511 511 55

I am using ITK version 4.9.1 on an Ubuntu 14.04 for dicomdev. I explicitly enabled ITK in cmake (GATE_USE_ITK set to ON).

[ThomasDeschler]

I'm glad that dicom-dev solve this issue.

@dsarrut : Did I need to make a pull request of the current version of dicom-dev to close this issue ?

[djboersma]

I have the suspicion that the rounding errors originate from the rather antique MetaIO copy we are using (copied from some 3.x release of ITK). If that suspicion is correct, then another way to fix the bug is to update our MetaIO copy; I think there were some plans to do that at some point anyway. I think that avoiding the copy altogether and using a system install of ITK is in principle a better long term solution but it looks like the MetaIO copy was originally included in Gate to make things easier for the users, it avoids that we have to add to the "install requirements list" that they should install ITK version 4.9 or newer. This is not something you can easily do using the package manager on a vanilla linux distro from last year or older. (ITK is not very hard to compile from source though, not harder than compiling Gate itself, except that it takes a little longer because ITK has more code than Gate.)

Also: my test is extremely limited, it involves only one use case. I suppose that there must be many more verification steps (other use cases) for such a significant update.

[BishopWolf]

IMHO ITK is way too large and extremely big for the average user to worry about installing it, bundle just the piece that we need was a good idea, what is wrong is to keep it outdated. The end user shall not be forced to do more than he usually do but less, thats a principle in software development.

[brenthuisman]

Not sure why pulling in ITK would be necessary for MetaIO: it's a separate (small) library so merging from time to time (like CLHEP) shouldn't be difficult. I just updated MetaIO with HEAD from https://github.com/Kitware/MetaIO, and after keeping Gates CMakeLists it compiles and works fine. @dsarrut I pushed it to the metaio branch on our server, I could not push here on github.

Going to disagree regarding the defaults: when attaching a voxellized output to a voxellized volume, and no geometric options are set, I would consider the logical default to inherit the geometric properties of the volume.

[dsarrut]

metaio updated now, thanks @brenthuisman !

(I will close the issue once the dicom-dev branch has been merged).