- Email: [email protected]
Numerical T1 computation from NMR intensity ratios
073&725X/86$3.00+.00 Copyright 0 1986 Pergamon Journals Ltd.
Magneric ResonanceImaging. Vol.4pp. 311-319.1986 Printed in the USA. All rights reserved.
l Research Article NUMERICAL
T1 COMPUTATION
FROM NMR INTENSITY RATIOS
MAX S. LIN, JAMES W. FLETCHER, FRANCIS K. HERBIG, AND ROBERT M. DONATI Nuclear MedicineService (115 JC), VA Medical Center, St. Louis, Missouri 63 106, and Department of Internal Medicine, St. Louis University School of Medicine, St. Louis, Missouri 63104.
Two-point measurementof tissue T, from NMR intensity ratios consists of forming an a priori ratio function describing a T, dependence of the ratio R(T,) and computing T, from an observed ratio Q by numerically solving R(T,) - Q = 0 or an equivalent equation. Impact of R(T,) designs on the numerical computation and dependence of relative speeds of three numerical methods on desired computational precisions q and on other factors are examined. All three methods begin with computing a table of R(T,) entries in uniform T, steps (AT,). In two iterative methods, a step containing the T, root is looked up, and the precise Tl location within the step is pinpointed to within a q value by either linear-interpolative (LI) or Newton-Raphson (NR)iteration.The third method simply consists of computing a large table of AT, = q for a mere “look-up” with no iterative search. All three methods require a monotonous R(T,) for uniformly effective computation over wide T, ranges. Speeds of either iterative method for computing Tl images are expected to vary with AT, and q with unsharp speed maxima at AT, near 20,6, and 2 ms for q = lo-‘, lo-*, and 10W3 ms, respectively. Either iterative method is suitable for both low- and high-precision computations, the LI method being generally faster. The simple look-up is the fastest of the three for T, image computation to low precisions of q > 1 ms, is likely the slowest for that to q = 0.1 ms, and is impractical for that to q ( 0.01 ms. Keywords: T, computation, Numerical method, Magnetic resonance imaging.
INTRODUCTION Tissue T, can be measured in commercial
subsequent iterative search tends to take less time. The pinpointing also depends on a computational precision q wanted of the ?, and on the iterative method used. Narrowing AT, to a width equal to a desired q value reduces the T, computation essentially to a simple look-up, dispensing with the iterative search. In this work, we explored proper designs of the ratio function from the standpoint of its suitability for numerical f, computation, which is an integral part of the two-point measurement. Dependence of the computing speeds of three numerical methods on various factors was investigated to seek guidelines for optimal use of these methods. One method is an iterative linear interpolation described elsewhere.7 Another is an adaptation of a standard Newton-Raphson iteration for solving transcendental equations.” The third is the noniterative simple look-up. Each method begins with computing a look-up table and then proceeds to a look-up.
‘H-NMR
imagers from ratios of two measured signals pixel by pixe1.6,‘0v’2A relation that equates an observed signal ratio Q to a R( T,) function describing an expected T, dependence of the ratio embodies the two-point method. A f, root of the equation, Q - R( f,) = 0, is the measured or calculated ?, as opposed to the true T,. The equation defies algebraic solution, but can be solved numerically in an iterative fashion.73” A “lookup” table of R vs T, is computed at uniform T, intervals (AT,) on a computer. From the table, the computer looks up a T, step over which the R interval contains a given Q. In suitably designed studies, this T, step will contain the f, sought and provide the initial conditions for pinpointing the f, by an iterative computation, wherein successive guesses converge on the f,. Tables with smaller T, steps take longer time to compute and also require longer look-up time, but the
RECEIVED l/26/85;
Address correspondence to Max S. Lin, M.D., Ph.D., Nuclear Medicine Service (115 JC), VA Medical Center, 9 I5 N. Grand Blvd., St. Louis, MO 63 106.
ACCEPTED l/24/86.
Acknowledgements-Author M.S.L. gratefully acknowledges spirited encouragements of Stephen Ayres, M.D., Department of Internal School of Medicine.
Medicine,
St.
Louis
University
311
312
Magnetic Resonance Imaging 0 Volume 4, Number 4, 1986 Table 1. Seven representative
STUDY DESIGNS The design of ratio functions We define R( T,) as
Sequence
ratio studies
parameters
(ms)
R(T,) notation
T,a
T,
TRI
TRZ
The S, is an inversion-recovery (IR) or a non-IR signal, and the S, always a non-IR one. We assume an ideal case,8 in which IR and spin echo (SE) signals normalized to their asymptote have the forms,
IRSE-A IRSE-B IRSE-C IRSE-D SESE-A SESE-B SESE-C
30 30 30 30 20 20 20
600 400 400 200
800 714 700 500 455 444 375
1700 1786 1500 1500 2045 1556 1125
IR: S = 1 - 2exp( - T,/T,)
“Values of TE are identical between S, and S, of Eq. I.
R(T,)
+ 2exp[-(T,
-
(1)
= &/S,.
- TE/2)/T,] + ew( - TRITI),
(3)
the S and TR being either S, and TR, or S2 and TR2 for Eq. 3, but always S, and TR2 for Eq. 2. Studies are referred to by the type of ratios formed to calculate ?,. Thus an IR/SE study implies a R(T,) in which SZ is given by Eq. 2, and S, by Eq. 3. Compared to this approach based on solving R(f,)
- Q = 0,
(4)
the reciprocal one based on R _, (?‘,) - Q, = 0. where R_, = 1/R and Q, = 1/Q, gives the same ?, root except for a computational failure in non-IR/lR studies (vide infra). Monotonicity
of a ratio
function
makes
its value
an
T,. The four IR/SE and three SE/SE R(T,) functions plotted in Fig. 1 are all monotonous and continuous in T,. Table 1 gives their sequence-parameter values. These seven studies were used to evaluate relative speeds of the three numerical indicator
of a unique
2
i
/
&
2750,150
i
/:
IR(Tn2vSEnn,)
27001200 \
___ _____2600, _____ 00
NT,)
\ 207
0
functions
‘..
\
\~\
_______________________ _
1,829
_ _ _ _ _ _ _ _ _ __ _ .
TdTFll TR( + T*x
R( T,)
2500 2500 2200 2000 2500 2000 1500
Look-up tables and look-up routines Due to noise errors, Eq. 4 predicts a ?, range much wider than the tissue T, range. Computable-?, range is limited to the table T, range, which should be at least
!’
monotonous
TRZ
methods and will be referred to by their identifying notations shown in Table 1. Stochastic simulation has shown that the IRSE-B study is superior to the other three IRSE studies for a wide coverage that minimizes f, noise errors overall across the 150- to 1500-ms T, band.’ The SESE-A is the best of the three SESE studies for the wide coverage, but is far inferior to the IRSE-B.9 In IR/SE studies, use of a large T,/T,, produces a nonmonotonous R(T,) (Fig. 2). For the same T,, TE, and TR, + TR2 as those of the IRSE-B in Fig. 1, an IR/SE R(T,) with TR, of 200 and 100 ms will exceed unity in values over the T, region, O-332 and O-481 ms, respectively. Each IRSE R(T,) in Figs. l-2 will give a reciprocal Rm,(T,) discontinuous at a T, where the R( T,) is null. The discontinuity causes a look-up failure.
3
Fig. 1. Plots of the seven identified in Table 1.
+
(2)
TE/2)1/TI - ev-T,dTJ
SE: S = 1 - 2exp[ -(TR
TRI
q 29cm
Fig. 2. Plots of five IR (T,,)/SE( TR,) ratio functions all of T, = 800 ms, TE = 20 ms, and TR, + TRz = 2900 ms. The TRz/TR, varies as indicated on individual lines. With increasing T,/ TR,, the nonmonotonous region (R > I ) widens.
Numerical T, computation 0 MAX S. LIN DAL.
as wide as the tissue T, range for a wide coverage. Setting f, values that are otherwise outside the table range to the nearer table cutoff value has a “bunching” effect. Use of a table range closer to the tissue range produces fewer computed unphysical f,, but compromises the contrast of constructed ?‘, images in the extreme T, regions. In this work, the look-up table contained [T,, R(T,)] entries evenly spaced at a AT, interval over the T, range, 45-1845 or 45-3045 ms. There were only two other entries, one at T, 40 ms and the other at 9000 ms for applications outside the present scope. In the R( T,) approach, a look-up process consists of a serial comparisons of a given Q against table R(T,) entries and eventually leads to a R(T,) step X,-Y, (Y,, > X0 by definition) satisfying X,, < Q < YO.
(5)
A monotonous R(T,) insures that the unknown T, is contained in x,-y,, the corresponding T, step with R(x,,) = X, and R(y,J = Y,. That is, x0 < f, < Yo
or
y, < T, < x0
search to converge on an erroneous f,. For Q_, that falls within a step straddling a discontinuity in the case of a discontinuous R_, (T,), there is no workable look-up routine. The R( T,) approach is preferable. iterative algorithms In the LI method, the ith interpolation
LI: ti =
_V_l -
i=
follows
[(yi-1 - Xi_,)
(K-1 - Q>/(yi-, - &,)I
(8)
1,2,3 ,...,
where Xi = R(xi) and Yi = R(yi) for each i.’ From the ith interpolation outcome ti, an iterative rule determines the operands for the (i + 1)th interpolation and forces either (Xi, xi) or ( yi, vi) to converge on (Q, f, ) .’ The same algorithm applies whether R(T,) is monotonously decreasing or increasing function of T,. In the NR method, the ith iteration follows
NR: ti = ti_, -
R(ti-,) - Q, i = 1 2 3 R’(ti_,) ’ ’ ““’
(9)
(6)
depending on rising or declining nature of the monotonicity. All look-up routines are based on a monotonicity of ratio functions. The serial comparison may begin at one end of the table and progress step by step toward the other end until the X,-Y, step of Eq. 5 is found. Alternatively, the comparison may involve a median R( T,) entry of a R(T,) region and narrow the Q whereabout by half the region. We refer to the first and second schemes as regular and fast, respectively. Look-up routines under each scheme can differ according to whether the routine is effective independent of increasing vs decreasing nature of the monotonicity. The inequality,
[Q - Ri(T,)I [Q - Rj(Tl)l < 0,
313
(7)
reveals that a Ri( T,) - R,( T,) step or region contains the Q independent of the nature of the monotonicity. Regular and fast routines were devised, and usually the fast ones were used with the R( T,)-based tables. In all routines, the look-up for each Q began with comparing it against R(45) and R(1845) or R(3045) to ascertain that it was contained in the R(45)-R( 1845) or R(45)-R(3045) range. The notation R(45) denotes a table entry or a Q value equal to the R(T,) value at T, 45 ms. For a Q that falls within a nonmonotonous R(T,) region, either type of routines can lead to an erroneous table step and cause a subsequent iterative
where R’(t,-,) is the derivative dR/dT, evaluated at T, = ti_,.” A ti from the NR iteration gives the f, root of Eq. 4. The iteration requires an initial guess to, which was usually taken to be the midpoint of the T, step looked up, to = ]y, - x01/2. No approximation was used in evaluating the dRf dT, . In both LI and NR methods, the iteration terminates when a desired precision is met: Iti - ti-11
5 9.
(10)
The termination after rth iteration gives t, as the f, sought. The precision requirement constitutes the criterion for convergence on the ??, root and is checked for i = 2, 3, 4, . . . intheLImethodandi= 1,2,3,...in the NR one. Regardless of the q value, the number of iterations required for the convergence (r) is at least two and one in the LI and NR methods, respectively. Measurement of computing times Precomputation of Q’s. Measurements of look-up and iterative computing times were taken on a precomputed set of evenly spaced Q’s. The number of Q’s in such a set is denoted by N. Ranges of Q values of a set were either narrow or wide. All Q’s in a narrowly ranging set required the same r iterations for convergence. Wide Q ranges were either R( lOO)-R( 1800) or R( 1OO)-R(3000). Tables covering R(45)-R( 1845) and R(45)-R(3045) were used to look up Q’s ranging R( lOO)-R( 1800) and R( lOO)-R( 3000), respectively.
314
Magnetic Resonance Imaging 0 Volume 4, Number 4, 1986
Component computing times. A LI or a NR computing run consisted of four stages: first, precomputation of a set of Q’s; second, calculation of a lookup table; third, looking up a step for each Q; and finally, iterative computation for the set of Q’s. The programs were flagged for timing on an Apple IIe microcomputer to give the second, third, and last stage times, respectively, as estimated table computing time (TT), look-up time (LT), and interpolation time (IPT) or Newton-Raphson time (NRT). The NRT value included the time to calculate t,. In runs on a set of widely ranging Q’s, the time estimate represented a speed averaging over the Q range. Single-block operations. In runs on 201 Q’s averaging over R( lOO)-R( 1800) and R( lOO)-R( 3000), respectively, the table AT, had to be at least 1.8 and 3.0 ms for the microcomputer memory to hold the table in its entiety. Use of larger tables with Al’, below these limits would require a piecemeal look-up operations for N 1 201. Such piecemeal operations on the microcomputer were not done to keep the measurement sufficiently simple for analytical and comparative purposes and to avoid arduously long computing time. Correction of estimates. The component times obtained by flagging computing stages represented an overestimate. The second through last stages each involved a looped computation, which did not actually commence immediately following the flagging. For a precise extrapolation of computing times from the AT, and N regions where measurements were taken to expected values in regions of smaller AT, and/or larger N, individual corrections to TT, LT, IPT, and NRT estimates were found and subtracted. To find the correction for TT, the TT estimate was obtained for varying AT, in runs averaging over wide Q ranges. Plotting the TT estimate against AT, for AT, < 200 ms gave a straight line fitting TT = k/AT, (k being a proportionality constant) except for a finite intercept, which arose from R(40) and R(9000) calculations preceding a looped computation of table R( T,) entries and from a delayed initiation of the looped computation. The intercept was the correction and amounted to 4.8-5.0 s among the various studies of Table 1. To find corrections for LT, IPT, and NRT, the estimates were taken in runs on varying number (N) of Q’s that all required the same r iterations for convergence. Plots of estimated LT. IPT, and NRT against N, respectively, gave straight lines with intercepts, 1.8, 1.6, and 2.8 s. which were the respective corrections. Averaging and scaling. Measurements averaging over wide Q ranges were usually taken using 201 Q’s. In LI and NR runs on Q sets ranging 10 l-601 in N and ranging R( lOO)-R(3000) in Q values for each N, the corrected LT, IPT, and NRT were all found to be
directly proportional to N. Clinical f, images derive from a number of nonuniformly distributed Q’s in the order of 128 x 256. Averaging with uniformly distributed Q’s provided a representative framework for comparisons. For evaluation of relative speeds at an N > 201, the LT, IPT, and NRT measured for 201 Q’s were linearly scaled up by the factor N/201. The scaling supposes a single-block operation. For convenience, the ?, computation for N r 64 x 64 is referred to at times as T, image computation. Total computing times on a minicomputer. To assess the order of magnitude of F,-image computing speeds of image processors, limited direct measurements were obtained on a Data General C-305 minicomputer “standing alone” to the exclusion of other service functions. Single-block operations on the minicomputer were possible for N values up to about 1K and 2K, respectively, in double- and single-precision computations when tables of 6-ms AT, over 3000-ms span were used. A look-up in tables of the said size for a set of Q’s and an iterative computation on that Q set were both repeated sufficient number of times for an N equivalent of 128 x 256. The sum of the time taken in calculating the table once plus that needed in the repetitive look-up and iterative pinpointing was the total computing time for N = 128 x 256 in a virtual single-block operation with no intrarun interaction between the central processing unit and the disc. No correction was applied to the measurement. RELATIVE
COMPUTING
SPEEDS
Iterative convergence and time per iteration Over the 50- to 3000-ms ?, range studied, the LI process always converged, and the NR one almost always. The LI convergence could be slow in short-f, regions when a wide AT, was used (Table 2). The NR process diverged when the initial guess was not sufficiently close to the f, sought. With AT, = 100 ms and y = 1 ms, the NR process converged on all F,‘s in the 50- to 1 14-ms region in all studies (Table 1) when the initial guess was 95 ms, the midpoint of the step looked up. When the initial guess was programed to be 45 ms however, the same process diverged for the same 4 for all f,‘s in the 90- to 114-ms region in all studies except IRSE-D. Divergence at 4 = 1 ms implies the same for all smaller y values. In computing runs on N evenly spaced Q’s that all required the same r iterations for convergence, the iteration time (IPT or NRT) divided by rNgave a time value, which tended to decrease with increasing r (Table 2, N = 101). At higher r values, IPT/rN and NRT/rN fell to a plateau value representing the arithmetic time of Eq. 8 or 9 per iteration (TPI) on the microcomputer. For IRSE and SESE studies, the TPI
Numerical T, computation 0 MAX S. LIN ETAL. Table 2. Interpolation
315
AT, = 120 ms
IPT (r)
60
345 322 176 104 103 103 103 80
90 120 200 500 1000 2000 3000
SESE-A NRT (r)
(14) (13) (7) (4) (4) (4) (4) (3)
344 232 232 177 175 177 176 121
1018 317 172 111 70 91 90 69
IRSE-B IPT (r)
NRT (r)
IPT (r)
(6) (4) (4) (3) (3) (3) (3) (2)
0.01ms using
AT, =4ms
IRSE-B i; ms
to q =
time (IPT) or Newton-Raphson time (NRT) in computing ?,on a microcomputer table steps, 120 or 4 ms, and the number of iterations required for convergencea
(49) (15) (8) (5) (3) (4) (4) (3)
301 204 203 156 107 156 155 107
(6) (4) (4) (3) (2) (3) (3) (2)
80 80 80 56 55 55 55 55
SESE-A IPT (r)
NRT (r) 177 120 121 121 120 121 120 122
(3) (3) (3) (2) (2) (2) (2) (2)
(3) (2) (2) (2) (2) (2) (2) (2)
69 69 69 49 49 49 48 49
NRT (r) 156 106 155 107 107 107 107 107
(3) (3) (3) (2) (2) (2) (2) (2)
(3) (2) (3) (2) (2) (2) (2) (2)
“The iterative IPT or NRT given is seconds taken to find 101 f;s that are clustered within 1 ms of the indicated f ,and require the same r iterations (given in parenthesis) for convergence. Table T, cutoffs are 45 and 3045 ms. Definitions of IRSE-B and SESE-A are given in Table 1.
of the Ll process was 0.24 and 0.21 s, respectively, and that of the NR process 0.55 and 0.48 s, respectively. A single LI iteration was about twice as fast as a single NR iteration. Table computing time and look-up time In LI and NR runs, the TT was significantly longer for IRSE studies compared to SESE, but the LT was essentially independent of the R(T,) type. Among various studies (Table 1) averaging over R(lOO)25r,,,,,,,,,,,,
R( 1800), there were a total of 22 LI or NR runs for each AT, value studied. Deviations of the 22 individual LT values from their mean were all under 3% of the mean for each of 17 AT, values ranging 2-200 ms, and differences among the TT values were similarly small. Figure 3 shows a fit of the mean TT value to log(TT) = logk - log(AT,). The line is representative of either the four IRSE or the three SESE studies. Figure 4 shows the mean LT value to be a linear function of log(AT,). Here both lines represent all seven studies. The linearity permitted a back-extrapolation to find TT (independent of N) and LT (N = 201) for AT, 5 2 ms.
125
\
LT = 727 . 332!q7(ATJ 100
75
NT,)
“S
0
-051 0
’
3045
l&SE
615
cu
3045
SESE
506
1845
IRSE
370
1845
SESE
305
’ 0.5
’
’ 1.0
’
LT
II
-
’ 1.5
LT = 120 - 332bgCAT)
50
’
’ 2.0
’
’ 2.5
’ 3.0
‘a9 (AT,)
Fig. 3. Plots of the common logarithm of table computing times (TT in seconds) against that of table step widths (AT, in milliseconds) for IRSE or SESE studies of Table 1. The table spans 1800 or 3000 ms and is identified by its upper bounds (UB). The k constant is proportional to the span. The TT was taken on a microcomputer.
01
’
’
0.5
’
’
l.0
’
’
1.5
’
’
20
’
’
2.5
’
30
IWCAT,l
Fig 4. Plots of “look-up” times (LT in seconds) for 201 ratios on a microcomputer against the logarithm of AT, in milliseconds. The look-up table covers the T,range, 45-l 845 or 45-3045 ms, identified by the upper bounds (Us).
Magnetic Resonance Imaging l Volume 4, Number 4, 1986
316
Speed dependence on A T, and 4 in LI and NR methods Tables 3 and 4 illustrate observations generally true (Table 1) in the AT, region above about 2 ms for all 9 values ranging 10-3-100 ms in both averagings using 201 Q’s, except as specifically noted. At all AT,, IPT and LT + IPT were shorter than NRT and LT + NRT, respectively (Table 3). Both IPT and NRT tended to decrease with narrowing AT, before reaching an absolute minimum fixed by the minimum possible r (Table 3). While not generally negligible relative to LT + IPT or LT + NRT for N = 201, the TT was negligible relative to either sum scaled up for N 2 64 x 64 (Table 3). Accordingly, the total computing time for p, images should be shorter with the Ll method than with the NR one. With decreasing AT,, one or more shallow minima occurred in LT + IPT or LT + NRT and so in the image computing time (Table 3). The shortest or nearly the shortest LT + NRT was about 1.5 times the LT + IPT counterpart (Table 4). The AT, that produced the shortest look-up plus iteration time was near 60, 20, 6, and 2 ms, respectively, for q = loo, lo-‘, 1O-‘, and 1O-3 ms in both LI and NR computations except for the NR one to I -ms q (Table 4). When LT + IPT or LT + NRT fell to nearly the shortest value with narrowing AT,, the r value for an overwhelming majority of Q’s had just fallen to two in all computations except the NR one to I-ms q (Table 4). In this exception, the r value was nearly exclusively one. Extending computations to the AT, region below 2 ms would produce LT lengthening that could not be compensated by any further IPT or
NRT shortening except possibly in NR computations to q = 10e3- 10-l ms (Tables 3 and 4). Minicomputing. The total computing time (TCT) for an equivalent of 128 x 256 IRSE-B ratios were measured at AT, = 6 ms and q = 1O-2 ms on the minicomputer in a single-block manner as described earlier. Five Q sets of differing N ranging 201-601 all averaging over R( lOO)-R(3000) were used in both LI and NR runs. The LI and NR TCT values were found, respectively, to range 78-79 s and 126-l 28 s in doubleprecision computations and 65-66 s and 101-103 s in single-precision ones. The NR TCT values were about 1.6 times the LI TCT. The relative LI vs NR speed at AT, above 2 ms on the microcomputer (Tables 3 and 4) translated into a practically significant difference on the minicomputer speed scale.
Extending considered AT, down to 0.03 ms To see if NR computations to q = 10-3- 10-l ms would improve in speeds at a small AT, below 2 ms, the n(r) distribution, i.e., the number of ?,‘s (among the 201) requiring r iterations, was determined in NR runs at a AT, -C2 ms on an IBM 370 mainframe computer. For the 201 f,‘s as a whole, the total number of iterations (TNI) required is Z,m(r). The NRT expected for a small AT, should be shorter than that measured at 3-ms AT, by an amount equal to the TPI times the TN1 difference. That is, NRT (AT,) = NRT(3) To account
- TPI[TNI(3)
for to computing
time, the NRT (AT,) was
Table 3. Table computing time (TT), look-up time (LT), and iterative interpolation time (IPT) or Newton-Raphson (NRT) all vs table step in computing T, to q = 0.01 ms from 201 ratios on a microcomputer” IRSE-B
(11)
- TNI(AT,)].
time
SESE-A
A T, ms
TT
LT
IPT
NRT
LT + IPT
LT + NRT
TT
LT
IPT
NRT
LT + IPT
LT + NRT
1000 600 315 200 100 50 25 15 IO 7.5 6 4 3
0.6 1.o 1.7 3.1 6.2 12.3 24.6 40.9 61.4 82.0 102.6 154.0 205.3
26.6 34.9 40.9 51.3 61.8 71.6 81.2 88.4 94.2 98.4 102.0 107.8 111.9
299 254 242 221 200 169 159 152 138 122 114 113 113
dvgb 437 402 373 342 336 316 295 264 247 245 241 241
326 289 283 272 262 240 240 240 233 220 216 221 224
dvgb 472 443 424 403 407 398 384 358 345 347 349 353
0.5 0.9 1.4 2.6 5.2 10.3 20.5 34.2 51.2 68.5 85.7 128.5 171.6
27.6 35.2 42.3 51.6 61.5 71.2 80.7 87.7 93.8 97.9 101.3 107.4 111.6
267 241 228 188 164 140 135 125 106 99 98 98 97
379 362 346 311 301 289 269 242 222 217 213 213 213
295 282 270 240 225 211 216 213 200 196 199 205 209
409 397 388 363 362 360 349 330 316 315 315 320 324
“Measurements are given in seconds. The ratios are evenly spaced over R( lOO)-R(3000). bNot measurable
due to divergence
in a short-T, region.
317
Numerical T, computation 0 MAX S. LIN ETAL.
Table 4. The shortest look-up plus pinpointing time found within the 1.8- to 900-ms A T, region and the A TI that produced
the said shortest time in computing 201 T,‘s to q = 10-3-100 ms from 201 ratios ranging R(lOO)-R(1800) linear-interpolative vs Newton-Raphson method on a microcomputer’ Linear interpolation
by
Newton-Raphson
R(T,)
loo
10-l
10-Z
10-3
loo
10-l
1O-2
lo-3
IRSE-A IRSE-B IRSE-C IRSE-D
179 176 177 179
193 191 193 195
209 208 208 211
227 225 226 230
242 241 240 240
323 321 320 320
339 338 338 336
356 356 353 357
SESE-A SESE-B
159 159
175 175
191 191
207 207
226 226
291 291
309 308
325 325
SESE-C
160
176
193
209
225
291
309
324
(185) (188) (186) (191)
20 20 20 15
(185) (190) (189) (195)
7.5 6 6 5
2 2 1.8 1.8
(176) (191) (192) (192)
(184) (191) (192) (177)
2 2 2 2
(0) (0) (0) (0)
20 20 20 15
(192) (193) (190) (197)
7.5 6 6 6
(188) (192) (191) (192)
2 2 2 1.8
(192) (189) (194) (195)
IRSE-A IRSE-B IRSE-C IRSE-D
60 60 60 50
SESE-A SESE-B
75 (184) 75 (182)
25 (184) 25 (183)
7.5 (188) 7.5 (182)
2.5 (187) 2.5 (188)
2 (0) 2 (0)
30 (182) 20 (191)
7.5 (190) 7.5 (191)
2 (193) 2 (194)
SESE-C
60 (190)
20 (186)
6 (187)
2 (187)
2 (1)
20 (191)
6 (190)
2 (194)
“Computational precisions are shown across the top. The upper half of the tables gives LT + IPT or LT + NRT in seconds. The lower half gives the said A T, in milliseconds followed in parenthesis by the number of f,‘s (among the 201) taking two iterations (r = 2) to find at that A T,. LT = look-up time. IPT = interpolation time. NRT = Newton-Raphson time.
not directly calculated as TPI x TN1 (AT,). The equations of Figs. 3 and 4 were used to find extrapolated TT and LT for N = 201 at the small AT,, and relative speeds were assessed in terms of the TCT for a specified N: TCT = TTt
[(LT + NRT)N/201].
(12)
No computation at AT, 5 2 ms brought the NR method up to the LI speed. Narrowing AT, to any value below 2 ms generally lengthened the NR TCT. These findings are illustrated with N = 128 x 256 in Table 5 using the IRSE-B study, for which the TPI was 0.55 s. The exceptional slight improvement in the NR speed for O.l-ms q and image-size N at AT, = 0.2 ms compared to AT, = 20 ms has little or no practical significance in view of the single-block assumption, which grossly underestimates the TCT at 0.2-ms AT,. At individual AT, optima, to recapitulate, the LI method is faster than the NR one for f,-image computations to any q in the range 10e3 - 10’ ms. For those to q = 10-3, 10-2, and 10-i ms by either method, the optima are near 2,6, and 20 ms, respectively. The simple look-up method Relative speeds of the simple look-up at AT, = q vs the LI method at AT, optima were assessed in terms of a TCT on the microcomputer. For the simple look-up, extrapolated TT and LT for N = 201 at AT, =
10m3- 10’ ms were calculated from the equations of Figs. 3-4, and the TCT was calculated as TT + (N/201)LT for a larger N. The LI TCT was calculated for the same N per Eq. 12, but replacing the NRT term with the IPT. Table 6 illustrates general findings for N ranging 64 x 64 to 256 x 256. In contrast to either iterative method, the TCT with the simple look-up increased rapidly with decreasing q. For computations to 1-ms q, the simple look-up was faster than the LI method for N exceeding about 2000, and slower for smaller N. The similar TCT values for O.l-ms q between simple look-up and LI methods under the single-block assumption (Table 6) suggest that either iterative method should be faster for O.l-ms q in practice. For computations to q 5 0.01 ms, the simple look-up is impractical.
DISCUSSION Signs of the IR signal in Eq. 2 were preserved in this work. The sign preservation assumes a phase-sensitive detection.4 The assumption is necessary to avoid a nonmonotonous (R(T,) 1 of the IR/non-IR variety (Figs. l-2) and to make an optimized IR/non-IR technique the best two-point method for observing T, in a wide-coverage’ or sharp-targeting’ fashion. The three numerical methods are readily extendable to two-point techniques that involve a multiple spin-echo (MSE) sequence intended for calculations of f2 as well
Magnetic Resonance Imaging l Volume 4, Number 4, 1986
318
Table 5. Newton-Raphson f,computation to q = 10-3-100 ms from IRSE-B ratios ranging R(lOO)-R(3000): n(r) distribution and component times in computing 201 fi’s and total computing times (TCT) per equation 12 for N = 128 x 256”
n(r)
TT
LT
(I)
S
S
S
0
62 42 0
139 159 201
205 246 308
112 114 117
164 154 130
754 731 678
20 3 I 0.5 0.3 0.2
6 0 0 0 0 0
I95 187 154 132 72 0
0 I4 47 69 129 201
31 205 615 1230 2050 3080
84 112 127 137 I45 150
244 233 215 203 170 130
894 942 941 945 889 815
10-2
7.5 3 I 0.3 0.1 0.05 0.03
10 0 0 0 0 0 0
191 200 194 184 164 125 66
0 I 7 17 37 76 135
82 205 615 2050 6150 12300 20510
98 112 127 145 160 170 I78
247 241 238 232 221 200 167
938 962 1000 1060 1140 1210 1280
10-r
3 2 I 0.1 0.05 0.03
32 5 2 0 0 0
169 193 199 197 194 190
0 3 0 4 7 11
205 308 615 6150 12300 20510
112 II7 127 160 170 I78
259 242 242 239 237 235
1010 982 1010 1190 1310 1460
4 ms
A TI ms
(3)
IO0
3 2.5 2
0 0
10-l
(2)
NRT
TCT min
“Under n(r) is given the number of ?,‘s (among 201) taking three, two, or one iteration to find. Time measurements on a microcomputer speed. TT = table computing time. LT = look-up time. NRT = Newton-Raphson time.
as ?,.I%’ As an excellent approximation for the purpose of calculating fI, full MSE signal expressions for the first echo can be simplified to a saturation-recovery form using an effective TR in place of the actual T 5.9.12 R.
Tissue T, ranges differ with the field strength.’ The table ranges of 45-1845 and 45-3045 ms with their somewhat arbitrary cutoffs are meant to be wider than a tissue T, range. Table entries were uniformly spaced in T, mainly for analytical purposes (Figs. 3-4). In a
Table 6. Computing
large-T,/ T, region of R( T,) where /dR/dT, 1 is very small, iterative convergence can be slow or impossible when the table step is wide (Tables 2-3). To facilitate convergence in such regions, table steps there can be made smaller than elsewhere. Little difference was found among the seven studies of Table 1 or between the two averaging ranges in AT, optima for f, image computation or in other comparative speed findings. The distribution of Q values in actual intensity-ratio images is diverse. A uniform Q
i$ to precision = q from 64 x 64 or 256 x 256 ratios by simple look-up vs iterative interpolation Relative total computing times (TCT) on a microcomputer speed scalea Interpolative
Simple look-up A T, (4) ms
TT .-___IRSE SESE
LT 201
IO0 10-1 10-z lO-3
10 103 1030 10250
2.12 2.67 3.23 3.78
8 85 848 8480
TCT( IRSE) ____~~___ 64 x 64 256 x 256 53 I57 1090 10330
702 975 2080 II480
TCT(SESE) ..___~__~~_ 64 x 64 256 x 256 52 139 914 8560
are based
700 957 1900 9710
q ms IO0 10-l lo-* lo-3
TCT vs q
IRSE-B ________ 64 x 64 256 x 256 62 68 75 87
997 1080 II80 1340
method:
SESE-A 64 x 64
256 x 256
56 61 68 76
896 979 1070 1170
“The TCT is calculated from measurements for 201 ratios ranging R( lOO)-R(3000). Also shown only for the simple look-up method (A T, = q) are LT = look-up time for 201 ratios and TT = table computing time. In the interpolation, the A T, used is near optimum. All measurements are given in minutes.
Numerical T, computation 0 MAX S. LINETAL.
distribution was used as a simple device to give a representative finding. Averaging over uniformly distributed Q’s instead of uniformly distributed T,‘s gives less weight to the asymptotic R(T,) regions and more weight to the more useful intermediate R(T,) region (Fig. l).’ Typically, stochastic accuracy of a two-point T, measurement is directly proportional to the square of S/N,’ and doubling image pixel resolutions can increase the inaccuracy by a factor about 16, a penalty that can not be reduced by averaging pixel T,.’ For this reason, the number of Q’s was limited to at most 256 x 512 in the comparative speed evaluation. There are multiple determinants of S/N other than the pixel resolution and multiple sources of p, errors other than the stochastic one. A reasonable percision to which f, should be computed in clinical studies depends on magnitudes of these noncomputational errors. Apart from the application to calculations of region-ofinterest p,‘s or whole FL maps from clinical intensity images, the numerical T, computation is a necessary
319
tool in simulation work analyzing stochastic or systematic f, errors in an optimization contextV3”*’To render purely computational errors negligible in clinical and analytical studies, respectively, computations to O.lms and 0.0 1-ms q levels probably suffice. With either the LI or the NR method, the absolute speed of ?‘, computation at 6-ms AT, from IRSE-B ratios to O.Ol-ms q on the minicomputer was about 450 times as fast as that on the microcomputer. On both computers, the LI method was 1.6 times as fast as the NR method. On a single-block basis, about the same relative speed is expected to apply to other computer systems. Similarly, details of algorithms used affected the absolute speed, but not the relative one in this work. Clearly, the comparative findings among the three methods studied will not quantitatively apply to computations on a piecemeal basis. Individual computing circumstances can be highly complex and need to be considered in using the present findings as a guide to an optimal application of the numerical methods studied.
REFERENCES 1. Bakker, C.J.G.; de Graaf, C.N.; van Dijk, P. Derivative of quantitative information in NMR imaging: a phantom study. Phys. Med. Biol. 29: 15 1 l-l 525; 1984. R.E.; Pfeif2. Bottomley, P.A.; Foster, T.H.; Argersinger, er, L.M. A review of normal tissue hydrogen NMR relaxation times and relaxation mechanisms from l-100 MHz: Dependence on tissue type, NMR frequency, temperature, species, excision, and age. Med. Phys. 11:425-448, 1984. 3. Hardy, C.J.; Edelstein, W.A.; Vatis, D. Calculated T, images derived from a partial saturation-inversion recovery pulse sequences. Scientific Program of the Society of Magnetic Resonance in Medicine, Third Annual Meeting (Abstract), pp. 299-300, 1984. R.E.; Nelson, T.R.; Hendee, W.R. Phase 4. Hendrick, detection and contrast loss in magnetic resonance imaging. Magn. Reson. Imag. 2:279-283; 1984. 5. Kurland, R.J. Strategies and tactics in NMR imaging relaxation time measurements. I. Minimizing relaxation time errors due to image noise-The ideal case. Magn. Reson. Med. 2: 136-l 58; 1985. of spin-lattice relaxation times 6. Lin, M.S. Measurement
in double spin-echo imaging. Magn. Reson. Med. 1:361369; 1984. 7. Lin, M.S. Interpolative computation of spin-lattice relaxation times from signal ratios. Magn. Reson. Med. 2~234-244; 1985. 8. Lin, MS. Accuracy of proton T, calculated by approximations from image signals. J. Nucl. Med. 26:54-58; 1985. 9. Lin, M.S.; Fletcher, J.W.; Donati, R.M. Two-point T, measurement: wide-coverage optimizations by stochastic simulations. Magn. Reson. Med. In press. 10. Pykett, I.L.; Rosen, B.R.; Buonanno, F.S.; Brady, T.J. Measurement of spin-lattice relaxation times in nuclear magnetic resonance imaging. Phys. Med. Biol. 28:723729; 1983. J.B. Numerical 11. Scarborough, 6th edition. Baltimore, Johns 213,1966.
mathematical analysis. Hopkins Press, pp. 201-
12. Technicare Corporation. Teslacon OTI documentation, preliminary version. Technicare Corporation, Solon, May 24, 1984.
Magneric ResonanceImaging. Vol.4pp. 311-319.1986 Printed in the USA. All rights reserved.
l Research Article NUMERICAL
T1 COMPUTATION
FROM NMR INTENSITY RATIOS
MAX S. LIN, JAMES W. FLETCHER, FRANCIS K. HERBIG, AND ROBERT M. DONATI Nuclear MedicineService (115 JC), VA Medical Center, St. Louis, Missouri 63 106, and Department of Internal Medicine, St. Louis University School of Medicine, St. Louis, Missouri 63104.
Two-point measurementof tissue T, from NMR intensity ratios consists of forming an a priori ratio function describing a T, dependence of the ratio R(T,) and computing T, from an observed ratio Q by numerically solving R(T,) - Q = 0 or an equivalent equation. Impact of R(T,) designs on the numerical computation and dependence of relative speeds of three numerical methods on desired computational precisions q and on other factors are examined. All three methods begin with computing a table of R(T,) entries in uniform T, steps (AT,). In two iterative methods, a step containing the T, root is looked up, and the precise Tl location within the step is pinpointed to within a q value by either linear-interpolative (LI) or Newton-Raphson (NR)iteration.The third method simply consists of computing a large table of AT, = q for a mere “look-up” with no iterative search. All three methods require a monotonous R(T,) for uniformly effective computation over wide T, ranges. Speeds of either iterative method for computing Tl images are expected to vary with AT, and q with unsharp speed maxima at AT, near 20,6, and 2 ms for q = lo-‘, lo-*, and 10W3 ms, respectively. Either iterative method is suitable for both low- and high-precision computations, the LI method being generally faster. The simple look-up is the fastest of the three for T, image computation to low precisions of q > 1 ms, is likely the slowest for that to q = 0.1 ms, and is impractical for that to q ( 0.01 ms. Keywords: T, computation, Numerical method, Magnetic resonance imaging.
INTRODUCTION Tissue T, can be measured in commercial
subsequent iterative search tends to take less time. The pinpointing also depends on a computational precision q wanted of the ?, and on the iterative method used. Narrowing AT, to a width equal to a desired q value reduces the T, computation essentially to a simple look-up, dispensing with the iterative search. In this work, we explored proper designs of the ratio function from the standpoint of its suitability for numerical f, computation, which is an integral part of the two-point measurement. Dependence of the computing speeds of three numerical methods on various factors was investigated to seek guidelines for optimal use of these methods. One method is an iterative linear interpolation described elsewhere.7 Another is an adaptation of a standard Newton-Raphson iteration for solving transcendental equations.” The third is the noniterative simple look-up. Each method begins with computing a look-up table and then proceeds to a look-up.
‘H-NMR
imagers from ratios of two measured signals pixel by pixe1.6,‘0v’2A relation that equates an observed signal ratio Q to a R( T,) function describing an expected T, dependence of the ratio embodies the two-point method. A f, root of the equation, Q - R( f,) = 0, is the measured or calculated ?, as opposed to the true T,. The equation defies algebraic solution, but can be solved numerically in an iterative fashion.73” A “lookup” table of R vs T, is computed at uniform T, intervals (AT,) on a computer. From the table, the computer looks up a T, step over which the R interval contains a given Q. In suitably designed studies, this T, step will contain the f, sought and provide the initial conditions for pinpointing the f, by an iterative computation, wherein successive guesses converge on the f,. Tables with smaller T, steps take longer time to compute and also require longer look-up time, but the
RECEIVED l/26/85;
Address correspondence to Max S. Lin, M.D., Ph.D., Nuclear Medicine Service (115 JC), VA Medical Center, 9 I5 N. Grand Blvd., St. Louis, MO 63 106.
ACCEPTED l/24/86.
Acknowledgements-Author M.S.L. gratefully acknowledges spirited encouragements of Stephen Ayres, M.D., Department of Internal School of Medicine.
Medicine,
St.
Louis
University
311
312
Magnetic Resonance Imaging 0 Volume 4, Number 4, 1986 Table 1. Seven representative
STUDY DESIGNS The design of ratio functions We define R( T,) as
Sequence
ratio studies
parameters
(ms)
R(T,) notation
T,a
T,
TRI
TRZ
The S, is an inversion-recovery (IR) or a non-IR signal, and the S, always a non-IR one. We assume an ideal case,8 in which IR and spin echo (SE) signals normalized to their asymptote have the forms,
IRSE-A IRSE-B IRSE-C IRSE-D SESE-A SESE-B SESE-C
30 30 30 30 20 20 20
600 400 400 200
800 714 700 500 455 444 375
1700 1786 1500 1500 2045 1556 1125
IR: S = 1 - 2exp( - T,/T,)
“Values of TE are identical between S, and S, of Eq. I.
R(T,)
+ 2exp[-(T,
-
(1)
= &/S,.
- TE/2)/T,] + ew( - TRITI),
(3)
the S and TR being either S, and TR, or S2 and TR2 for Eq. 3, but always S, and TR2 for Eq. 2. Studies are referred to by the type of ratios formed to calculate ?,. Thus an IR/SE study implies a R(T,) in which SZ is given by Eq. 2, and S, by Eq. 3. Compared to this approach based on solving R(f,)
- Q = 0,
(4)
the reciprocal one based on R _, (?‘,) - Q, = 0. where R_, = 1/R and Q, = 1/Q, gives the same ?, root except for a computational failure in non-IR/lR studies (vide infra). Monotonicity
of a ratio
function
makes
its value
an
T,. The four IR/SE and three SE/SE R(T,) functions plotted in Fig. 1 are all monotonous and continuous in T,. Table 1 gives their sequence-parameter values. These seven studies were used to evaluate relative speeds of the three numerical indicator
of a unique
2
i
/
&
2750,150
i
/:
IR(Tn2vSEnn,)
27001200 \
___ _____2600, _____ 00
NT,)
\ 207
0
functions
‘..
\
\~\
_______________________ _
1,829
_ _ _ _ _ _ _ _ _ __ _ .
TdTFll TR( + T*x
R( T,)
2500 2500 2200 2000 2500 2000 1500
Look-up tables and look-up routines Due to noise errors, Eq. 4 predicts a ?, range much wider than the tissue T, range. Computable-?, range is limited to the table T, range, which should be at least
!’
monotonous
TRZ
methods and will be referred to by their identifying notations shown in Table 1. Stochastic simulation has shown that the IRSE-B study is superior to the other three IRSE studies for a wide coverage that minimizes f, noise errors overall across the 150- to 1500-ms T, band.’ The SESE-A is the best of the three SESE studies for the wide coverage, but is far inferior to the IRSE-B.9 In IR/SE studies, use of a large T,/T,, produces a nonmonotonous R(T,) (Fig. 2). For the same T,, TE, and TR, + TR2 as those of the IRSE-B in Fig. 1, an IR/SE R(T,) with TR, of 200 and 100 ms will exceed unity in values over the T, region, O-332 and O-481 ms, respectively. Each IRSE R(T,) in Figs. l-2 will give a reciprocal Rm,(T,) discontinuous at a T, where the R( T,) is null. The discontinuity causes a look-up failure.
3
Fig. 1. Plots of the seven identified in Table 1.
+
(2)
TE/2)1/TI - ev-T,dTJ
SE: S = 1 - 2exp[ -(TR
TRI
q 29cm
Fig. 2. Plots of five IR (T,,)/SE( TR,) ratio functions all of T, = 800 ms, TE = 20 ms, and TR, + TRz = 2900 ms. The TRz/TR, varies as indicated on individual lines. With increasing T,/ TR,, the nonmonotonous region (R > I ) widens.
Numerical T, computation 0 MAX S. LIN DAL.
as wide as the tissue T, range for a wide coverage. Setting f, values that are otherwise outside the table range to the nearer table cutoff value has a “bunching” effect. Use of a table range closer to the tissue range produces fewer computed unphysical f,, but compromises the contrast of constructed ?‘, images in the extreme T, regions. In this work, the look-up table contained [T,, R(T,)] entries evenly spaced at a AT, interval over the T, range, 45-1845 or 45-3045 ms. There were only two other entries, one at T, 40 ms and the other at 9000 ms for applications outside the present scope. In the R( T,) approach, a look-up process consists of a serial comparisons of a given Q against table R(T,) entries and eventually leads to a R(T,) step X,-Y, (Y,, > X0 by definition) satisfying X,, < Q < YO.
(5)
A monotonous R(T,) insures that the unknown T, is contained in x,-y,, the corresponding T, step with R(x,,) = X, and R(y,J = Y,. That is, x0 < f, < Yo
or
y, < T, < x0
search to converge on an erroneous f,. For Q_, that falls within a step straddling a discontinuity in the case of a discontinuous R_, (T,), there is no workable look-up routine. The R( T,) approach is preferable. iterative algorithms In the LI method, the ith interpolation
LI: ti =
_V_l -
i=
follows
[(yi-1 - Xi_,)
(K-1 - Q>/(yi-, - &,)I
(8)
1,2,3 ,...,
where Xi = R(xi) and Yi = R(yi) for each i.’ From the ith interpolation outcome ti, an iterative rule determines the operands for the (i + 1)th interpolation and forces either (Xi, xi) or ( yi, vi) to converge on (Q, f, ) .’ The same algorithm applies whether R(T,) is monotonously decreasing or increasing function of T,. In the NR method, the ith iteration follows
NR: ti = ti_, -
R(ti-,) - Q, i = 1 2 3 R’(ti_,) ’ ’ ““’
(9)
(6)
depending on rising or declining nature of the monotonicity. All look-up routines are based on a monotonicity of ratio functions. The serial comparison may begin at one end of the table and progress step by step toward the other end until the X,-Y, step of Eq. 5 is found. Alternatively, the comparison may involve a median R( T,) entry of a R(T,) region and narrow the Q whereabout by half the region. We refer to the first and second schemes as regular and fast, respectively. Look-up routines under each scheme can differ according to whether the routine is effective independent of increasing vs decreasing nature of the monotonicity. The inequality,
[Q - Ri(T,)I [Q - Rj(Tl)l < 0,
313
(7)
reveals that a Ri( T,) - R,( T,) step or region contains the Q independent of the nature of the monotonicity. Regular and fast routines were devised, and usually the fast ones were used with the R( T,)-based tables. In all routines, the look-up for each Q began with comparing it against R(45) and R(1845) or R(3045) to ascertain that it was contained in the R(45)-R( 1845) or R(45)-R(3045) range. The notation R(45) denotes a table entry or a Q value equal to the R(T,) value at T, 45 ms. For a Q that falls within a nonmonotonous R(T,) region, either type of routines can lead to an erroneous table step and cause a subsequent iterative
where R’(t,-,) is the derivative dR/dT, evaluated at T, = ti_,.” A ti from the NR iteration gives the f, root of Eq. 4. The iteration requires an initial guess to, which was usually taken to be the midpoint of the T, step looked up, to = ]y, - x01/2. No approximation was used in evaluating the dRf dT, . In both LI and NR methods, the iteration terminates when a desired precision is met: Iti - ti-11
5 9.
(10)
The termination after rth iteration gives t, as the f, sought. The precision requirement constitutes the criterion for convergence on the ??, root and is checked for i = 2, 3, 4, . . . intheLImethodandi= 1,2,3,...in the NR one. Regardless of the q value, the number of iterations required for the convergence (r) is at least two and one in the LI and NR methods, respectively. Measurement of computing times Precomputation of Q’s. Measurements of look-up and iterative computing times were taken on a precomputed set of evenly spaced Q’s. The number of Q’s in such a set is denoted by N. Ranges of Q values of a set were either narrow or wide. All Q’s in a narrowly ranging set required the same r iterations for convergence. Wide Q ranges were either R( lOO)-R( 1800) or R( 1OO)-R(3000). Tables covering R(45)-R( 1845) and R(45)-R(3045) were used to look up Q’s ranging R( lOO)-R( 1800) and R( lOO)-R( 3000), respectively.
314
Magnetic Resonance Imaging 0 Volume 4, Number 4, 1986
Component computing times. A LI or a NR computing run consisted of four stages: first, precomputation of a set of Q’s; second, calculation of a lookup table; third, looking up a step for each Q; and finally, iterative computation for the set of Q’s. The programs were flagged for timing on an Apple IIe microcomputer to give the second, third, and last stage times, respectively, as estimated table computing time (TT), look-up time (LT), and interpolation time (IPT) or Newton-Raphson time (NRT). The NRT value included the time to calculate t,. In runs on a set of widely ranging Q’s, the time estimate represented a speed averaging over the Q range. Single-block operations. In runs on 201 Q’s averaging over R( lOO)-R( 1800) and R( lOO)-R( 3000), respectively, the table AT, had to be at least 1.8 and 3.0 ms for the microcomputer memory to hold the table in its entiety. Use of larger tables with Al’, below these limits would require a piecemeal look-up operations for N 1 201. Such piecemeal operations on the microcomputer were not done to keep the measurement sufficiently simple for analytical and comparative purposes and to avoid arduously long computing time. Correction of estimates. The component times obtained by flagging computing stages represented an overestimate. The second through last stages each involved a looped computation, which did not actually commence immediately following the flagging. For a precise extrapolation of computing times from the AT, and N regions where measurements were taken to expected values in regions of smaller AT, and/or larger N, individual corrections to TT, LT, IPT, and NRT estimates were found and subtracted. To find the correction for TT, the TT estimate was obtained for varying AT, in runs averaging over wide Q ranges. Plotting the TT estimate against AT, for AT, < 200 ms gave a straight line fitting TT = k/AT, (k being a proportionality constant) except for a finite intercept, which arose from R(40) and R(9000) calculations preceding a looped computation of table R( T,) entries and from a delayed initiation of the looped computation. The intercept was the correction and amounted to 4.8-5.0 s among the various studies of Table 1. To find corrections for LT, IPT, and NRT, the estimates were taken in runs on varying number (N) of Q’s that all required the same r iterations for convergence. Plots of estimated LT. IPT, and NRT against N, respectively, gave straight lines with intercepts, 1.8, 1.6, and 2.8 s. which were the respective corrections. Averaging and scaling. Measurements averaging over wide Q ranges were usually taken using 201 Q’s. In LI and NR runs on Q sets ranging 10 l-601 in N and ranging R( lOO)-R(3000) in Q values for each N, the corrected LT, IPT, and NRT were all found to be
directly proportional to N. Clinical f, images derive from a number of nonuniformly distributed Q’s in the order of 128 x 256. Averaging with uniformly distributed Q’s provided a representative framework for comparisons. For evaluation of relative speeds at an N > 201, the LT, IPT, and NRT measured for 201 Q’s were linearly scaled up by the factor N/201. The scaling supposes a single-block operation. For convenience, the ?, computation for N r 64 x 64 is referred to at times as T, image computation. Total computing times on a minicomputer. To assess the order of magnitude of F,-image computing speeds of image processors, limited direct measurements were obtained on a Data General C-305 minicomputer “standing alone” to the exclusion of other service functions. Single-block operations on the minicomputer were possible for N values up to about 1K and 2K, respectively, in double- and single-precision computations when tables of 6-ms AT, over 3000-ms span were used. A look-up in tables of the said size for a set of Q’s and an iterative computation on that Q set were both repeated sufficient number of times for an N equivalent of 128 x 256. The sum of the time taken in calculating the table once plus that needed in the repetitive look-up and iterative pinpointing was the total computing time for N = 128 x 256 in a virtual single-block operation with no intrarun interaction between the central processing unit and the disc. No correction was applied to the measurement. RELATIVE
COMPUTING
SPEEDS
Iterative convergence and time per iteration Over the 50- to 3000-ms ?, range studied, the LI process always converged, and the NR one almost always. The LI convergence could be slow in short-f, regions when a wide AT, was used (Table 2). The NR process diverged when the initial guess was not sufficiently close to the f, sought. With AT, = 100 ms and y = 1 ms, the NR process converged on all F,‘s in the 50- to 1 14-ms region in all studies (Table 1) when the initial guess was 95 ms, the midpoint of the step looked up. When the initial guess was programed to be 45 ms however, the same process diverged for the same 4 for all f,‘s in the 90- to 114-ms region in all studies except IRSE-D. Divergence at 4 = 1 ms implies the same for all smaller y values. In computing runs on N evenly spaced Q’s that all required the same r iterations for convergence, the iteration time (IPT or NRT) divided by rNgave a time value, which tended to decrease with increasing r (Table 2, N = 101). At higher r values, IPT/rN and NRT/rN fell to a plateau value representing the arithmetic time of Eq. 8 or 9 per iteration (TPI) on the microcomputer. For IRSE and SESE studies, the TPI
Numerical T, computation 0 MAX S. LIN ETAL. Table 2. Interpolation
315
AT, = 120 ms
IPT (r)
60
345 322 176 104 103 103 103 80
90 120 200 500 1000 2000 3000
SESE-A NRT (r)
(14) (13) (7) (4) (4) (4) (4) (3)
344 232 232 177 175 177 176 121
1018 317 172 111 70 91 90 69
IRSE-B IPT (r)
NRT (r)
IPT (r)
(6) (4) (4) (3) (3) (3) (3) (2)
0.01ms using
AT, =4ms
IRSE-B i; ms
to q =
time (IPT) or Newton-Raphson time (NRT) in computing ?,on a microcomputer table steps, 120 or 4 ms, and the number of iterations required for convergencea
(49) (15) (8) (5) (3) (4) (4) (3)
301 204 203 156 107 156 155 107
(6) (4) (4) (3) (2) (3) (3) (2)
80 80 80 56 55 55 55 55
SESE-A IPT (r)
NRT (r) 177 120 121 121 120 121 120 122
(3) (3) (3) (2) (2) (2) (2) (2)
(3) (2) (2) (2) (2) (2) (2) (2)
69 69 69 49 49 49 48 49
NRT (r) 156 106 155 107 107 107 107 107
(3) (3) (3) (2) (2) (2) (2) (2)
(3) (2) (3) (2) (2) (2) (2) (2)
“The iterative IPT or NRT given is seconds taken to find 101 f;s that are clustered within 1 ms of the indicated f ,and require the same r iterations (given in parenthesis) for convergence. Table T, cutoffs are 45 and 3045 ms. Definitions of IRSE-B and SESE-A are given in Table 1.
of the Ll process was 0.24 and 0.21 s, respectively, and that of the NR process 0.55 and 0.48 s, respectively. A single LI iteration was about twice as fast as a single NR iteration. Table computing time and look-up time In LI and NR runs, the TT was significantly longer for IRSE studies compared to SESE, but the LT was essentially independent of the R(T,) type. Among various studies (Table 1) averaging over R(lOO)25r,,,,,,,,,,,,
R( 1800), there were a total of 22 LI or NR runs for each AT, value studied. Deviations of the 22 individual LT values from their mean were all under 3% of the mean for each of 17 AT, values ranging 2-200 ms, and differences among the TT values were similarly small. Figure 3 shows a fit of the mean TT value to log(TT) = logk - log(AT,). The line is representative of either the four IRSE or the three SESE studies. Figure 4 shows the mean LT value to be a linear function of log(AT,). Here both lines represent all seven studies. The linearity permitted a back-extrapolation to find TT (independent of N) and LT (N = 201) for AT, 5 2 ms.
125
\
LT = 727 . 332!q7(ATJ 100
75
NT,)
“S
0
-051 0
’
3045
l&SE
615
cu
3045
SESE
506
1845
IRSE
370
1845
SESE
305
’ 0.5
’
’ 1.0
’
LT
II
-
’ 1.5
LT = 120 - 332bgCAT)
50
’
’ 2.0
’
’ 2.5
’ 3.0
‘a9 (AT,)
Fig. 3. Plots of the common logarithm of table computing times (TT in seconds) against that of table step widths (AT, in milliseconds) for IRSE or SESE studies of Table 1. The table spans 1800 or 3000 ms and is identified by its upper bounds (UB). The k constant is proportional to the span. The TT was taken on a microcomputer.
01
’
’
0.5
’
’
l.0
’
’
1.5
’
’
20
’
’
2.5
’
30
IWCAT,l
Fig 4. Plots of “look-up” times (LT in seconds) for 201 ratios on a microcomputer against the logarithm of AT, in milliseconds. The look-up table covers the T,range, 45-l 845 or 45-3045 ms, identified by the upper bounds (Us).
Magnetic Resonance Imaging l Volume 4, Number 4, 1986
316
Speed dependence on A T, and 4 in LI and NR methods Tables 3 and 4 illustrate observations generally true (Table 1) in the AT, region above about 2 ms for all 9 values ranging 10-3-100 ms in both averagings using 201 Q’s, except as specifically noted. At all AT,, IPT and LT + IPT were shorter than NRT and LT + NRT, respectively (Table 3). Both IPT and NRT tended to decrease with narrowing AT, before reaching an absolute minimum fixed by the minimum possible r (Table 3). While not generally negligible relative to LT + IPT or LT + NRT for N = 201, the TT was negligible relative to either sum scaled up for N 2 64 x 64 (Table 3). Accordingly, the total computing time for p, images should be shorter with the Ll method than with the NR one. With decreasing AT,, one or more shallow minima occurred in LT + IPT or LT + NRT and so in the image computing time (Table 3). The shortest or nearly the shortest LT + NRT was about 1.5 times the LT + IPT counterpart (Table 4). The AT, that produced the shortest look-up plus iteration time was near 60, 20, 6, and 2 ms, respectively, for q = loo, lo-‘, 1O-‘, and 1O-3 ms in both LI and NR computations except for the NR one to I -ms q (Table 4). When LT + IPT or LT + NRT fell to nearly the shortest value with narrowing AT,, the r value for an overwhelming majority of Q’s had just fallen to two in all computations except the NR one to I-ms q (Table 4). In this exception, the r value was nearly exclusively one. Extending computations to the AT, region below 2 ms would produce LT lengthening that could not be compensated by any further IPT or
NRT shortening except possibly in NR computations to q = 10e3- 10-l ms (Tables 3 and 4). Minicomputing. The total computing time (TCT) for an equivalent of 128 x 256 IRSE-B ratios were measured at AT, = 6 ms and q = 1O-2 ms on the minicomputer in a single-block manner as described earlier. Five Q sets of differing N ranging 201-601 all averaging over R( lOO)-R(3000) were used in both LI and NR runs. The LI and NR TCT values were found, respectively, to range 78-79 s and 126-l 28 s in doubleprecision computations and 65-66 s and 101-103 s in single-precision ones. The NR TCT values were about 1.6 times the LI TCT. The relative LI vs NR speed at AT, above 2 ms on the microcomputer (Tables 3 and 4) translated into a practically significant difference on the minicomputer speed scale.
Extending considered AT, down to 0.03 ms To see if NR computations to q = 10-3- 10-l ms would improve in speeds at a small AT, below 2 ms, the n(r) distribution, i.e., the number of ?,‘s (among the 201) requiring r iterations, was determined in NR runs at a AT, -C2 ms on an IBM 370 mainframe computer. For the 201 f,‘s as a whole, the total number of iterations (TNI) required is Z,m(r). The NRT expected for a small AT, should be shorter than that measured at 3-ms AT, by an amount equal to the TPI times the TN1 difference. That is, NRT (AT,) = NRT(3) To account
- TPI[TNI(3)
for to computing
time, the NRT (AT,) was
Table 3. Table computing time (TT), look-up time (LT), and iterative interpolation time (IPT) or Newton-Raphson (NRT) all vs table step in computing T, to q = 0.01 ms from 201 ratios on a microcomputer” IRSE-B
(11)
- TNI(AT,)].
time
SESE-A
A T, ms
TT
LT
IPT
NRT
LT + IPT
LT + NRT
TT
LT
IPT
NRT
LT + IPT
LT + NRT
1000 600 315 200 100 50 25 15 IO 7.5 6 4 3
0.6 1.o 1.7 3.1 6.2 12.3 24.6 40.9 61.4 82.0 102.6 154.0 205.3
26.6 34.9 40.9 51.3 61.8 71.6 81.2 88.4 94.2 98.4 102.0 107.8 111.9
299 254 242 221 200 169 159 152 138 122 114 113 113
dvgb 437 402 373 342 336 316 295 264 247 245 241 241
326 289 283 272 262 240 240 240 233 220 216 221 224
dvgb 472 443 424 403 407 398 384 358 345 347 349 353
0.5 0.9 1.4 2.6 5.2 10.3 20.5 34.2 51.2 68.5 85.7 128.5 171.6
27.6 35.2 42.3 51.6 61.5 71.2 80.7 87.7 93.8 97.9 101.3 107.4 111.6
267 241 228 188 164 140 135 125 106 99 98 98 97
379 362 346 311 301 289 269 242 222 217 213 213 213
295 282 270 240 225 211 216 213 200 196 199 205 209
409 397 388 363 362 360 349 330 316 315 315 320 324
“Measurements are given in seconds. The ratios are evenly spaced over R( lOO)-R(3000). bNot measurable
due to divergence
in a short-T, region.
317
Numerical T, computation 0 MAX S. LIN ETAL.
Table 4. The shortest look-up plus pinpointing time found within the 1.8- to 900-ms A T, region and the A TI that produced
the said shortest time in computing 201 T,‘s to q = 10-3-100 ms from 201 ratios ranging R(lOO)-R(1800) linear-interpolative vs Newton-Raphson method on a microcomputer’ Linear interpolation
by
Newton-Raphson
R(T,)
loo
10-l
10-Z
10-3
loo
10-l
1O-2
lo-3
IRSE-A IRSE-B IRSE-C IRSE-D
179 176 177 179
193 191 193 195
209 208 208 211
227 225 226 230
242 241 240 240
323 321 320 320
339 338 338 336
356 356 353 357
SESE-A SESE-B
159 159
175 175
191 191
207 207
226 226
291 291
309 308
325 325
SESE-C
160
176
193
209
225
291
309
324
(185) (188) (186) (191)
20 20 20 15
(185) (190) (189) (195)
7.5 6 6 5
2 2 1.8 1.8
(176) (191) (192) (192)
(184) (191) (192) (177)
2 2 2 2
(0) (0) (0) (0)
20 20 20 15
(192) (193) (190) (197)
7.5 6 6 6
(188) (192) (191) (192)
2 2 2 1.8
(192) (189) (194) (195)
IRSE-A IRSE-B IRSE-C IRSE-D
60 60 60 50
SESE-A SESE-B
75 (184) 75 (182)
25 (184) 25 (183)
7.5 (188) 7.5 (182)
2.5 (187) 2.5 (188)
2 (0) 2 (0)
30 (182) 20 (191)
7.5 (190) 7.5 (191)
2 (193) 2 (194)
SESE-C
60 (190)
20 (186)
6 (187)
2 (187)
2 (1)
20 (191)
6 (190)
2 (194)
“Computational precisions are shown across the top. The upper half of the tables gives LT + IPT or LT + NRT in seconds. The lower half gives the said A T, in milliseconds followed in parenthesis by the number of f,‘s (among the 201) taking two iterations (r = 2) to find at that A T,. LT = look-up time. IPT = interpolation time. NRT = Newton-Raphson time.
not directly calculated as TPI x TN1 (AT,). The equations of Figs. 3 and 4 were used to find extrapolated TT and LT for N = 201 at the small AT,, and relative speeds were assessed in terms of the TCT for a specified N: TCT = TTt
[(LT + NRT)N/201].
(12)
No computation at AT, 5 2 ms brought the NR method up to the LI speed. Narrowing AT, to any value below 2 ms generally lengthened the NR TCT. These findings are illustrated with N = 128 x 256 in Table 5 using the IRSE-B study, for which the TPI was 0.55 s. The exceptional slight improvement in the NR speed for O.l-ms q and image-size N at AT, = 0.2 ms compared to AT, = 20 ms has little or no practical significance in view of the single-block assumption, which grossly underestimates the TCT at 0.2-ms AT,. At individual AT, optima, to recapitulate, the LI method is faster than the NR one for f,-image computations to any q in the range 10e3 - 10’ ms. For those to q = 10-3, 10-2, and 10-i ms by either method, the optima are near 2,6, and 20 ms, respectively. The simple look-up method Relative speeds of the simple look-up at AT, = q vs the LI method at AT, optima were assessed in terms of a TCT on the microcomputer. For the simple look-up, extrapolated TT and LT for N = 201 at AT, =
10m3- 10’ ms were calculated from the equations of Figs. 3-4, and the TCT was calculated as TT + (N/201)LT for a larger N. The LI TCT was calculated for the same N per Eq. 12, but replacing the NRT term with the IPT. Table 6 illustrates general findings for N ranging 64 x 64 to 256 x 256. In contrast to either iterative method, the TCT with the simple look-up increased rapidly with decreasing q. For computations to 1-ms q, the simple look-up was faster than the LI method for N exceeding about 2000, and slower for smaller N. The similar TCT values for O.l-ms q between simple look-up and LI methods under the single-block assumption (Table 6) suggest that either iterative method should be faster for O.l-ms q in practice. For computations to q 5 0.01 ms, the simple look-up is impractical.
DISCUSSION Signs of the IR signal in Eq. 2 were preserved in this work. The sign preservation assumes a phase-sensitive detection.4 The assumption is necessary to avoid a nonmonotonous (R(T,) 1 of the IR/non-IR variety (Figs. l-2) and to make an optimized IR/non-IR technique the best two-point method for observing T, in a wide-coverage’ or sharp-targeting’ fashion. The three numerical methods are readily extendable to two-point techniques that involve a multiple spin-echo (MSE) sequence intended for calculations of f2 as well
Magnetic Resonance Imaging l Volume 4, Number 4, 1986
318
Table 5. Newton-Raphson f,computation to q = 10-3-100 ms from IRSE-B ratios ranging R(lOO)-R(3000): n(r) distribution and component times in computing 201 fi’s and total computing times (TCT) per equation 12 for N = 128 x 256”
n(r)
TT
LT
(I)
S
S
S
0
62 42 0
139 159 201
205 246 308
112 114 117
164 154 130
754 731 678
20 3 I 0.5 0.3 0.2
6 0 0 0 0 0
I95 187 154 132 72 0
0 I4 47 69 129 201
31 205 615 1230 2050 3080
84 112 127 137 I45 150
244 233 215 203 170 130
894 942 941 945 889 815
10-2
7.5 3 I 0.3 0.1 0.05 0.03
10 0 0 0 0 0 0
191 200 194 184 164 125 66
0 I 7 17 37 76 135
82 205 615 2050 6150 12300 20510
98 112 127 145 160 170 I78
247 241 238 232 221 200 167
938 962 1000 1060 1140 1210 1280
10-r
3 2 I 0.1 0.05 0.03
32 5 2 0 0 0
169 193 199 197 194 190
0 3 0 4 7 11
205 308 615 6150 12300 20510
112 II7 127 160 170 I78
259 242 242 239 237 235
1010 982 1010 1190 1310 1460
4 ms
A TI ms
(3)
IO0
3 2.5 2
0 0
10-l
(2)
NRT
TCT min
“Under n(r) is given the number of ?,‘s (among 201) taking three, two, or one iteration to find. Time measurements on a microcomputer speed. TT = table computing time. LT = look-up time. NRT = Newton-Raphson time.
as ?,.I%’ As an excellent approximation for the purpose of calculating fI, full MSE signal expressions for the first echo can be simplified to a saturation-recovery form using an effective TR in place of the actual T 5.9.12 R.
Tissue T, ranges differ with the field strength.’ The table ranges of 45-1845 and 45-3045 ms with their somewhat arbitrary cutoffs are meant to be wider than a tissue T, range. Table entries were uniformly spaced in T, mainly for analytical purposes (Figs. 3-4). In a
Table 6. Computing
large-T,/ T, region of R( T,) where /dR/dT, 1 is very small, iterative convergence can be slow or impossible when the table step is wide (Tables 2-3). To facilitate convergence in such regions, table steps there can be made smaller than elsewhere. Little difference was found among the seven studies of Table 1 or between the two averaging ranges in AT, optima for f, image computation or in other comparative speed findings. The distribution of Q values in actual intensity-ratio images is diverse. A uniform Q
i$ to precision = q from 64 x 64 or 256 x 256 ratios by simple look-up vs iterative interpolation Relative total computing times (TCT) on a microcomputer speed scalea Interpolative
Simple look-up A T, (4) ms
TT .-___IRSE SESE
LT 201
IO0 10-1 10-z lO-3
10 103 1030 10250
2.12 2.67 3.23 3.78
8 85 848 8480
TCT( IRSE) ____~~___ 64 x 64 256 x 256 53 I57 1090 10330
702 975 2080 II480
TCT(SESE) ..___~__~~_ 64 x 64 256 x 256 52 139 914 8560
are based
700 957 1900 9710
q ms IO0 10-l lo-* lo-3
TCT vs q
IRSE-B ________ 64 x 64 256 x 256 62 68 75 87
997 1080 II80 1340
method:
SESE-A 64 x 64
256 x 256
56 61 68 76
896 979 1070 1170
“The TCT is calculated from measurements for 201 ratios ranging R( lOO)-R(3000). Also shown only for the simple look-up method (A T, = q) are LT = look-up time for 201 ratios and TT = table computing time. In the interpolation, the A T, used is near optimum. All measurements are given in minutes.
Numerical T, computation 0 MAX S. LINETAL.
distribution was used as a simple device to give a representative finding. Averaging over uniformly distributed Q’s instead of uniformly distributed T,‘s gives less weight to the asymptotic R(T,) regions and more weight to the more useful intermediate R(T,) region (Fig. l).’ Typically, stochastic accuracy of a two-point T, measurement is directly proportional to the square of S/N,’ and doubling image pixel resolutions can increase the inaccuracy by a factor about 16, a penalty that can not be reduced by averaging pixel T,.’ For this reason, the number of Q’s was limited to at most 256 x 512 in the comparative speed evaluation. There are multiple determinants of S/N other than the pixel resolution and multiple sources of p, errors other than the stochastic one. A reasonable percision to which f, should be computed in clinical studies depends on magnitudes of these noncomputational errors. Apart from the application to calculations of region-ofinterest p,‘s or whole FL maps from clinical intensity images, the numerical T, computation is a necessary
319
tool in simulation work analyzing stochastic or systematic f, errors in an optimization contextV3”*’To render purely computational errors negligible in clinical and analytical studies, respectively, computations to O.lms and 0.0 1-ms q levels probably suffice. With either the LI or the NR method, the absolute speed of ?‘, computation at 6-ms AT, from IRSE-B ratios to O.Ol-ms q on the minicomputer was about 450 times as fast as that on the microcomputer. On both computers, the LI method was 1.6 times as fast as the NR method. On a single-block basis, about the same relative speed is expected to apply to other computer systems. Similarly, details of algorithms used affected the absolute speed, but not the relative one in this work. Clearly, the comparative findings among the three methods studied will not quantitatively apply to computations on a piecemeal basis. Individual computing circumstances can be highly complex and need to be considered in using the present findings as a guide to an optimal application of the numerical methods studied.
REFERENCES 1. Bakker, C.J.G.; de Graaf, C.N.; van Dijk, P. Derivative of quantitative information in NMR imaging: a phantom study. Phys. Med. Biol. 29: 15 1 l-l 525; 1984. R.E.; Pfeif2. Bottomley, P.A.; Foster, T.H.; Argersinger, er, L.M. A review of normal tissue hydrogen NMR relaxation times and relaxation mechanisms from l-100 MHz: Dependence on tissue type, NMR frequency, temperature, species, excision, and age. Med. Phys. 11:425-448, 1984. 3. Hardy, C.J.; Edelstein, W.A.; Vatis, D. Calculated T, images derived from a partial saturation-inversion recovery pulse sequences. Scientific Program of the Society of Magnetic Resonance in Medicine, Third Annual Meeting (Abstract), pp. 299-300, 1984. R.E.; Nelson, T.R.; Hendee, W.R. Phase 4. Hendrick, detection and contrast loss in magnetic resonance imaging. Magn. Reson. Imag. 2:279-283; 1984. 5. Kurland, R.J. Strategies and tactics in NMR imaging relaxation time measurements. I. Minimizing relaxation time errors due to image noise-The ideal case. Magn. Reson. Med. 2: 136-l 58; 1985. of spin-lattice relaxation times 6. Lin, M.S. Measurement
in double spin-echo imaging. Magn. Reson. Med. 1:361369; 1984. 7. Lin, M.S. Interpolative computation of spin-lattice relaxation times from signal ratios. Magn. Reson. Med. 2~234-244; 1985. 8. Lin, MS. Accuracy of proton T, calculated by approximations from image signals. J. Nucl. Med. 26:54-58; 1985. 9. Lin, M.S.; Fletcher, J.W.; Donati, R.M. Two-point T, measurement: wide-coverage optimizations by stochastic simulations. Magn. Reson. Med. In press. 10. Pykett, I.L.; Rosen, B.R.; Buonanno, F.S.; Brady, T.J. Measurement of spin-lattice relaxation times in nuclear magnetic resonance imaging. Phys. Med. Biol. 28:723729; 1983. J.B. Numerical 11. Scarborough, 6th edition. Baltimore, Johns 213,1966.
mathematical analysis. Hopkins Press, pp. 201-
12. Technicare Corporation. Teslacon OTI documentation, preliminary version. Technicare Corporation, Solon, May 24, 1984.