Automated medical image visualization using volume rendering with local histograms

[0001] The present invention relates to medical renderings of imaging data.

[0002] A portion of the disclosure of this patent document contains material to which a claim of copyright protection is made. The copyright owner has no objection to the facsimile or reproduction by anyone of the patent document or the patent disclosure, as it appears in the Patent and Trademark Office patent file or records, but reserves all other rights whatsoever.

[0003] Direct Volume Rendering (“DVR”) has been used in medical visualization research for a number of years. DVR can be generally described as rendering visual images directly from volume data without relying on graphic constructs of boundaries and surfaces thereby providing a fuller visualization of internal structures from 3-D data. DVR holds promise for its diagnostic potential in analyzing medical image volumes. Slice-by-slice viewing of medical data may be increasingly difficult for the large data sets now provided by imaging modalities raising issues of information and data overload and clinical feasibility with current radiology staffing levels. See, e.g., Adressing the Coming Radiology Crisis: The Society for Computer Applications in Radiology Transforming the Radiological Interpretation Process (TRIP™) Initiative, Andriole et al., at URL scarnet.net/trip/pdf/TRIP_White_Paper.pdf (November 2003). In some modalities, patient data sets can have large volumes, such as greater than 1 gigabyte, and can even commonly exceed 10's or 100's of gigabytes.

[0004] Despite its potential, DVR has not achieved widespread use for non-research medical imaging, particularly in computer network systems with visualization pipelines. This may be because DVR may need time-consuming manual adjustment using conventional transfer functions (TF) and/or editing tools. That is, the TF construction can be relatively complex and/or the tissue separation abilities may not be sufficient where dissimilar tissues have similar intensity values limiting the ability to generate diagnostic clinical renderings.

[0005] Embodiments of the present invention are directed to methods, systems and computer program products that employ local histogram data that is used for direct volume renderings.

[0006] Some methods include evaluating data associated with direct volume renderings by automatically electronically analyzing a plurality of local histograms derived from a volume rendering data set, and programmatically generating data that is used for at least one of (a) material classification or (b) material detection, of materials in the data set, including materials of similar intensity, based on the analyzing step.

[0007] The methods may include distinguishing different materials with overlapping intensity values using range weight data values obtained from the electronic analysis of the local histograms. In particular embodiments, the methods may include electronically employing a statistical model of classification probability to classify whether a voxel in the volume rendering data set belongs to a feature and/or material of interest. The statistical model of classification probability may be based on weighted inputs of confidence for selected neighborhood characterizing features.

[0008] Other methods are directed to visualizing images of volume data sets in visualization systems. The methods include iteratively electronically subdividing a global histogram of a respective volume data set using partial local histograms having associated value ranges and intensity values to electronically automatically identify different materials, including materials having overlapping image intensity values.

[0009] In some embodiments, the methods can include adapting at least one a priori transfer function to evaluate different volume rendering data sets of similar examination types using programmatically generated partial range histograms of the local histogram data.

[0010] Other embodiments are directed to methods for providing a tissue exploration tool to allow a physician to interactively analyze medical volume data sets in a visualization system. The methods include: (a) allowing a user to electronically select a partial range of interest in an intensity scale of voxel data to thereby allow the user to interactively investigate voxels in a volume rendering data set; (b) electronically generating a partial range histogram associated with the selected partial range; and (c) electronically fitting an adaptive trapezoid to the partial range histogram to render an image of material associated with the selected partial range.

[0011] In some embodiments, the methods may also include disconnecting the trapezoid from a transfer function while a physician investigates different voxels.

[0012] Still other embodiments are directed to methods for visualizing images of volume data sets. The methods include: (a) iteratively electronically subdividing a respective volume data set using local histograms of neighborhood voxel data to allocate the neighborhoods of data into partial range histograms, each partial range histogram having voxels with intensity values in a defined range; and (b) electronically automatically identifying different materials in the volume data set, including materials having overlapping image intensity values which may be distributed over a target volume of interest.

[0013] In particular embodiments, the local neighborhoods are selected, configured and sized so that the volume data set can be analyzed with non-overlapping subdivision of the voxel data.

[0014] Other embodiments are directed to systems for generating DVR medical images. The systems include a volume rendering medical image processor system configured to generate data for a diagnostic medical image of a target region of a patient by electronically subdividing a volume rendering data set using local histogram analysis to separate different tissues with overlapping image intensity values.

[0015] In some embodiments, the image processor system may be configured to automatically electronically: (a) identify range weight values of different tissue associated with the patient medical volume data set using the local histogram analysis: (b) generate partial range histograms based on the local histogram analysis: and (c) generate adaptive trapezoids using the partial range histograms, to thereby generate the diagnostic medical image.

[0016] Still other embodiments are directed to computer program products for providing physician interactive access to patient medical volume data for rendering diagnostic medical images using a computer network. The computer program product includes a computer readable storage medium having computer readable program code embodied in the medium. The computer-readable program code includes computer readable program code configured to generate partial range histograms having associated peak characteristics and intensity values to electronically identify different types of tissue having overlapping image intensity values to thereby provide data used to render a diagnostic medical image of a target region of a patient.

[0017] It is noted that any of the features claimed with respect to one type of claim, such as a system, apparatus, or computer program, may be claimed or carried out as any of the other types of claimed operations or features.

[0018] Further features, advantages and details of the present invention will be appreciated by those of ordinary skill in the art from a reading of the figures and the detailed description of the preferred embodiments that follow, such description being merely illustrative of the present invention.

[0019] FIG. 1 is a global histogram of number of voxels versus voxel values of a data set of MR renal arteries (full and magnified).

[0020] FIG. 2 is a block diagram of an electronic visualization pipeline according to embodiments of the present invention.

[0021] FIG. 3 is a block diagram of a rendering system with a DVR image processor system according to embodiments of the present invention.

[0022] FIG. 4 is a flow chart of operations that may be used to visualize volume data in diagnostic medical images according to embodiments of the present invention.

[0023] FIG. 5A is a flow chart of operations that may be used to visualize volume data in diagnostic medical images according to embodiments of the present invention.

[0024] FIG. 5B is a flow chart of operations that may be used to visualize volume data in diagnostic medical images according to embodiments of the present invention.

[0025] FIG. 5C is a flow chart of operations that may be used to provide an interactive electronic tissue evaluation tool according to embodiments of the present invention.

[0026] FIG. 5D is a flow chart of operations that may be used to selectively visualize a feature in an image according to embodiments of the present invention.

[0027] FIG. 6 is a schematic illustration of a data processing system according to embodiments of the present invention.

[0028] FIG. 7A is a full histogram and a partial range histogram of number of voxels versus voxel values with a Gaussian curve and an adaptive trapezoid having a transfer function with an opacity scale according to embodiments of the present invention.

[0029] FIG. 7B is a color rendering generated using the local histogram data shown in FIG. 7A.

[0030] FIG. 8 is a graph of classification probability shown as a function of range weight according to embodiments of the present invention.

[0031] FIG. 9 illustrates a 2-D transfer function (TF), where a conventional 1-D TF is extended with a material classification dimension according to embodiments of the present invention.

[0032] FIG. 10 is a graph of a full resolution histogram and a corresponding simplified histogram of number of voxels versus voxel values according to embodiments of the present invention.

[0033] FIG. 11A is a graph of an MRI bilary duct data set with a full histogram, partial range histograms and histograms of certain anatomical features according to embodiments of the present invention.

[0034] FIG. 11B is a graph of an MRI renal arteries data set with a full histogram, partial range histograms, and a histogram of a vessel of number of voxels versus voxel values according to embodiments of the present invention.

[0035] FIG. 12A is a medical image rendering illustrating spongy bone and vessels using a conventional 1-D TF. FIG. 12B is magnified view of the vessels and spongy bone shown in the image in FIG. 12A.

[0036] FIG. 13A is a medical image rendering illustrating spongy bone and vessels using voxel-centric local neighborhoods according to embodiments of the present invention.

[0037] FIG. 13B is magnified view of the vessels and spongy bone shown in the image in FIG. 13A.

[0038] FIG. 14A is a medical image rendering illustrating spongy bone and vessels using double-block local neighborhoods according to embodiments of the present invention.

[0039] FIG. 14B is magnified view of the vessels and spongy bone shown in the image in FIG. 14A.

[0040] FIG. 15A is a medical image rendering of the heart showing the coronary artery using a conventional 1-D TF.

[0041] FIG. 15B is a magnified view of the coronary artery shown in FIG. 15A.

[0042] FIG. 16A is a medical image rendering of the heart showing the coronary artery using voxel-centric local neighborhoods according to embodiments of the present invention.

[0043] FIG. 16B is magnified view of the vessels and spongy bone shown in the image in FIG. 16A.

[0044] FIG. 17A is a medical image rendering of the heart showing the coronary artery using double-block local neighborhoods according to embodiments of the present invention.

[0045] FIG. 17B is magnified view of the vessels and spongy bone shown in the image in FIG. 17A.

[0046] FIGS. 18A-18C are medical image renderings of a liver with a tumor. FIG. 18A was obtained using a conventional 1D TF rendering model.

[0047] FIG. 18B is an image that was altered to attempt to highlight the tumor also using a 1D TF. FIG. 18C illustrates the image rendered using a 2-D TF with local neighborhood data according to embodiments of the present invention.

[0048] FIG. 19 is a graph of opacity versus intensity overlaying a color gradation map illustrating an example of a transfer function where color and opacity mapping for an intensity value are decoupled, continuous functions.

[0049] FIG. 20A is a graph of voxels versus intensity of a global histogram with an overlying electronic user-selectable partial range tool according to embodiments of the present invention.

[0050] FIG. 20B is a graph of voxel opacity versus intensity of adaptive trapezoid transfer functions created based on the user selectable partial range tool shown in FIG. 20A according to embodiments of the present invention.

[0051] The present invention now is described more fully hereinafter with reference to the accompanying drawings, in which embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

[0052] Like numbers refer to like elements throughout. In the figures, the thickness of certain lines, layers, components, elements or features may be exaggerated for clarity. Broken lines illustrate optional features or operations unless specified otherwise. In the claims, the claimed methods are not limited to the order of any steps recited unless so stated thereat.

[0053] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. As used herein, the term “and/or” includes any and all combinations of one or more of the associated listed items. As used herein, phrases such as “between X and Y” and “between about X and Y” should be interpreted to include X and Y. As used herein, phrases such as “between about X and Y” mean “between about X and about Y.” As used herein, phrases such as “from about X to Y” mean “from about X to about Y.”

[0054] Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the specification and relevant art and should not be interpreted in an idealized or overly formal sense unless expressly so defined herein. Well-known functions or constructions may not be described in detail for brevity and/or clarity.

[0055] It will be understood that, although the terms first, second, etc. may be used herein to describe various elements, components, regions, layers and/or sections, these elements, components, regions, layers and/or sections should not be limited by these terms. These terms are only used to distinguish one element, component, region, layer or section from another region, layer or section. Thus, a first element, component, region, layer or section discussed below could be termed a second element, component, region, layer or section without departing from the teachings of the present invention. The sequence of operations (or steps) is not limited to the order presented in the claims or figures unless specifically indicated otherwise.

[0056] The term “Direct Volume Rendering” or DVR is well known to those of skill in the art. DVR comprises electronically rendering a medical image directly from volumetric data sets to thereby display color visualizations of internal structures using 3D data. In contrast to conventional iso-surface graphic constructs, DVR does not require the use of intermediate graphic constructs (such as polygons or triangles) to represent objects, surfaces and/or boundaries. However, DVR can use mathematical models to classify certain structures and can use graphic constructs.

[0057] The term “global histogram” means a main histogram of a volume rendering data set of numbers of voxels versus voxel values that represents a target region of an image, such as of a portion of a patient body, that can be subdivided into local histograms. The term “local histogram” is a histogram for a single local neighborhood. The term “local neighborhood” means a related region of voxel intensity data. The related region of data may be of any desired shape or size and number of voxels. The local neighborhoods may be configured as arbitrary non-overlapping subdivisions of voxel data from the volume rendering data set. The local neighborhoods may have the same shape and same size or different shapes and sizes. The local neighborhoods may be configured as arbitrary non-overlapping subdivisions of voxel data. In some particular embodiments, each local neighborhood may contain at least about 50 voxels. In some examples, the local neighborhood can be cubical blocks that are about at least 8 voxels wide. The term “partial range histograms” means histograms that are populated by one or more (the sum of) selected local histograms in the data set in a partial intensity range that may provide a relatively large foot print in the partial intensity range and may include blocks distributed through a volume (they are not required to have a common locality). A region (such as block) associated a local neighborhood may be selected to belong to a set of partial range histograms if the region has a sufficient number of voxels within a given range weight in a partial range. Although illustrated herein for discussion purposes as graphs when interrogating volume data, the histogram is not required to be generated in a graphic format and can be stored, manipulated or interrogated as a numerical data set.

[0058] The term “peak characteristic” describes one or more parameters associated with a peak lineshape (such as, but not limited to, height, width, average width, median location, volume, or other measure) associated with a peak or peaks in a histogram. The term “automatically” means that the operation can be substantially, and typically entirely, carried out without human or manual input, and is typically programmatically directed or carried out. The term “electronically” includes both wireless and wired connections between components. The term “adaptive” means that the shape and size of the geometric perimeter lineshape (typically a trapezoid) is automatically electronically fitted to a histogram based on selected parameters and/or characteristics of the local histogram. The term “tissue” means blood, cells, bone and the like. “Distinct or different tissue” or “distinct or different material” means tissue or material with dissimilar density or other structural or physically characteristic. For example, in medical images, different or distinct tissue or material can refer to tissue having biophysical characteristics different from other (local) tissue. Thus, a blood vessel and spongy bone may have overlapping intensity but are distinct tissue. In another example, a contrast agent can make tissue have a different density or appearance from blood or other tissue. The term “transfer function” means a mathematical conversion of volume data to image data. Visualization means to view in 3-D volume data representing features with different visual characteristics such as with differing opacity, color, texture and the like. The term “similar examination type” refers to corresponding anatomical regions or features in images having diagnostic or clinical interest in different data sets corresponding to different patients (or the same patient at a different time). For example, but not limited to, a coronary artery, organs, such as the liver, heart, kidneys, lungs, brain, and the like.

[0059] In the past, DVR has been carried out with the global histogram acting as the guiding tool used to find tissues of interest based on intensity scale. However, the target tissues of interest may be minor features, making their contribution almost, if not entirely, invisible in the global histogram. Generally stated, embodiments of the present invention can employ local histograms that can expose peaks, peak characteristics and/or value ranges that may be obscured or suppressed in the global histogram. The local histograms can be used to express properties in terms of histogram content in local neighborhoods. Embodiments of the present invention can employ data from local histogram analysis for either material (tissue) detection or material (tissue) classification, or both. Thus, the local histogram analysis can be employed for tissue classification without the tissue detection, tissue detection without tissue classification, or tissue detection and tissue classification together.

[0060] FIG. 1 is an intensity graph of an MR renal arteries data set showing a global (full) histogram (solid line with stippled fill), a magnified scale of the global histogram (dashed line with long segments separated by shorter segments), and arteries magnified (evenly broken line). The arteries are significantly, if not totally, obscured in the full histogram. As shown, in the global histogram, a single main peak that is associated with tissue/structure of lesser interest, dominates the minor peak of the arteries, the latter of which is the tissue of primary interest for certain evaluations.

[0061] Some embodiments of the present invention employ data derived from the local histograms (such as neighborhood characteristic data) to classify the voxels in an overlapping intensity range. The classification value can be used as a second attribute value for the voxels, which in connection with a multi-dimensional or multi-attribute (such as 2D) TF, can provide enhanced rendering. The second attribute or dimension may be constructed as a linear interpolation of tissue-specific trapezoids or other defined lineshape used to define one or more visual characteristics of the tissue (such as color-opacity), decreasing the complexity for a user. The TF's provided by embodiments of the present invention may be used to identify in vivo materials in an intensity scale, then apply a visual characteristic(s) of feature(s) to them (different color, opacity, density, etc . . . ). Embodiments of the invention in the specification will be discussed with respect to provide diagnostic medical images, but may be used in other non-medical DVR applications as well.

[0062] Some embodiments of the invention may be particularly suitable for processing MRI data (which may not have a calibrated intensity scale) wherein volume rendering may otherwise be difficult, particularly where there are different features of interest and/or different in vivo materials or tissues having overlapping intensity values.

[0063] As will be appreciated by one of skill in the art, embodiments of the invention may be embodied as a method, system, data processing system, or computer program product. Accordingly, the present invention may take the form of an entirely software embodiment or an embodiment combining software and hardware aspects, all generally referred to herein as a “circuit” or “module.” Furthermore, the present invention may take the form of a computer program product on a computer-usable storage medium having computer-usable program code embodied in the medium. Any suitable computer readable medium may be utilized including hard disks, CD-ROMs, optical storage devices, a transmission media such as those supporting the Internet or an intranet, or magnetic or other electronic storage devices.

[0064] Computer program code for carrying out operations of the present invention may be written in an object oriented programming language such as Java, Smalltalk or C++. However, the computer program code for carrying out operations of the present invention may also be written in conventional procedural programming languages, such as the “C” programming language or in a visually oriented programming environment, such as VisualBasic.

[0065] Certain of the program code may execute entirely on one or more of the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer. In the latter scenario, the remote computer may be connected to the user's computer through a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, some program code may execute on local computers and some program code may execute on one or more local and/or remote server. The communication can be done in real time or near real time or off-line using a volume data set provided from the imaging modality.

[0066] The invention is described in part below with reference to flowchart illustrations and/or block diagrams of methods, systems, computer program products and data and/or system architecture structures according to embodiments of the invention. It will be understood that each block of the illustrations, and/or combinations of blocks, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general-purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the block or blocks.

[0067] These computer program instructions may also be stored in a computer-readable memory or storage that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory or storage produce an article of manufacture including instruction means which implement the function/act specified in the block or blocks.

[0068] The computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions/acts specified in the block or blocks.

[0069] Turning now to FIG. 2, a visualization pipeline 10 is illustrated. As known to those of skill in the art and as shown by the broken line box 10b, the pipeline 10 can include a data capture circuit 15, a compression circuit 18, a storage circuit 20, a decompression circuit 22 and a rendering system 25. The visualization pipeline 10 can be in communication with an imaging modality 50 that electronically obtains respective volume data sets of patients and can electronically transfer the data sets to the data capture circuit 15. The imaging modality 50 can be any desirable modality such as, but not limited to, MRI, X-ray of any type, including, for example, CT (computed tomography) and fluoroscopy, ultrasound, and the like. The visualization system 10 may also operate to render images using data sets from more than one of these modalities. That is, the visualization system 10 may be configured to render images irrespective of the imaging modality data type (i.e., a common system may render images for both CT and MRI volume image data). In some embodiments, the system 10 may optionally combine image data sets generated from different imaging modalities 50 to generate a combination image for a patient.

[0070] As shown, the rendering system 22 may be in communication with a physician workstation 30 to allow user input and interactive collaboration of image rendering to give the physician the image views of the desired features in generally real time. The rendering system 25 can be configured to zoom, rotate, and otherwise translate to give the physician visualization of the patient data in numerous views, such as section, front, back, top, bottom, and perspective views. The rendering system 22 may be wholly or partially incorporated into the physician workstation 30, but is typically a remote or local module, component or circuit that can communicate with a plurality of physician workstations (not shown). The visualization system can employ a computer network and may be particularly suitable for clinical data exchange/transmission over an intranet.

[0071] FIG. 3 illustrates that the rendering system 25 can include a DVR image processor system 29. The image processor system 29 can include a digital signal processor and other circuit components that allow for collaborative user input 30i as discussed above. Thus, in operation, the image processor system 29 renders the visualization of the medical image using the medical image volume data, typically on a display at the physician workstation 30.

[0072] FIG. 4 illustrates exemplary operations that can be carried out to render images according to embodiments of the present invention. A plurality of local histograms derived from volume rendering data set can be automatically electronically analyzed (block 100). Data used for one or more of material classification or material detection, of materials in the data set including materials of similar intensity, can be programmatically generated (block 110).

[0073] A direct volume rendering diagnostic medical image can be programmatically generated based on the interrogating step (block 106). As shown, the (medical) volume data set can optionally comprise MRI data (block 105) or CT X-ray data (block 107).

[0074] In some embodiments, the analyzing local histograms can include automatically electronically detecting peak characteristics and/or range weight values of voxels associated with different materials in the local histograms. Each local histogram may be generated to have peak characteristics and an associated intensity value range associated therewith. The volume data set can be automatically electronically interrogated (analyzed) using iterative subdivision of the data set using local histograms to detect substantially all voxels with similar intensity and/or peaks and/or peak characteristics of distinct materials.

[0075] In some embodiments, user input can be obtained to define a priori domain data regarding at least one feature and/or material of interest in the image (block 108). In some embodiments, a statistical model of classification probability can be electronically employed to classify as to which feature or material type a voxel in the medical data set belongs, where there are overlapping intensity values (block 109).

[0076] FIG. 5A is a flow chart of operations that can be used for automatic or semi-automatic tissue-specific detection to obtain visualization of volumetric data (render images) according to embodiments of the present invention. Peak characteristics associated with a peak (typically the highest peak) of interest in the primary or global histogram is identified (block 120). A Gaussian curve can be fit to the shape of the histogram lineshape (block 122). The curve may be fit using a midpoint, height, and deviation value of the histogram. Local partial range histograms can be created based on data range weight values for a selected portion of the Gaussian curve (block 124) thereby revealing, detecting and/or generating a plurality of more narrow peaks associated with tissues in this portion of the data. The selected portion of the global histogram that is used to focus the local histogram analysis can be a medial portion. The partial range histogram blocks can be removed from the global histogram (block 126). That is, in some embodiments, as each local histogram is analyzed, the blocks associated with the respective local histogram can be electronically removed from the main histogram and a new peak for the adjusted main histogram can be generated. Operations describe in blocks 120, 122, 124 and 126 can be repeated until the global or main histogram is substantially, if not entirely, empty (block 128).

[0077] The term “Partial Range Histogram” (“PRH”) refers to a composite set of histograms for a set of neighborhoods that are typical for a given intensity range as discussed above. The neighborhoods can have an arbitrary shape, but typically fill the volume and are non-overlapping. The partial range histogram can be created in two ways, either through automated tissue detection (see, e.g., block 170 of FIG. 5B) or by allowing an operator or physician to manually define a partial range.

[0078] In some embodiments, PRHs can be based on cubic block neighborhoods. In this example, to select the blocks to be part of a PRH, a mathematical equation for a range weight w_r can be used to measure the size of the neighborhood footprint in the partial range as expressed in Equation (1). $\begin{array}{cc}{w}_r\left(\Phi ,N\right)=\frac{N\bigcap V_{\Phi }}{N}& \left(1\right)\end{array}$
N is an arbitrary voxel neighborhood, V_Φ is the set of voxels within a range Φ. |V| denotes the number of voxels in a set V. A block is added to the PRH if the range weight is sufficiently high, typically defined as where: w_r≧ε. The data set characteristics determine the threshold. In some embodiments, one can use ε=0.95 to capture narrow features and ε=0.5 to pick up inhomogeneous regions. In particular embodiments, a block size of 8³ can be used. The block size choice is a trade-off between the simplicity of larger blocks and the sensitivity for narrow features associated with smaller blocks. Note that a PRH does not need to find all blocks containing a material, since the position of the PRH is the primary input to TF construction (not the height).

[0079] Referring again to FIG. 5A, the partial range histograms can be compared and if peak characteristics of the partial range histograms (typically using Gaussian curves) are determined to be similar, such as the shapes have similar deviation and their means are relatively similar, the similar partial range histograms can be merged (block 132). A new Gaussian curve can be generated for merged curves to allow for multiple merges. In some embodiments, Gaussian curves can be generated for each PRH (whether merged or unmerged) (block 134). As such, a partial range histogram and associated Gaussian curve can be generated for each peak and/or peak characteristic that corresponds to different/distinct tissue material and/or feature. A transfer function can be generated using an adaptive selected lineshape associated with at least one visual characteristic can be generated for each partial range histogram using the respective Gaussian curve (block 136). The adaptive selected lineshape may be a color-opacity trapezoid.

[0080] In some particular embodiments, a fully automatic analysis of the global histogram using iterative subdivision of the histogram data (e.g., the volume data set) to detect all peaks corresponding to different materials using partial range histograms. The automated analysis steps may include:

[0081] 1. Find the highest peak of the main histogram.

[0082] 2. Fit a Gaussian to the main histogram using data associated with the peak.

[0083] 3. Create a PRH for a mid-portion of the Gaussian.

[0084] 4. Remove the PRH blocks from the main histogram.

[0085] 5. Run steps 1-4 until the main histogram is empty.

[0086] 6. Merge PRH's with peak characteristics (or peaks) that are similar.

[0087] Thus, in some embodiments, the analysis can start with volume data associated with the global or main histogram range. The highest peak is identified and a Gaussian curve is fitted to its shape. The Gaussian may be described by the midpoint μ, height ĥ, and deviation σ. The fitting process reduces the accumulated height difference between the histogram and the medial portion of the Gaussian. The error summation range is μ±ασ. In some embodiments, α=1. This small range may be selected because typically only the tip of the histogram peak is visible.

[0088] The next step is to create a PRH for a desired range. For example, the range can be selected as: μ±σ. This range choice yields many narrow, accurate peaks. Since the merging step follows, an abundance of peaks is not a problem. The blocks of the first PRH can be stored and/or removed from the main histogram. This exposes a new peak in the “main” histogram and the analysis algorithm repeats. If a PRH is empty, the ε value can be lowered and step 3 can be performed again. In this structured manner, all blocks can become part of a PRH. Then, PRH's having similar peak characteristics (e.g., similar peaks) can be merged.

[0089] To describe the PRH peak characteristics, a Gaussian can be fit to each PRH with the same method as described above. For the PRH Gaussian, the error summation may be adjusted. For example, α=2 may be employed because the entire peak is typically exposed. PRH's can be merged or joined if they have similar peak characteristic or peak deviation. In some embodiments, “similar deviation” can be defined as when σ_max/σ_min≦4, while their means are close. Their means may be considered close when, μ_max−μ_min≦σ_min·max(1.0, 2.0−σ_min/40.0). By defining the similar deviation in this manner, the second criterion is less strict for narrow peaks. For merged PRH's, a new Gaussian can be fitted to allow for multiple mergers. As used herein, the term PRH also refers to the merged PRHs.

[0090] It is noted that there are many suitable ways to construct a transfer function as known to those of skill in the art. In some examples herein, trapezoids are used to generate the transfer function that provides a graduated image intensity scale. The trapezoidal transfer functions may be particularly suitable for embodiments generating interactive tissue exploration tools for a physician as will be discussed further below. However, in some embodiments, particularly with respect to tissue detection, local histogram analysis can be used to generate (or adapt, see below) transfer functions of other types.

[0091] In some embodiments, the tissue detection can be seen as a calibration tool, independent of a transfer function concept, and can be compared to other calibration/standardization efforts. See e.g., On standardizing the MR Image Intensity Scale, Nyul & Udupa, Magnetic Resonance in Medicine, 1999 (which proposes a method based on global histogram comparisons). An example of another type of transfer function is when the color and opacity mapping for an intensity value are decoupled and generate continuous functions as shown graphically in FIG. 19. FIG. 19 illustrates opacity versus intensity, with colors varying in a general continuum from left to right across the graph (white, yellow, dark yellow, orange, red, pink, purple, blue, etc. with the colors intensifying from right to left within color segments).

[0092] In particular embodiments, the tissue detection methods of the present invention can be used to adapt, improve, calibrate, evaluate and/or extend existing manual transfer functions, rather than creating new ones “from scratch” as described in more detail herein under Example 3 for an MR angiography embodiment. In operation, manual or a priori transfer functions can be created or established with a minor training set (such as about three-ten cases of exemplary volume rendering data sets of a target region with similar materials and configurations of features, such as images of the same target region of a body). The automated tissue detection described herein can be applied on the training sets, and a relation between the detected tissues using the automated tissue detection methods and the manual and/or a priori transfer function can be derived. Then, the tissue detection can be applied to other, unknown, data sets of the same type. Knowing the relation, an appropriate transfer function (one with increased reliability) can be established for each individual anatomical feature or material of interest in an image.

[0093] In some embodiments, an adaptive trapezoid can be used to create a transfer function from a partial range histogram. This can be achieved via the fitting of a Gaussian to the partial range histogram as described above. However, as noted above, in other embodiments, other types of transfer functions can be employed.

[0094] FIG. 5B is a flow chart of operations that may be carried out to semi-automatically (with manual input) or automatically generate transfer functions to visualize the image volume data. Peak characteristics associated with material and/or features in the volume data set can be automatically detected using the iterative evaluation of the volume rendering data set using local histograms (block 170). Alternatively, or additionally, a partial range histogram can be created by manually defining a partial range of interest, along with or instead of automatically generating partial range histograms.

[0095] In some embodiments, the automated or manually assisted tissue detection (which may include peak characteristic detection) can employ one or more of those operations described in FIG. 5A. An adaptive trapezoid (or other selected transfer function shape) can be generated as a component that adapts its center, width and shape to the Gaussian approximation of a PRH as shown for example in FIG. 7A. FIG. 7A illustrates a PRH which is a portion of the global histogram data with a Gaussian curve fit over the PRH. A trapezoid can be fit over the Gaussian (defined by parameters of the Gaussian).

[0096] Referring again to FIG. 5B, the adaptive trapezoid can optionally be selectively passed over viewing areas of an image and then locked or affixed to engage selected tissue (block 173). A 1-D TF can be generated using the adaptive trapezoids (block 174). An image using a plurality of the generated trapezoids, each associated with distinct tissue of interest can be generated (block 176). The color and opacity of one or more trapezoids can be adjusted (block 177). The adjustment can be automatic or manual based on input from the physician. FIG. 7B is a rendering resulting from the local histogram data shown in FIG. 7A. A low opacity gray ramp was added to the TF for clarity.

[0097] In some cases, the neighborhood histogram data can be sufficient to generate the rendered image in that different surroundings of equal-intensity voxels can provide diversifying information. That is, local histogram analysis can be used to classify voxels in a type of morphological filtering and/or fuzzy classification. In some embodiments, a classification model can be used to define a second attribute for a 2-D TF (block 178) as will be discussed further below.

[0098] In some embodiments, adaptive trapezoids can allow an efficient workflow when manually defining or modifying a TF “on the fly” using the following operations, which allow for some physician input:

• • [0099] 1. Apply automatic tissue detection.
  • [0100] 2. Activate an adaptive trapezoid, during an initial stage of the analysis, the TF can consist of this single trapezoid.
  • [0101] 3. Browse the detected PRHs until the adapted trapezoid highlights a tissue or feature of interest.
  • [0102] 4. Fixate the trapezoid at this position, making it static and temporarily disconnected from the TF.
  • [0103] 5. Perform steps 2-4 until all tissues of interest are found.
  • [0104] 6. View the full TF, i.e. activate all defined trapezoids. Adjust the color and opacity for each trapezoid as desired.

[0105] In addition, adaptive trapezoids can be used as a semi-automatic or manual electronic tissue exploration tool for the user. In so doing, step 1 above can be omitted and step 3 can be replaced by manually sliding the center of a partial range across the intensity scale. For an example of operations that may be used for the physician interactive tool, see FIG. 5C and operations 180-188. For example, a partial range histogram can be defined (block 180), and an adaptive trapezoid associated with or connected to a partial range of interest can be generated (block 182). The display can be configured to allow a user to electronically investigate voxels of interest in the volume data set by sliding the partial range over an intensity scale (block 185). A rendered image can be viewed using a plurality of the adaptive trapezoids (which are active transfer functions), each typically associated with a distinct tissue and/or feature of interest (block 188).

[0106] In some embodiments, a user can electronically fixate or overlay a trapezoid over a tissue of interest in the display (block 186). In addition, the color and/or opacity of the trapezoid (the partial range histogram voxels) can be adjusted based on user input (block 187). In addition, as discussed above for operations shown in FIG. 5B, a 1-D transfer function can be generated using the adaptive trapezoids (block 183) and a classification model may be used to define a second attribute for a 2-D transfer function (block 184).

[0107] In some embodiments, selective application of adaptive trapezoids can be used to selectively view neighborhoods of data having particular interest as shown in FIG. 5D. In operation, an adaptive trapezoid can be connected to a partial range histogram (manually or automatically) (block 190). Typically the trapezoid/transfer function is applied globally to generally the whole volume of the rendered image. However, it is also possible to apply the trapezoid selectively as shown in block 192 to thereby enhance visualization of features of interest. This is because irrelevant blocks can be hidden from view in the rendered image (or visually differentiated or suppressed) even if they have voxel intensity contents in the interesting range.

[0108] FIGS. 20A and 20B illustrate graphically how a user can define a partial range of interest 193 from data associated with the voxel volume rendering data set (shown as the global histogram) 196. In this example, the user-defined range 193 is represented by the tall thin box with the medial vertical line overlying the global histogram 196. A partial range histogram 194 is electronically calculated for the user defined partial range 193 and an adaptive trapezoid 195 (FIG. 20B) corresponding to the partial range histogram 194 can be generated. The corresponding material, tissue or feature can be electronically automatically highlighted, visually enhanced otherwise visually distinguished in the rendering. The process can be repeated, i.e., the user can slide the partial range 193 to a different partial range of interest. A second partial range histogram can be calculated 197 (FIG. 20A). A second adaptive trapezoid 198 (FIG. 20B) associated with the second partial range histogram 197 can be generated. Again, the corresponding material, tissue and/or feature can be highlighted or otherwise visually enhanced or distinguished in the rendering.

[0109] Some embodiments of the present invention are configured to separate overlapping tissue ranges and provide for a truer visualization in the rendered image. That is, in the past, a common problem in medical DVR is in rendering tissues with overlapping intensity ranges. A typical example is CT angiographies, where vessels having a contrast agent have the same intensity as spongy bone. A conventional TF based only on intensity has difficulty in generating a rendering that separates the two tissues as shown in FIGS. 12A and 12B. In contrast, using the local histogram analysis can extend the DVR capabilities in this respect. Such an application of neighborhood operators to classify voxels can be considered to be a type of morphological filtering. See generally Maragos et al., Morphological systems for multi-dimensional signal processing, Proceedings of the IEEE 78, 4 (1990), 690-710.

[0110] In some embodiments, when the intensity data alone does not sufficiently separate materials, neighborhood histogram analysis data can be used to do so. This is because different surroundings of equal-intensity voxels can provide diversifying (classifying) information. One analysis type that can be used to classify the distinct tissue can be based on range weights as described above with respect to PRH.

[0111] In other embodiments, a classification probability model can be used with the local histogram analysis data to separate distinct tissue having substantially similar intensity voxels. This embodiment can be described referring to the example of a heart data set shown in FIGS. 15A, 15B, 16A, 16B and 17A, 17B. The classification model goal is to separate the coronary artery from the large vessels. The images shown in FIGS. 15A and 15B were obtained using a 1D TF. The brightness of the large vessels can impede the visualization of the coronary artery. The image shown in FIGS. 16A and 16B were obtained with a classifying 2D TF with a voxel-centric neighborhood. Additional discussion of exemplary “neighborhood” types is provided later in the application. FIG. 16B is a magnified version of FIG. 16A (the image shown was obtained using w_A1=0.3, w_B1=0.1, w_A2=0.2, w_B2=0.4). In this image, the coronary artery is more prominent. FIGS. 17A and 17B illustrate images rendered using a classifying 2D TF for a double block neighborhood (the image shown was obtained using w_A1=0.5, w_B1=0.2, w_A2=0.2, w_B2=0.4 with ranges: Φ₁=[1300,1500], Φ₂=[900,1150]). FIGS. 17A and 17B also highlight the coronary artery (but, as shown, introduce speckle artifacts at the large vessel boundaries).

[0112] The classification model can be configured to use domain knowledge of a physician to define feature or tissue specific criteria. This information can be defined at an OEM facility or at a use facility, then used for future analysis of similar regions of interest across different patients or within the same patient over time.

[0113] The classification model can accept input to establish one or more parameters that can specify a characteristic about the feature or tissue of interest. With respect to the coronary artery, a first criterion can be established as: “The coronary artery is more narrow than other vessels”. As such, one can define Φ₁ as the mutual range of all vessels. The narrowness can then be described as the range weight w_r(Φ₁). The system can be configured to allow a physician to define the range weights where the classification is confident. For instance, for a range weight of 0.3, a voxel can confidently be assumed to be a part of a large vessel, whereas a weight of 0.1 confidently implies a coronary artery voxel. These confidence levels can be denoted w_A1 and w_B1, respectively.

[0114] A measure p₁, defined as the classification certainty from the narrowness criterion can be determined. To do so, a signed representation can be employed, where the sign denotes the most likely tissue (in this example: negative for large vessels) and the magnitude is the confidence of the classification. This form can be achieved as described in Equations 2 and 3, where C₂ is a function that clamps to the [−2.0,2.0] interval. The resulting function of classification probability p₁ is shown in FIG. 8. The probability p₁ is a function of the range weight, w_r(Φ₁). The transformation can be defined by the confidence levels w_A1 and w_B1.
μ₁(w_A1+w_B1)/2, δ₁=(w_B1−w_A1)/2  (2)
p_i=C₂((w_r(Φ₁)−μ₁)/δ₁)  (3)
A second criterion of domain knowledge that may be used for this example can be associated with relative anatomical or spatial information about the feature of interest. For example, “the coronary artery is close to the heart muscle.” As such, Φ₂ can be defined as the range of heart muscle tissue. This proximity criterion can then be described as w_r(Φ₂). As for the narrowness, the physician can define two confidence levels w_A2 and w_B2 corresponding to confident classification as either material. The classification certainty p₂ is derived in the same way as p₁, see Equation 3. The total classification certainty P can be defined as a weighted sum of p₁ and p₂.

[0115] The general definition of P, valid for the separation of any two materials A and B, is given in Equation 4. There are n materials or criterion that can be used for diversifying criteria, defined by ranges Φ₁ through Φ_n. Each p_i can be derived in accordance with Equation 3. It is possible to steer the contribution from each criterion with the weights λ_i, but equal weighting is used for the examples shown and discussed herein. P is clamped to the [−1.0, 1.0] interval, where the extreme values correspond to more certain classification. Since the p_i components have a wider range, [−2.0,2.0], increased confidence from one component can dominate an uncertain classification in another. $\begin{array}{cc}P=C_1\left(\sum _{i=1}^n{\lambda }_i{\rho }_i\right),\mathrm{where}\sum _{i=1}^n{\lambda }_i=1& \left(4\right)\end{array}$
In summary, the methods can employ an a priori analysis of criterion to establish how neighborhood materials can achieve the wanted separation of A and B. The information for each material is the range Φ along with the reference levels w_A and w_B. Note that w_A and w_B can have “unnatural” values outside [0.0,1.0]. For example, if a range weight≧0.5 confidently implies tissue B, whereas 0.0 corresponds to equal probability of either tissue, then w_B=0.5 and w_A=−0.5.

[0116] For another data set of a different anatomical region, see FIGS. 12-14. FIGS. 12A and 12B were rendered using a 1D TF. In this rendering, spongy bone turns red. FIGS. 13A, 13B and 14A, 14B were obtained using a classifying 2D TF, in each image the vessels stand out more distinctly from the background. FIGS. 13A and 13B were rendered using voxel-centric local neighborhoods (w_A1=0.0, w_B1=0.1, w_A2=0.2, w_B2=−0.1). FIGS. 14A and 14B were rendered using double-block neighborhoods in which the material separation is achieved with minor artifacts (w_A1=0.0, w_B1=0.05, w_A2=0.2, w_B2=−0.2: ranges Φ₁=[1300,2000], Φ₂=[900,1150]).

[0117] In some embodiments, to separate the tissues in the rendering, the classification probability P from above can be used as a second attribute for each voxel in an overlapping range. This yields a multivariate data set, where the second value can be used to separate the overlapping tissues. In a general case there can be an arbitrary number of tissues, whose ranges may partly overlap. In particular embodiments, it may be assumed that there are no more than two overlapping tissues at any intensity value.

[0118] To render the multivariate data set, a 2D TF that is defined by tissue-centric trapezoids can be used. A trapezoid can be placed at either extreme of the classification dimension in a way that overlapping trapezoids are always on opposite sides as shown for example in FIG. 9. The second dimension can be filled by linear interpolation of the 1D TFs at the extremes. The interpolation is only necessary for overlapping intensity ranges, for other ranges the 1D trapezoid at one of the extremes can be used. FIG. 9 illustrates that a conventional 1D TF can be extended with a material classification dimension. Trapezoids for overlapping materials can be placed at extreme classification values and the 2D TF can be interpolated between them. The gray scale in FIG. 9 represents two intensity segments where the trapezoids overlap with a smooth transition from blue to red and blue to yellow, respectively. The red segment is marked by the upwardly extending rectangular post on the far left and top part of FIG. 9. The yellow segment is marked by the larger upwardly extending post on the right side of the figure. The blue transition zone is marked by the smaller upwardly extending post on the bottom part of the figure. The 1-D TF is extended with a material classification dimension thereby providing a 2-D TF with classifier attributes.

[0119] Advantageously, this approach has a number of desirable features. It is natural for the user to design a TF where one tissue at a time looks good. The adaptive trapezoid workflow is well suited to this task. The automatic interpolation of the second dimension then allows the power of a 2D TF without introducing complexity to the user. Another appealing benefit is that the classification dimension integrates a fuzzy segmentation into the rendering, reducing the impact of misclassification artifacts.

[0120] The implementation can be further simplified. In some embodiments, the used data sets have 12-bit precision, leaving 4 bits unused. The value of P, which does not need to be extremely precise, can be entered in these (top) 4 bits.

[0121] The local histograms are used to investigate neighborhoods as discussed above. The neighborhoods may be configured in a number of suitable manners as will be understood by those of skill in the art. For brevity of discussion, three different types of voxel surroundings are used in this specification: single block, voxel-centric, and double block neighborhoods. Single block neighborhoods are obtained by a static subdivision of non-overlapping blocks, typically cubes. A voxel-centric neighborhood can be more precise, since the region varies for each position in the volume. For further accuracy, the voxel itself can be defined to not be part of the neighborhood. In the voxel-centric examples of FIGS. 13A, 13B and 16A, 16B, a spherical neighborhood of radius 7 was employed, since this choice has proven to yield reliable statistics for many anatomical features.

[0122] Unfortunately, large voxel-centric neighborhoods can require extensive computation. Thus, an alternative complement, a double block neighborhood can be used that yields voxel-specific statistics while being fast to process. Two separate block subdivisions are applied to the volume, where one is shifted half a block size in all dimensions. Thus, each voxel is a member of exactly two blocks. Any neighborhood measure for a voxel is then simply derived as a linear combination of the values calculated for the two blocks it belongs to. The interpolation weights c₁ and c₂ are determined through two criteria: they are inversely proportional to distance to the center of each block (d₁ and d₂) and the sum is 1 for each voxel, see Equation 5. The block size of 8³ voxels was used for double block neighborhoods in the renderings shown in FIGS. 14A, 14B and 17A, 17B. $\begin{array}{cc}\{\begin{array}{c}{c}_1+c_2=1\\ c_1/c_2=d_2/d_1\end{array}\iff \{\begin{array}{c}{c}_1=d_2/\left(d_1+d_2\right)\\ c_2=d_1/\left(d_1+d_2\right)\end{array}& \left(5\right)\end{array}$

[0123] In some embodiments, an additional simplification can be implemented by avoiding the use of fully resolved block histograms. For double block neighborhoods, one can employ the histogram simplification approach described in Ljung et al., Transfer Function Based Adaptive Decompression for Volume Rendering of Large Medical Data Sets, Proceedings IEEE Volume Visualization and Graphics Symposium (2004), pp. 25-32. The simplification is based on piece-wise constant segments as shown for example in FIG. 10. This method may efficiently preserve the shape of the histogram, especially the peaks and/or peak characteristics. In this reference, the simplification consists of 12 segments with a minimum width of 10.

[0124] As illustrated in FIG. 6, embodiments of the invention may be configured as a data processing system, which can be used to carry out or direct operations of the rendering, and can include a processor circuit 200, a memory 236 and input/output circuits 246. The data processing system may be incorporated in, for example, one or more of a personal computer, workstation, server, router or the like. The processor 200 communicates with the memory 236 via an address/data bus 248 and communicates with the input/output circuits 246 via an address/data bus 249. The input/output circuits 246 can be used to transfer information between the memory (memory and/or storage media) 236 and another computer system or a network using, for example, an Internet protocol (IP) connection. These components may be conventional components such as those used in many conventional data processing systems, which may be configured to operate as described herein.

[0125] In particular, the processor 200 can be commercially available or custom microprocessor, microcontroller, digital signal processor or the like. The memory 236 may include any memory devices and/or storage media containing the software and data used to implement the functionality circuits or modules used in accordance with embodiments of the present invention. The memory 236 can include, but is not limited to, the following types of devices: ROM, PROM, EPROM, EEPROM, flash memory, SRAM, DRAM and magnetic disk. In some embodiments of the present invention, the memory 236 may be a content addressable memory (CAM).

[0126] As further illustrated in FIG. 6, the memory (and/or storage media) 236 may include several categories of software and data used in the data processing system: an operating system 252; application programs 254; input/output device drivers 258; and data 256. As will be appreciated by those of skill in the art, the operating system 252 may be any operating system suitable for use with a data processing system, such as IBM®, OS/2®, AIX®D or zOS® operating systems or Microsoft® Windows®95, Windows98, Windows2000 or WindowsXP operating systems Unix or Linux™. IBM, OS/2, AIX and zOS are trademarks of International Business Machines Corporation in the United States, other countries, or both while Linux is a trademark of Linus Torvalds in the United States, other countries, or both. Microsoft and Windows are trademarks of Microsoft Corporation in the United States, other countries, or both. The input/output device drivers 258 typically include software routines accessed through the operating system 252 by the application programs 254 to communicate with devices such as the input/output circuits 246 and certain memory 236 components. The application programs 254 are illustrative of the programs that implement the various features of the circuits and modules according to some embodiments of the present invention. Finally, the data 256 represents the static and dynamic data used by the application programs 254 the operating system 252 the input/output device drivers 258 and other software programs that may reside in the memory 236.

[0127] The data 256 may include a probability classification model 226 that provides an attribute for the transfer function module 224. As further illustrated in FIG. 6, according to some embodiments of the present invention application programs 254 include one or more of: a DVR Module 220, a (automatic or semi-automatic) Distinct Material Detection Module Using Local Histograms 224, and may include an Adaptive Trapezoid Transfer Function Module 225. The application programs 220, 224, 225 (and/or Probability Model 226) may be located in a local server (or processor) and/or database or a remote server (or processor) and/or database, or combinations of local and remote databases and/or servers.

[0128] While the present invention is illustrated with reference to the application programs 254, 220, 224, 225 in FIG. 6, as will be appreciated by those of skill in the art, other configurations fall within the scope of the present invention. For example, rather than being application programs 254 these circuits and modules may also be incorporated into the operating system 252 or other such logical division of the data processing system. Furthermore, while the application programs 220, 224, 225 are illustrated in a single data processing system, as will be appreciated by those of skill in the art, such functionality may be distributed across one or more data processing systems. Thus, the present invention should not be construed as limited to the configurations illustrated in FIG. 6, but may be provided by other arrangements and/or divisions of functions between data processing systems. For example, although FIG. 6 is illustrated as having various circuits and modules, one or more of these circuits or modules may be combined or separated without departing from the scope of the present invention.

[0129] The present invention is explained in greater detail in the following non-limiting Examples.

[0130] The automatic tissue detection system as well as the 2D TF classification system has been tested on a number of data sets from actual clinical use at CMIV (“Center for Medical Image Science and Visualization,” Linkoping University, Sweden). Results of the automatic tissue detection confirmed that the correct peaks were defined by manual segmentation for each distinct tissue. Some peak detection results are presented in FIGS. 11A and 11B. It is noted that very large and very small peaks were left out of these figures for clarity. The tests resulted in a total number of between about 5-10 peaks for each data set.

[0131] The first data set evaluation shown in FIG. 11A is an MR examination of the biliary duct, where the automatic peak detection reveals multiple PRHs where the distinct tissues are the liver, kidney, spleen and tumor. The spleen and the tumor are treated as one tissue below, since their histograms fully coincide. A value of ε=0.5 was used. The automatic peak detection performs well in finding the position of the tissue peaks. The height estimation varies in accuracy, depending on the extent of other soft tissue in the same intensity range. Deviations in height are a minor concern, since they have little effect on subsequent TF construction. The biliary duct is not found. Its histogram shows diffuse spikes below height 20, hardly characterizable as a peak.

[0132] The second data set shown in FIG. 11B is an MR angiography of the renal arteries (see also FIGS. 7A and 7B). The single interesting peak corresponds to the vessels with contrast agent. Despite the extreme difference in height compared to the main peak (4.2k vs 520k voxels), the tissue detection scheme accurately locates the vessel peak; the true value being 126, μ is found at 124. The height is slightly underestimated, since the narrow vessels do not dominate all their neighborhoods. A value of ε=0.95 was used. Thus, the detected peak corresponds with the true peak from segmentation.

[0133] The 2D TF approach described above has been tested on a CT pelvis angiography data set. The neighborhoods tested are of both voxel-centric spherical and double block type. The parameters of P were derived by manually probing the volume. Vessel voxels were characterized by neighborhoods with much soft tissue and little dense bone. The resulting renderings are shown in FIGS. 13A, 13B and 14A, 14B. The methods very well separate the two issues. The diagnostic value is enhanced since the vessels are not disturbed by red bone in the background. Even the difficult task of preserving thin vessels close to bone is successful. The fuzzy segmentation shows uncertain regions in pink, letting the user make the ultimate classification decision. The spherical neighborhood achieves a smooth classification. The double block approach yields faint speckle artifacts in the bone, but performs well in spite of its simplicity.

[0134] The separation can also be applied to identical tissues. An example for a CT angiography of the heart is given in FIGS. 15-17. The 1D TF (FIGS. 15A, 15B) makes all vessels and heart ventricles bright, since they are all filled with contrast agent. However, the primary objective of the examination is to study the coronary arteries. With P as a second attribute, it is possible to highlight only narrow vessels close to muscle tissue. Hence, the visually disturbing large bright regions are removed. A voxel-centric spherical surrounding provides high image quality (FIGS. 16A, 16B). The double block neighborhood (FIGS. 17A, 17B) causes artifacts for large vessel boundaries, being mistaken for narrow vessels, whereas the coronary artery shows few artifacts.

[0135] Another example on tissue separation is the MR biliary duct volume of FIG. 1A (see FIGS. 18A-18C). In this example, there is a large liver tumor that is seen as a slightly darker region when rendering the liver only as shown in FIG. 18A. A 1D TF attempting to highlight the tumor is of little help, since it becomes obscured by unimportant tissue as shown in FIG. 18B. When introducing the classification dimension to emphasize homogeneous regions with few low-intensity neighbors, the tumor as well as the spleen (top right) are clearly rendered as shown in FIG. 18C. The image shown in FIG. 18C used a voxel-centric neighborhood configuration (w_A1=0.1, w_B1=0.6, w_A2=0.2, w_B2=0.1, Φ₁=[150,200], Φ₂=0,140]).

[0136] As expected, the performance for the large spherical neighborhood is dependent on the number of voxels in the overlapping range. The classification technique may use about 0.02 ms per overlapping voxel, running a 1.8 GHz Intel Pentium M laptop. This amounts to about 32 s for the pelvis angiography data set and about 350 s for the heart. With the double block approach, these times are reduced to about 0.65 s and 1.85 s, respectively, of which about 60% is spent adapting the classification value to the 4-bit format.

[0137] Automatically adjusted parameters for volume rendering (VR) of MR angiography (MRA) data sets were used to analyze contrast-enhanced three-dimensional MRA from 18 patients with abdominal aortic aneurysm (AAA). An automatic algorithm based on partial range histograms was designed to detect the intensity range of the contrast agent in each data set. Three cases were used to derive a relation between the detected range and the optimal VR parameters. The individualized automatic parameters were then applied to each of the remaining 15 patients' data sets. The aortic diameter was measured proximal and distal to the AAA by three independent radiologists. The results were compared to free manual adjustment of VR parameters. Digital subtraction angiography (DSA) was used as the gold standard.

[0138] Results of the aorta measurements for VR with the automatic parameters were significantly more precise than with free manual adjustment. The average measurement error decreased from 1.5 mm to 0.7 mm, and the ratio of large errors (>2.5 mm) decreased from 13% to 0%. The interobserver agreement on the measurements was also significantly higher with the automatic parameters. Thus, automatically adjusted parameters for volume rendering of MR angiography data sets have been successfully been generated. Unexpectedly, automatically adjusted parameters were associated with smaller difference to DSA than time-consuming manual adjustment. Altogether, the proposed automatic algorithm makes VR suitable for MRA in clinical practice.

[0139] In summary, the use of local histogram analysis can be used to simplify and/or to extend TF design in VR. The simplification achieved by automatic tissue detection and adaptive trapezoids can facilitate a streamlined workflow of medical volume rendering. Preliminary results show that having detected all tissue peaks, it is feasible to automatically adapt TFs between data sets of similar examination types. Peaks that resemble each other, mainly in terms of ĥ and σ, can be identified as the same tissue. A TF connected to the tissues can then easily be repositioned.

[0140] For the more advanced user, voxel neighborhood properties are shown to extend the TF capabilities. We believe that the framework can be used for many complex rendering tasks.

[0141] The foregoing is illustrative of the present invention and is not to be construed as limiting thereof. Although a few exemplary embodiments of this invention have been described, those skilled in the art will readily appreciate that many modifications are possible in the exemplary embodiments without materially departing from the novel teachings and advantages of this invention. Accordingly, all such modifications are intended to be included within the scope of this invention as defined in the claims. The invention is defined by the following claims, with equivalents of the claims to be included therein.