Package gov.nih.mipav.model.algorithms

Class AlgorithmSingleMRIImageSNR

  • All Implemented Interfaces:
    de.jtem.numericalMethods.calculus.function.RealFunctionOfOneVariable, java.awt.event.ActionListener, java.awt.event.WindowListener, java.lang.Runnable, java.util.EventListener
    public class AlgorithmSingleMRIImageSNR
extends AlgorithmBase
implements de.jtem.numericalMethods.calculus.function.RealFunctionOfOneVariable
    This program uses a mandatory signal 1 VOI, an optional signal 2 VOI, a mandatory noise background VOI, and the number of NMR receivers to calculate the signal to noise ratio for the signal VOI(s). The contrast to noise ratio is simply the signal to noise ratio for signal VOI 1 minus the signal to noise ratio for signal VOI 2. This program requires a single MRI magnitude image. The program mathematics assumes a General Rician or noncentral chi probability distribution characteristic of a MRI magnitude image. For the case of a single receiver system, the General Rician distribution reduces to the Rician distribution.

    The program calculates the SNR by 2 different methods. First, the signal to noise ratio is found using a equation for the first moment of Generalized Rice distribution. The first method follows the general 3 step approach outlined in the Constantinides reference: 1.) The noise standard deviation equals the square root((sum over i for background VOI of pixel(i) * pixel(i))/ (2 * number of background pixels * number of NMR receivers)) 2.) For each signal VOI the mean is calculated. 3.) The the mean is divided by the noise standard deviation The signal to noise ratio is calculated from the VOI mean to noise ratio using an equation (mean/standard deviation) = 1 * 3 * 5 * ... (2n - 1) * (1/((2**(n-1))*((n - 1)!))) * sqrt(PI/2) * 1F1(-1/2, n, -SNR*SNR/2) where n = number of NMR receivers in the phase array 1F1 is the confluent hypergeometric function of the first kind For the case of 1 NMR receiver the confluent hypergometric function can be replaced by modified Bessel functions as follows: (mean/standard deviation) = sqrt(PI/2) * exp(-SNR*SNR/4)((1 + (SNR*SNR/2))*I0(SNR*SNR/4) + (SNR*SNR/2)*I1(SNR*SNR/4)) However, since this Bessel function equivalence can only be used in the case of 1 receiver, a port of the ACM Algorithm 707 conhyp routine is used to find the confluent hypergeometric function for values of SNR * SNR/2 <= 3000. The port of the conhyp routine failed for values of SNR*SNR/2 > 3055, so for SNR*SNR/2 > 3000 the asymptotic limiting formula of the confluent hypergeometric function used in Appendix C of the Sato reference was used.

    The second method uses the maximum likelihood approach outlined in the Sijbers thesis. This approach finds the signal value which maximizes the value of the log of the maximum likelihood function (Equation 5.38): Let pixelNumber = number of pixels in the signal VOI log(L) equals approximately -pixelNumber*(n - 1)*log(signal) -pixelNumber*signal**2/(2*backgroundVariance) + sum from i = 1 to pixelNumber of (log(BESSEL_I of order (n - 1) of (pixel[i]*signal/backgroundVariance)) Note that this method does not work for large values of signal to noise ratio because the BESSEL_I function cannot handle large arguments. For example, the UNSCALED BESSEL_I of order 0 of x can only handle values of x up to about 699. For x >= 700, the Bessel UNSCALED function overflows. Since we require the log of the Bessel function, the range is extended by using the SCALED Bessel function which returns the BESSEL function * exp(-realArg), so we obtain the log(Bessel) from the log(returned function) + realArg. For the scaled Bessel function errors only start for x >= 32768. In finding the signal value that minimizes the negative of the maximum likelihood equation, a one dimensional search method called Brent's method that uses inverse parabolic interpolation calls the log of the maximum likelihood function. Brent's method is started with 3 points, 0.75 * first moment calculated signal, the first moment calculated signal, and 1.25 * first moment calculated signal. This assumes that the first moment result differs from the maximum likelihood moment by less than 25%.

    Testing the validity of this program is easy enough. Each receiver in the phase array has 2 standard deviations - one associated with the real part of the complex image and one associated with the imaginary part of the complex image. signal with noise = square root((signal without noise)**2 + sum from k = 1 to n of (sigmak real)**2 + (simgak imaginary)**2) So to simulate noisy data at each pixel we simply use 2*n Gaussian distributions in the above fashion. Running testAlgorithm with 1000 counts of a 100.0 signal and 1000 counts of a varying background for the 1 receiver case gave: background Calculated SNR Ideal SNR from first moment 0.5 198.0 200.0 1.0 101.0 100.0 2.0 48.9 50.0 5.0 19.7 20.0 20.0 5.10 5.00 50.0 2.28 2.00 100.0 1.09 1.00

    A second run of testAlgorithm showed good agreement between the 2 methods of finding the SNR: background First moment SNR Maximum likelihood SNR 100.0 1.21 1.16 50.0 2.20 2.21 20.0 5.09 5.09 10.0 10.0 10.0 5.0 20.0 20.0

    For 1F1(-1/2, 1, x) tested with Shanjie Zhang and Jianming Jin Computation of Special Functions CHGM routine and ACM Algorithm 707 conhyp routine by Mark Nardin, W. F. Perger, and Atul Bhalla the results were: From x = -706 to x = +184 the results of the 2 routines matched to at least 1 part in 1E5.

    For x = -3055 to x = +184 the ACM routine seemed to give valid results. For x <= -3056 the ACM routine results oscillated wildly and at x = -3200 the result was frozen at 1.0E75. For x >= +185 the ACM result was frozen at 1.0E75. The log result was also seen to oscillate for x <= -3056.

    For x = -706 to x = +781 the CHGM routine seemed to give valid results. For x = -707 to x = -745 the CHGM routine gave infinity and for x <= -746 the CHGM routine gave NaN. For x >= +783 the CHGM routine gave a result of -infinity.

    So for x = -(signal/noise)*(signal/noise)/2, the CHGM routine can handle a maximum signal/noise value of 37.5 and in standard output the ACM routine can handle a maximum signal/noise value of 78.166.

    For 1F1(-1/2, 2, x) the results were very similar: From x = -706 to x = +189 the results of the 2 routines matched to at least 1 part in 1E5.

    From x = -3058 to x = +189 the ACM routine seemed to give valid results. At x <= -3059 the ACM results oscillated wildly. For x > +189 the ACM result was frozen at 1.0E75.

    For x = -706 to x = +791 the CHGM routine seemed to give valid results. For x = -707 to x = -745 the CHGM routine gave infinity and for x <= -746 the CHGM routine gave NaN. For x >= +792 the CHGM routine gave a result of -infinity.

    In this module use the ACM routine for x >= -3000. For large negative x, use the formula: As |x| approaches infinity, 1F1(a, b, x) approaches (gamma(b)/gamma(b-a))*((-x)**(-a))*[1 + Order(|x|**-1)} for the real part of x < 0. More specifically from the Sato reference, the asymptotic behavior of the confluent hypergeometric function is: 1F1(a, b, x) approaches (gamma(a)/gamma(b-a)*((-x)**(-a))*G(a, a-b+1, -x) where G is an asymptotic series G(a, b, x) = 1 + (a*b)/(1!*x) + (a*(a+1)*b*(b+1))/(2!*x**2) + ...

    References: 1.) PhD. Thesis Signal and Noise Estimation From Magnetic Resonance Images by Jan Sijbers, Universiteit Antwerpen, Department Natuurkunde. 2.) Signal-to-noise Measurements in Magnitude Images from NMR Phased Arrays by Chris D. Constantinides, Ergin Atalar, Elliot R. McVeigh, Magnetic Resonance in Medicine, November, 1997, 38(5), pp. 852-857. Erratum in Magnetic Resonance in Medicine, July, 2004, 52(1), p. 219. 3.) I. S. Gradshteyn and I. M. Rhyzik, Tables of integrals, series, and products, Sixth edition. 4.) Tohru Sato, Liviu F. Chibotaru, and Arnout Ceulemans, The Exe dynamic Jahn-Teller problem: A new insight from the strong coupling limit, The Journal of Chemical Physics, Vol. 122, 054104, 2005, Appendix C.

    Derivation of Bessel function formula when number of receivers == 1 E[M**v] = (2*sigma**2)**(v/2) * gamma(1 + v/2) * 1F1(-v/2, 1, -A**2/(2*sigma**2)) E[M] = sigma * sqrt(2) * gamma(3/2) * 1F1(-1/2, 1, -A**2/(2*sigma**2)) gamma(1/2) = sqrt(PI) gamma(n+1) = n*gamma(n) gamma(3/2) = sqrt(pi)/2 E[M] = sigma * sqrt(PI/2) * 1F1(-1/2, 1, -A**2/(2*sigma**2)) 1F1(alpha, gamma, z) = exp(z) * 1F1(gamma - alpha, gamma, -z) Let pow = A**2/(2*sigma**2) 1F1(-1/2, 1, -pow) = exp(-pow) * 1F1(3/2, 1, pow) = exp(-pow) * [1/2 * 1F1(3/2, 1, pow) + 1/2 * 1F1(3/2, 1, pow)]

    alpha*1F1(alpha+1, gamma, z) = (z + 2*alpha - gamma)*1F1(alpha, gamma, z) + (gamma - alpha) * 1F1(alpha - 1, gamma, z) 1F1(-1/2, 1, -pow) = exp(-pow) * [0.5 * 1F1(3/2, 1, pow) + pow * 1F1(1/2, 1, pow) + 0.5 * 1F1(-1/2, 1, pow)] = exp(-pow) * [(1 + pow) * 1F1(1/2, 1, pow) + 0.5*(1F1(3/2, 1, pow) - 1F1(1/2, 1, pow)) - 0.5*(1F1(1/2, 1, pow) - 1F1(-1/2, 1, pow))]

    (z/gamma)*1F1(alpha+1, gamma+1, z) = 1F1(alpha+1, gamma, z) - 1F1(alpha, gamma, z)

    1F1(-1/2, 1, -pow) = exp(-pow) * [(1 + pow)* 1F1(1/2, 1, pow) + (pow/2) * 1F1(3/2, 2, pow) - (pow/2) * 1F1(1/2, 2, pow)] = exp(-pow) * [(1 + pow) * 1F1(1/2, 1, pow) + ((pow/2)**2)*1F1(3/2, 3, pow)] = exp(-pow/2) * [(1 + pow) * exp(-pow/2) * 1F1(1/2, 1, pow) + ((pow/2)**2)* exp(-pow/2) * 1F1(3/2, 3, pow)]

    Iv(x) = ((2**-v)/gamma(v+1)) * (x**v) * exp(-x) * 1F1(1/2 + v, 1 + 2*v, 2*x) I0(x) = exp(-x) * 1F1(1/2, 1, 2*x) I1(x) = (x * exp(-x)/2) * 1F1(3/2, 3, 2*x)

    1F1(-1/2, 1, -pow) = exp(-pow/2) * [(1 + pow) * I0(pow/2) + pow * I1(pow/2)]

    • Field Detail

      • automatic

        private boolean automatic

      • backgroundIndex

        private int backgroundIndex
        The index of a rerquired noise background VOI.

      • backgroundVariance

        private double backgroundVariance
        DOCUMENT ME!

      • nf

        private java.text.DecimalFormat nf
        DOCUMENT ME!

      • numReceivers

        private int numReceivers
        The number of NMR receivers.

      • signal2Index

        private int signal2Index
        The index of a second optional signal VOI. >= 0 if present

      • signalBuffer

        private float[] signalBuffer
        DOCUMENT ME!

      • signalIndex

        private int signalIndex
        The index of a required signal VOI.

      • useMaxLikelihood

        private boolean useMaxLikelihood
        DOCUMENT ME!

    • Constructor Detail

      • AlgorithmSingleMRIImageSNR

        public AlgorithmSingleMRIImageSNR​(ModelImage srcImg,
                                  boolean automatic,
                                  int signalIndex,
                                  int signal2Index,
                                  int backgroundIndex,
                                  int numReceivers)
        Creates a new AlgorithmSingleMRIImageSNR object.
        Parameters:
        srcImg - source image
        automatic - If true, use fuzzy c means to separate all pixels into signal and background
        signalIndex - the index of the signal VOI
        signal2Index - the index of the signal2 VOI if >= 0
        backgroundIndex - the index of the background VOI
        numReceivers - the number of NMR receivers

    • Method Detail

      • finalize

        public void finalize()
        Prepares this class for destruction.
        Overrides:
        finalize in class AlgorithmBase

      • funcC

        private double funcC​(double meanDivStdDev,
                     boolean signal)
        This routine iterates to find the solution to the first moment equation for the Generalized Rice distribution.
        Parameters:
        meanDivStdDev - DOCUMENT ME!
        signal - DOCUMENT ME!
        Returns:
        DOCUMENT ME!

      • receiverMaximumLikelihood

        private double receiverMaximumLikelihood​(double signal)
        For an input signal value this routine returns the log of the maximimum likelihood function.
        Parameters:
        signal - DOCUMENT ME!
        Returns:
        DOCUMENT ME!

      • testAlgorithm

        private void testAlgorithm()
        DOCUMENT ME!

      • eval

        public double eval​(double x)
        Specified by:
        eval in interface de.jtem.numericalMethods.calculus.function.RealFunctionOfOneVariable