I survey the present status of three subfields of reconstruction tomography, X-ray CT, emission CT, and magnetic resonance CT, and mention some new results and insights, as well as open problems. This is for my Frontier in Imaging Science lecture at the IEEE Nuclear Science International Workshop on Physics and Engineering of Computerized Multidimensional Imaging and Processing, April 2-4, 1986, Irvine, CA.

The double Fourier decomposition of the sinogram is obtained by first taking the Fourier transform of each parallel-ray projection and then calculating the coefficients of a Fourier series with respect to angle for each frequency component of the transformed projections. The values of these coefficients may be plotted on a two-dimensional map whose coordinates are spatial frequency w (continuous) and angular harmonic number n (discrete). For |w| large, the Fourier coefficients on the line n=kw of slope k through the origin of the coefficient space are found to depend strongly on the contributions to the projection data that, for each view, come from a certain distance to the detector plane, where the distance is a linear function of k. The values of these coefficients depend only weakly on contributions from other distances from the detector. The theoretical basis of this property is presented in this paper and a potential application to emission computerized tomography is discussed.

Iterative 3-d reconstruction of images from a few projections for application in Digital Subtraction Angiography of cerebral blood vessels is investigated. The method reconstructs a 3-d vascular network in a series of 2-d sections by filtered back-projection and summation of the projections followed by non-linear iterative deconvolution of the filtered-summated image. The reconstruction method is tested on a 3-d wire phantom resembling vasculature. Four to eight projections are used as input, where each projection is separated by rotation and is assumed to represent a projected stack of mutually parallel axial sections. The results show that given at least 4 projections there is sufficient detail in the reconstructed phantom to determine the location of major structures.

Computed tomography (CT) can gain a significant increase in the signal-to-noise ratio using constraints to better define regions of known intensity. This is accomplished in positron emission tomography (PET) by reducing the region of unknown activity with time-coincidence circuitry. Mathematically, constraints can be implemented into reconstruction algorithms using a priori information, such as the use of an x-ray CT image to define regions of radionuclide uptake in single photon emission computed tomography (SPECT) or PET imaging. In addition to regional constraints, intensity constraints can also be included into the model equations. This may be especiallly useful where the data has a high degree of contrast as in DSA limited-angle tomography. The maximum entropy reconstruction problem with constraints is reduced to a dual optimization program using duality principles. The solution to the reconstruction problem is determined by solving for the optimum solution to this dual program. Inequality con-straints complicates the problem in that the derivatives of the function to be optimized may not exist at the boundary points. Using penalty function techniques, the problem can be structured as an unconstrained optimization program. This way, solutions can be determined using gradient type of algorithms which require an essentially smooth function as feasible solutions approach boundary points.

We have designed a HIgh Spatial resolution Positron Emission Tomograph with a 2π solid angle coverage, based on the use of MultiWire Proportional Chambers with dense drift space converters. The imaging capabilities of HISPET have been studied with a fully 3-D Monte Carlo simulation of the tomograph and various source/phantom configurations. By using a simple 3-D filtered back-projection algorithm a spatial resolution of 4 mm (FWHM) has been obtained. The imaging results for several radioisotope distributions are presented.

This paper proposes a series expansion method for the direct reconstruction of a three dimensional image from its line integrals. We clarify a relation between line integrals and plane integrals of a three dimensional image. The relation yields a pair of basis functions in the image space and the projection space. The basis functions are similar to the singular functions of the Radon transformation. Since projection data are practically obtainable along a finite number of beams, we also propose a reconstruction formula which is expressed by a finite sum of products of sampled radiographs and a set of functions. The present reconstruction method is applicable to any incomplete projections and any scanning modes.

The reconstruction of a three dimensional representation of arbitrary vascular beds from multiple projections is discussed. Operator notation is used to represent algorithms used in the reconstruction process. These operators provide the interconversion of the various data structures. The flow of the reconstruction process is reviewed and examples of reconstruction are given. Applications to fluid dynamic analysis of vascular function are presented.

Steepest descent (SD) and conjugate gradient (CG) methods have been applied to limited-angle, 3D density reconstructions of longitudinal, circular tomography. Both variable and constant stepsize have been considered with and without positivity constraint. Simulational as well as clinical studies are reported. From the simulation studies, CG with constant steplength and with positivty constraint is found to be most effective.

Rotate-rotate fan beam configuration commonly used in the third generation x-ray computed tomography scanners is susceptible to ring artifacts resulting from non-uniform detector response. A method which incorporates the relative variation of the angular velocities of the x-ray tube and the detector array is proposed. Simulation results utilizing the new technique are presented and compared with the third and fourth generation scanners.

Digital imaging techniques are being developed for analysis of gaseous flowfields. Methods for acquiring single frames ("snapshots"), sequences of frames ("movies") and three-dimensional images are described. The data can be processed to yield distributions of chemical species, density, pressure, temperature and velocity.

Diagnostic radiology is based on interpreting projection images of the scalar quantity known as the x-ray total linear attenuation coefficient. While this information is adequate in many cases, there are situations where more detailed knowledge (e.g. of the chemical composition of an organ) is desirable. Such information is precluded in principle when only transmission radiation is measured. It is shown that measurement of the angular variation of the coherent and Compton radiation scattered from a small sample allows the chemical composition of the sample to be determined. Using Hubbell's compilation of atomic form factors and incoherent scatter functions for the six commonest elements of the human body, scatter data have been simulated with realistic noise components for some representative compounds. Good agreement is obtained between the amounts of each element derived from fitting the scatter data and those used to generate the data. A scatter system for chemical imaging is proposed, based on a monochromatic pencil x-ray beam and detector arc, and incorporating translation and rotation movements as in first generation transmission CT. An iterative technique is described, analogous to those developed for SPECT, which allows the spatial and angular variation of the coherent and Compton scatter strengths to be reconstructed. Chemical imaging with x-ray scatter (CIXS) appears feasible as a technique to provide information for diagnostic radiology which is unobtainable by other means.

A technique, based upon the interrogation of x-ray scatter, has been used to construct precise animated displays of the three-dimensional surface of the heart throughout the cardiac cycle. With the selection of motion amplification, viewing orientation, beat rate, and repetitive playbacks of isolated segments of the cardiac cycle, these displays are used to directly visualize epicardial surface velocity and displacement patterns, to construct regional maps of old or new myocardial infarction, and to visualize diastolic stiffening of the ventricle associated with acute ischemia. The procedure is non-invasive. Cut-downs or injections are not required.

The theoretical basis of an approach towards the inversion of scatter-data is presented. The method is applicable to 3-D imaging, for a variety of radiations, but is specifically developed here for the case of acoustic wave scattering from an inhomogeneous medium consisting of velocity fluctuations. The inversion procedure is exact, and explicitly takes into account all orders of multiple scattering. An interesting feature is that progressively higher resolution of the recovered image is obtained as higher order scattering is progressively incorporated into the inversion procedure. The method represents a hitherto unexplored approach towards the problem of reconstructing an image of an object from the measurement of its scattering amplitude.

In a conventional ultrasound digital sector scanner, the sampling space is in the polar coordinates while the display space is in the Cartesian coordinates, which necessitates a coordinate transformation process. As a result a number of artifacts appear in displayed image. Recently, a new data space, to be called the uniform-ladder space (ULS), in which the sampling space is the same as the display space, has been introduced. In this approach, sampled points are uniformly distributed in each horizontal line as well as in each radial ray so that the data are displayed as they are acquired without any interpolation process required. In this paper mathematical models of various interpolation schems are devised to compare their performance.

The human brain produces a measurable magnetic field in response to specific stimuli. This magnetic field, called the neuromagnetic field, can be localized to regions of the brain that actively process the stimulus. Reconstruction of functional brain images from the neuromagnetic field, a process that we call neuromagnetic imaging, represents a new non-invasive modality dependent upon brain structure and activity. The basics of neuromagnetic imaging including problems inherent to this approach, models invoked to solve some of the problems, computer simulation studies and preliminary experimental results obtained with a line source are presented.

Although historically the Born or Rytov linear approximations have received a great deal of attention, it is now more apparent that only a full nonlinear formulation of the inverse scattering problem, such as those we have developed, provide the accuracy for quantitative clinical ultrasound imaging. Our inverse scattering solutions have been developed to reconstruct quantitative images of speed of sound, density, and absorption using the exact Helmholtz wave equation without perturbation approximations. We have developed fast algorithms which are based upon Galerkin or moment discretizations and use various iterative solution techniques such as back propagation and descent methods. In order to reconstruct images with reflection-only scanner geometries we have extended our algorithms to include multiple frequency data. We have demonstrated a procedure for imaging inhomogeneous density distributions. We also discuss the significance and potential applications of these new methods.

Three-dimensional (3D) surface display is an alternate way of presenting to the physician information that is available in a sequence of two-dimensional CT or MRI scans. The aim is to present organs (or parts of organs) as they would appear if they were removed from the body, possibly cut open, and viewed from user-selected directions. In recent years there have been a number of papers discussing the clinical utility of this approach. In nearly all these papers the presentation of the surface consists of single images of the objects of interest. In these monoscopic images, depth perception is conveyed by the differential shading that is computed as if light were shining on the surface. This is augmented by the silhouette of the external features. Since shading is dependent on the distance from the light source and the angle of the surface to the light rays, these two effects may oppose each other, especially with perception of details in depths of cavities. In addition, the detail inside a cavity cannot be silhouetted at any viewing angle. Since many anatomical surfaces have significant information in the depths of cavities (e.g., orbits, neural foramina, cardiac cavities), the addition of stereoscopic depth perception and motion should be clinically useful. In a recent article, we presented 3D surface displays in stereo, thereby providing an important additional cue for correct 3D depth perception. Here we discuss issues of software and hardware for dynamic stereo display and for 3D interaction with the stereoscopic images.

Adaptive histogram equalization (ahe) is a contrast enhancement method designed to be broadly applicable and having demonstrated effectiveness [Zimmerman, 1985]. However, slow speed and the overenhancement of noise it produces in relatively homogeneous regions are two problems. We summarize algorithms designed to overcome these and other concerns. These algorithms include interpolated ahe, to speed up the method on general purpose computers; a version of interpolated ahe designed to run in a few seconds on feedback processors; a version of full ahe designed to run in under one second on custom VLSI hardware; and clipped ahe, designed to overcome the problem of overenhancement of noise contrast. We conclude that clipped ahe should become a method of choice in medical imaging and probably also in other areas of digital imaging, and that clipped ahe can be made adequately fast to be routinely applied in the normal display sequence.

A fully interactive physician's workstation which supports the display and manipulation of 3-D medical datasets has been constructed. Based on the "Voxel Processor" architecture, the system uses special-purpose computer hardware to provide a variety of display modalities including shaded surface display, multi-planar reconstruction, and 3-D display of dynamically changing objects - all with real-time update rates (15 frames/second). Grey-scale voxel densities are retained throughout the processing so that windowing and thresholding can be accomplished interactively in true real-time. Other real-time features include rotation, scaling, and slice plane control. This paper describes the 3-D physician's workstation and summarizes the major requirements for future 3-D workstations as they relate to specific medical applications. The general architecture and functional characteristics of advanced 3-D display processors supporting these requirements are presented.

Three dimensional surface reconstruction and visualization are well recognized to be useful in the area of pre-operative surgical planning and radiation treatment therapy planning. This function has been integrated into the capabilities of the comprehensive PACS system of the Harbor-UCLA Diagnostic Imaging Center. This paper describes how the user can transfer and view CT or MR images from the scanners on the radiation treatment node thru SYNERCOM* or microwave. Afterwards, an interactive segmentation algorithm is used to isolate an arbitrary region of interest and gray level thresholding is applied to generate an intermediate bin'ary image. A ray tracing algorithm is used for visible surface formation and a shading value is then calculated. Different views of the object can be generated by interactive manipulation of the rotational angles. The final 256 x 256 images are shown in cine-mode on a color Deanza 5512 system.

At the Medical Image Processing Group, we have been developing for the past several years software packages for the display and analysis of 3D medical objects. 3D98 is the most recent among them. It has been designed to work on a General Electric 9800 CT scanner or on its independent console. This paper describes the processing features of 3D98 and outlines some of the algorithms implemented in it. From the information processing point of view, 3D98 goes through four distinct steps: house keeping, dialogue, data processing, and termination. The paper describes these steps and illustrates the available options using real CT data.

"Help me communicate more quickly and more effectively with referring clinicians". This request was the driving force behind the installation of the AT&T CommViewm System at Duke. The CommView System is a type of Digital Image Management System and Picture Archival Communication System whose chief purpose is to deliver interpreted diagnostic images to referring clinicians and attending physicians. The system acquires electronic images from modal-ities in a diagnostic imaging facility, stores these images in computer managed patient files and distributes these on demand over fiber optic cable to Display Consoles. The CommView System was designed at AT&T Bell Labs; it uses fiber optic ribbon cable between buildings fused to multistrand lightguide building cables to distribute images, typically around a medical center or campus at data transfer rates of 40 Mbps. This paper gives the rationale used in designing a start-up network and placing the initial equipment for a field trial of the AT&T CommView System in the Radiology Department of Duke University Medical Center.

There is a growing interest in the use of three-dimensional (3D) cranial reconstruction from CT scans for surgical planning. A low-cost imaging system has been developed, which provides pseudo-3D images which may be manipulated to reveal the craniofacial skeleton as a whole or any particular component region. The contrast between congenital (hydrocephalic), normocephalic and acquired (carcinoma of the maxillary sinus) anomalous cranial forms demonstrates the potential of this system.

Single-photon emission computed tomography (SPECT) is an inverse problem in which one wants to determine the distribution of a radionuclide from a set of measured projections. Like most inverse problems, it is ill-posed and does not admit of an exact solution. In this paper we review various methods from the literature on inverse problems that are applicable to SPECT. Topics considered include the discrete representation of continuous objects, the intrinsic dimensionality of an object, null functions, the role of prior information, and various reconstruction principles, including maximum likelihood, least squares, and Bayesian methods.

In an effort to overcome limitations of conventional collimation for imaging high energy (>300 keV) single photon emitters, studies were conducted to measure the performance parameters of an electronically collimated gamma camera for imaging a Cs-137 source at 662 keV. In addition to theoretical computations, coincident counts from a point, line and disk source of Cs-137 were acquired with a prototype system to investigate the resolution and sensitivity characteristics. The results show that, unlike conventional collimation, the performance of an electronically collimated system improves with energy, thereby enabling efficient imaging at high energies.

A post reconstruction method of attenuation compensation for Single Photon Emission Computed Tomography (SPECT) has been investigated that offers a new approach to the problem of quantitation. A modified correction matrix is generated for attenuation compensation in which the Linear Attenuation Coefficient (LAC) for each pixel is assigned a value depending on the radial distance of the pixel from the true section boundary. Attenuation compensation of transverse section images of small and large volume sources of Tc-99m in phantoms using this modified matrix indicated that a known quantity of radionuclide could be determined to better than 10%. The scatter fraction was estimated as the difference in the corrected section images using a multiplicative matrix generated with a constant LAC for each pixel and the modified matrix proposed in this report.

A simple and economically practical method of improving the sensitivity of camera-based SPECT was developed using converging (cone-beam) collimation. This geometry is particularly advantageous for SPECT devices using large field-of-view cameras in imaging smaller, centrally located activity distributions. Geometric sensitivities, spatial resolutions, and fields-of-view of a cone-beam collimator having a focal length of 48 cm and a similarly designed parallel hole collimator were compared analytically. At 15 cm from the collimator surface the point-source sensitivity of the cone-beam collimator was 2.4 times the sensitivity of the parallel-hole collimator. SPECT projection data (simulated using Monte Carlo methodology) were reconstructed using a 3-D filtered backprojection algorithm. Cone-beam emission CT (CE-SPECT) seems potentially useful for animal investigations, pediatric studies, and for brain imaging.

The ability of focused collimators to achieve both high sensitivity and high resolution at depth makes them highly desirable for SPECT imaging. Tomographic systems with stationary detector rings are known to have the advantages of simplicity, stability and tolerance to variations in detector response. We have combined these two desirable features in the design of a multidetector ring SPECT system with a unique modular focused collimator for brain imaging. In a manner similar to that employed by 4th generation CT scanners, each detector unit acquires a set of fan beam projection data while the collimator assembly makes a complete rotation. Our investigation shows that the new tomograph can be expected to achieve a spatial resolution of 8 mm and a system sensitivity of 2,500 cps/pCi/cc for a 13 mm thick image slice. The high performance characteristics of the new SPECT system should make it a useful tool in functional imaging of the brain, especially in perfusion studies using the current 1-123 or Tc-99m labeled agents.

Improved SPECT reconstructions have been obtained by the use of an Inverse Monte Carlo (IMOC) reconstruction algorithm. Monte Carlo solution to the photon transport equation for SPECT generates an accurate system response model for use in an iterative Maximum Likelihood estimation (MLE) algorithm. Reconstruction from simulated projections demonstrate that the MLE algorithm continues to converge at 1000 iterations for noise free data. Reconstruction from experimentally acquired projections reveals superior RMS noise and resolution recovery when compared with filtered backprojection for simple phantoms. Compensation for scatter and attenuation is provided by including these effects in the Monte Carlo simulation.

The performance of large diameter NaI (T1) scintillation crystals applied to the detection of nuclear images is bounded by specific physical limitations. The results presented are analytical expressions for the detection efficiency and the point spread function (PSF) of scintillation detectors in the absence of multiple scattering. The PSF is determined both for a perpendicular incident photon beam and for oblique incidence. The scintillation detection efficiency is determined only for a co-axial point source. An expression for the efficiency for oblique incidence is given in integral form.

This paper presents recent detector developments and perspectives for positron emission tomography (PET) instrumentation used for medical research, as well as the physical processes in positron annihilation, photon scattering and detection, tomograph design considerations, and the potentials for new advances in detectors.

It is well known that scattered radiation in PET reduces the contrast in images. There is about the same fraction of scattered events in both emission and transmission scans. The technique presented here reduces the scatter fraction in transmission scans. An orbiting rod source is used in place of the ring source normally used for transmission scans. The source moves in steps rather than in a continuous motion. At each point only coincidences between detector pairs for which the source and detectors are co-linear are accepted. In this way all wide angle scatter is rejected. The attenuation coefficient for water at 511 keV was measured as .095 cm-1 with the orbiting source compared to .086 cm-1 with a conventional source and the value due to the electron density of water .097 cm-1 The recovery coefficient for hollow cylinders 23 mm in diameter increased from 77% for a conventional to 93% for a rod source transmission scan. If the source is collimated to a flat fan beam, a much higher activity could be used, and transmission scans could be performed independent of activity in the patient section being scanned.

Different designs of high-resolution positron camera systems have been investigated. The design goal is an instrument that can measure the whole brain with a spatial resolution of 5 mm FWHM in all directions. An overall spatial resolution of 5 mm requires crystal dimensions of 6 x 6 x L mm3, or less, L being the length of the crystal. Timing and energy requirements necessitate high-performance photomultipliers. The small size scintillator crystals can currently only be read in schemes based on Anger techniques or by utilizing surface-sensitive PMTs. In the future, photodiodes can replace photomultipliers. In the present work, two different approaches have been investigated: Anger techniques and position-sensitive PMTs. Preliminary results of the two detector designs are reported.

The use of detector arrays in PET has been proposed and studied as a economical means of obtaining very high resolution images. The properties of detector arrays for use in high resolution PET scanners were investigated. Analytical approximations and Monte Carlo simulations were used to design detector arrays consisting of six to eight crystals coupled to two photomultipliers to allow identification of the individual crystals. Arrays of 2.85 mm thick crystals of Bismuth Germanate (BGO), Gadolinium Orthosilicate (GSO), and Barium Fluoride (BaF2) were examined. The effect of interdetector materials such as lead, and plastic on positioning accuracy was tested. Assembled arrays of six 2.85 mm thick BG0 crystals yielded line spread function FWHMs of 2.4 to 3.2 millimeters. The limiting resolution of detector arrays was found to be defined by the scintillation light yield of the crystals, the light gathering efficiency of the detector arrays, and the fraction of interdetector scatter.

The Maximum Likelihood Estimator (MLE) algorithm for tomographic image reconstruction is being investigated in substantial detail by a number of research groups, as it appears to promise images with very low noise and increased sharpness when compared with filtered backprojection techniques. Recently, however, it has been found that the reconstruction of data from uniform activity distributions exhibits strong peaks and valleys when the number of iterations increases toward a maximum in the likelihood function. This problem has now been investigated with our Positron Emitter Beam Analyzer (PEBA) camera, which, because of its small size and favorable geometry, has allowed an analysis with enough detail to find the origin of that apparent instability. The findings can be summarized as follows: 1) The very low noise of the MLE reconstructions comes about by the ability of the Poisson-based MLE algorithm to generate an image which favors the matching of experimental data (detector pairs) containing few counts. 2) The image instability at a high number of iterations is a direct consequence of the above characteristic. 3) The matrix of probability elements needed for the MLE reconstruction provides the link between the two above observed phenomena. It appears that, by proper system design, it is possible to obtain the favorable low noise characteristic without the instability. The applicability of the above findings to true tomography (PEBA does not carry out a true tomographic reconstruction) seems direct, but confirmation should be obtained by further research on the question.

A computerized brain atlas data base has been developed, primarily for use in positron emission tomography. The underlying information was derived from a digitized cryosectioned cadaver brain. The atlas can be individually adjusted to fit a wide range of patients with reasonable accuracy. The necessary transformations are chosen so that the atlas will fit an initial set of CT or NMR images of the patient. The individualized brain atlas opens several new possibilities in the quantification and evaluation of PET data. It can - be used to select a suitable patient orientation during the PET study, - improve the attenuation and scatter corrections, - supply external information to be used in the image reconstruction, such as proper 3-dimensional regions of interest, - serve as a vehicle in the comparison of different examinations of the same patient, thus reducing the need of reproducible fixation systems, - facilitate the merging and comparison of results from different individuals or groups of individuals (i.e. serve as a reference atlas). The brain atlas can also be used with other imaging modalities as well as in stereo-taxic surgery.

The capabilities of superconducting magnets for magnetic resonance imaging and spectroscopy are described. Emphasis is given to recent developments directly affecting the imaging process. Improved homogeneity and higher field magnets are described.

A method for the multidimensional measurement of flow fields in liquids by nuclear magnetic resonance (NMR) is discussed. It is assumed that the flow field at a given point in space is given by the superposition of two components. The first component, characterizing the average bulk flow, is called the coherent velocity. The incoherent flow field super-imposed on top of the coherent one is assumed to be random as one may find in turbulent flow. Limiting cases are discussed. It is shown that coherent flow introduces a phase into NMR images whereas the incoherent one causes a reduction of signal intensity similar to the one given by T2 or diffusion processes.

The parameters for designing surface coils for NMR imaging have been identified by introducing an equivalent circuit for it. Also a method of evaluation and testing of surface coils has been discussed in the lab environment.

Magnetic resonance imaging may be extended by various orders of modulation-encodings to acquire simultaneously FID data-sets representing very high dimensionality. For example, we may encode three spatial-position dimensions, three molecular-velocity dimensions, a chemical shift dimension, three perfusive flow dimensions, and so on. Careful observation on even the common 2D-MRI scans, however, yield some puzzling results about the behavior of contrast/detail detectability. The multi-dimensional options and capabilities of MRI arise because the physical system inherently is an interferometric data-acquisition modality. There is a "hidden" dimension in the data-acquisition, namely, the resolution scale of digitization of the FID-data values. This dimension implies requirements on the dynamic-range of the data-acquisition, just as the Nyquist-limit implies requirements on the fineness of the temporal sampling mesh. One must account for the limitations so imposed. We propose here a new quantitative imaging descriptor, the Detail-in-Image to Totality-in-object Ratio (DITOR), whose values impose limitations upon the system dynamic-range required. These limitations become more severe as the imaging dimensionality becomes larger, and MRI system behavior in this regard is very different from ordinary x-ray medical imaging. This paper shows several examples of phantom and human studies which may well illustrate the limitations implied by DITOR, and develops the role of stochastic noise contributions to the DITOR requirements.

An equation which relates the signal-to-noise ratio (SNR) in a two dimensional NMR tomographic image to several basic image parameters is derived. These variables are the number of FID's averaged per view, number of views, pixel size and slice thickness. The theoretical equation was studied experimentally by varying one parameter at a time. Good agreement between experimental results and theoretical prediction was found.

Moving spins have a significant effect upon the received MRI signal, which is seen, depending upon the pulse sequence utilized, as a modulation in signal intensity and/or phase relative to that of stationary spins. This can be used in MRI to distinguish between blood vessels and stationary anatomic structures. Monoplanar time-of-flight techniques [1-9] use pulse sequences which modulate the signal intensity from spins in a vessel flowing perpendicularly through a slice. This modulation is the result of the fact that the signal is composed of signal strengths from different spin populations whose relative ratio will depend upon the thickness of the slice, the velocity of fluid perpendicular to the slice, and pulse sequence parameters. If two different pulse sequences are used with identical signals from stationary anatomy, one producing increased signal in the region of flowing spins and the other producing decreased signal, then a subtraction image can be formed depicting only regions of flow [2]. Here we develop expressions for the signal intensity of two pulse sequences which are expected to give optimum contrast for imaging flowing blood free of overlying non-vascular anatomy. This imaging technique would provide a means for imaging blood vessels without the use of ionizing radiation or contrast injection; furthermore, it provides information about the presence of flow in vessels. As such, it may be promising as a method to evaluate vessels which are not accessible to standard angiographic imaging.

Since August Krogh's original work in 1919, organ perfusion has been modelled extensively after what is known as "single capillary" model. In these models, the perfused organ is conceptualized as composed of identical building blocks, each consisting of either a circular or a hexagonal tissue cylinder supplied by a single capillary. A vast literature has grown up on such models, the most noteworthy of which are listed in references. Krogh's original model, shown in Fig. 1, consists of identical parallel capillaries each supplying it's own circular tissue cylinder. A variant of Krogh cylinder model for a well-perfused organ with parallel capillaries with equal flow and supplying hexagonal axisymmetric tissue columns of equal lengths is that of Bassingthwaighte and is shown in Fig. 2. Another approach to perfusion modelling is that of compartmental analysis where the perfused organ is divided into a blood and a tissue compartment. Various assumptions are then made regarding the nature of these two compartments. The tissue compartment is usually considered to be well stirred and for the blood compartment, a linearly distributed capillary bed is assumed. This is equivalent to a single capillary with an axial variation of concentration down its length. Usually the mixing in the radial direction is assumed to take place infinitely rapidly. Most important of these models are those of Renkinll , Johnson and Wilson, and Levitt. The model presented here is shown in Fig. 3, and assumes that the blood compartment, where the microcirculation takes place, is a random network of interconnected microchannels, the network being in contact with the tissue compartment. Blood from the arterial side enters this capillary network in which microcirculation takes place with or without exchange with the tissue compartment, and finally the perfused blood exits at the venus end. This clearly departs from the single capillary models since no assumption is made as to the orientation, distribution and identity in structure of the microcirculatory capillaries. For the same reason, this also differs from the linearly or any other deterministically distributed capillary bed of the compartmental analysis. Moreover, the concentration distribution in such a network is not inferred or prescribed, but is calculated from the basic hydrodynamics of the medium. This leads to a concentration distribution in both parallel and perpendicular (to mean flow) directions in terms of the dispersion of statistical hydrodynamics.