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 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 = +185 the ACM result was frozen at 1.0E75. The log result was also seen to oscillate for x

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 = +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 +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 = +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

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)]