OpenGATE / GateTools

Python tools and functions for Gate
GNU Lesser General Public License v3.0
30 stars 22 forks source link

Method roi_utils.get_mask() does not work for HFP patient set up #22

Open Golpette opened 4 years ago

Golpette commented 4 years ago

Hi,

I tried to use this tool in a radiotherapy simulation in which the patient was in head first prone position (HFP). The mhd image I am using has the header:

ObjectType = Image NDims = 3 BinaryData = True BinaryDataByteOrderMSB = False CompressedData = False TransformMatrix = -1 0 0 0 -1 0 0 0 1 Offset = 250 250 -100 CenterOfRotation = 0 0 0 AnatomicalOrientation = LPI ElementSpacing = 0.97656200000000004 0.97656200000000004 2.5 DimSize = 512 512 164 ElementType = MET_SHORT ElementDataFile = ct.raw

Seems the issue is that the get_mask method is expecting the Origin to represent the centre of the voxel with the minimum coordinates, which is not the case for this HFP image.

I'm not comfortable enough with ITK / mhd / dicom to suggest a robust fix to this. But my hack solution for now is to determine a "direction" for each axis using GetDirection() (corresponding to the mhd TransformMatrix). I'm then using these +/-1 values to determine if the Origin field is representing a max or min position of the image. So I have:

(1) Replaced lines 517-518 of roi_utils.py to find the centre of the minimal voxel and the direction of the axes:

        orig = roimask.GetOrigin()  
        space = roimask.GetSpacing()   

        # directions +/- 1 for each axis       
        directions = list( roimask.GetDirection()*(1,1,1) )         

        # Find minimum corner (voxel centre)
        orig_min_corner = []
        for o,dim,s,direction in zip(orig,dims,space,directions):
            if direction>0:
                orig_min_corner.append(o)
            else:
                orig_min_corner.append(o-(dim-1)*s)

(2) Give this orig_min_corner to check "contained" in line 527:

for o,s,d,rmin,rmax in zip(orig_min_corner,space,dims,self.bb.mincorner,self.bb.maxcorner):

(3) Replaced 537-538 in case this min/max orientation issue ever happens with the z-axis too:

        if directions[2]>0:
            zmin = orig[2] - 0.5*space[2]
            zmax = orig[2] + (dims[2]-0.5)*space[2]
        else:
            zmin = orig[2] + (0.5-dims[2])*space[2]
            zmax = orig[2] + 0.5*space[2]

(4) Replaced 557-558 to give the xpoints and ypoints the correct ordering, using the 'directions' list above, before feeding them into the meshgrid:

        xpoints=np.linspace(orig[0],orig[0]+directions[0]*space[0]*dims[0],dims[0],False)
        ypoints=np.linspace(orig[1],orig[1]+directions[1]*space[1]*dims[1],dims[1],False)

The above seems to have solved my problem but I'm not comfortable enough with ITK/mhd/dicom images and their properties, or the rest of roi_utils.py, to know if this hack will generate problems elsewhere.

djboersma commented 4 years ago

Hi @Golpette , thanks for the report. I must admit that I only tested this code with HFS images and I agree that we should make sure that it also works with other anatomical orientations. I think I saw in some of the other gate tools that @tbaudier managed to ensure that they work with other orientations. The roi utils code is more messy than the other tools, unfortunately. We'll have look at your proposed fix. I'm currently very busy with other things, so I cannot promise that this is fixed very soon, but hopefully within a few weeks this should at least partly be addressed.

Golpette commented 4 years ago

Great, I'll also update the thread if I manage to do anything extra. I think using GetDirection() as above is fine if the matrix is always diagonal with values of +/-1 (which it should be for any standard image for radiotherapy). I initially wanted to just take the AnatomicalOrientation tag from the mhd file ("RAI", "LPI" etc) and generate the directions list from that, but as far as I can see ITK totally ignores that information.