WO1995030914A1

WO1995030914A1 - Post-processing of mri images

Info

Publication number: WO1995030914A1
Application number: PCT/US1995/005673
Authority: WO
Inventors: John Listerud; James F. Kohli; Scott Atlas
Original assignee: The Trustees Of The University Of Pennsylvania
Priority date: 1994-05-03
Filing date: 1995-05-03
Publication date: 1995-11-16
Also published as: AU2436195A

Abstract

Post-processing of angiographic images formed by, e.g., a Magnetic Resonance Imaging (MRI) system or a Computed Axial Tomography (CAT) system for producing highly resolved images is presented. The system preferably contains a data gathering device (20) for creating a data set, a computer processor (10) interfaced with the data gathering device and adapted to receive the data set for projecting an image of the object and thereby creating a projected data set, and an extracting device (100) interfaced with a processor (10) and adapted to receive the projected data set for extracting image features therefrom by thresholding the projected data set to at least two predetermined values. Post-processing techniques applied in accordance with the present invention produce highly resolved and useful clinical images which are particularly suited for the diagnosis of diseases and maladies which afflict the body and are particularly useful for magnetic resonance angiography or applications wherein fluid flow in the body must be monitored.

Description

POST-PROCESSING OF MRI IMAGES

Field of the Invention

This invention relates generally to methods and apparatus for obtaining images. More specifically, this invention relates to imaging objects of interest and extracting features from the images of objects of interest to obtain high quality clinical images.

Background of the Invention

Imaging for medical purposes dates from the earliest days of X-ray imaging to the modern CAT scans and magnetic resonance imaging (MRI) devices which use sophisticated machines and advanced software to provide images of the human body. Today, MRI has become an invaluable tool for the clinical and diagnostic study of virtually all parts of the human body. Generally, MRI entails subjecting a body part to an alternating radio frequency magnetic field in the presence of a large static field, thereby causing the hydrogen nuclei in the body part to absorb and emit radio waves. The radio waves emitted by the hydrogen nuclei are useful in characterizing the body part and can yield significant information about the state of the body part, including information about some diseases which afflict it. Radiologists or other clinicians can then analyze this information to investigate and diagnose the maladies and afflictions which affect the patient.

There are several MRI techniques which have been developed that are tailored to perform analyses on different parts of the human body. One such technique is Magnetic Resonance Angiography (MRA) which utilizes the characteristics of flowing material to separate moving blood from stationary tissue, thereby obtaining images of blood vessels which may be afflicted with arteriosclerosis, aneurysms, blood clots, and other cardiovascular diseases. MRA is a useful improvement over conventional angiography, since conventional angiography requires that a small catheter be threaded into a vessel so that X-ray opaque dye can be injected therethrough while X-ray film is being exposed to produce images of the vessel.

Conventional angiography is an invasive procedure, due to the insertion of the catheter and the exposure of the patient to X-ray opaque dye. Although complications are infrequent, when they occur they can be severe. In contrast, MRA is non- invasive, as it can be performed without insertion of a catheter and usually without the use of a dye. Therefore MRA is much safer. However, cross-sectional imaging techniques tend to produce large numbers of images, and because good quality MRA requires tracking vessels with high resolution over long segments, MRA produces some of the largest numbers of images. At the University of Pennsylvania for example, a routinely performed MRA exam of the vessels of the leg for peripheral vascular disease acquires multiple 2.0 mm thick axial cross- sectional images starting in the pelvis and covering the entire length of the leg, producing over 400 "slices". Direct examination of these "slices" is tedious, and more importantly, cannot always guarantee that the reviewer will be able to appreciate the three dimensional relationships of structures represented in these images.

It is common to project these data onto lateral projection images in order to reduce the number of images that it will be necessary to review, and to better present the three- dimensional relationships of the imaged structure. Figures 1A and IB illustrate prior techniques of projecting onto lateral projection images. Figure 1A illustrates ray tracing spatial assignment of three-dimensional data on a lateral projection image, while Figure IB illustrates prior maximum intensity projection (MIP) feature extraction which assigns the maximum value found along the ray of Figure 1A to the corresponding pixel in the lateral projection image. The inventors of the subject matter herein claimed and disclosed have recognized a long-felt need in the art for techniques to improve MRA and to obtain sharp, clear, highly resolved lateral projection images of blood vessels which would be useful for clinical diagnosis. This long-felt need has not been solved by current MRA techniques, or other MRI procedures generally. In MRA, as with other MRI techniques, it is necessary to obtain three-dimensional data to produce lateral projection images. In order to achieve such images as discussed above, the body part of interest is subjected to a large, static magnetic field and radio frequency pulses so that the hydrogen nuclei will radiate radio signals that carry information regarding the affliction of the body part. The radio signal is gathered by a radio antenna, usually referred to as a "coil", and electrical signals are induced therein which correspond to data carried by the radiated energy about the body part. These electrical signals are usually digitized and stored in a data storage device of a computer that is associated with the MRI machine that has scanned the body part. The data is then processed or "reconstructed" to obtain a three dimensional data set formed by a "stack" of cross-sectional images of the body part. Figure 2 illustrates prior sequential cross-sectional images forming a 3-dimensional data set of the lung, for example.

Other cross-sectional imaging methods, such as the X- ray Computed Axial Tomography (CAT) scan, produce similar three dimensional data sets. Recently, reports in the scientific literature describe the use of CAT scans to produce angiographic images as well. In this technique, X-ray opaque dye is injected into a superficial vein, then the body part of interest and its vascular anatomy is scanned to produce a stack of cross- sectional axial images. See, for example, S.A. Napel, M.P. Marks, et.al. "Vascular Imaging with Spiral CT and Maximum Intensity Projection" RSNA Abstracts, #271 (1992) and H. Yoshida, K. Tsuchiya et.al. "3D CT Angioimaging of Cerebral Aneurysms with Spiral Scanning", RSNA Abstracts, #400 (1992).

In order to construct the cross-sectional image of the body part and to display the image for clinical diagnosis, cross-sectional imaging machines are typically supplied with programmable, special purpose, high speed computational devices lcnown as array processors to process the data and form an image. Figure 3 illustrates a prior imaging system and the data flow scheme typically implemented therein to form the image. For example, a General Electric Company (Milwaukee, Wisconsin) Signa 3.x MRI machine having a Floating Point Systems array processor, shown generally at 10, has in the past been utilized to implement MRI. Because additional processing of imaging data is computationally intensive, it is desirable to use the array processor to further manipulate the data and obtain magnetic resonance images.

The system further usually comprises a detector or coil 20 in which object 30 is placed for imaging. Raw data storage means 40, such as a hard disk or tape drive, receives raw data 50 which has first been digitized by an analog-to- digital converter 60. The array processor 10 interfaced with the raw data storage means 40 creates an image of object 30 which is then stored on an image storage device shown generally at 70. Image storage device 70 is similarly a disk drive, tape drive, or other type of mass memory unit which can store large volumes of data. Cross-sectional image data 80 and lateral projection images 90 are produced by array processor 10. The system is usually controlled by a master computer 100 which is interfaced to one or more control and communication busses 110 and data busses 120. The data which has been digitized is typically stored as a stack of axial images which are then "post-processed" into lateral projections in order to produce images which appear like conventional X-rays. This "post-processing" step may conceptually be broken down into two separate aspects: "spatial assignment" and "feature extraction". Spatial assignment techniques include "ray tracing" and "surface rendering" and the term pertains to the internal representation of data within the computer. Feature extraction techniques include "MIP" and "3D depth cues" and the term pertains to the appearance of objects on the lateral projection image.

In order to further explain "spatial assignment", it is necessary to understand the steps involved in ray tracing. In order to accomplish ray tracing, data projection is performed by passing rays parallel to a line-of-sight through the volume. See Figure 1A again. Each ray projects onto a single picture element or "pixel" on the lateral projection image, and the intensity along the ray within the three dimensional volume is then analyzed in some fashion to determine the projected pixel's intensity.

The locations of intensities within the volume are represented by a lattice of nodes, wherein each node corresponds to a small volume element or "voxel" in the three dimensional stack of cross-sectional images. These nodes do not necessarily fall on the projection ray passed through the volume, and some algorithm must be first used to prescribe intensity along the ray based on lattice nodes in close proximity to the ray. This step can be computationally intensive.

For example, in "trilinear interpolation" as shown in Figure 4, to compute the intensity at an arbitrary "test" point along a ray, it is necessary to calculate a weighted average of the eight lattice nodes, shown generally at 130, which then form the corners of the voxel containing the ray point. This requires identifying the locations of the eight lattice nodes nearest the test point which involves at least eight machine operations. Fetching the eight node intensities requires another eight machine operations. Computation of a forward and backward weight per axis requires two additional machine operations for trilinear interpolation, and yet additional machine operations for other more elaborate algorithms.

Thus, as many as twenty-one operations are required for the trilinear weighted average, and about forty-five machine operations in total are performed per test point in the voxel. When one considers the volume of data which is typical for a ray tracing projection, performing a ray trace through a 256 x 256 x 64 voxels requires that the test point be rastered through roughly four million locations. In this scenario, about 180 million computer operations are performed in the spatial assignment phase even before feature extraction techniques are applied to resolve the image. The number of machine operations grows proportionally to the third power of the desired digital resolution of the image. This is a computationally expensive and time consuming process which often must be accomplished either "off-line" in batch rather than in real time, or by a high speed specialized computer like an array processor.

Spatial assignment techniques such as "surface rendering" have been utilized to increase efficiency by reducing the amount of data, optimally preserving only that information required to project the image with the features of interest. Such increases in the efficiency of spatial assignment by selectively reducing data may result in loss of flexibility in the spatial assignment's ability to portray features of clinical interest which are unrecognized at that time. There has not heretofore been a solution to this problem and it would be useful to design data processing and post-processing techniques which introduce efficiency in this endeavor, but which also ensure that high quality clinical images are obtained, both now and in the future.

Once three-dimensional data has been acquired and spatial assignment is performed, it is generally necessary to assign an intensity to the projection image pixel. This post¬ processing step will be referred to here as "feature extraction" . In the case of ray tracing, feature extraction describes the algorithm which, given a plot of intensities along a ray, assigns a pixel intensity to the pixel in the lateral projection image onto which the ray falls. The industry standard in feature extraction is called "maximum intensity projection" ("MIP"). See L. Dumoulin, and H.R. Hart, Jr., Radiology, 161, 717-720 (1986); P.N. Ruggiere, G.A.. Laub, T.J. Masaryk, M.T. Modic, Radiology, 171, 785-791 (1989), the teachings of which are specifically incorporated herein by reference. In MIP, magnitudes of intensity of points along a particular ray trace are compared and the maximum (i.e., brightest) value amongst them is assigned to the pixel in the lateral projection.

An important alternative to MIP is "summation" wherein the intensity of the lateral projection pixel is the sum of intensities along the ray. This mimics the effects of an X-ray beam which is attenuated by X-ray dense structures as it passes through the three dimensional volume, and where the final intensity on the film represents the total accumulated intensity. In order to further understand the properties of various feature extraction methods, it is useful to examine the common characteristics of several prior cross-sectional angiography techniques. The acquisition parameters are optimized to increase the contrast between stationary tissue and flow. It should be remembered that vasculature constitutes a very small volume of the total organ tissue. For example, a large vessel in the head may be 3 or 4 mm in diameter, while the head may perhaps be 15 cm in diameter. Thus, the volume of the vessel is on the order of 0.1 % of the total head volume. In a typical MR image of 256 pixels on a side, out of a total of 65,536 pixels approximately a few hundred belong to vessels.

An MRA sequence is optimized to obtain high contrast (i.e., signal difference) between flow and stationary tissue in order to separate the two on the basis of image intensity. It is less widely recognized that the range of signal intensities associated with flowing signal constitute from 50% to 90% of the total range of signal intensities across the image. In Figures 5A and 5B, histograms of a type of MRA are shown which demonstrate this phenomenon. A histogram is a plot of the "pixel count" (i.e., the number of pixels or which contain a particular intensity) versus the intensity value. Note the dynamic range of the flow signal is so large and that the spatial volume of vasculature is so comparatively small that in almost all cases, the range of intensities over which flow is spread is much greater than the number of flow pixels. Consequently, it can be expected that the "average bin occupancy" in the flow range ought to be significantly less than one. Because tissues usually have a characteristic average intensity, one may often discern peaks in a typical histogram. For example, if a body part contains fat and muscle, one may see two possibly overlapping peaks on the histogram, each corresponding to their respective tissue types. In Figures 5A and 5B, no such peak may be associated with flowing signal. Rather, the pixel count is fairly constant across the range of flow intensities, or in other words, there is great variation in flowing signal intensities. There are many reasons that intra- and intervessel signal variations exist in MRA. Greater fluid flow related enhancement at the center of a flow stream, and greater intravoxel dephasing at the lumen of large vessels where the shear strain rate (that is, the velocity gradient) is greater both contribute to signal variations. With these features of the MRA histogram in mind, the strengths and weaknesses of common feature extraction methods can be discerned, starting with Maximum Intensity Projection. MIP is a popular feature extraction technique because it is computationally simple and since practically no computer operation is as efficient as the comparison of magnitudes accomplished in MIP. Furthermore, MIP works well to display delicate structures which would often be ignored by some more sophisticated feature extraction methods. To be seen in the lateral projection image with MIP, a structure must be slightly brighter than the maximum signal along the ray from stationary tissue plus any noise signals. The dynamic range of these structures may be on the order of 5% of the total dynamic range of image intensities. By comparison, in the summation technique even though a thin structure is brighter than dimmer surrounding tissues, the sum of a few bright pixels may not be greater than the sum of many dim pixels.

A desirable property of MIP is that no threshold need be chosen by the human operator as is common to more sophisticated feature extraction methods such as summation. The maximum signal from stationary tissue forms a de facto minimum threshold for that ray. In comparison, summation techniques usually require that stationary contributions be "screened out", typically by "thresholding" the pixel intensities. A threshold value must be chosen, typically by a human operator, which divides the range of intensities.

A pixel with intensity below the threshold is assumed by the subsequent computer processing to be stationary signal

(plus noise) and is set to zero, while a pixel with intensity above the threshold is assumed to be flow within a vessel and is left unaltered. Subsequent summation of all the ray intensities will accumulate the unaltered values from the vascular structure, while the zeroed values from stationary tissue contribute nothing to the total value of the sum. The requirement of having to select a threshold is an important drawback of any technique depending on thresholding, as it constitutes a potential opportunity for human error.

Another problematic aspect for feature extraction in MRA is that the dynamic range of flowing signal is so large. For example, because intra- and intervessel signal variations can be quite large when MIP is used in MRA, an aesthetic adjustment of the computer monitor "Window" and "Level" controls required for this kind of feature extraction technique can easily eliminate fainter structures and result in narrowing the apparent vessel diameters in the resultant image. This is an inherent disadvantage with MIP which contributes to difficulty in appropriately filming the lateral projection images. Film is the overwhelmingly dominant medium on which these images are presented for reading to the radiologist, and on which the images are stored.

Variability in flow also contributes to geometric distortion of structures in MIP. For example, with MIP a bright vascular structure will always be superimposed on a dimmer vascular structure, even when the dim structure lies between the viewer and the bright structure, producing deceptive depth cues. To overcome the confusion associated with this aspect of MIP, an increased number of views at different view angles are required in order that the radiologist may rebuild the true three- dimensional character of the vasculature.

In order to further clarify the geometry of complicated vascular structures, it is standard practice to "segment" the three-dimensional volume of interest by restricting the lateral projection to a smaller sub-volume containing a particular vascular structure of interest and excluding other structures. It is further routine practice to perform multiple segmentations on each vascular structure of interest contained within the total three-dimensional volume. Both segmentation and multiple view angles increase the number of images which must be produced to overcome the false geometric cues caused by applying MIP to the large flow related signal variation seen in MRA.

In the case of summation, the large, variable dynamic range of fluid flow reduces the radiologist's ability to make an^' analogy between the summation image and conventional X-ray angiography, as in the interpretation of the "double density" sign. In conventional X-ray angiography of cerebral vasculature, a normal vascular loop which is imaged end on can be distinguished from a small berry aneurysm which is a pathological condition. This interpretation depends on a feature extraction in which the intensity of a pixel on the lateral projection image is proportional to the length of the segment of the corresponding ray describing the intersection of the ray and the vascular structure. Conventional angiography demonstrates this effect because an X-ray beam is doubly attenuated as it passes through the vascular loop. Summation feature extraction in MRA cross-sectional images as described above cannot reliably show this effect because of the existence of variability in the intensity of flowing signal.

In summary, because variability in flow causes widely variable signal intensity in MRA, because there are complications involved in picking a threshold, and because there is computational intensity due to processing cross-sectional images in order to obtain lateral projection images, there is a long-felt need in the art for feature extraction methods which can improve on the current, industry wide standard of MIP.

The inventors of the subject matter herein claimed and disclosed have realized that feature extraction and ray tracing techniques have heretofore been inadequate in cross-sectional imaging generally, and specifically for MRA applications. The prior feature extraction and ray tracing techniques cannot produce highly resolved images resembling X-rays with appropriate depth cues, appropriate recognition of "double density" situations, and minimal operator interaction in the choice of thresholds. Therefore, an unresolved long-felt need exists in the art for post-processing methods and apparatus to produce high quality images in MRA.

Stannary of the Invention

Methods and apparatus in accordance with the present invention solve the aforementioned problems and fulfill long- felt needs in the art for efficient imaging of objects and for obtaining high caliber images for clinical and diagnostic purposes. In accordance with the invention, methods of imaging an object comprise the steps of gathering data from the object and storing the data in a storage device, performing ray tracing on a data set which is representative of an object that has been scanned to obtain a substantially three-dimensional image of the object and a three dimensional data set corresponding to the three-dimensional image of the object, and thresholding the three-dimensional data set according to at least two levels to obtain a thresholded data set and to project the three- dimensional image, thereby extracting features of interest from the three-dimensional image.

In a further preferred embodiment, a magnetic resonance imaging system is provided in accordance with the present invention. Preferably, the system comprises gathering means for gathering data from an object and creating a data set that can be further manipulated by the system to obtain useful information concerning the object, processing means interfaced to the gathering means and adapted to receive the data set for projecting an image of the object based on the data set and creating a projected data set which can be further manipulated to obtain additional useful information concerning the object, and extracting means interfaced to the processing means, and adapted to receive the projected data set for extracting image features from the projected data set by thresholding the projected data set to at least two predetermined levels.

Post processing techniques employed in accordance with the invention produce highly resolved useful clinical images which are particularly suited to diagnosing diseases and maladies which afflict the body. In a more further preferred embodiment, MRA is performed by methods and apparatus provided in accordance with the present invention so that highly resolved images of blood vessels are obtained to diagnose cardiac and pulmonary afflictions. Such results have not heretofore been achieved in the art.

The methods and apparatus disclosed and claimed herein are implementable with conventional imaging equipment. The invention can thus be practiced on standard MRI machines having either serial or array processors for data manipulation. Post processing techniques employed in accordance with the invention produce high conspicuity of vessels, and other bright objects in clinical images which are particularly suited to the diagnosis of diseases and maladies which afflict the body. In a further preferred embodiment, MRA is performed by methods and apparatus provided in accordance with the present invention so that highly conspicuous images of blood vessels are obtained to diagnose cardiac, cerebral, pulmonary, and peripheral vascular afflictions of the limbs. Such results have also not heretofore been achieved in the art. The present invention will be better understood by reading the following detailed description of preferred embodiments in conjunction with the drawings.

Brief Description of the Drawings Figures 1A and IB are schematic diagrams of image projection using a ray tracing spatial assignment step and a maximum intensity feature extraction step wherein three- dimensional data is projected onto a two-dimensional lateral projection image. Figure 2 is a schematic showing the manner in which cross-sectional images, or "slices", are collected sequentially in order to represent a body part (in this case, a lung) three- dimensionally.

Figure 3 is a block diagram of a computer system for controlling the acquisition of data and for processing such data to produce cross-sectional images implementing methods provided in accordance with the present invention.

Figure 4 is a diagram of the relationship of a "test" point on a ray and its relationship to its. nearest neighbors in the lattice of nodes.

Figures 5A and 5B are histograms of an MR Angiogram. Figures 6A and 6B are block diagrams of an MRI system and a CAT scan system respectively showing the relationship of the computer system in Figure 3 to these respective imaging system architectures.

Figure 7 illustrates a sigmoidal curve adequate for the implementation of soft thresholding in accordance with the invention Figure 8 is a "typical histogram" computed for the actual histograms of Figures 5A and 5B, and demonstrating the desired properties of the typical histogram.

Figures 9A - 9E present data supporting the use of an automated algorithm for choosing Window and the Level parameters in accordance with the present invention.

Figure 10A and 10B are schematics showing the result of applying soft thresholding to a plot of intensities plotted along a ray passing through a three-dimensional, cross-sectional data set.

Figures 11A - 11D illustrate the enhancement provided by the additional post-processing techniques of depth cuing and of summation after soft thresholding in accordance with the invention.

Detailed Description of Preferred Embodiments

Referring now to the drawings wherein like reference numerals refer to like elements as discussed above with respect to Figure 3, gathering means 20 is provided for gathering data 50 from an object 30 and for creating a three-dimensional image data set 80 which can be further manipulated by the system to obtain useful information 90 concerning object 30. Gathering means 20 could simply be an antenna of appropriate design for receiving reflected or stimulated energy from object 30. It will be further recognized by those with skill in the art that the electromagnetic energy received by gathering means 20 is continuous or "analog" in nature, and therefore in a preferred embodiment, an analog-to-digital (A/D) converter 60 is interfaced to gathering means 20 to digitize the gathered analog electromagnetic energy and to create digitized data for further use by the system to obtain an image. Other embodiments may utilize similar computer configurations with alternative detector methods for collection of signals of other types. For^' example, a Computed Axial Tomography (CAT) Scanner uses an X-Ray source and detector system to digitize data from which cross- sectional images will be constructed based on X-ray density with an almost identical computer.

After the data is digitized, the raw data 50 is saved in intermediate storage means 70 until it is reconstructed by an array processor 10 into a set of images 80 which are then stored on some final image storage device 70. Images stored on device 70 can then be: (1) viewed on a computer monitor, (2) archived on some computer medium such as optical disk or magnetic tape, or (3) filmed for review by the radiologist.

In a preferred embodiment of this invention, the images 80 in storage device 70 are directed to a high speed computer subsystem such as array processor 10. The image data is then reprocessed and redirected to the image storage device 70 from where it can be viewed, archived and/or filmed. In a further preferred embodiment, the array processing means 10 performs ray tracing on the data set as heretofore discussed. While ray tracing is preferably used in accordance with the invention, virtually any type of graphics acceleration imaging technique can be utilized to project an image of object 30, and to further create a projected data set 90 which can be manipulated to obtain useful information about the object.

In still a further preferred embodiment, spatial assignment and feature extraction as described with respect to Figures 1A and IB are performed by the array processing means 10 to extract image features from the projected data set. This is preferably accomplished by estimating thresholds from analysis of histogram data in sub-volumes of the total three-dimensional volume, and by thresholding the projected data to at least two predetermined levels. Further features are extracted from the data set after the data set has been thresholded to at least two predetermined values. The data can be output as an image 90 for storage and display. The image 90 can then be viewed by a radiologist or other clinician to diagnose the problems that affect object 30.

Referring to Figures 6A and 6B, both MRI and CAT X-ray systems in accordance with the present invention are particular¬ ly useful for performing non-invasive angiography. In each case, a detector system is interfaced to digitization elec¬ tronics which are adapted to receive data from the blood vessels when a patient is placed inside the detector system. For the MRI system of Figure 6A, a magnet 140 is interfaced to pulse programming electronics 150 which initiates a sequence to control the magnet 140 and to set the radiofrequency fields and magnetic field gradients which will induce the blood vessels and fluid flowing therein to emit radio signals for analysis. Digital data from receiver and A/D electronics block 160 as well as data from pulse programmer 150 are then output along a data communications highway or bus 170 so that the data can be communicated to a computer system which contains the appropriate software to control the process that takes place in the MRI system of Figure 3, for example.

In the case of an X-ray CAT scanner, an X-ray source 170 and an array 180 of X-ray detectors are positioned on opposite sides of the object of interest 30. The array 180 is mounted on a motorized, rotatable ring apparatus 190 which allows rotation of the X-ray source 170 and array 180 about the object of interest. The signal from the detectors is digitized by an analog-to-digital system 200 and then transported across a data communications highway or bus 210 so that the data can be communicated to a computer system which contains the appropriate software to control the process that takes place in the X-ray CAT system of Figure 3.

In a further preferred embodiment, the data received from the magnet 140 through receiver and A/D electronics block 160 which has been digitized can then be stored on a storage means 40 which is adapted to store the digitized data. Storage means 40 is, for example, a hard disk or a large Random Access Memory (RAM) bank which can store many runs of data from several different patients.

The data storage device 40 is further preferably interfaced to array processor 10 which is also interfaced with the master computer 100. The array processor 10 in accordance with the invention provides image projection and feature extraction of the raw data so that images of the blood vessels can be obtained. An array processor is a high speed, vector processing computer which is particularly well suited to process graphics data to obtain images. The spatial assignment and feature extraction functions performed by array processor 10 are preferably accomplished by software packages which program the array processor to perform the spatial assignment and feature extraction tasks. While it is possible to perform both image projection and feature extraction in accordance with the present invention with a scalar processor, the array processor is a preferred device for processing and post-processing MRA images with methods provided in accordance with the present invention. In particular, the Floating Point Systems array processor on a General Electric Signa 3.x MRI machine has been found to be an excellent tool for performing feature extraction in accordance with the invention. The array processor is designed to act on lists of numbers, and as such is well suited to manipulating line segment geometries such as are associated with ray tracing. Because all commercial scanners have array processors, algorithms adapted to these high speed processors are particularly appropriate to current scanner architectures. As discussed above, in prior generic thresholding operations for feature extraction, thresholding assigns a "zero" or a "one" for every pixel intensity below or above the threshold value respectively. This is a harsh criterion, and the inventors of the subject matter herein claimed and disclosed have "softened" this criterion to obtain "soft thresholding" for feature extraction.

With soft thresholding as taught herein, at least two intensity ranges are created such that the dynamic values of high and low intensity bands corresponding to pixel intensities are compressed in these ranges. In yet a further preferred embodiment of soft thresholding in accordance with the invention, an intermediate range is also created which is not characterized on the basis of signal intensity alone. The intermediate range is defined herein as an "interpretable" intensity range in which the data is expanded for analysis purposes. By removing the gross inhomogeneities in signal intensities, the fine contrast necessary to distinguish delicate structures in MRA is accomplished. The data compression performed by soft thresholding in accordance with the invention is preferably accomplished by applying a monotonically increasing function to the data which is approximately linear over a range of "interpretable" intensities corresponding to the fine structures observed in MRA, and suppressing differences outside of this range. The monotonically increasing function is more preferably applied by the feature extraction block techniques described on the data received after spatial assignment.

More preferably, soft thresholding in accordance with the present invention is accomplished by applying a sigmoidal- shaped curve to the data. In a more preferred embodiment, three ranges of signal intensity are defined, either interactively by the operator or by a threshold prediction algorithm described below. These are a signal intensity corresponding to definite flow, a signal intensity corresponding to an interpretable range wherein flow may or may not be occurring within the vessel, and a signal intensity corresponding to definite, stationary tissue. In this way, the controlling parameters in the sigmoid function are analogous to "Window" and "Level" controls on a standard gray scale monitor. The ranges of definite signal interpretation are thus compressed in soft thresholding since it can be assured that flow is or is not occurring in these areas, while the intermediate "interpretable" range is expanded to determine the fine structure associated with these intermediate^' signal intensities.

In principle, any sigmoid-type curve can be used. Preferably, an arctangent function is applied to the projected data to extract the information. The arctangent function follows the form:

outpu = arctan { (I^_put - Level) / Window }+ 1/2 .

In the arctangent equation expressed above, the signal intensity output, I_output , is a function of the input intensity, li_nut_* and "Level" and "Window" constants. These constants are chosen automatically by an algorithm described below, or by the technologist, diagnostician or clinician in accordance with the particular contrast desired for the lateral projection MRA images.

Referring to Figure 7, the arctangent intensity output function most preferably used in accordance with the invention to perform soft thresholding is shown. The width of the linear region is defined by the Window value. When I_input is greater than the Level + Window/2 level value, that is, when there is high intensity corresponding to definite flow, I_output approaches 1 and these high intensity values are compressed. Similarly, when li_nu is less than the Level - Windows/2, that is a stationary tissue outputs low intensity energy, I_output is close to zero and low intensity values are compressed. When I^_put is in the neighborhood of the Level, the linear range of the arctangent sigmoidal function is reached and the interpretable range is expanded, thereby increasing the conspicuity of any features with intensity outputs falling within this interpretable range.

In a further preferred embodiment, the Window and Level parameters may be chosen by the following algorithm, which is based on the characteristics of the histogram of the cross- sectional MRA image described previously. Referring to Figure 8, a "typical histogram" is computed for the intersection of a cross-sectional image with the segmented sub-volume. In a preferred embodiment, a histogram is computed for the intersection of each cross-sectional image with the segmented sub-volume, and then an average histogram is computed over all such individual histograms. Finally, to produce an adequate "typical histogram", this average histogram is smoothed by convolution with a simple low pass filter of width approximately 4% of the total dynamic range of intensities in the average histogram.

This "typical histogram" plots the average bin occupancy in a slice intersecting the sub-volume, and has the following properties: (1) it is smooth, (2) it is monotonically decreasing in the range of intensities between definite stationary and definite flow, (3) the average bin occupancy in the region of definite flow is less than one, and (4) the average bin occupancy in the region of definite stationary signal is greater than one. In a further preferred embodiment, for a choice of Window equal to 2% to 4% of the dynamic range of intensities in the "typical histogram", the Level is chosen to be that intensity at which the average bin occupancy plotted in the "typical histogram" exactly equals one. In Figure 8, a series of simple triangular, low-pass filters with half-width at half-maximum (HWHM) varying between 10, 40, 80, 120 and 160 intensity units was convolved with the histograms of Figures 5A and 5B to output the series of "average bin occupancy" curves. Relative to the total dynamic range of the histogram, this range of HWHM represents about 0.5%, 2%, 4%, 8% and 12% of the maximum image intensity range respectively.

Referring to Figures 9A - 9E, it can be seen that this automatic algorithm is an excellent predictor of the optimal Level parameter as chosen by three experienced readers of lateral projections of MRA. The quantities noted at the bottom of each of the graphs of Figures 9A - 9E are the Window widths and range of optimum values for each reader given as percentages of the maximum pixel value within the volume of interest. Three experienced readers viewed a series of images projected by three-dimensional MRA data. Each series of images consisted of multiple projections of multiple views of the data set in which the Level parameter was varied. The readers were asked to identify the maximum Level for which the lateral projection image contained too much noise (the low bound of acceptable range of Levels) and the minimum Level for which the lateral projection image contained loss of important detail (the upper bound of acceptable range of Levels) .

In two cases, the estimate fell within the optimal range of at least one reader, in another two cases within the optimal range of two readers, and in the last case within the optimal range of all three readers. It is notable that the optimal ranges for the readers themselves does not always intersect.

Soft thresholding in this fashion suppresses intravessel variation, thereby increasing the accuracy of perception of the diameter of the signal present within the lumen of the blood vessel. Thus, the apparent diameter of a vessel will be increased both on the display means, and during and after photography. Referring to Figures 10A and 10B, graphs of signal intensities before and after soft thresholding respectively illustrate how soft thresholding in accordance with the present invention enhances the accuracy of perception of a vessel and extracts features of interest regarding the vessel while suppressing unwanted information. Figure 10A shows a large vessel 220, a small vessel 230, and an exemplary noise spike 240. Stationary tissue is also shown generally at 250. As can be seen from Figure 10A, in the definite flow range, that is the upper limit of the soft thresholded range, the large vessel 220 has an intense pixel brightness. The small vessel 230 lies in both the interpretable range and the stationary tissue range, that is the range of medium intensity and low intensity, so that it is not clear from the non-soft thresholded data whether these features are of importance. The noise signal 240 is also in the interpretable range, while the stationary tissue 250 lies completely in the low intensity range.

Figure 10B shows the results of applying soft thresholding. Large vessel signal variation within the definite flow range has been compressed 260 to an almost uniform intensity. The stationary tissue within the low intensity range has also been compressed at 270 to values almost equal to zero. The small vessel now lies in an expanded interpretable range which allows features 280 in the small vessel to be more easily appreciated. The noise signals 290 are also expanded in the interpretable range, where they can be more closely examined for diagnostic significance.

In this manner, soft thresholding in accordance with the invention involving setting two pixel intensity levels results in the expansion of the mid-range pixel intensity so that intensities which exist in the mid-range can then be more accurately surveyed to examine important features in these structures. This amounts to a "stretching" of the mid-range so that dynamic range expansion of mid-range pixel intensities produces an image which allows a clinician to more closely examine features in the mid-range, wherein known flow and known stationary pixels are uniformly bright or dark, respectively.

As a preliminary step, soft thresholding may be used to enhance the other features extraction methods and techniques to obtain good clinical images. Referring now to Figures 11A - 11D, once soft thresholding has been accomplished in accordance with the present invention, it is preferred to obtain images with "three-dimensional depth cues" by first applying standard artificial lighting followed by MIP. Because flow structures are of almost uniform weighting after soft thresholding, weighting by an artificial light will cause structures closer to the viewer to be brighter. A subsequent MIP will cause near objects to be superimposed on top of more distant, darker objects as shown in Figure 11A. A similar procedure performed on non-thresholded data will not reliably superimpose near structures on far structures, due to the flow induced variation in signal intensity as shown in Figure 11B. The attainment of good quality depth cues reduces the need for computation of multiple views, and for segmentation of the original three dimensional volume into several sub-volumes.

Alternatively, a preferred implementation of "X-ray density"-like images may be obtained after soft thresholding by integration or summation of intensity along the projected rays. Summation of intensity after soft thresholding (hereinafter referred to as "SIST") has advantages over direct summation of vessel intensities in MRA because intra-vessel variation is reduced, thereby making the value of the summation more reflective of vessel diameter as illustrated in Figure llC. Before soft thresholding in Figure 11D, it is unclear whether summation of intensities would produce meaningful density based on the length of a ray's intersection with attenuating structures. Consequently, SIST produces an improvement in the "double-density sign" as described previously. Soft thresholding in accordance with the invention provides superior feature extraction for MRA specifically, and for MRI generally. When coupled with other feature extraction techniques, soft thresholding ensures that three-dimensional depth cues are efficiently and accurately obtained. These results have not been achieved with prior feature extraction techniques. Thus, soft thresholding in accordance with the invention provides an excellent clinical and diagnostic tool to obtain highly resolved images in MRA.

There have thus been described certain preferred embodiments of post-processing of MRI images provided in accordance with the present invention. While preferred embodiments have been described and disclosed, it will be recognized by those with skill in the art that modifications are within the true spirit and scope of the invention. The appended claims are intended to cover all such modifications.

Claims

- 3 CT -CLAIMSWhat is claimed is:

1. A method of imaging an object comprising the steps of: gathering data from the object and storing the data in a storage device; performing ray tracing on a data set which is representative of an object that has been scanned to obtain a substantially three-dimensional image of the object and a three-dimensional data set corresponding to the three dimensional image of the object; and thresholding the three-dimensional data set according to at least two levels to obtain a thresholded data set and to project the three- dimensional image, thereby extracting features of interest from the three-dimensional image.

2. The method recited in claim 1 wherein the thresholding step further comprises: further thresholding the three-dimensional data set to an intermediate range which is an interpretable range that expands the three-dimensional data set so that features of interest from the three-dimensional image can be extracted.

3. The method recited in claim 2 wherein the thresholding step comprises applying a monotonically increasing function to the three-dimensional data set which is approximately linear over a defined range.

4. The method recited in claim 3 wherein the monotonically increasing function is a sigmoidal shaped curved.

5. The method recited in claim 4 wherein the sigmoidal function is an arctangent function having the form:

•∑outut = arctan { (I^p_ut - Level)/Window },

wherein I_output is an output intensity, the "Level" is a constant corresponding to an initial intensity, the "Window" corresponds to an intensity range, and I_input is an input intensity.

6. The method recited in claim 5 wherein the at least two levels are signal ranges wherein the three-dimensional set is compressed, and wherein the interpretable range is a range in which signal intensity is expanded so that intermediate signal intensities can be extracted.

7. The method recited in claim 6 further comprising the step of further extracting features of interest from the three-dimensional image by applying artificial lighting to the thresholded data.

8. The method recited in claim 6 wherein the thresholded data is summed to extract features of interest after the data has been thresholded.

9. The method recited in claim 6 further comprising the step of applying maximum intensity projection to the thresholded data after the data has been thresholded.

10. A magnetic resonance imaging system comprising: gathering means for gathering data from an object and creating a data set that can be further manipulated by the system to obtain useful information concerning the object; processing means interfaced to the gathering means and adapted to receive the data set for projecting an image of the object based on the data set and creating a projected data set which can be further manipulated to obtain additional useful information concerning the object; and extracting means interfaced to the processing means and adapted to receive the projected data set for extracting image features from the projected data set by thresholding the projected data set to at least two predetermined levels.

11. The magnetic resonance imaging system recited in claim 10 further comprising output means interfaced to the extracting means for outputting an image of the object.

12. The magnetic resonance imaging system recited in claim 11 wherein the extracting means further thresholds the data to an intermediate range in which the projected data is expanded so that features can be extracted from the projected data set.

13. The magnetic resonance imaging system recited in claim 12 wherein the thresholded data is obtained by applying a monotonically increasing function to the data set which is approximately linear over a range of intensities, and suppressing differences outside of the range.

14. The magnetic resonance imaging system recited in claim 13 wherein the monotonically increasing function is a sigmoidal shaped curve.

15. The magnetic resonance imaging system recited in claim 14 wherein the sigmoidal shaped curve is an arctangent function having the form:

l_outp_ut = arctan { (I_input - Level)/Window },

wherein I_output is an output intensity, the "Level" is a constant corresponding to an initial intensity, the "Window" corresponds to an intensity range, and Imp_ut. s an input intensity.

16. The magnetic resonance imaging system recited in claim 15 wherein the processing means is an array processor.

17. The magnetic resonance imaging system recited in claim 16 further comprising a magnet which emits a sinusoidal magnetic field that stimulates the object to emit radiation so that the data set can be gathered.

18. A method of obtaining images of an object with a magnetic resonance imaging machine comprising the steps of: gathering data about the object and storing the data in a memory device for manipulation by a processing block in the magnetic resonance imaging machine; manipulating the data stored in the memory to obtain a projected image of the object and a projected image of the data set of the object; storing the projected image data set in a memory for further processing of the projected image data set; extracting features from the projected image data set by thresholding the projecting image data set to at least two predetermined values such that data in the projected image data set which are linear are retained and non-linear data are suppressed, thereby obtaining an extracted data set; and outputting a substantially three-dimensional image of the object by manipulating the extracted data set to form the three-dimensional image.

19. The method recited in claim 18 wherein the extracting step further comprises the step of: further thresholding the projected image data set to a third value, thereby creating an interpretable intensity range which is expanded for further extraction of the data set.

20. The method recited in claim 19 wherein the extracting step comprises the step of applying the monotonically increasing function to the projected image data to expand the projected image data which is in a linear range, and to suppress non-linear data which are outside of the linear range.

21. The method recited in claim 20 wherein the monotonically increasing function is a sigmoidal shaped curve.

22. The method recited in claim 21 wherein the sigmoidal shaped curve is an arctangent curve.

23. The method recited in claim 22 further comprising the steps of applying graphics extraction to the extracted data set to obtain images of the object.