Abstract:
A cross-time inspection method for medical image diagnosis. A first set of medical images of a subject is accessed wherein the first set is captured at a first time period. A second set of medical images of the subject is accessed, wherein the second set is captured at a second time period. The first and second sets are each comprised of a plurality of medical image. Image registration is performed by mapping the plurality of medical images of the first and second sets to predetermined spatial coordinates. A cross-time image mapping is performed of the first and second sets. Means are provided for interactive cross-time medical image analysis.

Description:
CROSS REFERENCE TO RELATED APPLICATIONS  
       [0001]     Reference is made to, and priority is claimed from, U.S. Provisional Patent Application No. 60/755,156, titled “CROSS-TIME INSPECTION FOR MEDICAL IMAGE DIAGNOSIS” in the names of Chen et al., provisionally filed on Dec. 29, 2005.  
         [0002]     Reference is made to U.S. Provisional Patent Application No. 60/754,884, titled “CROSS-TIME AND CROSS-MODALITY INSPECTION FOR MEDICAL IMAGE DIAGNOSIS” in the names of Chen et al., provisionally filed on Dec. 29, 2005. 
     
    
     FIELD OF THE INVENTION  
       [0003]     The present invention relates to a digital image processing/computer vision method for image analysis and, in particular, to cross-time inspection of tissues of different properties in medical image as a time function (cross-time image sequences).  
       BACKGROUND OF THE INVENTION  
       [0004]     Digital imaging techniques in medicine were implemented in the 1970&#39;s with the first clinical use and acceptance of the Computed Tomography or CT scanner. Later, extensive use of x-ray imaging (CT) and the advent of the digital computer and new imaging modalities like ultrasound and magnetic resonance imaging (MRI) have combined to create an explosion of diagnostic imaging techniques in the past three decades.  
         [0005]     There are benefit to using digital medical imaging technology in health care. For example, angiographic procedures for looking at the blood vessels in the brain, kidneys, arms and legs, and heart all have benefited from the adaptation of digital medical imaging and image processing technologies.  
         [0006]     With digital images, computerized multi-dimensional (e.g., spatial and temporal) image analysis becomes possible. Multi-dimensional image analysis can be used in applications such as automatic quantification of changes (anatomical or functional) in serial image volume scans of body parts, foreign objects localization, consistent diagnostic rendering, and the like.  
         [0007]     Also, different medical imaging modalities produce images providing different view of human body function and anatomy that have the potential of enhancing diagnostic accuracy dramatically with the help of the right medical image processing software and visualization tools. For example, X-ray computed tomography (CT) and magnetic resonance imaging (MRI) demonstrate brain anatomy but provide little functional information. Positron emission tomography (PET) and single photon emission computed tomography (SPECT) scans display aspects of brain function and allow metabolic measurements but poorly delineate anatomy. Furthermore, CT and MRI images describe complementary morphologic features. For example, bone and calcifications are best seen on CT images, while soft-tissue structures are better differentiated by MRI. Modalities such as MRI and CT usually provide a stack of images for certain body parts.  
         [0008]     It is known that the information gained from different dimensions (spatial and temporal) or modalities is often of a difference or complementary nature. Within the current clinical setting, this difference or complementary image information is a component of a large number of applications in clinical diagnostics settings, and also in the area of planning and evaluation of surgical and radiotherapeutical procedures.  
         [0009]     In order to effectively use the difference or complementary information, image features from different dimensions or different modalities had to be superimposed to each other by physicians using a visual alignment system. Unfortunately, such a coordination of multiple images with respect to each other is extremely difficult and even highly trained medical personnel, such as experienced radiologists, have difficulty in consistently and properly interpreting a series of medical images so that a treatment regime can be instituted which best fits the patient&#39;s current medical condition.  
         [0010]     Another problem encountered by medical personnel today is the large amount of data and numerous images that are obtained from current medical imaging devices. The number of images collected in a standard scan can be in excess of 100 and frequently numbers in the many hundreds. In order for medical personnel to properly review each image takes a great deal of time and, with the many images that current medical technology provides, a great amount of time is required to thoroughly examine all the data.  
         [0011]     Accordingly, there exists a need for an efficient approach that uses image processing/computer vision techniques to automatically detect/diagnose diseases.  
         [0012]     U.S. Publication No. 2004/0064037 (Smith), incorporated herein by reference, is directed to a technique that applies pre-programmed rules that specify the manner in which medical image data is to be classified or otherwise processed.  
         [0013]     U.S. Publication No. 2003/0095147 (Daw), incorporated herein by reference, relates to a computerized method of medical image processing and visualization.  
         [0014]     It is known that malignant breast tumors begin to grow their own blood supply network once they reach a certain size; this is the way the cancer can continue to grow. In a breast MRI scan, a contrast agent injected into the bloodstream can provide information about blood supply to the breast tissues; the agent “lights up” a tumor by highlighting its blood vessel network. Usually, several scans are taken: one before the contrast agent is injected and at least one after. The pre-contrast and post-contrast images are compared and areas of difference are highlighted. It should be recognized that if the patient moves even slightly between the two scans, the shape or size of the image may be distorted—a big loss of information.  
         [0015]     An contrast agent for MRI is Gadolinium or gadodiamide, and provides contrast between normal tissue and abnormal tissue in the brain and body.  
         [0016]     Gadolinium looks clear like water and is non-radioactive. After it is injected into a vein, Gadolinium accumulates in the abnormal tissue that may be affecting the body or head. Gadolinium causes these abnormal areas to become bright (enhanced) on the MRI. This makes it easy to see. Gadolinium is then cleared from the body by the kidneys. Gadolinium allows the MRI to define abnormal tissue with greater clarity. Tumors enhance after Gadolinium is given. The exact size of the tumor and location is important in treatment planning and follow up. Gadolinium is also helpful in finding small tumors by making them bright and easy to see.  
         [0017]     Dynamic contrast enhanced MRI is used for breast cancer imaging; in particular for those situation that have an inconclusive diagnosis based on x-ray mammography. The MRI study involves intravenous injection of a contrast agent (typically gadopentetate dimeglumine) immediately prior to acquiring a set of T1-weighted MR volumes with a temporal resolution of around a minute. The presence of contrast agent within an imaging voxel results in an increased signal that can be observed over the time course of the experiment.  
         [0018]     Study of these signal-time curves enables identification of different tissue types due to their differential contrast uptake properties as illustrated in  FIG. 1 . Typically, cancerous tissue shows a high and fast uptake due to a proliferation of “leaky” angiogenic microvessels, while normal and fatty tissues show little uptake. The uptake (dynamic) curves have often been fitted using a pharmacokinetic model to give a physiologically relevant parameterisation of the curve (refer to P. S. Tofts, B. Berkowitz, M. Schnall, “Quantitative analysis of dynamic Gd-DTPA enhancement in breast tumours using a permeability model”,  Magn Reson Med  33, pp 564-568, 1995).  
         [0019]     U.S. Pat. No. 6,353,803 (Degani, Hadassa), incorporated herein by reference, is directed to an apparatus and method for monitoring a system in which a fluid flows and which is characterized by a change in the system with time in space. A preselected place in the system is monitored to collect data at two or more time points correlated to a system event. The data is indicative of a system parameter that varies with time as a function of at least two variables related to system wash-in and wash-out behavior.  
         [0020]     Study of these curves/parameters has been used clinically to identify and characterize tumors into malignant or benign classes, although the success has been variable with generally good sensitivity but often very poor specificity (refer to S. C. Rankin “MRI of the breast”,  Br. J. Radiol  73, pp 806-818, 2000).  
         [0021]     While such systems may have achieved certain degrees of success in their particular applications, there is a need for an improved digital image processing method for medical image analysis that overcomes the problems set forth above and addresses the utilitarian needs set forth above.  
         [0022]     The present invention provides a method for image analysis and, in particular, for cross-time inspection of tissues of different properties in medical image as a time function.  
       SUMMARY OF THE INVENTION  
       [0023]     An object of the present invention is to provide a method for cross-time inspection of tissues of different properties (for example, abnormal and normal tissues) in medical image as a time function (cross-time image sequences).  
         [0024]     Any objects provided are given only by way of illustrative example, and such objects may be exemplary of one or more embodiments of the invention. Other desirable objectives and advantages inherently achieved by the disclosed invention may occur or become apparent to those skilled in the art. The invention is defined by the appended claims.  
         [0025]     The present invention provides a pattern recognition method for cross-time inspection of tissues of different properties using contrast enhanced MRI images augmented with other physical or non-physical factors. The method includes the steps of acquiring a plurality of medical image (e.g. MRI images before and after the injection of contrast enhancement agent) cross-time sequences; performing intra-registration of the plurality of medical image cross-time sequences with respect to spatial coordinates; performing inter-registration of the plurality of medical image cross-time sequences with respect to spatial coordinates; classifying tissues of different properties for the registered plurality of medical image cross-time sequences; and presenting the classification results for cross-time inspection.  
         [0026]     According to one aspect of the invention, there is provided a method for automatic abnormal tissue detection and differentiation using contrast enhanced MRI images augmented with other physical or non-physical factors. The method includes the steps of acquiring a plurality of MRI breast image sets; aligning the plurality of MRI breast images with respect to spatial coordinates; differencing the plurality of MRI breast image sets with a reference MRI image set, producing a plurality of difference image sets; segmenting the plurality of difference image sets, producing a plurality of MRI breast images with segmented intensity pixels; applying dynamic system identification to the segmented intensity pixels, producing a plurality of dynamic system parameters; and classifying the plurality of system parameters augmented with other physical or non-physical factors into different classes.  
         [0027]     According to another aspect of the invention, there is provided a method for automatic material classification. The method includes the steps of: acquiring a plurality of image sets of an object sequentially in time; aligning the plurality of image sets with respect to spatial coordinates; differencing the plurality of image sets with a reference image set to produce a plurality of difference image sets; segmenting the plurality of difference image sets to produce a plurality of images with segmented intensity pixels; applying dynamic system identification to the segmented intensity pixels of the plurality of images to produce a plurality of dynamic system parameters; and classifying the plurality of system parameters into different classes.  
         [0028]     According to still another aspect of the invention, there is provided a method for abnormal tissue detection using contrast enhanced MRI images. The method includes the steps of: acquiring a plurality of MRI breast image sets sequentially in time; aligning the plurality of MRI breast image sets with respect to spatial coordinates; differencing the plurality of MRI breast image sets with a reference MRI image set to produce a plurality of difference image sets; segmenting the plurality of difference image sets to produce a plurality of MRI breast image sets with segmented intensity pixels; applying a dynamic system identification to the segmented intensity pixels of the plurality of MRI breast image sets to produce a plurality of dynamic system parameters; and classifying the plurality of system parameters into different classes to detect abnormal tissue. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0029]     The foregoing and other objects, features, and advantages of the invention will be apparent from the following more particular description of the embodiments of the invention, as illustrated in the accompanying drawings. The elements of the drawings are not necessarily to scale relative to each other.  
         [0030]      FIG. 1  is a graph illustrating dynamic contrast uptake properties (curves) for different breast tissues.  
         [0031]      FIG. 2  is a schematic diagram of an image processing system useful in practicing the method in accordance with present invention.  
         [0032]      FIG. 3  is a flow chart illustrating one embodiment of the automatic abnormal tissue detection method in accordance with the present invention.  
         [0033]      FIG. 4  is a graph illustrating dynamic contrast uptake properties (curves) for malignant and benign tumor tissues.  
         [0034]      FIG. 5  is a schematic diagram illustrating the concept of step function response and system identification.  
         [0035]      FIG. 6  is a flowchart illustrating a method of system identification in accordance with the present invention.  
         [0036]      FIG. 7  is a graph illustrating two cross-time image sequences.  
         [0037]      FIG. 8  is a flowchart illustrating one embodiment of the cross-time tissue property inspection method in accordance with the present invention.  
         [0038]      FIG. 9  is a graph illustrating a method of cross-time tissue property inspection visualization presentations of the present invention.  
         [0039]      FIG. 10  is a flowchart illustrating a method of image registration in accordance with the present invention.  
         [0040]      FIG. 11  is a graph illustrating image registration concept. 
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0041]     The following is a detailed description of the preferred embodiments of the invention, reference being made to the drawings in which the same reference numerals identify the same elements of structure in each of the several figures.  
         [0042]      FIG. 2  shows an image processing system  10  useful in practicing the method in accordance with the present invention. System  10  includes a digital MRI image source  100 , for example, an MRI scanner, a digital image storage device (such as a compact disk drive), or the like. The digital image from digital MRI image source  100  is provided to an image processor  102 , for example, a programmable personal computer, or digital image processing work station such as a Sun Sparc workstation. Image processor  102  can be connected to a display  104  (such as a CRT display or other monitor), an operator interface such as a keyboard  106 , and a mouse  108  or other known input device. Image processor  102  is also connected to computer readable storage medium  107 . Image processor  102  transmits processed digital images to an output device  109 . Output device  109  can comprise a hard copy printer, a long-term image storage device, a connection to another processor, an image telecommunication device connected, for example, to the Internet, or the like.  
         [0043]     In the following description, a preferred embodiment of the present invention will be described as a method. However, in another preferred embodiment, the present invention comprises a computer program product for detecting abnormal tissues in a digital MRI image in accordance with the method described. In describing the present invention, it should be recognized that the computer program of the present invention can be utilized by any well-known computer system, such as the personal computer of the type shown in  FIG. 2 . However, other types of computer systems can be used to execute the computer program of the present invention. For example, the method of the present invention can be executed in the computer contained in a digital MRI machine or a PACS (picture archiving communication system). Consequently, the computer system will not be discussed in further detail herein.  
         [0044]     It will be further recognized that the computer program product of the present invention can make use of image manipulation algorithms and processes that are well known. Accordingly, the present description will be directed in particular to those algorithms and processes forming part of, or cooperating more directly with, the method of the present invention. Thus, it will be understood that the computer program product embodiment of the present invention may embody algorithms and processes not specifically shown or described herein that are useful for implementation. Such algorithms and processes are conventional and within the ordinary skill in such arts.  
         [0045]     Other aspects of such algorithms and systems, and hardware and/or software for producing and otherwise processing the images involved or co-operating with the computer program product of the present invention, are not specifically shown or described herein and can be selected from such algorithms, systems, hardware, components, and elements known in the art.  
         [0046]     A computer program for performing the method of the present invention can be stored in a computer readable storage medium. This medium may comprise, for example: magnetic storage media such as a magnetic disk (such as a hard drive or a floppy disk) or magnetic tape; optical storage media such as an optical disc, optical tape, or machine readable bar code; solid state electronic storage devices such as random access memory (RAM), or read only memory (ROM); or any other physical device or medium employed to store a computer program. The computer program for performing the method of the present invention may also be stored on computer readable storage medium that is connected to the image processor by way of the Internet or other communication medium. Those skilled in the art will readily recognize that the equivalent of such a computer program product may also be constructed in hardware.  
         [0047]     Turning now to  FIG. 8 , the method of cross-time inspection of tissues of different properties in medical image as a time function will be outlined.  FIG. 8  is a flow chart illustrating one embodiment of the method of the cross-time inspection of tissues of different properties in medical image of the present invention. In the embodiment shown in  FIG. 8 , a plurality of medical image cross-time sequences goes through a series of processes. Each of these processes performs a specific functionality such as intra-sequence registration, inter-sequence registration, dynamic curve classification, and visualization and diagnosis.  
         [0048]     Next, the concept of image registration is to be introduced. The method of curve classification will be discussed in depth later.  
         [0049]     Referring now to  FIG. 10 , the flow chart of the method of a generic image registration process is shown. The process of image registration is to determine a mapping between the coordinates in one space (a two dimensional image) and those in another (another two dimensional image), such that points in the two spaces that correspond to the same feature point of an object are mapped to each other. The process of determining a mapping between the coordinates of two images provides a horizontal displacement map and a vertical displacement map of corresponding points in the two images. The found vertical and horizontal displacement maps are then used to deform one of the involved images to minimize the misalignment between the two.  
         [0050]     In terms of image registration terminology the two images involved in registration process are referred as a source image  1020  and a reference image  1022 . Denote the source image and the reference image by I(x t ,y t ,t) and I(x t+1 ,y t+1 ,t+1) respectively. The notations x and y are the horizontal and vertical coordinates of the image coordinate system, and t is the image index (image  1 , image  2 , etc.). The origin, (x=0,y=0), of the image coordinate system is defined at the center of the image plane. It should be pointed that the image coordinates, x and y, are not necessarily integers.  
         [0051]     For the convenience of implementation, the image (or image pixel) is also indexed as I(i,j) where i and j are strictly integers and parameter t is ignored for simplicity. This representation aligns with indexing a matrix in the discrete domain. If the image (matrix) has a height of h and a width of w, the corresponding image plane coordinates, x and y, at location (i,j) can be computed as x=i−(w−1)/2.0, and y=(h−1)/2.0−j. The column index i runs from 0 to w−1. The row index j runs from 0 to h−1 .  
         [0052]     In general, the registration process is to find an optimal transformation function Φ t+1 (x t ,y t ) (see step  1002 ) such that 
 
[ x   t+1   ,y   t+1 ,1] T =Φ t+1 ( x   t   ,y   t )[ x   t   ,y   t ,1] T    (10-1) 
 
         [0053]     The transformation function of Equation (10-1) is a 3×3 matrix with elements shown in Equation (10-2).  
             Φ   =     [           ϕ   00           ϕ   01           ϕ   02               ϕ   10           ϕ   11           ϕ   12             0       0       1         ]             (     10   ⁢     -     ⁢   2     )             
 
         [0054]     In fact, the transformation matrix consists of two parts, a rotation sub-matrix  
         [           ϕ   00           ϕ   01               ϕ   10           ϕ   11           ]             
 
 and a translation vector  
         [           ϕ   02               ϕ   12           ]     .       
 
         [0055]     Noted that the transformation function Φ is either a global function or a local function. A global function Φ transforms every pixel in an image in a same way. A local function Φ transforms each pixel in an image differently based on the location of the pixel. For the task of image registration, the transformation function Φ could be a global function or a local function or a combination of the two.  
         [0056]     In practice, the transformation function Φ generates two displacement maps (step  1004 ), X(i,j), and Y(i,j), which contain the information that could bring pixels in the source image to new positions that align with the corresponding pixel positions in the reference image. In other words, the source image is to be spatially corrected in step  1008  and become a registered source image  1024 . For both displacement maps, X(i,j) and Y(i,j), the column index i runs from 0 to w−1 and the row index j runs from 0 to h−1.  
         [0057]     An exemplary result of misalignment correction is shown in  FIG. 11 . In  FIG. 11 , on the left is the source image  1102 ; on the right is the reference image  1106 . Clearly, there are varying vertical misalignments between the source image  1102  and the reference image  1106 . By applying the steps shown in  FIG. 10  to these two images, the vertical misalignment corrected source image is obtained (image  1104 ).  
         [0058]     Noted that the registration algorithm used in computing the image transformation function Φ could be a rigid registration algorithm, a non-rigid registration algorithm or a combination of the two. People skilled in the art understand that there are numerous registration algorithms that can carry out the task of finding the transformation function Φ that generates the needed displacement maps for the correction of the misalignment in two relevant images. Exemplary algorithms can be found in “Medical Visualization with ITK”, by Lydia Ng, et al. at http://www.itk.org. Also, people skilled in the art understand that spatially correcting an image with a displacement map could be realized by using any suitable image interpolation algorithms (see “Robot Vision” by Berthold Klaus Paul Horn, The MIT Press Cambridge, Mass.)  
         [0059]     For the present invention, the above discussed image registration process can be viewed as a black box  1000  with input terminal A ( 1032 ), input terminal B ( 1034 ) and output terminal D ( 1036 ). Box  1000  will be used in the following description of the present invention of cross-time inspection of tissues with different properties.  
         [0060]     Now turning back to  FIG. 8 , the processes of intra-sequence and inter-sequence registration will be described. Exemplary MRI image sequences for an object (a beast, for example) are depicted in  FIG. 7 . An MRI image sequence  704  contains an exemplary collection of MRI slice sets  706 , 708  and  710  for the same object (the breast). Each MRI slice set contains a number of slices that are images (cross-sections) of the object (the breast). Exemplary slices are slice (image)  712  for set  706 , slice (image)  714  for set  708 , and slice (image)  716  for set  710 . Purposely, MRI slice sets are taken at different time to capture functional changes of the object in time space when contrast enhancement agent is administrated. Exemplary time gap between the MRI slice sets could be 1 minute, 2 minutes, etc.  
         [0061]     For cross-time inspection of tissues with different properties, besides sequence  704 , one or more sequences of MRI image for the same object (the breast) are needed. An exemplary MRI sequence  724  is such a sequence. Sequence  724  is captured at a different time. Exemplary time gap between sequence  724  and sequence  704  could be several months.  
         [0062]     Similarly, sequence  724  contains an exemplary collection of MRI slice sets  726 , 728  and  730  for the same object (the breast). Each MRI slice set contains a number of slices that are images (cross-sections) of the object (the breast). Exemplary slices are slice (image)  732  for set  726 , slice (image)  734  for set  728 , and slice (image)  736  for set  730 . Purposely, MRI slice sets are taken at different time to capture functional changes of the object in time space. Exemplary time gap between the MRI slice sets could be 1 minute, 2 minutes, etc.  
         [0063]     An intra-sequence registration ( 804 ) is defined as registering slices (images) of the same cross-section of an object within a sequence of MRI image sets. Exemplary slices are slices (images)  712 ,  714 , and  716  for sequence  704 , and slices (images)  732 ,  734 , and  736  for sequence  724 . An embodiment of intra-sequence registration is discussed in the context of the method of tissue property inspection of a set of images, which acts as an independent entity, illustrated in  FIG. 3 . The need of intra-sequence registration stems from the fact that during the process of capturing MRI images, due the inevitable object (breast, for example) motion, images (for example,  712 ,  714  and  716 ) for the same cross-section of the object present misalignment. This misalignment may cause errors in the process of tissue property inspection.  
         [0064]     As stated previously, for cross-time inspection of tissues with different properties, two or more image sequences (such as sequences  704 , and  724 ) obtained at different times are required for the same object. Corresponding slices (such as  712  and  732 ) in different sequences are most likely misaligned and may have somehow different shapes. An inter-sequence registration ( 806 ) is thus needed and defined as registering slices (images) of the same cross-section of an object from different sequences. One embodiment of inter-sequence registration is pair-wise (2D) registration. Exemplary pairs of slices to be inter-registered are pairs  712  and  732 ,  714  and  734 , and  716  and  736 . Another embodiment of inter-sequence registration is volume-wise (3D) registration. In volume-wise (3D) registration, intra-registration is applied to individual sequences (e.g.  704  and  724 ) first. Then the intra-registered sequences are input to box  1000 .  
         [0065]     Turning now to  FIG. 3 , the method of tissue property inspection of a set of images (also step  808 , dynamic curve classification) will be outlined.  FIG. 3  is a flow chart illustrating one embodiment of the automatic abnormal tissue detection method of the present invention. Noted that the flow chart illustrated in  FIG. 3  serves as an independent entity that constitutes a self-contained process. Therefore, the flow chart illustrated in  FIG. 3  is not interpreted as an expansion of step  808 . Rather, step  808  and step  804  are explained using the steps shown in the flow chart in  FIG. 3 . In the embodiment shown in  FIG. 3 , a plurality of MRI breast images sets acquired before and after contrast agent injection go through a series of processes. Each of these processes performs a specific functionality such as alignment, subtraction, segmentation, system identification, and classification. In the present invention, abnormal tissue detection tasks are accomplished by a means of dynamic system parameter classification.  
         [0066]     In the embodiment shown in  FIG. 3 , a first step  202  (also step  802 ) is employed for acquiring a plurality of MRI breast image sets before and after an injection of contrast agent at one time. For cross-time inspection, step  202  repeats to acquire another plurality of MRI breast image sets before and after an injection of contrast agent at another time.  
         [0067]     Denote I 0 (x,y,z) as a set of MRI image for a breast with a number of images (slices) in a spatial order before an injection of contrast agent, where zε[1, . . . S] is the spatial order index, s is the number of images in the set, x and y are the horizontal and vertical indices respectively for an image where xε[1, . . . X] and yε[1, . . . Y]. After the administration of contrast agent, a plurality of MRI image sets is acquired with the same number (s) of images of the same breast for each set in the same spatial order z. The plurality of MRI image sets is taken with a temporal resolution, for example, of around one minute. This MRI image sets can be expressed by I k (x,y,z) where k is the temporal order index and kε[1, . . . K]; K is the number of sets. Exemplary sets are  706 ,  708  and  710 , (three set, K=3), or  726 ,  728  and  730 , (three sets, K=3). An exemplary slice I k (x,y,1) (at location  1 ) for set  706  (the first set for sequence  704 , k=1) is slice  712 .  
         [0068]     The presence of a contrast agent within an imaging voxel results in an increased signal that can be observed over the time course of the image acquisition process. Study of these signal-time curves enables identification of different tissue types due to their differential contrast uptake properties. For the purpose of automatic detection of abnormal tissues, the K sets of MRI images, I k (x,y,z), taken after the injection of contrast agent have to be spatially aligned (misalignment correction), in a step  204  (also step  804  intra-sequence registration), with a reference set of MRI images with respect to spatial coordinates x,y. In general, the reference set of MRI image is the set of MRI images, I 0 (x,y,z), taken before the injection of the contrast agent. The alignment process ensures that pixels belong to a same tissue region of the breast have the same x,y coordinates in all the K sets of images. The alignment process executes the following:  
                                                                                               for k = 1 : K                for   z = 1 : S                align(I k (x,y,z),I 0 (x,y,z))                end                end                      
 
         [0069]     Using the black box  1000 , I k (x,y,z) is input to terminal A ( 1032 ), I 0 (x,y,z) is input to terminal B ( 1034 ) and the registered image of I k (x,y,z) is obtained at output terminal D ( 1036 ). An exemplary method employable to realize the alignment function, align(A,B), is a non-rigid registration that aligns A with B and is widely used in medical imaging and remote sensing fields. The registration process (misalignment correction) has been discussed previously. Persons skilled in the art will recognize that other registration methods could also be used.  
         [0070]     As was shown in  FIG. 1 , after the injection of contrast agent, image pixel intensity increases differently for different breast tissues. This phenomenon indicates that subtracting the image taken before the injection from the image taken after the injection will provide radiologists with clearer information of locations of abnormal tissues in the image. This information can also be used to extract regions from the original MRI breast images for automatic abnormal tissue detection and differentiation. This information is obtained in step  206  in  FIG. 3  that carries out differencing the plurality of MRI breast image sets, I k (x,y,z), kε[1, . . . K] with a reference MRI image set to produce a plurality of  10  difference image sets, δI k (x,y,z), kε[1, . . . K]. The set of MRI images, I 0 (x,y,z), is selected as intensity reference images. The differencing process is executed as following:  
                                                                                                   for   k = 1 : K                for   z = 1 : S                δI k (x,y,z) = subtraction(I k (x,y,z),I 0 (x,y,z))                end                end                      
 
 wherein the function, subtraction(A,B), subtracts B from A. 
 
         [0071]     In  FIG. 3  at step  208 , the difference images, δI k (x,y,z), are subject to a segmentation process that first evaluates the plurality of difference image sets δI k (x,y,z), and produces a plurality of mask image sets, M k (x,y,z), kε[1, . . . K] obtained by executing:  
                                                                                                                                                   for   k = 1 : K                for   z = 1 : S                for   x = 1 : X                for   y = 1 : Y                if δI k (x,y,z) &gt; T            M k (x,y,z) = 1            end                 end                 end                 end                end                      
 
 wherein mask image sets, M k (x,y,z), kε[1, . . . K], are initialized with zeros, T is a statistical intensity threshold. An exemplary value of T is an empirical value 10. 
 
         [0072]     The segmentation process in step  208  segments the images in the plurality of MRI breast image sets, I k (x,y,z), according to the non-zero pixels in the mask images, M k (x,y,z), to obtain segmented intensity pixels in the images of the plurality of MRI breast image sets. Denote the resultant images by S k (x,y,z), kε[1, . . . K], the segmentation operation can be expressed as:  
                                                                                                                                                             for   k = 1 : K                for   z = 1 : S                for   x = 1 : X                for   y = 1 : Y                if M k (x,y,z) = 1                 S k (x,y,z) = I k (x,y,z)            end                 end                end                 end                end                      
 
         [0073]     wherein images, S k (x,y,z), are initialized as zeros. Persons skilled in the art will recognize that, in practical implementation, the stage of generating mask images can be omitted and the segmentation process can be realized by executing the following:  
                                                                                                                                                             for   k = 1 : K                for   z = 1 : S                for   x = 1 : X                for   y = 1 : Y                if δI k (x,y,z) &gt; T                 S k (x,y,z) = I k (x,y,z)            end                 end                end                 end                end                      
 
 wherein images, S k (x,y,z), are initialized as zeros. 
 
         [0074]     Referring now to  FIG. 4 , there is shown a chart that is a replica to the chart shown in  FIG. 1  except that  FIG. 4  includes the insertions of a step function, ∫(t), curve  302  and the removal of the normal and fat tissue curves.  
         [0075]     It is the intention of the present invention to detect abnormal tissues and more importantly to differentiate Malignant from Benign tissues. (Note: the step function, ∫(t), is defined as ∫(t&lt;0)=0; ∫(t≧0)=|λ|; λ≠0). Pixels that belong to normal and fat tissues are set to zeros in images S k (x,y,z) in the segmentation step  208 . The remaining pixels in images S k (x,y,z) belong to either malignant or benign tissues. It is practically difficult if not impossible to differentiate malignant tissue from benign tissue by just assessing the pixels brightness (intensity) in a static form, that is, in individual images. However, in a dynamic form, the brightness changes present a distinction between these two types of tissues. As shown in  FIG. 4 , starting from time zero, the brightness (contrast) curve  304 , m(t), of the malignant tissue rises quickly above the step function curve  302  and then asymptotically approaches the step function curve  302 ; while the brightness (contrast) curve  306 , b(t), of the benign tissue rises slowly underneath the step function curve  302  and then asymptotically approaches the step function curve, ∫(t),  302 .  
         [0076]     Persons skilled in the art can recognize that the brightness (contrast) curve  304 , m(t), resembles a step response of an underdamped dynamic system, while the brightness (contrast) curve  306 , b(t), resembles a step response of an overdamped dynamic system.  
         [0077]     An exemplary generic approach to identifying a dynamic system behavior is generally depicted in  FIG. 5 . For an unknown dynamic system  404 , a step function  402  is used as an excitation. A response  406  to the step function  402  from the dynamic system  404  is fed to a system identification step  408  in order to estimate dynamic parameters of system  404 .  
         [0078]     An exemplary realization of dynamic system modeling  212  (of  FIG. 3 ) is shown in  FIG. 6  where it is shown an ARX (autoregressive) model  500  (refer to “System identification Toolbox”, by Lennart Ljung, The Math Works).  
         [0079]     A general ARX model can be expressed as the equation: 
 
 y ( t )= G ( q )∫( t )+ H ( q )ε( t )   (1) 
 
 where G(q) ( 506 ) and H(q) ( 504 ) are the system transfer functions as shown in  FIG. 6 , u(t) ( 502 ) is the excitation, ε(t) ( 508 ) is the disturbance, and y(t) ( 510 ) is the system output. It is known that the transfer functions G(q) ( 506 ) and H(q) ( 504 ) can be specified in terms of rational functions of q −1  and specify the numerator and denominator coefficients in the forms:  
               G   ⁡     (   q   )       =       q     -   nk       ⁢       B   ⁡     (   q   )         A   ⁡     (   q   )                   (   2   )                 H   ⁡     (   q   )       =     1     A   ⁡     (   q   )                 (   3   )             
 
 wherein A and B are polynomials in the delay operator q −1 : 
 
 A ( q )=1 +a   1   q   1   + . . . +a   na   q   na    (4) 
 
 B ( q )= b   1   +b   2   q   −1   + . . . +a   nb   q   −nb+1    (5) 
 
 When A and B are polynomials, the ARX model of the system can be explicitly rewritten as: 
 
 y ( t )=−a 1   y ( t− 1)− . . . − a   na   y ( t−na )+ b   1   u ( t−nk )+ . . .  b   nb   u ( t−nk−nb+ 1)+ e ( t )   (6) 
 
 Equation (6) can be further rewritten as a regression as follows:  
                 y   ⁡     (   t   )       =         φ   ⁡     (   t   )       T     ⁢   θ       ⁢     
     ⁢       where   ⁢           ⁢     φ   ⁡     (   t   )         =         [           -     y   ⁡     (     t   -   1     )                 ⋮             -     y   ⁡     (     t   -   na     )                   u   ⁡     (     t   -   nk     )               ⋮             u   ⁡     (     t   -   nk   -   nb   +   1     )             ]     ⁢           ⁢   and   ⁢           ⁢   θ     =     [           a   1             ⋮             a   na               b   1             ⋮             b   nb           ]                 (   7   )             
 
 The system identification solution for the coefficient vector θ is  
               θ   ^     =         (       Φ   T     ⁢   Φ     )       -   1       ⁢     Φ   T     ⁢   Y             (   8   )             where                         Φ   =     [             φ   T     ⁡     (     t   0     )               ⋮               φ   T     ⁡     (       t   0     +     N   t     -   1     )             ]             (   9   )             and                         Y   =     [           y   ⁡     (     t   0     )               ⋮             y   ⁡     (       t   0     +     N   t     -   1     )             ]             (   10   )             
 
 In Equations (9) and (10), t 0  is the data sampling starting time and N t  is the number of samples. 
 
         [0080]     In relation to the brightness (contrast) curve m(t)  304 , and the brightness (contrast) curve b(t)  306 ,  
         φ   ⁡     (   t   )       =         [           -     m   ⁡     (     t   -   1     )                 ⋮             -     m   ⁡     (     t   -   na     )                   u   ⁡     (     t   -   nk     )               ⋮             u   ⁡     (     t   -   nk   -   nb   +   1     )             ]     ⁢           ⁢   and   ⁢           ⁢     φ   ⁡     (   t   )         =     [           -     b   ⁡     (     t   -   1     )                 ⋮             -     b   ⁡     (     t   -   na     )                   u   ⁡     (     t   -   nk     )               ⋮             u   ⁡     (     t   -   nk   -   nb   +   1     )             ]           
 
 respectively. 
 
         [0081]     In this particular case, u(t) is a step function. And the corresponding solutions are {circumflex over (θ)} m  and {circumflex over (θ)} b . The computation of {circumflex over (θ)} realizes the step of Dynamic system identification  210  (also step  408 ).  
         [0082]     Referring again to  FIG. 3 , in order to classify (step  214 ) a region with high contrast brightness in MRI images as benign or malignant tumor, a supervised learning step  218  is needed.  
         [0083]     A supervised learning is defined as a learning process in which the exemplar set consists of pairs of inputs and desired outputs. In this MRI image breast tissue classification case, the exemplar inputs are {circumflex over (θ)} m  and {circumflex over (θ)} b  (or the known curves), the exemplar desired outputs are indicators O m  and O b  for malignant and benign tumors respectively. In  FIG. 3 , step  218  receives M sample breast MRI dynamic curves with known characteristics (benign or malignant) from step  216 . An exemplary value for M could be 100. Within the M curves, there are M m  curves belong to malignant tumors and M b  curves belong to benign tumors. Exemplary values for M m  and M b  could be 50 and 50. In step  218 , applying Equation (8) to all the sample curves generates M coefficient vectors {circumflex over (θ)} among which, M m  coefficient vectors (denoted by {circumflex over (θ)} m   i , i=1 . . . M m ) represent malignant tumor with indicator O m , and M b  coefficient vectors (denoted by {circumflex over (θ)} b   i , i=1 . . . M b ) represent benign tumor with indicator O b . These learned coefficient vectors {circumflex over (θ)} m   i  and {circumflex over (θ)} b   i  are used to train a classifier that in turn is used to classify a dynamic contrast curve in a detection or diagnosis process.  
         [0084]     To increase the specificity (accuracy in differentiating benign tumors from malignant tumors) other factors (step  220 ) can be incorporated into the training (learning) and classification process. It is known that factors such as the speed of administration of the contrast agent, timing of contrast administration with imaging, acquisition time and slice thickness (refer to “Contrast-enhanced breast MRI: factors affecting sensitivity and specificity”, by C. W. Piccoli, Eur. Radiol. 7 (Suppl. 5), S281-S288 (1997)).  
         [0085]     Denote the speed of administration of the contrast agent by α, the timing of contrast administration with imaging by β, the acquisition time by γ and slice thickness by δ. These exemplary factors are to be used in conjunction with the coefficient vectors {circumflex over (θ)} m   i  and {circumflex over (θ)} b   i  to train the classify that that in turn is used to classify a region in the MRI breast image into malignant or benign tumor classes. Noted that these exemplary factors should be quantified in a range comparable to that of the coefficient vectors {circumflex over (θ)} m   i  and {circumflex over (θ)} b   i .  
         [0086]     For the learning and training purpose, construct the training data set 
 
{ p   j τ j   },j= 1  . . . l,τ   j ={−1,1 },p   j    ε R   d    (11) 
 
 wherein τ j  are the class labels. 
 
         [0087]     For example, if the tumor is malignant, τ j =1, otherwise, τ j =−1. The vector p j =[{circumflex over (θ)},α,β,γ,δ] is traditionally called feature vector in computer vision literature. The notion R d  represents a domain, d is the domain dimension. For this exemplary case, assume that the coefficient vector θ has five elements, then d=9. The data format in Equation (11) is used in leaning step  218  as well as in classification step  214 . Persons skilled in the art understand that the data vector p j  can be constructed in a different manner and augmented with different physical or non-physical numerical elements (factors) other than the ones aforementioned.  
         [0088]     There are known types of classifiers that can be used to accomplish the task of differentiating malignant tumors from benign tumors with the use of dynamic contrast curves along with other physical or non-physical factors. An exemplary classifier is an SVM (support vector machine) (refer to “A Tutorial on Support Vector Machines for Pattern Recognition”, by C. Burges,  Data Mining and Knowledge Discovery,  2(2), 1-47, 1998, Kluwer Academic Publisher, Boston, with information available at the website: http://aya.technion.ac.il/karniel/CMCC/SVM-tutonial.pdf).  
         [0089]     An example case of an SVM classifier would be training and classification of data representing two classes that are separable by a hyper-plane. A hyper-plane that separates the data satisfies 
 
 w·p+σ= 0   (12) 
 
 where · is a dot product. 
 
         [0090]     The goal of training the SVM is to determine the free parameter w and σ. A scaling can always be applied to the scale of w and σ such that all the data obey the paired inequalities: 
 
τ j ( w·p   j +σ)−1≧0, ∀j    (13) 
 
 Equation (13) can be solved by minimizing a Lagrangian function  
               L   ⁡     (     w   ,   ξ     )       =         1   2     ⁢          w        2       -       ∑     j   =   1     l     ⁢       ξ   j     ⁡     (       τ   j     ⁡     (       w   ·     p   j       +   σ     )       )                   (   14   )             
 
 with respect to the parameter w, and maximize it with respect to the undetermined multipliers ξ j ≧0. 
 
         [0091]     After the optimization problem has been solved, the expression for w in Equation (13) can be rewritten in terms of the support vectors with non-zero coefficients, and plugged into the equation for the classifying hyper-plane to give the SVM decision function:  
               Ψ   ⁡     (     p   new     )       =       (       w   ·     p   new       +   σ     )     =         ∑     j   =   1       l   s       ⁢       τ   j     ⁢     ξ   j     ⁢       p   j     ·     p   new           +   σ               (   15   )             
 
 wherein l s  is the number of support vectors. Classification of a new vector p new  into one of the two classes (malignant and benign) is based on the sign of the decision function. Persons skilled in the art will recognize that in non-separable case, non-linear SVMs can be used. 
 
         [0092]     The above described method of tissue property inspection of a set of images (also steps  804  and  808 ) is applied to all the cross-time image sequences such  704  and  724  for cross-time tissue property inspection. It is understood that in the present invention, the cross-time image sequence go through the steps of intra-registration and inter-registration before entering step  808 . One exemplary execution procedure of the steps of intra-registration and inter-registration for the exemplary sequences is applying intra-registration to sequence  704  first, then applying inter-registration to sequences  704  and  724 . People skilled in the art should know that the roles of sequences  704  and  724  are exchangeable.  
         [0093]     For intra-registering sequence  704  for this particular exemplary execution procedure, select arbitrarily a set of images as the reference image set, e.g. set  706 . Images of set  706  are input to terminal B ( 1034 ), other image sets ( 708  and  710 ) are input to terminal A ( 1032 ). The registered images of image sets ( 708  and  710 ) are obtained at terminal D ( 1036 ).  
         [0094]     For inter-registration for this particular exemplary execution procedure, images of sequence  724  are input to terminal A ( 1032 ), images of sequence  704  are input to terminal B ( 1034 ) and the registered images of sequence  724  are obtained at output terminal D ( 1036 ).  
         [0095]     Upon the completion of step  808 , multiple dynamic curves (two curves in the current exemplary case) are generated reflecting tissue properties captured in multiple cross-time image sequences (two sequences  704  and  724  for the current exemplary case) at multiple time instances (two for the current exemplary case). It is well known that these dynamic curves provide the medical professionals with valuable information regarding disease conditions (or progressions) for patients. In step  810 , visualization tools are employed for medical professional to examine concerned regions of the object (regions of interest in the images) for better diagnosis. One embodiment of such visualization facility is illustrated in  FIG. 9 .  
         [0096]     There is shown in  FIG. 9 a  computer monitor screen  900  (also  104  in  FIG. 2 ) hooked up to an image processor ( 102 ) that executes previously described steps. On screen  900 , two representative image slices  712  and  732  are shown on the left. For example, slice  712  is the first image of I k (x,y,1)|kε[1,2,3] across three sets ( 706 ,  708  and  710 ) at spatial location  1 ; slice  732  is the first image of I k (x,y,1)|kε[1,2,3] across three sets ( 726 ,  728  and  730 ) at spatial location  1 . Breast images  902  and  912  are shown in slices  712  and  732 . Breast images  902  and  912  are the images of a same cross-section of a breast. In operation, a medical professional moves a computer mouse  906  (as a user interface) over a location  908  in slice  712 . Simultaneously, a ghost mouse  916  appears at the same spatial location  918  in slice  732  as  908  in slice  712 . The user also can move a computer mouse  916  (as a user interface) over a location  918  in slice  732 . Simultaneously, a ghost mouse  906  appears at the same spatial location  908  in slice  712  as  918  in slice  732 . In either case, two dynamic curses ( 924  solid and  926  dashed) appear at the left side ( 922 ) of the screen. The exemplary curves  924  and  926  reflect different tissue properties for the same spot of a breast at two different times. For example, the image sequence containing slice  712  may be taken  6  months prior capturing the sequence containing slice  732 . The medical professional can move the mouse to other locations to examine the change of the tissue properties over time (6 months). With this visualization facility, disease progression can be readily analyzed. People skilled in the art should understand that tissue properties could be represented by other means besides the dynamic curve plots  924  and  926 . For example, tissue properties could be represented by colored angiogenesis maps. People skilled in the art also understand that multiple cross-time image sequences can be processed by the method of the current invention and multiple dynamic curves can be displayed simultaneously for medical diagnosis.  
         [0097]     The subject matter of the present invention relates to digital image processing and computer vision technologies, which is understood to mean technologies that digitally process a digital image to recognize and thereby assign useful meaning to human understandable objects, attributes or conditions, and then to utilize the results obtained in the further processing of the digital image.  
         [0098]     The invention has been described in detail with particular reference to a presently preferred embodiment, but it will be understood that variations and modifications can be effected within the spirit and scope of the invention. The presently disclosed embodiments are therefore considered in all respects to be illustrative and not restrictive. The scope of the invention is indicated by the appended claims, and all changes that come within the meaning and range of equivalents thereof are intended to be embraced therein.