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

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.