Library for calculating gamma index? (Preferrably in R or Python, but any language is OK)

Question

In physics--especially medical physics--the gamma index is a criterion for comparing data from two particle detectors. More abstractly, the gamma index takes two 2D arrays (let's say array1 and array2) and compares each element of array1 with spatially-nearby elements of array2.

There are hundreds of academic papers that use the gamma index in their analysis sections. These papers don't seem to mention what tools/libraries they use to calculate the gamma index. It's possible the authors implement their own gamma index calculations (it's not that hard). However, I'm guessing that there are libraries/extensions/tools for calculating a gamma index.

Can anyone suggest a gamma index library to use in R or Python? (Other languages would be ok if there's nothing off-the-shelf for Python or R.)

Answer 1

I found basic MATLAB implementation of the 2D gamma index in Appendix A of this thesis.

I copy/pasted the following code from the thesis, and I made a couple of simplifications for readability. I talked to the author and confirmed that my version of the code (below) is correct. Recently, I have been using this code in the analysis portion of a medical physics study that I'll be publishing soon.

The inputs A1 and A2 are 2D arrays (which, in practice, are dose maps or fluence maps). We let A1 serve as the reference data, and A2 is the data that is being evaluated. If we use a typical 2%, 2mm acceptance criterion, then we set distance to agreement as DTA=2mm, and we set the dose threshold dosed=0.02, which is 2%.

In this simple implementation, we assume that the array indices are spaced in 1mm increments. If your data doesn't use 1mm increments, then scale your dosed value accordingly (e.g. if your A1 and A2 are in 0.5mm increments, then use DTA=4 to get a 2mm criterion).

The output, G, is a 2D array of gamma values.

function G = gamma2d (A1, A2, DTA, dosed)
    size1=size (A1) ;
    size2=size (A2) ;
    dosed = dosed *  max(A1 ( : ) ) ; %scale dosed as a percent of the maximum dose

    G=zeros ( size1 ) ; %this will be the output
    Ga=zeros ( size1 ) ;
    if size1 == size2
        for i = 1 : size1( 1 )
            for j = 1 : size1( 2 )
                for k = 1 : size1( 1 )
                    for l = 1 : size1( 2 )
                        r2 = ( i - k )^2 + (j - l) ^2 ; %distance (radius) squared
                        d2 = ( A1( i , j ) - A2( k , l ) )^2 ; %difference squared
                        Ga( k , l ) = sqrt(r2 / (DTA^2) + d2/ dosed ^ 2);
                    end
                end
                G( i , j )=min(min(Ga)) ;
            end
        end
    else
        fprintf=('matrices A1 and A2 are do not share the same dimensions! \n')
    end
end

To see an explanation of the gamma index in math notation, I recommend looking at this blog post.