In earlier work, we proposed a super-resolution (SR) method that required the availability of two low resolution (LR) sequences corresponding to two different sampling rates, where images from one sequence were used as a basis to represent the polyphase components (PPCs) of the high resolution (HR) image, while the other LR sequences provided the reference LR image (to be super-resolved). The (simple) algorithm implemented by Salem and Yagle is only applicable when the scene is static. In this paper, we recast our approach to SR as a two-stage example-based algorithm to process dynamic scenes. We employ feature selection to create, from the LR frames, local LR dictionaries to represent PPCs of HR patches. To enforce sparsity, we implement Gaussian generative models as an efficient alternative to L1-norm minimization. Estimation errors are further reduced using what we refer to as the anchors, which are based on the relationship between PPCs corresponding to different sampling rates. In the second stage, we revert to simple single frame SR (applied to each frame), using HR dictionaries extracted from the super-resolved sequence of the previous stage. The second stage is thus a reiteration of the sparsity coding scheme, using only one LR sequence, and without involving PPCs. The ability of the modified algorithm to super-resolve challenging LR sequences reintroduces sampling rate diversity as a prerequisite of robust multiframe SR.

Multiframe super-resolution is the problem of reconstructing a single high-resolution (HR) image from several low-resolution (LR) versions of it. We assume that the original HR image undergoes different linear transforms, where each transform can be approximated as a set of linear shift-invariant transforms over different subregions of the HR image. The linearly transformed versions of the HR image are then downsampled, resulting in different LR images. Under the assumption of linearity, these LR images can form a basis that spans the set of the polyphase components (PPCs) of the HR image. We propose sampling rate diversity, where a secondary LR image, acquired by a secondary sensor of different (lower) sampling rate, is used as a reference to make known portions (subpolyphase components) of the PPCs of the reconstructed HR image. This setup allows for non-parametric reconstruction of the PPCs, where no knowledge of the underlying transforms is required, by solving for the expansion coefficients of the PPCs, in terms of the LR basis.

The problem of reconstructing an image from irregular samples of its 2-D DTFT arises in synthetic aperture radar (SAR), magnetic resonance imaging (MRI), limited angle tomography, and 2-D filter design. Since there is no 2-D Lagrange interpolation, sufficient conditions for the uniqueness and conditioning of the reconstruction problem are both not apparent. The Good-Thomas FFT is used to unwrap the 2-D problem into a 1-D problem, from which uniqueness results and, more importantly, insights into the problem conditioning are available. We propose the variance of distances between adjacent frequency locations as a sensitivity measure, which aids in determining a well-conditioned configuration of frequency values. The sensitivity measure is analyzed on its accuracy of estimating the conditioning and on its computational speed. The image is then reconstructed by solving the problem using the conjugate gradient method.

Toeplitz-like matrices are matrices that are Toeplitz, block Toeplitz with Toeplitz blocks, and mosaic Toeplitz (Toeplitz blocks of different sizes). Toeplitz-like systems of equations arise in 2-D interpolation, 2-D linear prediction, and 2-D least-squares deconvolution. In this paper, we use the Woodbury formula to reformulate a Toeplitz-like system of equations into a circulant-plus-diagonal system, and use the discrete Fourier transform to transform this into a banded system which can be solved quickly and which requires little storage. We propose the use of this as a preconditioner for Toeplitz-like systems, as an alternative to the circulant-block-circulant preconditioner commonly used for Toeplitz-block-Toeplitz systems. Using condition number as a figure-of-merit, this preconditioner seems to work much better for Toeplitz matrices than the usual circulant preconditioner.

The 2D discrete phase retrieval problem is to reconstruct an image defined at integer coordinates and having a known finite spatial extent from the magnitude of its discrete Fourier transform. Most methods for solving this problem are iterative but not POCS, and they tend to stagnate. Recently, we developed a new approach that unwrapped the 2D problem into a 1D problem with bands of zeros in it, using the Good-Thomas FFT (see Petroudi, S. and Yagle, A.E., Proc. ICASSP, 2000). However, this approach reconstructed the even part of the image much better than the odd part, and it was sensitive to the zero locations. This paper presents a modification of this approach. New features include: (1) an overdetermined problem less sensitive to the zero locations; (2) the solution of a Toeplitz-block-Toeplitz-plus-Hankel-block-Hankel linear system; and (3) details of characteristics of images for which the approach works best.

A Toeplitz-block-Toeplitz (TBT) matrix is block Toeplitz with Toeplitz blocks. TBT systems of equations arise in 2D interpolation, 2D linear prediction and 2D least-squares deconvolution problems. Although the doubly Toeplitz structure should be exploitable in a fast algorithm, existing fast algorithms only exploit the block Toeplitz structure, not the Toeplitz structure of the blocks. Iterative algorithms can employ the 2D FFT (fast Fourier transform), but usually take thousands of iterations to converge. We develop a new fast algorithm that assumes a smoothness constraint on the matrix entries. For an M/sup 2//spl times/M/sup 2/ TBT matrix with M M/spl times/M Toeplitz blocks along each edge, the algorithm requires only O(6M/sup 3/) operations to solve an M/sup 2//spl times/M/sup 2/ linear system of equations; parallel computing on 2M processors can be performed on the algorithm as given. Two examples show the operation and performance of the algorithm.

The convergence speed of the steepest descent (SD) adaptive algorithm is determined by the eigenvalue spread of the correlation matrix. When the correlation matrix is singular, the algorithm will take a long time to converge. The problem can be regularized by adding a small constant to the diagonal, but this comes at the price of introducing bias. This paper introduces a new adaptive algorithm that uses a family of steepest descent iterations which converge quickly while making the bias arbitrarily small. A numerical example discusses in some detail the operation of the new algorithm.

The 2-D discrete phase retrieval problem is to reconstruct an image defined at integer coordinates and having known finite spatial extent from the magnitude of its discrete Fourier transform. Most methods for solving this problem are iterative but not POCS, and they tend to stagnate. We have shown that a subband decomposition can be used to divide a large phase retrieval problem into smaller problems defined over subbands. However, these numerous smaller phase retrieval problems must be solved efficiently. This paper uses the Good-Thomas FFT to map the 2-D problem into a 1-D problem of reconstructing the signal formed by concatenating rows of the image with zero bands. The z-transform zeros of this 1-D signal are close enough to the unit circle to be found from local minima of the (known) DFT magnitude. The 1-D signal can then be reconstructed in closed form.

Computation of the inverse filter s(n) to a signal r(n) (so the convolution of r(n) and s(n) is an impulse) can require a surprisingly large amount of computation. The obvious approach of taking the reciprocal of the discrete Fourier transform (DFT) of r(n) may result in aliasing unless a large order DFT is used. O(Nlog(N)) operations are required, where s(N) is negligible. We develop a new algorithm which requires O(M/sup 2/log(N)) operations, where M is the length of r(n), and so is faster if N>M/sup 2/, as is often the case. The algorithm uses successive polynomial division to identify vectors orthogonal to the stable part of 1/R(z); it requires the solution of a resultant-like Toeplitz-plus-Hankel linear system of equations, which requires O(M/sup 2/) operations. The inverse filter can then be used to compute the complex cepstrum and minimum phase spectral factor in O(M/sup 2/) operations.

Cardiac muscle tissue consists of multiple, electrically-active layers, whose electrical potentials can be reconstructed from surface measurements as a volume conductor inverse problem. These are useful in diagnosis following myocardial infarction. The measurements include an unknown point-spread function, which must be blindly deconvolved from the bioelectrical potentials. We present a new divide-and-conquer approach to this problem: (1) unwrap the 2-D problem into a huge 1-D problem; (2) perform a multiresolution (subband) decomposition of the huge 1-D problem into small 1-D problems; (3) solve the small 1-D problems by finding the null vector of a Toeplitz-plus-Hankel matrix. Depending on the noise level, the polynomial Euclidian algorithm or structured total least squares may be used.

The problem of reconstructing an image from irregular (not on a cartesian grid) samples of its 2-D DTFT arises in synthetic aperture radar (SAR) and magnetic resonance imaging (MRI), in which the 2-D DTFT is known only on part of a polar raster. It also arises in limited angle tomography, in which the 2-D DTFT is known in a bowtie region. Reconstruction requires either nearest-neighbor interpolation or solution of a large linear system of equations, both of which are computationally intensive and can lead to errors due to poor conditioning of the problem. An explicit formula does not seem to exist, since there is no 2-D Lagrange interpolation formula. This paper uses the Good-Thomas FFT to unwrap the 2-D problem into a 1-D problem, to which the 1-D Lagrange interpolation formula can be applied, either directly or recursively (the latter is much faster). This approach results in a sufficient condition to ensure a unique reconstruction.

The 2-D blind deconvolution problem is to reconstruct an image having known finite spatial extent from its 2-D convolution with an also-unknown point-spread function. This is significantly more difficult than the typical image restoration problem of deconvolving a known blurring function. Many methods for solving this problem are iterative but not POCS, and they tend to stagnate. We present a completely novel approach that assumes the unknown point-spread function varies slowly in the 2-D Z-transform domain, as would be the case if the point-spread function had small spatial extent. Sampling the problem along one axis in the 2-D Z-transform domain, and assuming the image transform varies between samples while the point-spread function transform does not, results in an approximate polynomial greatest common divisor problem. This can be solved by finding the eigenvector associated with the minimum eigenvalue of a resultant matrix. Repeating for several samples allows the point-spread function to be reconstructed using Lagrange interpolation.

We present a finite-state memory model for parametric system identification of the lip seal friction process in a hydraulic actuator. The performance of the finite-state memory model is compared with two Hammerstein type models using experimental results.

We derive two new algorithms for reconstructing a discrete-time 1-D signal from the phase of its discrete-time Fourier transform (DTFT) at irregular frequencies. Previous algorithms for this problem have either required the computation of a matrix nullspace, requiring O(N/sup 3/) computations, or have been iterative in nature; for the latter, the irregularity of the frequency samples precludes use of the fast Fourier transform. Our first algorithm requires only O(N/sup 2/) computations (O(N log/sup 3/ N) asymptotically). In the special case of equally-spaced frequency samples, it is related to a previous algorithm. The second algorithm is recursive-at each recursion a meaningful magnitude retrieval problem is solved. This is useful for updating a solution; it also allows checking of the result at each recursion, avoiding any errors due to computational roundoff error and ill-conditioning of the problem.

The phase retrieval problem arises when a signal must be reconstructed from only the magnitude of its Fourier transform; if the phase information were also available, the signal could simply be synthesized using the inverse Fourier transform. In continuous phase retrieval, most previous solutions rely on discretizing the problem and then employing an iterative algorithm. We avoid this approximation by using wavelet expansions to transform this uncountably infinite problem into a linear system of equations. The wavelet bases permit a solution by incorporating a priori signal information and they provide a structured system of equations which results in a fast algorithm. Our solutions obviate the stagnation problems associated with iterative algorithms, they are computationally simpler and more stable than previous non-iterative algorithms, and they can accommodate noisy Fourier magnitude information. This paper develops our 1-D continuous, non-minimum phase retrieval algorithm and illustrates its effectiveness with numerical examples.

We present two Hammerstein-type models for parametric system identification of the lip seal friction process in a hydraulic actuator. Adaptive algorithms with least squares criteria are derived, and the performances of the two models are evaluated using experimental results.

We consider the problem of computing the maximum likelihood estimates of the reflection coefficients of a discrete 1-D layered medium from noisy observations of its impulse reflection response. We have side information in that a known subset of the reflection coefficients are known to be zero; this knowledge could come from either a priori knowledge of a homogeneous subregion inside the scattering medium, or from a thresholding operation in which noisy reconstructed reflection coefficients with absolute values below a threshold are known to be zero. Our procedure is simple, noniterative, and requires only solutions of systems of linear equations. Numerical examples are provided which demonstrate not only the operation of the algorithm, but also that the side information improves the reconstruction of unconstrained reflection coefficients as well as constrained ones, due to the nonlinearity of the problem.

Sahiner and Yagle (1992, 1993, 1995) developed a procedure for image reconstruction from projections under wavelet-domain image constraints (WC). We constrain the wavelet transform (WT) of the noisy image to zero in WT regions determined by thresholding the WT of the noisy image. The least-squares image satisfying these constraints shows improved reconstruction not just in the constrained region, but over the entire image. This amounts to image-constrained localized low pass filtering. We now present new results on iterative implementation of this technique. At each iteration we use the image from the previous iteration as the source of WC for the next iteration. This can be viewed as iterative spatially-varying image-constrained Wiener filtering. At the third iteration the RMS error is reduced by more than twice as much as in our previous results, while preserving high-resolution image features.

We present an explicit formula for reconstructing a finite-support object defined on a lattice of points and taking on integer values from a finite number of its discrete projections over a limited range of angles. We make extensive use of the discrete Fourier transform in doing so. Our approach computes the object sample values directly as a linear combination of the projections sample values. The well-known ill-posedness of the limited angle tomography problem manifests itself in some very large coefficients in these linear combinations; these coefficients (which are computed off-line) provide a direct sensitivity measure of the reconstruction samples to the projections samples. The discrete nature of the problem implies that the projections must also take on integer values; this means noise can be rejected. This makes the formula practical.

Discrete tomography is the problem of reconstructing a binary image defined on a discrete lattice of points from its projections at only a few angles. It has applications in X-ray crystallography, in which the projections are the number of atoms in the crystal along a given line, and nondestructive testing. The 2-D version of this problem is fairly well understood, and several algorithms for solving it are known, most of which involve discrete mathematics or network theory. However, the 3-D problem is much harder to solve. This paper shows how the problem can be recast in a purely algebraic form. This results in: (1) new insight into the number of projection angles needed for an almost surely unique solution; (2) non-obvious dependencies in projection data; and (3) new algorithms for solving.

The 2-D discrete phase retrieval problem is to reconstruct a discrete-time signal whose support is known and compact from the magnitude of its discrete Fourier transform. We show how a subband decomposition of this problem can be performed using filter banks. The result is a series of smaller 2-D phase retrieval problems whose solutions are the phases of the original problem in different frequency bands. The 2-D phase retrieval problem is first mapped to a 1-D phase retrieval problem and the subband decomposition is applied to this problem to decompose it into smaller 1-D problems, each of which can be viewed as a smaller 2-D problem. This is more flexible than performing the decomposition directly in 2-D. The filter used is a truncated Gaussian. While these results can also perform subband decomposition of any 1-D phase retrieval problem, the general lack of uniqueness for 1-D phase retrieval creates problems in reassembling the solutions to the smaller 1-D problems.

We apply a maximum-likelihood stochastic system identification technique to the problem of identifying parameters of a model for friction in a hydraulic actuator. This is the first application of stochastic modelling and identification techniques to the problem of modelling friction, which is a very complex physical process. The deterministic part of the model characterizes energy dissipation mechanisms of friction and the associated transient responses. The stochastic part of the (ARMAX) model characterizes unmodelled dynamics due to process disturbances and measurement noise. The model and identification algorithm are validated by comparison with experimental data.

The Krein integral equation of one-dimensional inverse scattering, which has a symmetric Toeplitz kernel, is transformed using wavelets into a "block-slanted Toeplitz" system of equations. The kernel of the integral equation does not satisfy the Calderon-Zygmund conditions and as a result, application of the wavelet transform to the integral equation does not yield a sparse system matrix. There is therefore a need for a fast algorithm which directly exploits the (symmetric block-slanted-Toeplitz) structure of the system matrix and does not rely on sparsity. The first such O(N/sup 2/) algorithm is presented.

The discrete phase retrieval problem is to reconstruct a discrete-time signal whose support is known and compact from the magnitude of its discrete Fourier transform. We assume knowledge of some values of the signal, e.g. bands of zeros or known values, and solve the problem by solving linear systems of equations. No rooting of polynomials or tracking zero curves of algebraic functions (both very expensive computationally and unstable numerically) is required. Not only is our method much simpler computationally than previous methods, but it also allows the use of total least squares type techniques to be applied to the linear systems of equations, so that noisy data can be handled.

We define two forms of stochastic tomography. In global tomography, the goal is to reconstruct an object from noisy observations of all of its projections. In region-of-interest (ROI) tomography, the goal is to reconstruct a small portion of an object (an ROI) from noisy observations of its projections densely sampled in and near the ROI and sparsely sampled away from the ROI. We solve both problems by expanding the object and its projections in a circular harmonic (Fourier) series in the angular variable so that the Radon transform becomes Abel transforms of integer orders applied to the harmonics. The algorithm has three major components. First, we fit state-space models to each order of Abel transform and thus represent the Radon transform operation as a parallel bank of systems, each of which computes the appropriate Abel transform of a circular harmonic. A variable transformation here allows either the global or ROI problem to be solved. Second, the object harmonics are modeled as a Brownian branch. This is a two-point boundary value system, which is Markovianized into a form suitable for the Kalman filter. Finally, a parallel bank of Kalman smoothing filters independently estimates each circular harmonic from the noisy projection data. Numerical examples illustrate the proposed procedure.

Provides a comparison between one-sided linear prediction (OSP) and two-sided linear prediction (TSP) with respect to prediction error, relationships to AR modeling and to two-sided AR modeling, and the application to time series interpolation, linear-phase filter design, and spectral estimation. New contributions of the paper include: (1) proof that TSP produces smaller, non-white residuals than OSP, extending previous results; (2) specification of the frequency-domain error criterion minimized by TSP, and comparison with the analogous OSP criterion; (3) demonstration that TSP and two-sided AR modeling are different problems, unlike OSP; (4) interpretation of performance of TSP interference-rejection filters.

Combines several ideas, including nonuniform sampling and circular harmonic expansions, into a new procedure for reconstructing a small region of interest (ROI) of an image from a set of its projections that are densely sampled in the ROI and coarsely sampled outside the ROI. Specifically, the radial sampling density of both the projections and the reconstructed image decreases exponentially with increasing distance from the ROI. The problem and data are reminiscent of the local tomography problem; however, the present algorithm reconstructs the ROI of the image itself, not the filtered version of it obtained using local tomography. The new algorithm has the added advantages of speed (it can be implemented entirely using the FFT) and parallelizability (each image harmonic is independent).

We provide a novel approach to the design of fast algorithms for matrix multiplication. The operation of matrix multiplication is reformulated as a convolution, which is implemented using pseudo-number-theoretic transforms. Writing the convolution as multiplication of polynomials evaluated off the unit circle reduces the number of multiplications without producing any error, since the (integer) elements of the product matrix are known to be bounded. The new algorithms are somewhat analogous to the arbitrary precision approximation (APA) algorithms, but have the following advantages: (1) a simple design procedure is specified for them; (2) they do not suffer from roundoff error; and (3) reasons for their existence is clear. The new algorithms are also non-commutative, so that they may be applied recursively to block matrix multiplication. This work establishes a link between matrix multiplication and fast convolution algorithms and so opens another line of enquiry for the fast matrix multiplication problem. Some numerical examples illustrate the operation of the new proposed algorithms.

The 2-D blind deconvolution problem is to reconstruct an image (2-D signal) from the 2-D convolution of the image with another, unknown 2-D signal which represents some unknown distortion. The only known information is the result of the convolution and the fact that both 2-D signals have compact support. We solve the 2-D discrete blind deconvolution problem by partitioning it into a mostly-decoupled set of 1-D blind deconvolution problems. We define discrete and modulated Radon transforms to formulate two coupled 1-D problems, the solution to which then specifies solutions to the other decoupled 1-D problems. The latter may in turn be solved in parallel; however, using the solution to one problem as input to a neighboring problem reduces the computation significantly for serial computers. Unlike other exact 2-D blind deconvolution methods which rely on tracking zero curves of algebraic functions or equivalent operations, no continuous-function-based methods are used here. This makes the procedure more robust numerically. By coupling more than two initial problems, the method can be made to work on data with small amounts of noise.

Some authors have shown that wavelets can be used to localize the Radon transform, i.e., to reconstruct a small region of interest (ROI) from projections densely sampled in the ROI and sparsely sampled outside it. A different approach using Fourier expansions of the image and its projections shows that this localization is a property of the Radon transform itself, not just the wavelet representation of it, provided that a wavelet-reminiscent exponential sampling scheme is used. This suggests that the Radon transform naturally lends itself to a multiresolution representation. We present a new procedure for reconstructing the wavelet transforms of the angular harmonics of an image from the wavelet transforms of the angular harmonics of its projections. The relation is a convolution over scale, rather than translation, and preservation of compact support of the wavelet basis functions since a derivative-Hilbert-transform is not required. It generalizes a procedure using exponential radial sampling, and decouples over angular harmonics.

Motivated by the radioglaciology inverse problem, we present new algorithms for reconstructing a lossy, stratified dielectric from its impulsive plane wave reflection responses at two different angles of incidence. Novel features of these algorithms include: 1) a digital signal processing formulation that does not require discretization of continuous equations; 2) use of the asymmetric Levinson algorithm for rapid solution of the forward and inverse problems; 3) a novel use of layer stripping ideas, featuring iteration between the forward and inverse problems, with each iteration recursively determining another layer of the medium; and 4) another recursive algorithm for determining the bottom lossy half-space from reflection data only. Numerical examples illustrate the new algorithms on the reconstruction of a synthetic but realistic layered ice shelf from both noiseless and noisy radar reflection data.

The authors combine several ideas, including nonuniform sampling and circular harmonic expansions, into a new procedure for reconstructing a small region of interest (ROI) of an image from a set of its projections that are densely sampled in the ROI and coarsely sampled outside the ROI. Specifically, the radial sampling density of both the projections and the reconstructed image decreases exponentially with increasing distance from the ROI. The problem and data are reminiscent of the recently formulated local tomography problem; however, the authors' algorithm reconstructs the ROI of the image itself, not the filtered version of it obtained using local tomography. The new algorithm has the added advantages of speed (it can be implemented entirely using the FFT) and parallelizability (each image harmonic is computed independently). Numerical examples compare the new algorithm to filtered backprojection.

The paper presents a multidimensional nonlinear edge-preserving filter for restoration and enhancement of magnetic resonance images (MRI). The filter uses both interframe (parametric or temporal) and intraframe (spatial) information to filter the additive noise from an MRI scene sequence. It combines the approximate maximum likelihood (equivalently, least squares) estimate of the interframe pixels, using MRI signal models, with a trimmed spatial smoothing algorithm, using a Euclidean distance discriminator to preserve partial volume and edge information. (Partial volume information is generated from voxels containing a mixture of different tissues.) Since the filter's structure is parallel, its implementation on a parallel processing computer is straightforward. Details of the filter implementation for a sequence of four multiple spin-echo images is explained, and the effects of filter parameters (neighborhood size and threshold value) on the computation time and performance of the filter is discussed. The filter is applied to MRI simulation and brain studies, serving as a preprocessing procedure for the eigenimage filter. (The eigenimage filter generates a composite image in which a feature of interest is segmented from the surrounding interfering features.) It outperforms conventional pre and post-processing filters, including spatial smoothing, low-pass filtering with a Gaussian kernel, median filtering, and combined vector median with average filtering.

Low-pass filtering computed tomography (CT) images to reduce noise may smooth or modify image features which are very important to the physician. Image features are often more easily identified and processed in the time-frequency plane. The authors use time-frequency distributions for spatially varying filtering of noisy CT images, constraining time-frequency representation coefficients of the projection data or of the reconstructed image to be zero in certain regions of the time-frequency plane. The authors consider two different applications: 1) filtering the projection data and then performing image reconstruction; and 2) filtering the reconstructed image directly. Criteria minimized, subject to constraints, may be either a deterministic minimum weighted perturbation of the given projection data or a stochastic minimum mean-square error in colored Gaussian noise. Results show improvement over processing the image with a linear spatially invariant filter.

We present a new procedure for reconstructing the wavelet transforms of the angular harmonics of an image from the wavelet transforms of the angular harmonics of its projections. The relation is a convolution over scale, rather than over translation, and preservation of compact support of the wavelet basis functions after a derivative-Hilbert transform is not required. It generalizes a recently proposed procedure using exponential radial sampling and decouples over angular harmonics.>

The eigenimage filter generates a composite image in which a desired feature is segmented from interfering features. The signal-to-noise ratio (SNR) of the eigenimage equals its contrast-to-noise ratio (CNR) and is directly proportional to the dissimilarity between the desired and interfering features. Since image gray levels are analytical functions of magnetic resonance imaging (MRI) parameters, it is possible to maximize this dissimilarity by optimizing these parameters. For optimization, the authors consider four MRI pulse sequences: multiple spin-echo (MSE); spin-echo (SE); inversion recovery (IR); and gradient-echo (GE). The authors use the mathematical expressions for MRI signals along with intrinsic tissue parameters to express the objective function (normalized SNR of the eigenimage) in terms of MRI parameters. The objective function along with a set of diagnostic or instrumental constraints define a multidimensional nonlinear constrained optimization problem, which the authors solve by the fixed point approach. The optimization technique is demonstrated through its application to phantom and brain images. The authors show that the optimal pulse sequence parameters for a sequence of four MSE and one IR images almost doubles the smallest normalized SNR of the brain eigenimages, as compared to the conventional brain protocol.

In using filtered backprojection to compute the inverse Radon transform, the ramp filter amplifies noise. Spatially invariant noise filters reduce resolution. It is desirable to filter noise where projections have no local high-frequency components. Using the short-time Fourier transform, the authors apply a time-frequency mask filter that zeroes out projections where local signal energy is below a threshold. Results show improvement over reconstructions using spatially invariant smoothing filters.

In the filtered backprojection procedure for image reconstruction from projections, backprojection dominates the computation time. A simple algorithm that reduces the number of multiplications in linear interpolation and backprojection stage by 50%, with a small increase in the number of additions, is proposed. The algorithm performs the interpolation and backprojection of four views together. Examples of implementation are given and extension to interpolation of more than four views is discussed.

The generalized Landweber iteration with a variable shaping matrix is used to solve the large linear system of equations arising in the image reconstruction problem of emission tomography. The method is based on the property that once a spatial frequency image component is almost recovered within in in the generalized Landweber iteration, this component will still stay within in during subsequent iterations with a different shaping matrix, as long as this shaping matrix satisfies the convergence criterion for the component. Two different shaping matrices are used: the first recovers low-frequency image components; and the second may be used either to accelerate the reconstruction of high-frequency image components, or to attenuate these components to filter the image. The variable shaping matrix gives results similar to truncated inverse filtering, but requires much less computation and memory, since it does not rely on the singular value decomposition.

A new algorithm for adaptive filtering of autoregressive processes is presented. The new algorithm is based on a three-term recurrence that replaces, but is mathematically equivalent to, the lattice equations of the gradient adaptive lattice (GAL) algorithm of Griffiths (1977). This three-term recurrence requires only half as many multiplications per recursion as do the lattice equations. A single-stage convergence analysis of the new algorithm is presented, and some numerical examples comparing it to the GAL are given. The new algorithm performs similarly to the GAL, but requires only 2/3 as many multiplications.>

A procedure that speeds up convergence during the initial stage (the first 100 forward and backward projections) of Landweber-type algorithms, for iterative image reconstruction for positron emission tomography (PET), which include the Landweber, generalized Landweber, and steepest descent algorithms, is discussed. The procedure first identifies the singular vector associated with the maximum singular value of the PET system matrix, and then suppresses projection of the data on this singular vector after a single Landweber iteration. It is shown that typical PET system matrices have a significant gap between their two largest singular values; hence, this suppression allows larger gains in subsequent iterations, speeding up convergence by roughly a factor of three.

The performance of the eigenimage filter is compared with those of several other filters as applied to magnetic resonance image (MRI) scene sequences for image enhancement and segmentation. Comparisons are made with principal component analysis, matched, modified-matched, maximum contrast, target point, ratio, log-ratio, and angle image filters. Signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), segmentation of a desired feature (SDF), and correction for partial volume averaging effects (CPV) are used as performance measures. For comparison, analytical expressions for SNRs and CNRs of filtered images are derived, and CPV by a linear filter is studied. Properties of filters are illustrated through their applications to simulated and acquired MRI sequences of a phantom study and a clinical case; advantages and weaknesses are discussed. The conclusion is that the eigenimage filter is the optimal linear filter that achieves SDF and CPV simultaneously.

An isotropic random field is expanded into its circular harmonics. Computation of its Radon transform is equivalent to computation of the n/sup th/-order Abel transform of the n/sup th/ order circular harmonic. A state-space model is fitted to the Abel transform of each order and augmented with a state-space model describing the random field circular harmonic. The latter is derived using backward Markovianization of a two-point boundary value model. The tomographic problem of computing the inverse Radon transform is then solved by using a bank of Kalman filters to estimate each random field harmonic separately. Combining these gives the linear least-squares estimate of the random field. The authors also consider a simpler Wiener process model of the circular harmonics.>

An estimate is made of the number and locations of planar bioelectric sources in a horizontally layered volume conductor, from noisy surface measurements of electrical potential on a plane. The effects of the volume conductor are modeled by a medium filter, a 2-D filter which relates the potential on one plane to the potential on another plane; this permits simple regularization. The location estimation problem can then be formulated as a multiple hypothesis minimum error probability generalized likelihood ratio test. Numerical results are presented. The Cramer-Rao bound on the variance of any unbiased depth estimator is discussed to interpret the results.>

Novel fast algorithm for solving arbitrary Toeplitz-plus-Hankel systems of equations are presented. The algorithms are analogues of the split Levinson and Schur algorithms, although the more general Toeplitz-plus-Hankel structure requires that the algorithms be based on a four-term recurrence; relations with previous split algorithms are noted. The algorithms require roughly half as many multiplications as previous fast algorithms for Toeplitz-plus-Hankel systems.>

Discrete generalized split Levinson and Schur algorithms for the two-dimensional linear least-squares prediction problem on a polar raster are derived. The algorithms compute the prediction filter for estimating a random field at the edge of a disk from noisy observations inside the disk. The covariance functions of the random field is assumed to have a Toeplitz-plus-Hankel structure for its radial part and its transverse part. This assumption can be shown to be closely related with some types of random fields, such as isotropic random fields. The algorithms generalized the split Levinson and Schur algorithms in two ways: (1) to two dimensions; and (2) to Toeplitz-plus-Hankel covariances.>

The exponential distribution (ED) developed by H.-I. Choi and W.J. Williams (see IEEE Trans. Acoust. Speech, Sig. Proc., ASSP-37, no.6, p.862-71, 1989) is applied to the problem of simultaneous determination of the velocity and spin density of a moving object in magnetic resonance imaging (MRI). The object velocity is assumed to be an affine function of location (reasonable for selectively imaged thin slice) and is determined from the trajectory of the maxima of the time-frequency distribution. Once the velocity is computed, its effects can be deconvolved from the MRI signal. This improves on previous results which require prior knowledge of either the velocity or discontinuities in the spin density. Numerical examples are presented.>

The phase retrieval problem, i.e. computing a function known to have finite spatial support from knowledge of its Fourier magnitude, is considered. By writing the support constraint as a state equation in the Fourier domain (dual to the system used in state-space harmonic retrieval), an extended Kalman filter is used to compute the Fourier magnitudes as observations. The filter is initialized by solving a linear system of equations, followed by a nonlinear set of trigonometric equations.>

The performance of two fast algorithms for estimating the power spectral density of an autoregressive process is analyzed. The algorithms perform similarly to the Burg algorithm, but require only two-thirds as many multiplications as the most efficient implementation of the Burg algorithm. This allows the high resolution associated with the Burg alogorithml to be obtained using many fewer computations. One algorithm is the deterministic form of the split lattice algorithm adjoined to the split Levinson recursions; however, its resolution is relatively poor. The other algorithm corrects a bias in the first algorithm and has resolution similar to the Burg algorithm.>

New multichannel split or immittance three-term recurrence algorithms are derived for non-Hermitian block Toeplitz matrices, extending previous work to the most general case. Since matrix inversions require much more computation than matrix multiplications (a distinction not present in the scalar case), so-called unnormalized versions of the new algorithms are given; these require only two matrix inversions per recursion. Reductions to the block-Hermitian, block-Toeplitz, Hermitian Toeplitz, asymmetric real Toeplitz, and symmetric real Toeplitz cases are noted; this allows previous work on these special cases to be put in a general context.>

The blind deconvolution problem is to reconstruct a sequence (r/sub k/) from noisy observations y/sub k/=S/sub k/*r/sub k/+n/sub k/ where (r/sub k/) and (S/sub k/) are both unknown. It is assumed that (S/sub k/) is deterministic and minimum phase, (r/sub k/) is a realization of a white process, and (n/sub k/) is white Gaussian noise. The novelty of the approach is the use of lattice-based algorithms for all parts of the problem. First, the Burg method is used to estimate the magnitude of S/sub s/(f)=F (s/sub k/). Second, the Schur algorithm is used to compute the spectral factor of mod S/sub s/(f) mod , yielding (s/sub k/). Finally, the Levinson algorithm is used to solve the Toeplitz equations for the estimate of (r/sub k/). Results of numerical simulations are presented.>

Fast algorithms are derived for the scalar and tangential directional Nevanlinna-Pick interpolation problems. These fast algorithms require fewer multiplications than standard algorithms for these problems. The scalar interpolation problem, which arises in SISO (single input single output) H infinity optimization, is solved using a three-term recurrence similar to that of the so-called split algorithms for linear prediction. The tangential directional interpolation problem, which arises in MIMO (multiple input multiple output) H infinity optimization, is solved using an algorithm similar to a series of Householder transformations. Applications to H infinity optimization are noted.>

Several fast algorithms for the recovery of the impedance and reflection coefficients of a discrete asymmetric lossy medium are presented. The data required are the reflection response to an impulsive input on the surface, and the transmission response to an impulsive input incident from below. It is shown that these efficient algorithms are related to the generalized Levinson algorithm, the Schur algorithm, and the QR lattice algorithm for solutions of symmetric and asymmetric Toeplitz matrices equations. These unified, conceptual connections aid in the understanding of the symmetric and asymmetric Toeplitz matrix algorithms in the linear prediction literature.>

Fast algorithms for computing the linear least-squares estimate of a three-dimensional random field from noisy observations inside a sphere are derived. The algorithms can be viewed as generalized split Levinson and fast Cholesky algorithms, since they exploit the (assumed) Toeplitz structure of the double Radon transform of the random field covariance, and therefore they require fewer computations than would solution of the multidimensional Wiener-Hopf equation. Unlike previous generalized Levinson algorithms, no quarter-plane or asymmetric half-plane support assumptions for the filter are necessary; nor is the multidimensional filtering problem treated as a multichannel (vector) filtering problem.>

A layer-recursive algorithm is derived for reconstructing an equal-travel-time layered acoustic medium from its transmission response to an impulsive plane wave incident from below. The medium is bounded below by an infinite half-space from which it is probed, and it is bounded above by a free surface at which its transmission response is measured. This problem arises in the nondestructive evaluation of layered composite materials, in which the transmission response is measured. The algorithm is related to the lattice algorithm for computing the QR factorization of a matrix. The algorithm is related to the Levinson, fast Cholesky, and lattice algorithms, all of which are known fast algorithms for signal processing.>

The author shows that by interpreting the lattice filter as a discrete transmission line, the problem of determining the reflection coefficients associated with a Gaussian autocorrelation can be solved easily using the Schur algorithm. These reflection coefficients have been shown to be in geometric progression; it is claimed that this has been done in a much simpler and more enlightening manner than in the presentation of G. Jacoriti and G. Scarano (see ibid., vol.75, no.7, p.960-961, 1987). The geometric progression of reflection coefficients leads to a stationarity property of the discrete transmission line, which accounts for the striking simplicity of the expressions for the waves traveling in the line.>

This paper presents an algorithm for estimating a three-dimensional homogeneous random field from noisy observations inside a sphere of finite radius. It thus constitutes an alternative to solving a multi-dimensional Wiener-Hopf equation. The algorithm is fast in that it exploits the structure of the integral equation kernel to reduce the computation required to construct the optimal filter. It is thus an extension of similar fast algorithms that have been obtained for the one-dimensional and isotropic random field estimation problems. The problem of estimating an isotropic field at the center of a sphere of observations is treated as a special case.