Add common PET lession evaluation metrics

[AnderBiguri]

In PET there are many image evaluation metrics that are common in the clinic and it would be good for SIRF to be able to provide some of these. As @ALEXJAZZ008008 has already pointed out, some of these are either badly standardized or simply not standardized, and thus providing the actual implementation may help current and future users of SIRF a common code to use for this.

We can add, among others:

• SUVmax
• SUVmean
• SUVpeak (various)

Any other? Is this a good idea?

[AnderBiguri]

I have adapted @ALEXJAZZ008008's code to take SIRF inputs for the EANM definition of SUVpeak:

def SUVpeak(img,loc,tumour_diameter=None):
    """ Computes the SUVpeak by EANM convention[1]

        SUVpeak is the maximum mean value of an uniform sphere of 1.2cm of diameter around the tumour area. 

        Input: 
            - img, a SIRF.ImageData object
            - loc, a location around where you want to search the tumour
            - tumour_diameter, the expected diameter of the tumour in loc. If None (default), 12 mm is used.
        Output:
            - SUVpeak of that tumour
    [1]Boellaard, R., Delgado-Bolton, R., Oyen, W.J., Giammarile, F., Tatsch, K., Eschner, W., Verzijlbergen, F.J., Barrington, S.F., Pike, L.C., Weber, W.A. and Stroobants, S., 2015. 
    FDG PET/CT: EANM procedure guidelines for tumour imaging: version 2.0. 
    European journal of nuclear medicine and molecular imaging, 42(2), pp.328-354.
    """

    if isinstance(img, SIRF.ImageData):
        voxel_sizes=numpy.array(img.voxel_sizes())
        img=img.as_array()
    else:
        raise ValueError("First input image is not of type SIRF.ImageData, is of type " + type(img) + " instead.")
    
    if tumour_diameter is None:
        tumour_diameter = 12 # SUVpeak sphere diameter

    # It should be 1.2 cm sphere
    radius = 6

    # lets crop a VOI around the tumour location.
    max_off = np.ceil(tumour_diameter * 2 / voxel_sizes)
    roi = img[numpy.max([int(numpy.round(loc[0])-max_off[0]),0]):int(numpy.round(loc[0])+max_off[0]), 
              numpy.max([int(numpy.round(loc[1])-max_off[1]),0]):int(numpy.round(loc[1])+max_off[1]), 
              numpy.max([int(numpy.round(loc[2])-max_off[2]),0]):int(numpy.round(loc[2])+max_off[2])]

    # Create a filter
        # make sure its odd
    conv_size = np.int16(numpy.ceil(radius*2 / voxel_sizes) // 2 * 2 + 1)
    conv_filter = numpy.zeros((conv_size[0], conv_size[1], conv_size[2]))
        # mask
    X, Y, Z = numpy.ogrid[:conv_size[0], :conv_size[1], :conv_size[2]]
    dist_from_center = numpy.sqrt(voxel_sizes[0]**2 * (X - conv_size[0]/2)**2 + \
                                  voxel_sizes[1]**2 * (Y - conv_size[1]/2)**2 + \
                                  voxel_sizes[2]**2 * (Z - conv_size[2]/2)**2)
    conv_filter = dist_from_center <= radius # Should we have interpolation, instead of sharp filter boundary? it also fits the definition
    conv_filter = conv_filter.astype(float)
        # normalize
    conv_filter = conv_filter / numpy.sum(conv_filter)
  
    # Convolution with the image
    conv_output = scipy.ndimage.convolve(roi, conv_filter, mode='constant')
    # maximum value
    suv_peak = numpy.nanmax(conv_output)
    
    return suv_peak

[KrisThielemans]

I certainly think this is a good idea, but ideally we check carefully if our implementation fits the EANM definition (which I didn't do!).

A few things that I wonder about now that I look at it

• 'loc' would have to be in LPS coordinates, not indices. Probably you'd want to pass an ROI as opposed to location/diameter.
• I would suspect that we need to use floor/ceiling in the ranges as opposed to 2 rounds
• Not sure what they recommend for convolution and cropping (e.g. does it need to be cuboid cropping with constant boundary conditions?)

[AnderBiguri]

Good points!

• loc definetly needs to be in LPS. 100% agree.
• I used round because max_off is already wide enough, but yeah, agree that its tidier with floor/ceil.
• I have not seen any recommendation. There is a reference to a paper thrown there where they mention the VOI, but have not seen any explicit mention on the document itself. I can check again.

In order to get this out of a "discussion" and make a draft PR with it, one thing needs to be discussed first: where do we put it? as a method of ImageData, or as some new/existing class/method on utilities?