**MEDICAL IMAGING PROJECTS**

Edge-flat-grey scale image enhancement by convex functional minimization.

50 percent noise | TV filtered | EFG filtered |

These images are related to a baseline image constructed as follows. The baseline image is flat on a centered square where it assumes its maximum value, it is flat in the background where it assumes its minimum value, and otherwise it slopes linearly at the left between these two values. The left-most image shown is the baseline image plus 50 percent noise, i.e., the noise is normally distributed with zero mean and standard deviation equal to half the range of image values. The center image shown is the result of applying TV filtering to the noisy image; see [Rudin et al., 1992]. Note the classic staircasing in the region where the baseline image is linearly sloped. The right-most image shown is the result of applying an edge-flat-grey (EFG) filtering approach to the noisy image; see [Ito & Kunisch, 1999] and [Keeling, 2003]. Here, edge, flat, and grey scales, i.e., scales where the gradient is large, small, and intermediate, respectively, are clearly captured well.

The TV and EFG filters are defined in terms of minimizers for the functional:

for which the necessary optimality condition is expressed in the steady state for the following descent evolution equation:
(1)

This diffusion process can be further elucidated by the orthogonal decomposition:
(2)

where
(3)
*n* denotes the derivative in the direction of grad *u*. Thus, with respect to level sets of *u*, the first term captures normal diffusion and the second captures tangential diffusion. Here, *g*(*s*) = *f*'(*s*)/*s*, and *f*''(*s*) is globally smooth and non-negative. Specifically:

For simplicity, suppose
(4)
*c* = 1 and observe that the EFG penalty function takes the form on [0,1]:

Thus, the EFG penalty corresponds roughly to an
(5)
*L _{p}* penalty on |grad

Also, shown on the right is a plot of *f*''(*s*) on the interval [0,1]. Note further that *p*(0) = *p*(1) = 1 and that *f*'''(*s*) is bounded. In light of Eqn. (3), the plot of *f*''(*s*) illustrates that diffusion normal to level sets of *u* is strictly positive for a range of gradient values, but drops to zero for sufficiently large and for vanishing gradients.

The numerical implementation of these filters involves discretization of the evolution equation in Eqn. (2). Specifically, the spatial discretization is built from Cartesian central differences among cells, where the nonlinear diffusivity *g*(|grad *u*|) is evaluated on cell interfaces. Also, for the TV filter, the function *g*(*s*) = 1/*s* is treated numerically as *g*(max{*m*,*s*}) for *m* equal to machine zero. The temporal discretization is achieved with implicit Euler time stepping, where the time stepping equations are solved with a Jacobi preconditioned conjugate gradient scheme. Full details and given in [Keeling, 2003].

Multiscale edge enhancement by nonlinear anisotropic diffusion.

highly blurred | PM filtered | BFB filtered |

The left-most image shown consists in a sum of six functions of the form *a*(1+*bx*^{2}+*cy*^{2})^{-d} chosen to create ever weaker edges while approaching the center of the image. The center image shown is the result of applying a routinely used Perona-Malik (PM) filter to the blurred image; see [Perona & Malik, 1987] and [Gerig et al., 1992]. Note that for the particular PM contrast parameter chosen here, the innermost edge is flattened, the next adjacent edge is sharpened, and other edges are essentially unaffected. The right-most image shown is the result of applying a balanced-forward-backward (BFB) diffusion filter to the blurred image; see [Keeling & Stollberger] for details. Here, all edges are sharpened equally well.

The PM and BFB filters are defined in terms of the nonlinear anisotropic diffusion equation:

where the diffusivities are given by:
(6)

Comparing the two evolution equations in Eqns. (2) and (6) reveals the variational connection with
(7)
*g*(*s*) = *f*'(*s*)/*s* and:

Thus, the variational penalty function for PM is first convex and then concave, while for BFB it is globally concave. The diffusion decomposition in Eqn. (3) shows that both filters create forward diffusion tangent to level sets of
(8)
*u*. However, for PM the diffusion normal to level sets is first forward and then backward. For BFB, not only is normal diffusion always backward, but its coefficient *f*'' has the same magnitude as that for the tangential diffusion *g*. This balance is the basis of the multiscale edge enhancement capability of the balanced-forward-backward diffusion filter.

The numerical implementation of these filters involves the discretization of the evolution equation shown in Eqn. (6). Specifically, the spatial discretization is built from diagonal central differences among cells, where the nonlinear diffusivity *g*(|grad *u*|) is evaluated at cell corners. Also, for the BFB filter, the function *g*(*s*) = 1/*s*^{2} is treated numerically as *g*(max{*m*,*s*}) for *m* equal to machine zero. For the PM filter, temporal discretization is accomplished with explicit Euler time stepping. For the BFB filter, temporal discretization is achieved with implicit Euler time stepping, where the time stepping equations are solved only approximately with a Jacobi preconditioned conjugate gradient scheme. Full details are given in [Keeling & Stollberger, 2000].

Application to contrast enhancement of medical images.

abdominal MRI | Gauss filtered | Gauss-TV filtered |

The left-most image shown is a moderately noisy MRI of the abdomen. The center image shown is the result of applying a routinely used Gauss filter to the MRI. Note the evident blur in the Gauss filtered image. The right-most image shown is the result of applying a Gauss-TV diffusion filter to the MRI. Here, there is a high contrast level along with a marked reduction in the noise level.

The Gauss filtered image was computed pixelwise merely by performing a weighted average of values in neighboring pixels from the MRI; see [Gonzalez & Woods]. On the other hand, the Gauss-TV filtered image was obtained using the diffusion approach in Eqn. (2) with the diffusivity:

which has roughly the same bell shape as that associated with Perona-Malik diffusivities:
(9)

However, here the diffusivity *g*(*s*) = 1 corresponds to Gaussian filtering in the sense that the image evolves according to the linear heat equation. Also, the diffusivity *g*(*s*) = *c*/*s* corresponds to TV filtering as seen earlier in the comparison with EFG filters. The combination of the Gaussian and TV filtering in Eqn. (9) gives:

Therefore, because of Eqn. (3), diffusion normal to level sets is never backward as seen with the edge enhancing filters. As a result, the diffusivity in Eqn. (9) is edge-preserving for edges with |grad
(10)
*u*| > *c*. Also, it is isotropically dissipative for noise with |grad *u*| < *c*.

In terms of EFG filters, the Gauss-TV filter treats images as consisting of two scales: edge and grey. This combination is quite effective for contrast enhancement of medical images, and can be particularly useful to increase lesion conspicuity for example. On the other hand, the additional treatment of flat scales with EFG filters has been found to excessively flatten small slope regions in these images. Specifically, this flattening effect causes a loss of a natural three dimensional appearance. On the other hand, the flattening resulting from EFG, PM, and BFB filters can be used to advantage as a preprocessor for a segmentation procedure. For further details, see [Keeling, 2003].

Application to magnetic resonance diffusion tensor imaging.

FA brain map | Gauss filtered | Gauss-TV filtered |

Magnetic resonance diffusion tensor imaging is a procedure used for mapping the diffusion properties of water in tissues; see [Bammer, Augustin et al., 2000]. Such diffusion is considerably anisotropic and is influenced by the underlying tissue organization and damage thereto. As a result, this imaging technique has naturally advanced into everyday routine clinical applications. The work shown here is aimed at the contrast enhancement and noise reduction in such images, using the Gauss-TV approach outlined above.

The water diffusion properties of interest can be detailed as follows. The flux of water molecules is described locally by * J* = -

(11) |

In particular, imaging the rotational invariants of *D* has been used to advantage for diagnostic purposes. For instance, the trace of *D* provides an orientationally averaged apparent diffusivity:

and is used to increase lesion conspicuity by suppressing nearby anisotropic effects. On the other hand, when anisotropic effects are of primary concern, as in the evaluation of nerve fiber tracts, the fractional anisotropy (FA) can be imaged:
(12)

It is this fractional anisotropy which is mapped in the three images appearing above.
(13)

The left-most image shown above is a noisy FA map of the brain, obtained using unfiltered diffusion-weighted images, from which the diffusion tensor images were derived and used in Eqn. (13). The next two images are the results of using the filters discussed in detail earlier for the abdominal MRI. Specifically, the center image shown is the result of applying the Gauss filter to the diffusion tensor images before computing the FA map. As before, note the evident blur in the Gauss filtered image. The right-most image shown is the result of applying the Gauss-TV diffusion filter to the diffusion tensor images before computing the FA map. Again, there is a high contrast level along with a marked reduction in the noise level. See [Keeling, Bammer, Fazekas & Stollberger, 2000].

Restoration of images corrupted by noise and background variations.

nonuniform background | uniform background |

The purpose of this project is to investigate the enhancement of magnetic resonance images in which broad-band clinical information coexists not only with high gradient noise but also with a slowly varying background variation often seen in MR images. Specifically, this background variation occurs because of the difficulty in generating an ideal, uniform magnetic field for the production of an MR image; see [Hornak, 1996]. For instance, the left image shown above is an MRI of the pelvis. Here, the wide field of view frustrates the task of creating a uniform magnetic field over the entire image region. As a result, the nonuniform background is manifested with the image intensity much higher near the center than toward the periphery. On the other hand, the right image shown is an MRI obtained from an *in vitro* arthrosclerotic vessel. Here, the field of view is much narrower and there is less difficulty with generating a uniform magnetic field over the sample.

An accepted model for this background variation involves approximating a measured intensity *u*^{~} with a product *mu* where *m* > 0 is a smooth multiplier field adjusting for the nonhomogeneous magnetic field strength so that *u* corresponds to the intensity that would be measured in a uniform magnetic field; see [Pham & Prince, 1998]. In fact, the objective is to determine *m* and *u*, where *u*^{~}/*m* is free of the background variation appearing in *u*^{~}, and *u* is free of the noise still present in *u*^{~}/*m*.

To achieve this objective, the pair (*m*,*u*) is determined as a minimizer for the functional:

Here,
(14)
*f* is a convex function, as in Eqn. (9), used to penalize noisy oscillations in *u* without smoothing intensity discontinuities at tissue boundaries. The remaining penalty terms control the noisy oscillations and the curvature accumulations in *m*, depending upon regularization parameters *a* > 0 and *b* > 0, respectively. The optimality system for *J*(*m*,*u*) has been derived and novel numerical techniques have been formulated for the determination of *m* and *u*. For further details, see [Keeling & Bammer, 2004] and [Keeling & Haase, 2007].

Phase unwrapping by regularized phase construction in a region growing setting.

wrapped noisy phase | unwrapped noisy phase |

The raw data measured for an MR image are understood in terms of both a magnitude and a phase of a complex signal on the basis of the following model; see [Hornak, 1996]. To create the MR signal, a uniform magnetic field is first applied to align the magnetization of the hydogen proton spins in the body. Then, a sufficient secondary magnetic field pulse is applied to excite the protons at a resonant frequency so that the local net magnetization is pushed away from the alignment axis and into the orthogonal plane as it begins to precess about the axis at the resonant frequency. The magnetization precessing in this plane is modeled as a complex magnetization, with real and imaginary parts corresponding to projections onto orthogonal axes therein. The magnitude of this complex magnetization is related to the local hydrogen density, the property normally imaged with MRI. On the other hand, the phase of the magnetization can be related to local conditions which alter the local resonant frequency of hydrogen protons. For instance, the hydrogen protons in water have a slightly different resonant frequency than the hydrogen protons in fat.

However, given a noisy complex magnetization *m*^{~} = *r*^{~}exp(*ip*^{~}), determining the underlying phase *p*^{~} is nontrivial; see [Song, 1995] and [Liang, 1996]. For instance, the derived phase tan^{-1}(Im{*m*^{~}}/Re{*m*^{~}}) is necessarily wrapped into the range (-Pi/2,Pi/2]. As a result, phase jumps of Pi can appear as seen in the wrapped phase image shown above on the left but not in the unwrapped phase image on the right. This phase map arises in the course of using the Dixon method to distinguish fat and water around the knee; see [Flevaris, 2000]. A detailed examination of these images shows that the presence of image noise can frustrate the computation of an accurately unwrapped phase. Furthermore, steep gradients in the phase must evidently be allowed at tissue boundaries. Therefore, it is required not only to avoid phase wrappings but also to avoid smoothing natural discontinuities while removing noise to obtain an appropriately filtered phase *p* ~= *p*^{~}.

It is proposed that the objectives can be achieved with the following natural variational approach:

Here,
(15)
*r* is an estimate of the signal magnitude. Also, *f*_{1} and *f*_{2} are convex functions, as in Eqn. (9), designed to penalize oscillations in *r* and *p* while allowing them to exhibit natural discontinuities. Preliminary computations on model problems demonstrate that this variational approach can be performed in a region-growing fashion to achieve highly accurate phase distributions.

Well-posedness analysis of new image enhancement schemes.

This area of research is aimed most immediately at establishing well-posedness results for the diffusion filters in Eqns. (2) and (6), where the diffusivity is unbounded. For example, Eqns. (4) and (7) show that the TV, EFG, and BFB diffusivities are unbounded for vanishingly small gradients. As stated earlier, such an unbounded diffusivity *g*(*s*) is treated numerically as *g*(max{*m*,*s*}) where *m* is equal to machine zero. Nevertheless, computational findings demonstrate that results are insensitive to variations in a sufficiently small *m*. This observation suggests the following continuous dependence property.

First, let * U*(

(16) |

Observations about the desired continuous dependence property include the following. First, an *m* independent *C*^{1} bound on the solutions to Eqn. (16) provide subsequential convergence in *C*^{0} by compactness. The challenge is to obtain a unique characterization of subsequential limits. Studies of the fixed points of:

show that there exist examples where, for
(17)
*m* = 0, it is possible that this equation can only be solved nonuniquely by discontinuous functions. However, the corresponding limit of fixed points as *m* -> 0 is unique and continuous. Thus, the characterization of **U**^{*} must not involve equations of the form Eqn. (16) or (17).

Image registration by optical flow with maximal rigidity.

breast MR at t_{1} |
breast MR at t_{2} |

The diagnostic use of medical image sets in a clinical setting implicitly requires a point by point correspondence between the same tissue sites in separate images. For example, the two images shown above were taken before and after the injection of a contrast agent during a mammography examination. To monitor the uptake of the contrast agent at specific, potentially pathological sites, it is necessary to track the site of interest through the tissue motion occuring during the examination. Thus, what is needed finally is an explicit coordinate transformation that will map any point in one image to its corresponding point in the other. With such a mapping, images are said to be registered. See [Thirion, 1998] and [Rueckert et al., 1999].

In this work, two given images *I*_{0}, *I*_{1}: *P* (in * R^{N}*) ->

(18) |

(19) |

The registration given by (*I*,* F*) is determined by minimizing the functional:

(20) |

Kernel estimation, for

perfusion imaging, and

tracer exchange kinetics imaging,

by discontinuity-preserving spatio-temporal regularization.

Pathophysiological properties of breast tumors are evaluated clinically from an examination in which the tumor is imaged repeatedly to measure the uptake of an injected contrast agent. The amount of injected tracer and the amount of absorbed tracer are related by a convolution whose kernel provides the desired diagnostic information. Specifically, let *C*_{a}(*t*) be the tracer concentration at time *t* in an artery feeding the tumor. Also, let *C*_{r}(* x*,

(21) |

Note that the measurement *C*_{r}(* x*,

(22) |

Coil sensitivity estimation for high-speed parallel magnetic resonance imaging strategies.

sensitivity weighted | coil sensitivity |

The continuing development of magnetic resonance imaging has driven the need for increased image acquisition speed, for instance, to increase temporal resolution. Sensitivity encoding (SENSE) is a technique for reducing the scan time by employing multiple receiver coils which measure the magnetic signal in parallel at a reduced sampling rate; see [Bammer, Stollberger et al., 2000]. The sampling rate is reduced by a factor corresponding to the number of coils. Also, since the reduced sampling rate is below that required to capture the full information content in the signal, each measurement is aliased.

Specifically, suppose the sampling rate is only reduced in the phase-encoding direction, which corresponds to the vertical or say *y* direction in the image after Fourier transforming the signal; see [Hornak, 1996] for details. Suppose there are *N* coils each sampling in the phase-encode direction at a rate reduced by a factor of *N*. Now let M_{i} be the image that would be measured by the *i*th coil with full sampling, and let *M _{i}* be the image measured with reduced sampling. If the vertical dimension of the image M

(23) |

(24) |

(25) |

To achieve the objective, *C* is determined as a minimizer for the following cost:

Here,
(26)
*a* and *b* are positive weights for penalty terms respectively controlling the noisy oscillations and the curvature accumulations in *C*. Furthermore, these parameters control the trade-off between their regularizations and the reduction of the residual |*CM*_{body} - *M*_{surf}|. The optimality system for *J*(*C*) has been derived and novel numerical techniques have been formulated for the determination of *C*. For further details, see [Keeling & Bammer, 2004] and [Keeling & Haase, 2007]

References.

- R. Bammer, M. Augustin, T. Seifert, S. Strasser-Fuchs, R. Stollberger, H.-P. Hartung and F. Fazekas,
*High Resolution Diffusion Tensor Imaging using Interleaved Echo-Planar Imaging*, SPIE International Symposium on Medical Imaging, San Diego, CA, February 12-18, 2000. - R. Bammer, R. Stollberger, M. Augustin, T. Seifert, S. Strasser-Fuchs, H.P. Hartung and F. Fezekas,
*Parallel Imaging Strategies for High-Speed Magnetic Resonance Diffusion Imaging*, SPIE International Symposium on Medical Imaging, San Diego, CA, February 12-18, 2000. - E. Flevaris,
*Phase Unwrapping füer die Zwei-Punkt-Dixon-Methode bei Gradientenechosequenzen*, Diplomarbeit, Institut fuer Elektro- und Biomedizinische Technik, Technische Universitaet Graz, 2000. - G. Gerig, O. Kuebler, R. Kikinis and F.A. Jolesz,
*Nonlinear Anisotropic Filtering of MRI Data*, IEEE Transactions on Medical Imaging, Vol. 11, No. 2, pp. 221 - 232, June 1992. - R. Gonzalez and R. Woods,
*Digital Image Processing*, Addison-Wesley, Reading, MA, 1992. - J.P. Hornak,
*The Basics of MRI*,`http://www.cis.rit.edu/htbooks/mri`, 1996. - K. Ito and K. Kunisch,
*BV-type Regularization Methods for Convoluted Objects with Edge-Flat-Grey Scales*,*Inverse Problems*, Vol. 65, pp. 273 - 298, 2000. - S.L. Keeling,
*Total Variation Based Convex Filters for Medical Imaging*, Applied Mathematics and Computation, Vol. 139, pp. 101-119, 2003. - S.L. Keeling and R. Bammer,
*A Variational Approach to Magnetic Resonance Coil Sensitivity Estimation*, Applied Mathematics and Computation, Vol. 158, No. 2, pp. 53-82, 2004. - S.L. Keeling, R. Bammer, F. Fazekas and R. Stollberger,
*Total Variation Denoising for Improved Diffusion Tensor Calculation*, ISMRM 2000, Denver, CO, April 1-7, 2000. - S.L. Keeling, R. Bammer and R. Stollberger,
*Revision of the Theory of Tracer Transport and the Convolution Model of Dynamic Contrast Enhanced Magnetic Resonance Imaging*, Journal of Mathematical Biology, Vol. 55, pp. 389-411, 2007. - S.L. Keeling and G. Haase,
*Geometric Multigrid for High-Order Regularizations of Early Vision Problems*, Applied Mathematics and Computation, Vol. 184, pp. 536-556, 2007. - S.L. Keeling, T. Kogler and R. Stollberger,
*Deconvolution for DCE-MRI Using an Exponential Approximation Basis*, Medical Image Analysis, Vol. 13, pp. 80-90, 2009. - S.L. Keeling and W. Ring,
*Medical Image Registration and Interpolation by Optical Flow with Maximal Rigidity*, Journal of Mathematical Imaging and Vision, Vol. 23, No. 1, pp. 47-65, July 2005. - S.L. Keeling and R. Stollberger,
*Nonlinear Anisotropic Diffusion Filters for Multiscale Edge Enhancement*, Inverse Problems, Vol. 18, pp. 175-190, January 2002. - Z. Liang,
*A Model-Based Method for Phase Unwrapping*, IEEE Trans. Med. Imag., Vol. 15, No. 6, pp. 667 - 676, December 1996. - L. Ostergaard, R.M. Weisskoff, D.A. Chesler, C. Glydensted and B.R. Rosen,
*High Resolution Measurement of Cerebral Blood Flow using Intravascular Tracer Bolus Passages. Part I: Mathematical Approach and Statistical Analysis*, Magn. Reson. Med., Vol. 36, pp. 715 - 725, 1996. - P. Perona and J. Malik,
*Scale Space and Edge Detection Using Anisotropic Diffusion*, Proc. IEEE Comp. Soc. Workshop on Computer Vision, Miami Beach, November 30 - December 2, 1987, IEEE Computer Society Press, Washington, pp. 16 - 22, 1987. - D.L. Pham and J.L. Prince,
*An Adaptive Fuzzy C-Means Algorithm for Image Segmentation in the Presence of Intensity Inhomogeneities*, Proceedings of the SPIE Conference on Image Processing, San Diego, CA, February 1998. - D. Rueckert, L.I. Sonoda, E. Denton, S. Rankin, C. Hayes, M.O. Leach, D. Hill and D.J. Hawkes,
*Evaluation of non-rigid registration using free-form deformations for breast MR images*, Medical Image Understanding and Analysis 99, Oxford, UK, July 19-20, 1999. - L.I. Rudin, S. Osher and E. Fatemi,
*Total Variation Based Noise Removal Algorithms*, Physica D, Vol. 60, pp. 259 - 269, 1992. - S.M. Song, S. Napel, N.J. Pelc, and G.H. Glover,
*Phase Unwrapping of MR Phase Images Using Poisson Equation*, IEEE Trans. Image Proc., Vol. 4, No. 5, pp. 667 - 676, May 1995. - J.-P. Thirion,
*Image Matching as a Diffusion Process: An Analogy with Maxwell's Demons*, Medical Image Analysis, Vol. 2, No. 3, pp. 243 - 260, 1998. - J. Weickert,
*Anisotropic Diffusion in Image Processing*, B. G. Teubner Stuttgart, 1998.