WO1999054705A2 - Iterative warping for temporal subtraction of radiographs - Google Patents

Iterative warping for temporal subtraction of radiographs Download PDF

Info

Publication number
WO1999054705A2
WO1999054705A2 PCT/US1999/004290 US9904290W WO9954705A2 WO 1999054705 A2 WO1999054705 A2 WO 1999054705A2 US 9904290 W US9904290 W US 9904290W WO 9954705 A2 WO9954705 A2 WO 9954705A2
Authority
WO
WIPO (PCT)
Prior art keywords
shift
pixels
histogram
image
selecting
Prior art date
Application number
PCT/US1999/004290
Other languages
French (fr)
Other versions
WO1999054705A3 (en
Inventor
Takayuki Ishida
Shigehiko Katsuragawa
Kunio Doi
Original Assignee
Arch Development Corporation
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Arch Development Corporation filed Critical Arch Development Corporation
Priority to EP99917300A priority Critical patent/EP1295244A4/en
Priority to AU35452/99A priority patent/AU3545299A/en
Priority to CA002326776A priority patent/CA2326776C/en
Priority to JP2000545001A priority patent/JP2003530722A/en
Publication of WO1999054705A2 publication Critical patent/WO1999054705A2/en
Publication of WO1999054705A3 publication Critical patent/WO1999054705A3/en

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G06T7/0014—Biomedical image inspection using an image reference approach
    • G06T7/0016—Biomedical image inspection using an image reference approach involving temporal comparison
    • G06T3/14—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/254—Analysis of motion involving subtraction of images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/337—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30061—Lung

Definitions

  • the present is related to automated techniques for automated detection of abnormalities in digital images, for example as disclosed in one or more of U.S. Patents 4,839,807; 4,841,555; 4,851,984; 4,875,165; 4,907,156; 4,918,534; 5,072,384; 5,133,020; 5,150,292; 5,224,177; 5,289,374; 5.319,549; 5,343,390; 5,359,513; 5,452,367; 5,463.548; 5,491,627; 5,537,485; 5.598,481 ; 5,622,171; 5,638,458; 5,657,362; 5,666,434; 5.673,332; 5,668,888; as well as U.S.
  • the present invention also relates to technologies referenced and described in the references identified in the appended APPENDIX and cross-referenced throughout the specification by reference to the number, in brackets, of the respective reference listed in the APPENDIX, the entire contents of which are also incorporated herein by reference.
  • Various of these publications may correspond to various of the cross-referenced patents and patent applications.
  • the present invention is related to temporal analysis of medical images and, in particular, to the analysis of chest radiograph images using automated temporal subtraction. Discussion of the Background
  • radiologists For the interpretation of chest radiographs, radiologists commonly compare a current film with previous films in order to facilitate the detection of abnormalities on chest radiographs, such as pulmonary nodules, interstitial infiltrates, pleural effusions, and cardiomegaly.
  • abnormalities on chest radiographs such as pulmonary nodules, interstitial infiltrates, pleural effusions, and cardiomegaly.
  • investigators at the University of Chicago Department of Radiology have developed a temporal subtraction scheme based on a nonlinear geometric warping technique.
  • Another object of this invention is to provide a novel image processing technique to overcome severe misregistration errors mainly due to differences in a subject's inclination and/or rotation in temporally sequential images.
  • a further object of this invention is to provide a novel temporal image subtraction technique to assist radiologists in the detection of interval changes on chest radiographs.
  • the method includes the steps of (a) determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image; (b) warping said second digital image based on the first shift values to obtain a warped image in which spatial locations of pixels are varied in relation to said first shift values; (c) determining second shift values between pixels of said first digital image and pixels of said warped image obtained in step (b); (d) wa ⁇ ing said warped image obtained in step (b) based on the second shift values to obtain an iteratively wa ⁇ ed image in which spatial locations of pixels of said wa ⁇ ed image obtained in step (b) are varied in relation to said second shift values; and (e) subtracting the iteratively wa ⁇ ed image from said first digital image.
  • initial shift values are obtained by cross- correlation techniques using a template and search regions of interest. Shift vectors and a histogram of shift vectors are obtained from each initial shift value. Based on the histogram of shift vectors, variations in various groups of the initial shift values are smoothed using two-dimensional fitting.
  • the present invention similarly includes a computer readable medium storing program instructions by which the method of the invention can be performed when the stored program instructions are appropriately loaded into a computer, and a system for implementing the method of the invention.
  • Figure 1 is an illustration of selecting different ROIs for analysis of histograms of subtraction images.
  • Figures 2(a)(A) and 2(a)(B) are photographs of subtraction images with ROI placement for analysis of histograms for pixel values, with the photograph of Fig. 2(a)(A) being an example of a poor subtraction image and the photograph of Fig. 2(a)(B) being an example of a good subtraction image.
  • Figure 2(b) is a graph of histograms of subtraction images of Figures 2(a)(A) and 2(a)(B).
  • Figure 3(a)-3(e) are illustrations of distributions of histogram widths in both lungs for cases rated as (a) 1 (very poor quality), (b) 2 (poor quality), (c) 3 (adequate quality), (d) 4 (good quality), and (e) 5 (excellent quality), respectively.
  • Figure 4(a) is a schematic diagram of the temporal subtraction method using an iterative wa ⁇ ing technique according to the present invention.
  • Figures 4(b)-4(d) are highly schematic diagrams of shift vector analysis, first wa ⁇ ing, and second wa ⁇ ing steps, respectively, included in the method of Figure 4(a).
  • Figures 4(e) 1 and 4(e)(2) illustrate in more detail than the highly schematic Figures 4(a)-4(d) the processing performed in the iterative process of the present invention.
  • Figures 5(a) and 5(b) are photographs of respective subtraction images obtained using different shift vectors, i.e., (a) good subtraction image and (b) a poor subtraction image, wherein shift vectors are overlaid on each image, dark and light vectors (short lines) are used for inside and outside, respectively, of approximate lung areas, and small white dots at one end of dark vectors indicate initial locations of template ROIs, which are used in the wa ⁇ ing technique
  • Figures 6(a) and 6(b) are graphs respectively illustrating a distribution of accumulated shift vector values as a function of shift vector orientation for (a) a good subtraction case and (b) a poor subtraction case.
  • Figures 7(a)-7(d) are photographs respectively illustrating (a) dominant shift vectors (dark lines) in a temporal image, (b) subtraction image; (c) remaining shift vectors (dark lines) and (d) a subtraction image.
  • Figures 8(a), 8(b) and 8(c) are photographs of subtraction images respectively obtained with (a) a single wa ⁇ ing technique, (b) an iterative wa ⁇ ing technique, and (c) an iterative wa ⁇ ing technique plus linear inte ⁇ olation of shift values.
  • Figures 9(a), 9(b), 9(c), (d), and 9(e) are graphs respectively illustrating distributions of histogram widths in both lungs for the 181 cases, obtained with the prior wa ⁇ ing technique (dots) and the wa ⁇ ing technique of the present invention (x), rated as 1 (very poor quality), 2 (poor quality), 3 (adequate quality), 4 (good quality), and 5 (excellent quality).
  • Figures 10(a) and 10(b) are photographs of subtraction images for comparison, respectively obtained (a) with the prior wa ⁇ ing technique and (b) with the wa ⁇ ing technique of the present invention.
  • Figure 11 is a bar graph illustrating the distribution of the degree of improvement, with the subjective rating method, in the quality of the subtraction images obtained with the wa ⁇ ing technique of the present invention.
  • Figure 12 is a schematic illustration of a general pu ⁇ ose computer 100 programmed according to the teachings of the present invention.
  • the image database used in the derivation of the present invention included 181 pairs of chest radiographs obtained from the Iwate Prefectural Hospital, Morioka, Japan. These 181 cases were obtained sequentially from chest screening images which were made with a Fuji Computed Radiography (FCR) system (Fuji Medical Systems Co., Ltd., Tokyo, Japan). The time interval between the current and previous images for all cases was 1 year. The pixel size and gray level of CR chest images were 0.2 mm and 1024, respectively.
  • the technique of the present invention was developed by use of a Silicon Graphics 02 workstation. Methodology 1. Subjective evaluation of the quality of subtraction images
  • the final rating for each case was determined in one of the five categories above based on the average (or the majority) of all observers' ratings.
  • Figs. 1(a)- 1(d) For determination of the histogram of pixel values in each lung field of the subtraction image, various ROIs of different size and shape in each lung were examined, as shown in Figs. 1(a)- 1(d). Initially, the ROI included all of the lung and mediastinum areas, as illustrated in Fig. 1(a). Second, we segmented the fight and left lungs were segmented by using the heart border information, as shown in Fig. 1(b), in order to evaluate the quality of the subtraction image in each lung field separately. However, it was found that the top and the bottom areas of the ROIs tended to include misregistration errors caused by the difference between the clavicle positions, and also the difference between different diaphragm levels, respectively.
  • the areas around the ribcage edge boundaries and the boundaries of the mediastinum were excluded from the ROIs, as shown in Fig. 1(c).
  • small ROIs as shown in Fig. 1(d) were employed, because the extent of the misregistration error for the ribs can be detected more sensitively with the ROI placed over the central area of the lung rather than with the other ROIs shown in Figs. 1(a), 1(b), and 1(c).
  • the size of the small ROIs for both lungs was 30 pixels (width) x 120 pixels (length).
  • the centers of the ROIs were determined by using the x,y coordinates of the ribcage edges and the cardiac edges.
  • the vertical locations of the two centers of the ROIs were the same and equal to the y-location between the top lung, which was obtained from the ribcage edge detection scheme, and the bottom of the right ribcage edge.
  • the horizontal centers of the ROIs were determined as the x-location in each lung between the ribcage edge and the cardiac edge, at the vertical center level of the ROIs.
  • Figure 2 (a) shows a comparison of a poor (A) and a good (B) subtraction images of the right lung.
  • the ROIs used for determination of the histogram width are illustrated by black lines.
  • the corresponding histograms of these images are shown in Figure 2 (b). It is seen that the width of the histogram for the poor subtraction image (A) is much wider than that for the good subtraction image (B).
  • the "misregistration" the difference between an abnormal current image and a normal previous image
  • due to the actual interval abnormal change would usually be localized and small compared with the misregistration errors due to a failure in image matching. Therefore, the quality of the subtraction image was evaluated by using the widths of the histograms in the ROI for each lung.
  • the width of the histogram at 10%) of the largest peak was obtained. It was found that the width of the histogram at a low level such as 10%> of the peak was more sensitive than the histogram width at a high level in distinguishing between good and poor subtraction images.
  • the distributions of the histogram widths for each of the five groups of different qualities of subtraction images, which were grouped based on subjective ratings, obtained by the previous technique are shown in Fig. 3 (a)-3(e). In these figures, the horizontal and vertical axes correspond to the width of the histogram in the right and the left lung, respectively.
  • the histogram widths for "very poor” and “poor” subtraction images tend to be large and are distributed in the upper right, as shown in Figs. 3 (a) and 3(b).
  • the histogram widths of "adequate” subtraction images are distributed in the intermediate range, as shown in Fig. 3 (c).
  • the histogram widths for "good” subtraction images are shifted toward the lower left, as shown in Fig. 3 (d).
  • the "excellent" subtraction images as represented by a subjective rating score of 5 are characterized by narrow histogram widths, as shown in Fig. 3 (e). 3.
  • a nonlinear density correction technique [1,7] was applied first for adjustment of the density and contrast in the two digitized images, i.e., the current and previous chest radiographs.
  • the difference in the lateral inclination between the two images was corrected by use of an image rotation technique.
  • the global shift values were then determined for the initial registration by use of the cross-correlation of a pair of blurred low-resolution images obtained from the current and the previous image.
  • a number of template ROIs 32 x 32 matrix
  • the corresponding search area ROIs 64 x 64 matrix
  • Shift values, ⁇ x and ⁇ y, for all pairs of selected ROIs were determined by using a cross-correlation technique to find the "best" matched areas in the current and previous images.
  • a two-dimensional surface fitting by use of polynomial functions was then applied to each set of mapped shift values, ⁇ x and ⁇ y, for conversion of the x, y coordinates of the previous image, i.e., for wa ⁇ ing of the image.
  • the wa ⁇ ed previous image was then subtracted from the current image.
  • the temporal subtraction technique is illustrated in overview in Fig. 4 (a).
  • the matrix size for the current and the previous chest image obtained with a CR system is reduced to 586 x 586 (one third of the original image matrix size).
  • a density correction technique was not applied because the chest images obtained with the CR system can maintain consistent density and contrast with the use of the exposure data recognition (EDR) system [4] included in most recent CR systems.
  • EDR exposure data recognition
  • global matching 420 is performed using an image rotation [2], determining the angle between the two images by comparison of the two midlines of these images [5], and applying an initial image matching technique for determination of the global shift value between the two images, which corresponds to the shift in the x, y coordinates of one image relative to the other.
  • the two images are globally registered by use of the initial image matching technique based on the cross-correlation of blurred low-resolution images.
  • an iterative image wa ⁇ ing technique is employed to improve local image matching and thus to reduce registration errors.
  • local segmentation 430 is performed in which the ribcage edges of the previous image are detected based on image profile analysis.
  • the cardiac edges are determined based on an edge detection technique for segmentation of the lungs of the current image.
  • the shift values ⁇ x and ⁇ y, for all of the selected ROIs are determined by a cross-correlation technique.
  • Shift vectors based on the shift values ⁇ x and ⁇ y are then determined, as discussed in more detail hereinafter.
  • a first image wa ⁇ ing (step 450), typically of the previous or older image, is performed to derive a first wa ⁇ ed image (step 460). Then a second wa ⁇ ing (step 470) is performed on the first wa ⁇ ed image, followed by obtaining of a subtraction image 480 for diagnosis.
  • Fig. 4(b) illustrates very schematically one embodiment of determination of fitted shift values for the first wa ⁇ ing step.
  • a shift vector orientation histogram is obtained (step 411) for each lung.
  • step 412 peaks in the histogram of each lung are determined, and utilized to identify dominant and remaining shift vectors (step 413).
  • a first wa ⁇ ing is performed.
  • step 451 optimized shift values are determined and then in step 452 subjected to surface fitting [1]. After fitting, a coordinate transformation is performed in step 453 and a wa ⁇ ed image obtained in step 454.
  • a second image wa ⁇ ing technique which is applied to the wa ⁇ ed previous image, is then performed.
  • the optimized shift values between the current image and the wa ⁇ ed previous image are determined by the cross-correlation technique and vectors are again determined in step 471 to derive optimized shift values, as discussed in more detail hereinafter.
  • a linear inte ⁇ olation is performed on the optimized shift values in step 472 to determine the final shift values on all of the x,y coordinates of the wa ⁇ ed previous image.
  • Another coordinate transformation is performed in step 473 to obtain a second wa ⁇ ed image in step 474.
  • the temporal subtraction image is obtained by subtraction of the second wa ⁇ ed previous image from the current chest image, as shown schematically in step 480 of Fig. 4(a). 4. Weighting factors for surface fitting of shift values
  • weighting factors determined based on the cross-correlation values were used for surface fitting to determine the local shift values. [1] In general, the larger the cross-correlation value, the larger the weighting factor. However, it was found that the cross-correlation values near the ribcage borders, in the mediastinum area, and below the diaphragm area were generally very large, and much greater than the values in the lung areas.
  • the fitted shift values would have been affected considerably by the shift values near the ribcage borders, in the mediastinum, and below the diaphragms area. This is not desirable, because accurate subtraction is commonly required in the lung fields rather than in other areas in chest images.
  • a lung segmentation method which includes the detection of ribcage edges and heart edges for identifying the shift values in the lung area.
  • a weighting factor of 1.0 is assigned for ROIs in the lung areas, and 0.25 for ROIs in the mediastinum and below the diaphragm, whereupon weighted polynomial fitting is performed. [1]
  • the fitted shift values are determined for the first wa ⁇ ing step. When the new weighting factors were used, the registration errors in the lung areas were reduced in some cases. 5. Shift- vector orientation analysis
  • the surface fitting technique of the shift values was generally effective in improving local registration around the poorly matched regions. [1] However, if some shift vectors with incorrect orientation are included within a relatively small region, the correct shift vectors may be affected significantly by the fitting. In the previous method [1.3], the orientation (or angle) of the shift vectors was not considered in the process for determining the surface fitting of the shift values. However, according to the present invention it has since been determined that the orientation of the shift vectors is an important factor for accurate surface fitting, and therefore two different components, i.e., the dominant and the remaining (non-dominant) shift vectors, are identified by using the shift vector orientation histogram, which is a histogram of accumulated shift values as a function of the shift vector orientation.
  • the ribcage edges and the cardiac edges are determined by using the edge detection technique [6] for segmentation of the lungs.
  • Many ROIs are selected in the lung field for both the previous and current images.
  • Shift values, ⁇ x and ⁇ y, are determined for each pair of ROIs based on the cross-correlation technique.
  • the orientation (or angle) of the shift vector for each ROI is determined by the arc tangent of ⁇ y/ ⁇ x. Then the shift- vector orientation histograms for the right and the left lung are obtained.
  • Figures 5 (a) and 5(b) illustrate the distributions of shift vectors obtained with a good and a poor subtraction image, respectively.
  • the orientations of the shift vectors for good subtraction images tend to be similar. This may be because the rib contrast in this chest image is relatively high, and thus the template ROI which includes a rib edge could easily match the corresponding rib edge in the search area ROI.
  • the orientations of the shift vectors for a poorly registered case were varied because of the low-contrast ribs.
  • the shift-vector orientation histograms of these cases, Figs. 5(a) and 5(b), for the right and left lungs are shown in Fig. 6(a) and 6(b). It is apparent that the shift-vector orientation histograms for a good registered case have a large peak in each lung. However, there is no obvious peak in the orientation histograms for a poor case, as shown in Fig. 6(b).
  • FIGS 7(a) and 7(b) show the selected dominant shift vectors and the corresponding subtraction image obtained by fitting with the dominant shift vectors, respectively.
  • Figures 7(c) and 7(d) illustrate non-dominant shift vectors and the corresponding subtraction image obtained by fitting with the non-dominant shift vectors, respectively. It is clear that the subtraction image obtained with the dominant shift vectors is superior to that obtained with the non-dominant shift vectors.
  • Figures 8(a) and 8(b) show the temporal subtraction images obtained with the prior technique and that of the present invention using shift- vector analysis, respectively. The quality of the subtraction image obtained according to the present invention is improved slightly.
  • the present invention selects for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second wa ⁇ ed images exhibiting the narrower histogram of pixel values, and then performs further two-dimensional fitting on the further selected fitted shift values, and then wa ⁇ s the previous image using said further fitted shift values to produce a wa ⁇ ed image for further iterative processing.
  • the image wa ⁇ ing technique is applied repeatedly and thus iteratively first on the previous image, next on the first wa ⁇ ed image, and then on the second wa ⁇ ed image, and so on, until the desired quality of the subtraction image is obtained.
  • the second and subsequent image wa ⁇ ing steps are employed for improvement of the quality of the temporal subtraction image.
  • the wa ⁇ ed previous image from the first wa ⁇ ing step is obtained.
  • the current image and the wa ⁇ ed previous image are then used for the second image wa ⁇ ing.
  • the shift values, ⁇ x and ⁇ y are determined by use of the cross-correlation technique for all of the selected ROIs in the current image and the wa ⁇ ed previous image.
  • shift vectors and a histogram of shift vectors are derived, angular ranges determined, and two sets of shift vectors again selected for production of interim first and second twice wa ⁇ ed images.
  • corresponding subtraction images are obtained, and the histograms of pixel values in each subtraction image are evaluated.
  • one subtraction image exhibits a narrower histogram of pixel values in both lung regions in comparison to that of the other subtraction image, then that subtraction image is considered the final subtraction image for diagnosis (assuming that no further wa ⁇ ing iterations are to be performed). If however that is not the case, then another shift value fitting is performed.
  • the present invention selects for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second wa ⁇ ed images exhibiting the narrower histogram of pixel values, and then performs further two-dimensional fitting on the further selected fitted shift values, and then wa ⁇ s the previously wa ⁇ ed image using the further fitted shift values to produce a secondly wa ⁇ ed image for further iterative processing, if desired, or for production of the final subtraction image for diagnosis.
  • the shift values of all x,y coordinates over the entire lung fields for the second image wa ⁇ ing are preferably obtained by a linear inte ⁇ olation technique instead of the surface fitting technique with polynomial functions used during the first wa ⁇ ing, since linear inte ⁇ olation was observed to result in a somewhat better subtraction image after the second image wa ⁇ ing.
  • the present invention can also be practiced using linear inte ⁇ olation for the first image wa ⁇ ing or polynomial fitting for the second image wa ⁇ ing.
  • Figure 8(c) shows a subtraction image obtained by application of the second image wa ⁇ ing technique.
  • the registration errors for the ribs in the temporal subtraction image are reduced substantially by use of the second image wa ⁇ ing technique.
  • a third image wa ⁇ ing technique which is basically the same as the second wa ⁇ ing process, can also be used.
  • the quality of the subtraction image was improved slightly in some cases by applying the third image wa ⁇ ing; however, at the expense of increased computational time for the iterative wa ⁇ ing process.
  • two iterations of iterative image wa ⁇ ing are preferably employed, i.e., the second image wa ⁇ ing technique is employed as a standard for local image matching.
  • Figures 4(e) 1 and 4(e)(2) illustrate in more detail than the highly schematic Figures 4(a)-4(d) the processing performed in the iterative process of the present invention.
  • the process begins with obtaining first and second temporal sequential digital images of a subject (step 4000). Then these first and second digital images are preprocessed in step 4010, as above discussed in connection with step 420 of Fig. 4(a).
  • step 4020 initial shift values are determined by selecting plural template regions of interest (ROIs) in the first digital image and corresponding search area ROIs in the second image, and determining shift values between pixels centered in each template ROI and pixels centered in a respective search area ROI and exhibiting highest cross-correlation with respect to said template ROI, [1,3]
  • step 4020 shift vectors are determined for each of said shift values; a first cumulative histogram of the shift vectors for pixels in one lung, and a second cumulative histogram of the shift vectors for pixels in the other lung are also produced.
  • the shift vector histograms are used to select two sets of shift values for two-dimensional fitting.
  • the process includes selecting, based on said first histogram, a first set of pixels in the one lung with shift vectors within a ⁇ 90° range of angles of a peak in the first histogram of shift vectors and selecting the second set of pixels as the remaining pixels in the one lung, and selecting, based on the second histogram, a third set of pixels in the other lung with shift vectors within a ⁇ 90° range of angles of a peak in the second histogram of shift vectors and selecting the fourth set of pixels as the remaining pixels in the other lung.
  • the shift values of the first and third sets of pixels are fitted using polynomial fitting to produce a first fitted set of shift values, and likewise the shift values of the second and fourth sets of pixels are fitted using polynomial fitting to produce a second fitted set of shift values.
  • the method may include selecting, based on the first histogram, a first set of pixels in the one lung with shift vectors within a ⁇ 90° range of angles of a peak in the first histogram of shift vectors and selecting the second set of pixels as the remaining pixels in the one lung (assuming a peak exists in the histogram of that lung), and selecting, based on the second histogram, a third set of pixels in the other lung with shift vectors within an arbitrary range, e.g., a 0-180° range, of angles in the second histogram of shift vectors and selecting the fourth set of pixels as the remaining pixels in the other lung (assuming a peak does not exist in the shift vector histogram for that lung), and then forming first and second sets of fitted shift values as previously noted. If no peak
  • step 4040 the first and second sets of fitted shift values derived in step 4030 are used to wa ⁇ one of the original sequentially temporal images to derive two wa ⁇ ed images.
  • the method then proceeds in step 4050 by producing first and second interim subtraction images between the first digital image and the two wa ⁇ ed image, respectively.
  • the two interim subtraction images are then used to derive a final first wa ⁇ ed image for further iterative processing.
  • the method proceeds in step 4060 to determine a histogram of pixel values in lung regions of the first and second subtraction images and to determine if one of the first and second subtraction images exhibits a narrower histogram of pixel values in both lung regions.
  • the method proceeds by selecting the wa ⁇ ed image from which the subtraction image having the narrower histograms was produced for further iterative wa ⁇ ing.
  • the method proceeds in step 4080 by selecting for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second wa ⁇ ed images exhibiting the narrower histogram of pixel values, performing further two- dimensional fitting in step 4090 on the further selected fitted shift values, and wa ⁇ ing, in step 4100, the second digital image using the further fitted shift values to produce a final first wa ⁇ ed image for further iterative wa ⁇ ing.
  • a further iteration of wa ⁇ ing commences between the final first wa ⁇ ed image and the other original digital image by determining shift values, shift vectors and shift vector histograms etc.
  • step 4120 a subtraction image for diagnosis is produced by subtraction of the final second wa ⁇ ed image from the non- wa ⁇ ed original digital image.
  • the subtraction images scored as 1 or 2 may be considered as poor and inadequate, and would need to be improved for clinical use.
  • the subtraction images scored as 3, 4, or 5 are good subtraction images, which would be adequate for clinical use.
  • the number of adequate subtraction images increased from 78.5% to 97.7% with the new scheme. It is clear thus that the performance of the temporal subtraction was improved substantially by use of the iterative image wa ⁇ ing method of the present invention.
  • the quality of the 181 subtraction images were evaluated by using the histogram widths.
  • the histogram widths for each subjective rating group are shown in Figures 9 (a)-9(e).
  • the distributions of the histogram widths of the subtraction images obtained with the previous and the new temporal subtraction scheme are plotted with dots and x, respectively.
  • the results show that the histogram widths of the subtraction images obtained with the iterative wa ⁇ ing method of the present invention tend to be small and distributed in the lower left on the graphs. This indicates that the misregistration errors in the subtraction images obtained with the new method are smaller than those obtained with the previous method.
  • Subtraction images obtained with the previous and the new method are shown in Figs. 10 (a) and 10(b), respectively. It is apparent that the registration errors in the subtraction image obtained with the previous method are decreased substantially with the new temporal subtraction method, as shown in Fig. 10 (b).
  • the relative change in the quality of temporal subtraction images was subjectively evaluated by comparing two subtraction images obtained with the previous and the new methods. Using a subjective rating scale from -2 to 2, the quality of the new subtraction image compared to the previous subtraction image is rated as +2, clearly improved;
  • This invention may be conveniently implemented using a conventional general pu ⁇ ose digital computer or micro-processor programmed according to the teachings of the present specification, as will be apparent to those skilled in the computer art.
  • Appropriate software coding can readily be prepared by skilled programmers based on the teachings of the present disclosure, as will be apparent to those skilled in the software art.
  • the present invention includes a computer program product which is a storage medium including instructions which can be used to program a computer to perform a process of the invention.
  • the storage medium can include, but is not limited to, any type of disk including floppy disks, optical discs, CD-ROMs, and magneto-optical disks, ROMs, RAMs, EPROMs, EEPROMs, magnetic or optical cards, or any type of media suitable for storing electronic instructions.
  • FIG 12 is a schematic illustration of a general pu ⁇ ose computer 100 programmed according to the teachings of the present invention.
  • the general pu ⁇ ose computer 100 includes a computer housing 102 having a motherboard 104 which contains a CPU 106 and memory 108.
  • the computer 100 also includes plural input devices, e.g. , a keyboard 122 and mouse 124, and a display card 110 for controlling monitor 120.
  • the computer system 100 further includes a floppy disk drive 114 and other removable media devices (e.g., tape, and removable magneto-optical media (not shown)), a hard disk 112, or other fixed, high density media drives, connected using an appropriate device bus, e.g., a SCSI bus or an Enhanced IDE bus.
  • the computer 100 may additionally include a compact disc reader/ writer 118 or a compact disc jukebox (not shown).
  • the present invention includes programming for controlling both the hardware of the computer 100 and for enabling the computer 100 to interact with a human user.
  • Such programming may include, but is not limited to, software for implementation of device drivers, operating systems, and user applications.
  • Such computer readable media further includes programming or software instructions to direct the general pu ⁇ ose computer 100 to perform tasks in accordance with the present invention.
  • general pu ⁇ ose computer 100 may include a software module for digitizing and storing PA radiographs obtained from an image acquisition device.
  • the present invention can also be implemented to process digital data derived from a PA radiograph elsewhere.
  • the invention may also be implemented by the preparation of application specific integrated circuits or by interconnecting an appropriate network of conventional component circuits, as will be readily apparent to those skilled in the art.

Abstract

A method of computerized analysis of temporally sequential digital images, including (a) determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image (440); (b) warping the second digital image based on the first shift values to obtain a first warped image in which spatial locations of pixels are varied in relation to the first shift values (450); (c) determining second shift values between pixels of the first digital image and pixels of the first warped image; and (d) warping the first warped image based on the second shift values to obtain a second warped image (470). Iterative warping enhances image for the first subtraction of the first digital image and the final warped image to produce a difference image (480).

Description

Method, System And Computer Readable Medium For Iterative Image Warping Prior to Temporal Subtraction of Chest Radiographs in The Detection of Interval Changes
CROSS-REFERENCE TO RELATED APPLICATIONS AND PUBLICATIONS
The present is related to automated techniques for automated detection of abnormalities in digital images, for example as disclosed in one or more of U.S. Patents 4,839,807; 4,841,555; 4,851,984; 4,875,165; 4,907,156; 4,918,534; 5,072,384; 5,133,020; 5,150,292; 5,224,177; 5,289,374; 5.319,549; 5,343,390; 5,359,513; 5,452,367; 5,463.548; 5,491,627; 5,537,485; 5.598,481 ; 5,622,171; 5,638,458; 5,657,362; 5,666,434; 5.673,332; 5,668,888; as well as U.S. applications 08/158,388; 08,173,935; 08/220,917; 08/398,307; 08/428,867; 08/523,210; 08/536,149; 08/536,450; 08/515,798; 08/562,188; 08/562,087; 08/757,611; 08/758,438; 08/900,188; 08/900,189; 08/900,191; 08/900,192; 08/900,361; 08/900,362; 08/979,623; 08/979,639; 08/982.282; 09/028,518; and 09/027,685. each of which are incorporated herein by reference in their entirety. Of these patents and applications, 5,319,549 and 08/900,362 are of particular interest.
The present invention also relates to technologies referenced and described in the references identified in the appended APPENDIX and cross-referenced throughout the specification by reference to the number, in brackets, of the respective reference listed in the APPENDIX, the entire contents of which are also incorporated herein by reference. Various of these publications may correspond to various of the cross-referenced patents and patent applications.
STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH
The present invention was made in part with U.S. Government support under grant numbers CA 62625 and CA 64370 (National Institutes of Health). The U.S. Government has certain rights in the invention.
BACKGROUND OF THE INVENTION Field of the Invention
The present invention is related to temporal analysis of medical images and, in particular, to the analysis of chest radiograph images using automated temporal subtraction. Discussion of the Background
For the interpretation of chest radiographs, radiologists commonly compare a current film with previous films in order to facilitate the detection of abnormalities on chest radiographs, such as pulmonary nodules, interstitial infiltrates, pleural effusions, and cardiomegaly. However, it is a difficult task for radiologists to identify subtle interval changes on chest radiographs because lesions can overlap with anatomic structures such as ribs, vessels, heart, and diaphragm. In order to assist radiologists in their evaluation of temporal changes in chest radiographs, investigators at the University of Chicago Department of Radiology have developed a temporal subtraction scheme based on a nonlinear geometric warping technique. [1] In this scheme, the subtle changes and/or newly developed abnormalities on chest radiographs were enhanced on the subtracted image. In fact, in an observer test, the detection accuracy of interval changes in chest radiographs was improved significantly by use of temporal subtraction images. [2] In a more recent study, the temporal subtraction scheme has been improved by an initial image matching technique which is based on the cross-correlation of a pair of blurred low-resolution images for determination of accurate global shift values. [3] However, misregistration errors caused by failure in local image matching were still observed in some cases.
SUMMARY OF THE INVENTION It is an object of the present invention to provide improved automated temporal subtraction of images derived from chest radiographs.
Another object of this invention is to provide a novel image processing technique to overcome severe misregistration errors mainly due to differences in a subject's inclination and/or rotation in temporally sequential images.
A further object of this invention is to provide a novel temporal image subtraction technique to assist radiologists in the detection of interval changes on chest radiographs.
These and other objects are achieved according to the present invention by providing a novel method, system and computer readable medium for computerized analysis of temporally sequential digital images. The method includes the steps of (a) determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image; (b) warping said second digital image based on the first shift values to obtain a warped image in which spatial locations of pixels are varied in relation to said first shift values; (c) determining second shift values between pixels of said first digital image and pixels of said warped image obtained in step (b); (d) waφing said warped image obtained in step (b) based on the second shift values to obtain an iteratively waφed image in which spatial locations of pixels of said waφed image obtained in step (b) are varied in relation to said second shift values; and (e) subtracting the iteratively waφed image from said first digital image. Additional iterations of image waφing are possible to enhance image registration between the first digital image and the waφed version of the second digital image, followed by image subtraction of the first digital image and the final waφed image to produce a difference image from which diagnosis of temporal changes ensues. In the production of shift values used in image waφing, initial shift values are obtained by cross- correlation techniques using a template and search regions of interest. Shift vectors and a histogram of shift vectors are obtained from each initial shift value. Based on the histogram of shift vectors, variations in various groups of the initial shift values are smoothed using two-dimensional fitting.
The present invention similarly includes a computer readable medium storing program instructions by which the method of the invention can be performed when the stored program instructions are appropriately loaded into a computer, and a system for implementing the method of the invention.
BRIEF DESCRIPTION OF THE DRAWINGS
A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:
Figure 1 is an illustration of selecting different ROIs for analysis of histograms of subtraction images.
Figures 2(a)(A) and 2(a)(B) are photographs of subtraction images with ROI placement for analysis of histograms for pixel values, with the photograph of Fig. 2(a)(A) being an example of a poor subtraction image and the photograph of Fig. 2(a)(B) being an example of a good subtraction image.
Figure 2(b) is a graph of histograms of subtraction images of Figures 2(a)(A) and 2(a)(B).
Figure 3(a)-3(e) are illustrations of distributions of histogram widths in both lungs for cases rated as (a) 1 (very poor quality), (b) 2 (poor quality), (c) 3 (adequate quality), (d) 4 (good quality), and (e) 5 (excellent quality), respectively.
Figure 4(a) is a schematic diagram of the temporal subtraction method using an iterative waφing technique according to the present invention.
Figures 4(b)-4(d) are highly schematic diagrams of shift vector analysis, first waφing, and second waφing steps, respectively, included in the method of Figure 4(a).
Figures 4(e) 1 and 4(e)(2) illustrate in more detail than the highly schematic Figures 4(a)-4(d) the processing performed in the iterative process of the present invention.
Figures 5(a) and 5(b) are photographs of respective subtraction images obtained using different shift vectors, i.e., (a) good subtraction image and (b) a poor subtraction image, wherein shift vectors are overlaid on each image, dark and light vectors (short lines) are used for inside and outside, respectively, of approximate lung areas, and small white dots at one end of dark vectors indicate initial locations of template ROIs, which are used in the waφing technique
Figures 6(a) and 6(b) are graphs respectively illustrating a distribution of accumulated shift vector values as a function of shift vector orientation for (a) a good subtraction case and (b) a poor subtraction case.
Figures 7(a)-7(d) are photographs respectively illustrating (a) dominant shift vectors (dark lines) in a temporal image, (b) subtraction image; (c) remaining shift vectors (dark lines) and (d) a subtraction image.
Figures 8(a), 8(b) and 8(c) are photographs of subtraction images respectively obtained with (a) a single waφing technique, (b) an iterative waφing technique, and (c) an iterative waφing technique plus linear inteφolation of shift values.
Figures 9(a), 9(b), 9(c), (d), and 9(e) are graphs respectively illustrating distributions of histogram widths in both lungs for the 181 cases, obtained with the prior waφing technique (dots) and the waφing technique of the present invention (x), rated as 1 (very poor quality), 2 (poor quality), 3 (adequate quality), 4 (good quality), and 5 (excellent quality).
Figures 10(a) and 10(b) are photographs of subtraction images for comparison, respectively obtained (a) with the prior waφing technique and (b) with the waφing technique of the present invention.
Figure 11 is a bar graph illustrating the distribution of the degree of improvement, with the subjective rating method, in the quality of the subtraction images obtained with the waφing technique of the present invention.
Figure 12 is a schematic illustration of a general puφose computer 100 programmed according to the teachings of the present invention.
DESCRIPTION OF THE PREFERRED EMBODIMENTS
In the derivation of the present invention, 181 pairs of chest radiographs with the current temporal subtraction scheme were examined. Initial image matching technique as disclosed in U.S. patent 5,319,549 was employed. One hundred forty-two (78.5%) of 181 subtraction images showed adequate, good, or excellent quality. However, the remaining 39 (21.5%) cases were totally or partly misregistered, mainly due to A-P inclination and/or insufficient rotation correction. In order to reduce these misregistration errors, we developed a new temporal subtraction technique, as disclosed hereinafter, was developed by applying an iterative image waφing technique. Materials and Methods Materials
The image database used in the derivation of the present invention included 181 pairs of chest radiographs obtained from the Iwate Prefectural Hospital, Morioka, Japan. These 181 cases were obtained sequentially from chest screening images which were made with a Fuji Computed Radiography (FCR) system (Fuji Medical Systems Co., Ltd., Tokyo, Japan). The time interval between the current and previous images for all cases was 1 year. The pixel size and gray level of CR chest images were 0.2 mm and 1024, respectively. The technique of the present invention was developed by use of a Silicon Graphics 02 workstation. Methodology 1. Subjective evaluation of the quality of subtraction images
To evaluate the quality of subtraction images obtained with the previous and the new subtraction scheme, subjective ratings were obtained from two chest radiologists and two physicists independently. A five-point rating scale was employed, that is,
1. (very poor): most ribs were not registered and appear in the entire intercostal space,
2 (poor): most ribs were not well registered and appear in half of the intercostal space,
3 (adequate): most ribs were well registered, with some minor misregistration error, 4 (good): most ribs were almost completely registered with very minor misregistrations partly, and
5 (excellent): all ribs were perfectly registered.
The final rating for each case was determined in one of the five categories above based on the average (or the majority) of all observers' ratings.
2. Objective evaluation of the quality of subtraction images
For objective evaluation of the quality of a subtraction image, we employed the width of the histogram of pixel values for the right and the left lung area of the subtraction image was employed. It was found that in the case of subtraction images without an actual interval change, a high-quality subtraction generally yields a low-contrast image, whereas the contrast of poor subtraction images tends to be high, because poor subtraction images partly contain localized mixtures of very dark and very light areas due to misregistration. Therefore, the histogram of pixel values in each lung field of the subtraction image was employed to evaluate the quality of subtraction images. [3]
For determination of the histogram of pixel values in each lung field of the subtraction image, various ROIs of different size and shape in each lung were examined, as shown in Figs. 1(a)- 1(d). Initially, the ROI included all of the lung and mediastinum areas, as illustrated in Fig. 1(a). Second, we segmented the fight and left lungs were segmented by using the heart border information, as shown in Fig. 1(b), in order to evaluate the quality of the subtraction image in each lung field separately. However, it was found that the top and the bottom areas of the ROIs tended to include misregistration errors caused by the difference between the clavicle positions, and also the difference between different diaphragm levels, respectively. Therefore, third, the areas around the ribcage edge boundaries and the boundaries of the mediastinum were excluded from the ROIs, as shown in Fig. 1(c). Finally, small ROIs as shown in Fig. 1(d) were employed, because the extent of the misregistration error for the ribs can be detected more sensitively with the ROI placed over the central area of the lung rather than with the other ROIs shown in Figs. 1(a), 1(b), and 1(c). The size of the small ROIs for both lungs was 30 pixels (width) x 120 pixels (length). The centers of the ROIs were determined by using the x,y coordinates of the ribcage edges and the cardiac edges. The vertical locations of the two centers of the ROIs were the same and equal to the y-location between the top lung, which was obtained from the ribcage edge detection scheme, and the bottom of the right ribcage edge. The horizontal centers of the ROIs were determined as the x-location in each lung between the ribcage edge and the cardiac edge, at the vertical center level of the ROIs.
Figure 2 (a) shows a comparison of a poor (A) and a good (B) subtraction images of the right lung. The ROIs used for determination of the histogram width are illustrated by black lines. The corresponding histograms of these images are shown in Figure 2 (b). It is seen that the width of the histogram for the poor subtraction image (A) is much wider than that for the good subtraction image (B). Although an actual interval change in chest images can broaden the histogram of the subtraction image, the "misregistration" (the difference between an abnormal current image and a normal previous image) due to the actual interval abnormal change would usually be localized and small compared with the misregistration errors due to a failure in image matching. Therefore, the quality of the subtraction image was evaluated by using the widths of the histograms in the ROI for each lung.
As a measure for evaluation of the quality of subtraction images, the width of the histogram at 10%) of the largest peak was obtained. It was found that the width of the histogram at a low level such as 10%> of the peak was more sensitive than the histogram width at a high level in distinguishing between good and poor subtraction images. The distributions of the histogram widths for each of the five groups of different qualities of subtraction images, which were grouped based on subjective ratings, obtained by the previous technique are shown in Fig. 3 (a)-3(e). In these figures, the horizontal and vertical axes correspond to the width of the histogram in the right and the left lung, respectively. The histogram widths for "very poor" and "poor" subtraction images (subjective rating scores: 1 and 2) tend to be large and are distributed in the upper right, as shown in Figs. 3 (a) and 3(b). The histogram widths of "adequate" subtraction images (subjective rating: 3) are distributed in the intermediate range, as shown in Fig. 3 (c). The histogram widths for "good" subtraction images (subjective rating: 4) are shifted toward the lower left, as shown in Fig. 3 (d). The "excellent" subtraction images as represented by a subjective rating score of 5 are characterized by narrow histogram widths, as shown in Fig. 3 (e). 3. Overall scheme of the new temporal subtraction technique
In our previous temporal subtraction scheme, a nonlinear density correction technique [1,7] was applied first for adjustment of the density and contrast in the two digitized images, i.e., the current and previous chest radiographs. The difference in the lateral inclination between the two images was corrected by use of an image rotation technique. The global shift values were then determined for the initial registration by use of the cross-correlation of a pair of blurred low-resolution images obtained from the current and the previous image. [3] For local matching, a number of template ROIs (32 x 32 matrix) and the corresponding search area ROIs (64 x 64 matrix) were selected from the previous and the current image, respectively. The number of pairs of ROIs was approximately 300 per case. Shift values, Δx and Δy, for all pairs of selected ROIs were determined by using a cross-correlation technique to find the "best" matched areas in the current and previous images. A two-dimensional surface fitting by use of polynomial functions was then applied to each set of mapped shift values, Δx and Δy, for conversion of the x, y coordinates of the previous image, i.e., for waφing of the image. The waφed previous image was then subtracted from the current image. [1] These steps are preferably also employed in practice of the present invention.
The temporal subtraction technique is illustrated in overview in Fig. 4 (a). The matrix size for the current and the previous chest image obtained with a CR system is reduced to 586 x 586 (one third of the original image matrix size). In this study, a density correction technique was not applied because the chest images obtained with the CR system can maintain consistent density and contrast with the use of the exposure data recognition (EDR) system [4] included in most recent CR systems.
In order to correct for the lateral inclination caused by the variation in patient positioning between the current chest image 400 and the previous image 410 schematically shown in Fig. 4(a), global matching 420 is performed using an image rotation [2], determining the angle between the two images by comparison of the two midlines of these images [5], and applying an initial image matching technique for determination of the global shift value between the two images, which corresponds to the shift in the x, y coordinates of one image relative to the other. [3] The two images are globally registered by use of the initial image matching technique based on the cross-correlation of blurred low-resolution images.
In the method of the present invention, an iterative image waφing technique is employed to improve local image matching and thus to reduce registration errors. To select the ROIs for the local image matching, local segmentation 430 is performed in which the ribcage edges of the previous image are detected based on image profile analysis. [5] In addition, the cardiac edges are determined based on an edge detection technique for segmentation of the lungs of the current image. [6] In step 440, the shift values Δx and Δy, for all of the selected ROIs are determined by a cross-correlation technique. [1] Shift vectors based on the shift values Δx and Δy are then determined, as discussed in more detail hereinafter. A first image waφing (step 450), typically of the previous or older image, is performed to derive a first waφed image (step 460). Then a second waφing (step 470) is performed on the first waφed image, followed by obtaining of a subtraction image 480 for diagnosis.
Fig. 4(b) illustrates very schematically one embodiment of determination of fitted shift values for the first waφing step. First, once the initial values of shift values are obtained [1,3], a shift vector orientation histogram is obtained (step 411) for each lung. In step 412 peaks in the histogram of each lung are determined, and utilized to identify dominant and remaining shift vectors (step 413). After further processing to derive optimized shift values, a first waφing is performed.
The steps for the first and the second image waφing technique are very schematically shown in Figs. 4(b) and 4(d), respectively. With the first image waφing of step 450, in step 451 optimized shift values are determined and then in step 452 subjected to surface fitting [1]. After fitting, a coordinate transformation is performed in step 453 and a waφed image obtained in step 454.
As shown in Fig. 4(d), a second image waφing technique, which is applied to the waφed previous image, is then performed. The optimized shift values between the current image and the waφed previous image are determined by the cross-correlation technique and vectors are again determined in step 471 to derive optimized shift values, as discussed in more detail hereinafter. A linear inteφolation is performed on the optimized shift values in step 472 to determine the final shift values on all of the x,y coordinates of the waφed previous image. Another coordinate transformation is performed in step 473 to obtain a second waφed image in step 474. Finally, the temporal subtraction image is obtained by subtraction of the second waφed previous image from the current chest image, as shown schematically in step 480 of Fig. 4(a). 4. Weighting factors for surface fitting of shift values In the previous technique, weighting factors determined based on the cross-correlation values were used for surface fitting to determine the local shift values. [1] In general, the larger the cross-correlation value, the larger the weighting factor. However, it was found that the cross-correlation values near the ribcage borders, in the mediastinum area, and below the diaphragm area were generally very large, and much greater than the values in the lung areas. Therefore, the fitted shift values would have been affected considerably by the shift values near the ribcage borders, in the mediastinum, and below the diaphragms area. This is not desirable, because accurate subtraction is commonly required in the lung fields rather than in other areas in chest images. In order to avoid this problem, a lung segmentation method, which includes the detection of ribcage edges and heart edges for identifying the shift values in the lung area, is employed. A weighting factor of 1.0 is assigned for ROIs in the lung areas, and 0.25 for ROIs in the mediastinum and below the diaphragm, whereupon weighted polynomial fitting is performed. [1] In this way, the fitted shift values are determined for the first waφing step. When the new weighting factors were used, the registration errors in the lung areas were reduced in some cases. 5. Shift- vector orientation analysis
The surface fitting technique of the shift values was generally effective in improving local registration around the poorly matched regions. [1] However, if some shift vectors with incorrect orientation are included within a relatively small region, the correct shift vectors may be affected significantly by the fitting. In the previous method [1.3], the orientation (or angle) of the shift vectors was not considered in the process for determining the surface fitting of the shift values. However, according to the present invention it has since been determined that the orientation of the shift vectors is an important factor for accurate surface fitting, and therefore two different components, i.e., the dominant and the remaining (non-dominant) shift vectors, are identified by using the shift vector orientation histogram, which is a histogram of accumulated shift values as a function of the shift vector orientation.
In the method of the present invention, first, the ribcage edges and the cardiac edges are determined by using the edge detection technique [6] for segmentation of the lungs. Many ROIs are selected in the lung field for both the previous and current images. Shift values, Δx and Δy, are determined for each pair of ROIs based on the cross-correlation technique. [1] The orientation (or angle) of the shift vector for each ROI is determined by the arc tangent of Δy/Δx. Then the shift- vector orientation histograms for the right and the left lung are obtained.
Figures 5 (a) and 5(b) illustrate the distributions of shift vectors obtained with a good and a poor subtraction image, respectively. As shown in Fig 5 (a), the orientations of the shift vectors for good subtraction images tend to be similar. This may be because the rib contrast in this chest image is relatively high, and thus the template ROI which includes a rib edge could easily match the corresponding rib edge in the search area ROI. The orientations of the shift vectors for a poorly registered case were varied because of the low-contrast ribs. The shift-vector orientation histograms of these cases, Figs. 5(a) and 5(b), for the right and left lungs are shown in Fig. 6(a) and 6(b). It is apparent that the shift-vector orientation histograms for a good registered case have a large peak in each lung. However, there is no obvious peak in the orientation histograms for a poor case, as shown in Fig. 6(b).
If a shaφ peak exists in the shift- vector orientation histogram, and if the largest peak value is greater than twice the average of accumulated shift values, then the shift vectors within ± 90° of the peak angle are considered as the dominant orientation. The remaining vectors outside this angle range are considered as the non-dominant orientation. In that case, subsequently curve fitting of shift values (See U.S. patent 5,359,513) is performed based only on shift vectors within the ± 90° peak angular range, all other shift values being ignored during curve fitting. However, if a dominant orientation does not exist, then the shift vectors in the upper direction (from 0° to 180°) and lower direction (from 180° to 360°) are separately employed for curve fitting of the shift values. In other words, curve fitting of shift values for only those pixels in both lungs having shift vectors in the upper direction (from 0 ° to 180°) is performed separately of curve fitting of shift values for the remaining pixels of both lungs having shift vectors in the lower direction (from 181 ° to 360 °). Figures 7(a) and 7(b) show the selected dominant shift vectors and the corresponding subtraction image obtained by fitting with the dominant shift vectors, respectively. Figures 7(c) and 7(d) illustrate non-dominant shift vectors and the corresponding subtraction image obtained by fitting with the non-dominant shift vectors, respectively. It is clear that the subtraction image obtained with the dominant shift vectors is superior to that obtained with the non-dominant shift vectors. Figures 8(a) and 8(b) show the temporal subtraction images obtained with the prior technique and that of the present invention using shift- vector analysis, respectively. The quality of the subtraction image obtained according to the present invention is improved slightly.
To determine the efficacy of curve fitting of shift values based on the concepts of dominant and non-dominant shift vectors, as above described, subtraction images were obtained using a waφed image in which curve fitting limited to dominant shift vectors was performed for both lungs, and using a waφed image in which curve fitting using non- dominant shift vectors was performed. As measure of how image registration between two images, a histogram of pixel values of pixels in the lung regions was derived from each subtraction image. Subtraction images with narrow widths of histograms in the lung regions occur with good image registration since in that case relatively complete subtraction of rib edges, which contribute to wide histograms, occurs. Most of the good subtraction images with narrow pixel histogram widths, i.e, subtraction images derived from images with good registration and in which regions of the lungs have pixel values exhibiting a narrow histogram width, were obtained by use of curve fitting limited to fitting of dominant shift vectors for both lungs. However, in some cases, it is possible that the width of the pixel value histogram in one lung of the subtraction images obtained with the non-dominant shift vectors is narrower than that obtained with the dominant shift vectors because of a rotation error between the current and the previous image. Therefore, when the histogram width of lung pixel values of the subtraction image obtained with the dominant shift vectors is narrower in one lung, but wider in the other lung than the corresponding histogram width obtained with the non-dominant shift vectors, whichever shift values result in smaller pixel value histogram widths of lung pixels for a particular lung are used for a further curve fitting of shift values . In other words, a pair of waφed images of one of the digital images, such as the older in time or previous image, are obtained based on fitted dominant and fitted non-dominant shift values. Then a pair of subtraction images is produced between the current or recently obtained image and each of the waφed images of the previous image. These subtraction images are then evaluated to determine which of the two waφed images results in the better subtraction, i.e., reduced residual permanent structure such as occurs due to ribs, etc.. If a single subtraction image exhibits a narrower histogram of pixel values in both lungs, indicating a more uniform texture due to good subtraction of, e.g., ribs, than it is decided that no further improvement in shift values is necessary and the waφed image which resulted in the better subtraction is selected for further iterative processing.
If, however, one subtraction image exhibits a narrower histogram of pixel values in one lung and the other subtraction image exhibits a narrower histogram of pixel values in the other lung, i.e., one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, the present invention then selects for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, and then performs further two-dimensional fitting on the further selected fitted shift values, and then waφs the previous image using said further fitted shift values to produce a waφed image for further iterative processing. Thus, the waφed previous image is obtained by means of the fitted shift values based on the better subtraction image in each lung. Details of the technique of the invention are discussed further hereinafter in connection with Figures 4(e)(1) and 4(e)(2). 6. Iterative image waφing technique
According to the present invention, the image waφing technique is applied repeatedly and thus iteratively first on the previous image, next on the first waφed image, and then on the second waφed image, and so on, until the desired quality of the subtraction image is obtained. The second and subsequent image waφing steps are employed for improvement of the quality of the temporal subtraction image.
In order to perform the second image waφing, the waφed previous image from the first waφing step is obtained. The current image and the waφed previous image are then used for the second image waφing. The shift values, Δx and Δy, are determined by use of the cross-correlation technique for all of the selected ROIs in the current image and the waφed previous image. As in the first waφing step, shift vectors and a histogram of shift vectors are derived, angular ranges determined, and two sets of shift vectors again selected for production of interim first and second twice waφed images. Once again, corresponding subtraction images are obtained, and the histograms of pixel values in each subtraction image are evaluated. If one subtraction image exhibits a narrower histogram of pixel values in both lung regions in comparison to that of the other subtraction image, then that subtraction image is considered the final subtraction image for diagnosis (assuming that no further waφing iterations are to be performed). If however that is not the case, then another shift value fitting is performed. In that case, the present invention then selects for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, and then performs further two-dimensional fitting on the further selected fitted shift values, and then waφs the previously waφed image using the further fitted shift values to produce a secondly waφed image for further iterative processing, if desired, or for production of the final subtraction image for diagnosis.
The shift values of all x,y coordinates over the entire lung fields for the second image waφing are preferably obtained by a linear inteφolation technique instead of the surface fitting technique with polynomial functions used during the first waφing, since linear inteφolation was observed to result in a somewhat better subtraction image after the second image waφing. However, the present invention can also be practiced using linear inteφolation for the first image waφing or polynomial fitting for the second image waφing.
Figure 8(c) shows a subtraction image obtained by application of the second image waφing technique. The registration errors for the ribs in the temporal subtraction image are reduced substantially by use of the second image waφing technique.
A third image waφing technique, which is basically the same as the second waφing process, can also be used. The quality of the subtraction image was improved slightly in some cases by applying the third image waφing; however, at the expense of increased computational time for the iterative waφing process. Thus, to avoid lengthy computational times, two iterations of iterative image waφing are preferably employed, i.e., the second image waφing technique is employed as a standard for local image matching.
Figures 4(e) 1 and 4(e)(2) illustrate in more detail than the highly schematic Figures 4(a)-4(d) the processing performed in the iterative process of the present invention. The process begins with obtaining first and second temporal sequential digital images of a subject (step 4000). Then these first and second digital images are preprocessed in step 4010, as above discussed in connection with step 420 of Fig. 4(a). In step 4020 initial shift values are determined by selecting plural template regions of interest (ROIs) in the first digital image and corresponding search area ROIs in the second image, and determining shift values between pixels centered in each template ROI and pixels centered in a respective search area ROI and exhibiting highest cross-correlation with respect to said template ROI, [1,3] In step 4020, shift vectors are determined for each of said shift values; a first cumulative histogram of the shift vectors for pixels in one lung, and a second cumulative histogram of the shift vectors for pixels in the other lung are also produced. In step 4030, the shift vector histograms are used to select two sets of shift values for two-dimensional fitting. This may preferably be done based on a detection of peaks in the shift vector histograms of each lung, or in by arbitrarily selecting predetermined angular ranges of shift vectors. If the former approach is implemented, the process includes selecting, based on said first histogram, a first set of pixels in the one lung with shift vectors within a ± 90° range of angles of a peak in the first histogram of shift vectors and selecting the second set of pixels as the remaining pixels in the one lung, and selecting, based on the second histogram, a third set of pixels in the other lung with shift vectors within a ± 90° range of angles of a peak in the second histogram of shift vectors and selecting the fourth set of pixels as the remaining pixels in the other lung. Then, the shift values of the first and third sets of pixels are fitted using polynomial fitting to produce a first fitted set of shift values, and likewise the shift values of the second and fourth sets of pixels are fitted using polynomial fitting to produce a second fitted set of shift values.
It may occur that a clear peak does not exist in the histogram of shift vectors of one lung. In that case, an arbitrary range of angles may be selected for fitting with dominant and non-dominant shift values of the other lung. For example, in that case, the method may include selecting, based on the first histogram, a first set of pixels in the one lung with shift vectors within a ± 90° range of angles of a peak in the first histogram of shift vectors and selecting the second set of pixels as the remaining pixels in the one lung (assuming a peak exists in the histogram of that lung), and selecting, based on the second histogram, a third set of pixels in the other lung with shift vectors within an arbitrary range, e.g., a 0-180° range, of angles in the second histogram of shift vectors and selecting the fourth set of pixels as the remaining pixels in the other lung (assuming a peak does not exist in the shift vector histogram for that lung), and then forming first and second sets of fitted shift values as previously noted. If no peak exists in either shift vector histogram, than arbitrary ranges of angles can be selected in the histograms of shift vectors for both lungs and sets of shift values derived from both lungs then derived and fitted.
In step 4040, the first and second sets of fitted shift values derived in step 4030 are used to waφ one of the original sequentially temporal images to derive two waφed images. The method then proceeds in step 4050 by producing first and second interim subtraction images between the first digital image and the two waφed image, respectively. The two interim subtraction images are then used to derive a final first waφed image for further iterative processing.
In particular, the method proceeds in step 4060 to determine a histogram of pixel values in lung regions of the first and second subtraction images and to determine if one of the first and second subtraction images exhibits a narrower histogram of pixel values in both lung regions. When one of the subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, in step 4070 the method proceeds by selecting the waφed image from which the subtraction image having the narrower histograms was produced for further iterative waφing. When one of the subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions, the method proceeds in step 4080 by selecting for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two- dimensional fitting in step 4090 on the further selected fitted shift values, and waφing, in step 4100, the second digital image using the further fitted shift values to produce a final first waφed image for further iterative waφing. In step 4110 a further iteration of waφing commences between the final first waφed image and the other original digital image by determining shift values, shift vectors and shift vector histograms etc. in much the same way as above discussed in connection with above discussed steps 4020-4100, with the exception that linear inteφolation two-dimensional fitting is performed rather polynomial two- dimensional fitting as above noted, to derive a final second set of shift values for waφing of the final first waφed image in step 4120. In this way a second iteration of waφing is performed. Once the final second waφed image is obtained in step 4120, and assuming that no further waφing iterations are necessary, in step 4130 a subtraction image for diagnosis is produced by subtraction of the final second waφed image from the non- waφed original digital image. Results By using the subjective rating method described earlier, the quality of 181 subtraction images obtained independently with the previous temporal subtraction scheme and the iterative subtraction scheme were evaluated. The results of the subjective evaluation in terms of the number of cases in each rating scale are shown in Table 1 for subtraction images obtained with the previous image waφing technique and with the new iterative waφing technique of the present invention.
TABLE 1
Subjf active rating score
1 2 3 4 5
(very poor) (adequate) (excellent
Prior image 24 15 72 55 15 waφing method
(13.2%) (8.3%) (39.8%) (30.4%) (8.3%)
New iterative 1 3 14 102 61 image waφing
(0.6%) technique (1.7%) (7.7%) (56.4%) (33.7%)
The subtraction images scored as 1 or 2 may be considered as poor and inadequate, and would need to be improved for clinical use. The subtraction images scored as 3, 4, or 5 are good subtraction images, which would be adequate for clinical use. The number of adequate subtraction images increased from 78.5% to 97.7% with the new scheme. It is clear thus that the performance of the temporal subtraction was improved substantially by use of the iterative image waφing method of the present invention.
On the other hand, the quality of the 181 subtraction images were evaluated by using the histogram widths. The histogram widths for each subjective rating group are shown in Figures 9 (a)-9(e). The distributions of the histogram widths of the subtraction images obtained with the previous and the new temporal subtraction scheme are plotted with dots and x, respectively. The results show that the histogram widths of the subtraction images obtained with the iterative waφing method of the present invention tend to be small and distributed in the lower left on the graphs. This indicates that the misregistration errors in the subtraction images obtained with the new method are smaller than those obtained with the previous method. Subtraction images obtained with the previous and the new method are shown in Figs. 10 (a) and 10(b), respectively. It is apparent that the registration errors in the subtraction image obtained with the previous method are decreased substantially with the new temporal subtraction method, as shown in Fig. 10 (b).
In addition, the relative change in the quality of temporal subtraction images was subjectively evaluated by comparing two subtraction images obtained with the previous and the new methods. Using a subjective rating scale from -2 to 2, the quality of the new subtraction image compared to the previous subtraction image is rated as +2, clearly improved;
+2, clearly improved;
+1, moderately improved; 0, unchanged;
-1, moderately declined; and
-2, clearly declined. Figure 11 shows the distribution of the rating scores for all subtraction images. The results show that 156 (86.2%>) of the 181 cases were improved by use of the new temporal subtraction method. Thus it is concluded that the new temporal subtraction method based on an iterative waφing technique enhances many subtle changes in chest images and results in a substantially improved overall performance of the temporal subtraction.
This invention may be conveniently implemented using a conventional general puφose digital computer or micro-processor programmed according to the teachings of the present specification, as will be apparent to those skilled in the computer art. Appropriate software coding can readily be prepared by skilled programmers based on the teachings of the present disclosure, as will be apparent to those skilled in the software art.
The present invention includes a computer program product which is a storage medium including instructions which can be used to program a computer to perform a process of the invention. The storage medium can include, but is not limited to, any type of disk including floppy disks, optical discs, CD-ROMs, and magneto-optical disks, ROMs, RAMs, EPROMs, EEPROMs, magnetic or optical cards, or any type of media suitable for storing electronic instructions.
Figure 12 is a schematic illustration of a general puφose computer 100 programmed according to the teachings of the present invention. The general puφose computer 100 includes a computer housing 102 having a motherboard 104 which contains a CPU 106 and memory 108. The computer 100 also includes plural input devices, e.g. , a keyboard 122 and mouse 124, and a display card 110 for controlling monitor 120. In addition, the computer system 100 further includes a floppy disk drive 114 and other removable media devices (e.g., tape, and removable magneto-optical media (not shown)), a hard disk 112, or other fixed, high density media drives, connected using an appropriate device bus, e.g., a SCSI bus or an Enhanced IDE bus. Also connected to the same device bus or another device bus, the computer 100 may additionally include a compact disc reader/ writer 118 or a compact disc jukebox (not shown).
Stored on any one of the above described storage media (computer readable media), the present invention includes programming for controlling both the hardware of the computer 100 and for enabling the computer 100 to interact with a human user. Such programming may include, but is not limited to, software for implementation of device drivers, operating systems, and user applications. Such computer readable media further includes programming or software instructions to direct the general puφose computer 100 to perform tasks in accordance with the present invention.
The programming of general puφose computer 100 may include a software module for digitizing and storing PA radiographs obtained from an image acquisition device. Alternatively, it should be understood that the present invention can also be implemented to process digital data derived from a PA radiograph elsewhere.
The invention may also be implemented by the preparation of application specific integrated circuits or by interconnecting an appropriate network of conventional component circuits, as will be readily apparent to those skilled in the art.
Obviously, numerous modifications and variations of the present invention are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described herein. APPENDIX REFERENCES:
1 A. Kano, K. Doi, H. MacMahon, D. D. Hassell, M. L. Giger, "Digital image subtraction of temporally sequential chest images for detection of interval change," Med. Phys. 21 : 453-461 (1994). (See U.S. patent 5,359,513)
2 M. C. Difazio, H. MacMahon, XW Xu, P Tsai, J. Shiraishi, S. G. Armato III, K. Doi, "Digital chest radiography: Effect of temporal subtraction images on detection accuracy," Radiology 202, 447-452 (1977).
3 T. Ishida, K. Ashizawa, R. Engelman, S. Katsuragawa, H. MacMahon, K. Doi, "Application of temporal subtraction for detection of interval change in chest radiographs: Improvement of subtraction images using automated initial image matching," Submitted to Med. Phys. (See U.S. application serial number 08/900,362)
4 H. Takeo, N. Nakajima, M. Ishida, H. Kato, "Improved automatic adjustment of density and contrast in FCR system using neural network," Proc. SPIE, 2163, 98-109 (1994).
5 X. W. Xu, K. Doi, "Image feature analysis for computer-aided diagnosis: Accurate determination of ribcage boundary in chest radiographs," Med. Phys. 22: 617-626, (1995).
6 N. Nakamori, K. Doi, V. Sabeti, H. MacMahon, "Image feature analysis and computer-aided diagnosis in digital radiography: Automated analysis of sizes of heart and lung in digital chest images," Med. Phys. 17, 342-350 (1990).
7 H. Yoshimura, X. W. Xu, K. Doi, H. MacMahon, K. R. Hoffmann, M. L. Giger, S.M. Montner, "Development of a high quality film duplication system using a laser digitizer: Comparison with computed radiography," Med. Phys. 20, 179-186 (1993).

Claims

Claims:
1. A method for computerized analysis of temporally sequential digital images, comprising the steps of:
(a) determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image;
(b) waφing said second digital image based on the first shift values to obtain a waφed image in which spatial locations of pixels are varied in relation to said first shift values;
(c) determining second shift values between pixels of said first digital image and pixels of said waφed image obtained in step (b);
(d) waφing said waφed image obtained in step (b) based on the second shift values to obtain an iteratively waφed image in which spatial locations of pixels of said waφed image obtained in step (b) are varied in relation to said second shift values; and
(e) subtracting the iteratively waφed image from said first digital image.
2. The method of Claim 1, further comprising the steps of: determining third shift values between pixels of said first digital image and pixels of said iteratively waφed image; and waφing said iteratively waφed image based on the third shift values to obtain a further iteratively waφed image in which spatial locations of pixels of said iteratively waφed image are varied in relation to said third shift values.
3. The method of Claim 1, wherein said step (a) comprises: selecting plural template regions of interest (ROIs) in the first digital image and corresponding search area ROIs in the second image, determining shift values between pixels centered in each template ROI and pixels centered in a respective search area ROI and exhibiting highest cross-correlation with respect to said template ROI, determining shift vectors for each of said shift values, a first cumulative histogram of said shift vectors for pixels in one lung, and a second cumulative histogram of said shift vectors for pixels in the other lung, selecting, based on characteristics of said first and second cumulative histograms of shift vectors, plural of said shift values for derivation of fitted shift values to serve as said first shift values, and performing two-dimensional fitting on the shift values selected in the preceding step to derive said fitted shift values which serve as said first shift values; and wherein said waφed image obtained in said step (b) is obtained using said fitted shift values.
4. The method of Claim 3, wherein: said step of selecting plural shift values comprises, selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said step of performing two-dimensional fitting comprises, performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said step (b) comprises, using said first set of fitted shift values to waφ said second digital image to obtain a first waφed image, using said second set of fitted shift values to waφ said second digital image to obtain a second waφed image, producing first and second subtraction images between said first digital image and said first waφed image, and between said first digital image and said second waφed image, respectively, using said subtraction images to derive the waφed image for further processing in said steps (c) and (d).
5. The method of Claim 4, wherein the step of using said subtraction images to derive the waφed image for further processing in said steps (c) and (d) comprises: determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said steps (c), and (d), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting for further two- dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two-dimensional fitting on the further selected fitted shift values, and waφing said second digital image using said further fitted shift values to produce said waφed image for further processing in said steps (c) and (d).
6. The method of Claim 4, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
7. The method of Claim 4, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
8. The method of Claim 4, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
9. The method of Claim 6, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
10. The method of Claim 7, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
11. The method of Claim 3, wherein said step of performing two-dimensional fitting on the shift values comprises: using a two-dimensional nth order polynomial function.
12. The method of Claim 3, wherein said step (c) comprises: selecting plural template regions of interest (ROIs) in the first image and corresponding search area ROIs in the waφed image obtained in step (b) , determining shift values between pixels centered in each template ROI in said first image and pixels centered in a respective search area ROIin said waφed image obtained in step (b) and exhibiting highest cross-correlation with respect to said template ROI, determining shift vectors for each of said shift values determined in the preceding step, determining a first cumulative histogram of said shift vectors determined in the preceding step for pixels in one lung, and a second cumulative histogram of said shift vectors determined in the preceding step for pixels in the other lung , selecting plural of said shift values for derivation of fitted shift values based on characteristics of said first and second cumulative histograms of shift vectors, and performing two-dimensional fitting on the shift values selected in the preceding step to derive said fitted shift values; and wherein said waφed image obtained in said step (d) is obtained using said fitted shift values.
13. The method of Claim 12, wherein: said step of selecting plural shift values comprises, selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said step of performing two-dimensional fitting comprises, performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said step (d) comprises, using said first set of fitted shift values to waφ said waφed image obtained in said step (b) to obtain a first twice-waφed image, using said second set of fitted shift values to waφ said waφed image obtained in said step (b) to obtain a second twice-waφed image, producing first and second subtraction images between said first digital image and said first twice-waφed image and between said first digital image and said second twice- waφed image, respectively, using said subtraction images to derive the iteratively waφed image for further processing in said step (e).
14. The method of Claim 13, wherein the step of using said subtraction images to derive the iteratively waφed image for further processing in said step (e) comprises: determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said step (e), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting for further two- dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two-dimensional fitting on the further selected fitted shift values to derive said second shift values, and waφing said waφed image obtained in step (b) using said further fitted shift values serving as said second shift values to produce said iteratively waφed image for further processing in said step (e).
15. The method of Claim 13 , wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90 ┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
16. The method of Claim 13 , wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
17. The method of Claim 13, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
18. The method of Claim 15, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
19. The method of Claim 16, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
20. The method of Claim 12, wherein said step of performing two-dimensional fitting on the shift values in said step (c) comprises: performing linear inteφolation of shift values.
21. A computer readable medium storing computer instructions for computerized analysis of temporally sequential digital images, by performing the steps of:
(a) determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image;
(b) waφing said second digital image based on the first shift values to obtain a waφed image in which spatial locations of pixels are varied in relation to said first shift values;
(c) determining second shift values between pixels of said first digital image and pixels of said waφed image obtained in step (b);
(d) waφing said waφed image obtained in step (b) based on the second shift values to obtain an iteratively waφed image in which spatial locations of pixels of said waφed image obtained in step (b) are varied in relation to said second shift values; and
(e) subtracting the iteratively waφed image from said first digital image.
22. The computer readable medium of Claim 21, further storing computer instructions for performing the steps of: determining third shift values between pixels of said first digital image and pixels of said iteratively waφed image; and waφing said iteratively waφed image based on the third shift values to obtain a further iteratively waφed image in which spatial locations of pixels of said iteratively waφed image are varied in relation to said third shift values.
23. The computer readable medium of Claim 21, wherein said step (a) comprises: selecting plural template regions of interest (ROIs) in the first digital image and corresponding search area ROIs in the second image, determining shift values between pixels centered in each template ROI and pixels centered in a respective search area ROI and exhibiting highest cross-correlation with respect to said template ROI, determining shift vectors for each of said shift values, a first cumulative histogram of said shift vectors for pixels in one lung, and a second cumulative histogram of said shift vectors for pixels in the other lung, selecting, based on characteristics of said first and second cumulative histograms of shift vectors, plural of said shift values for derivation of fitted shift values to serve as said first shift values, and performing two-dimensional fitting on the shift values selected in the preceding step to derive said fitted shift values which serve as said first shift values; and wherein said waφed image obtained in said step (b) is obtained using said fitted shift values.
24. The computer readable medium of Claim 23, wherein: said step of selecting plural shift values comprises, selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said step of performing two-dimensional fitting comprises, performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said step (b) comprises, using said first set of fitted shift values to waφ said second digital image to obtain a first waφed image, using said second set of fitted shift values to waφ said second digital image to obtain a second waφed image, producing first and second subtraction images between said first digital image and said first waφed image, and between said first digital image and said second waφed image, respectively, using said subtraction images to derive the waφed image for further processing in said steps (c) and (d).
25. The computer readable medium of Claim 24, wherein the step of using said subtraction images to derive the waφed image for further processing in said steps (c) and (d) comprises: determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said steps (c), and (d), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting for further two- dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two-dimensional fitting on the further selected fitted shift values, and waφing said second digital image using said further fitted shift values to produce said waφed image for further- processing in said steps (c) and (d).
26. The computer readable medium of Claim 24, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
27. The computer readable medium of Claim 24, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
28. The computer readable medium of Claim 24, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
29. The computer readable medium of Claim 26, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
30. The computer readable medium of Claim 27, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
31. The computer readable medium of Claim 23, wherein said step of performing two-dimensional fitting on the shift values comprises: using a two-dimensional nth order polynomial function.
32. The computer readable medium of Claim 23, wherein said step (c) comprises: selecting plural template regions of interest (ROIs) in the first image and corresponding search area ROIs in the waφed image obtained in step (b) , determining shift values between pixels centered in each template ROI in said first image and pixels centered in a respective search area ROIin said waφed image obtained in step (b) and exhibiting highest cross-correlation with respect to said template ROI, determining shift vectors for each of said shift values determined in the preceding step, determining a first cumulative histogram of said shift vectors determined in the preceding step for pixels in one lung, and a second cumulative histogram of said shift vectors determined in the preceding step for pixels in the other lung , selecting plural of said shift values for derivation of fitted shift values based on characteristics of said first and second cumulative histograms of shift vectors, and performing two-dimensional fitting on the shift values selected in the preceding step to derive said fitted shift values; and wherein said waφed image obtained in said step (d) is obtained using said fitted shift values.
33. The computer readable medium of Claim 32, wherein: said step of selecting plural shift values comprises, selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said step of performing two-dimensional fitting comprises, performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said step (d) comprises, using said first set of fitted shift values to waφ said waφed image obtained in said step (b) to obtain a first twice-waφed image, using said second set of fitted shift values to waφ said waφed image obtained in said step (b) to obtain a second twice-waφed image, producing first and second subtraction images between said first digital image and said first twice-waφed image and between said first digital image and said second twice- waφed image, respectively, using said subtraction images to derive the iteratively waφed image for further processing in said step (e).
34. The computer readable medium of Claim 33, wherein the step of using said subtraction images to derive the iteratively waφed image for further processing in said step (e) comprises: determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said step (e), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, selecting for further two- dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two-dimensional fitting on the further selected fitted shift values to derive said second shift values, and waφing said waφed image obtained in step (b) using said further fitted shift values serving as said second shift values to produce said iteratively waφed image for further processing in said step (e).
35. The computer readable medium of Claim 33 , wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
36. The computer readable medium of Claim 33, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90 ┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
37. The computer readable medium of Claim 33, wherein said step of selecting plural shift values comprises, selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
38. The computer readable medium of Claim 35, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
39. The computer readable medium of Claim 36, wherein said step of selecting plural shift values comprises, determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
40. The computer readable medium of Claim 32, wherein said step of performing two-dimensional fitting on the shift values in said step (c) comprises: performing linear inteφolation of shift values.
41. A system for computerized analysis of temporally sequential digital images, comprising:
(a) means for determining first shift values between pixels of a first digital image and corresponding pixels of a second digital image;
(b) means for waφing said second digital image based on the first shift values to obtain a waφed image in which spatial locations of pixels are varied in relation to said first shift values;
(c) means for determining second shift values between pixels of said first digital image and pixels of said waφed image obtained in means (b);
(d) means for waφing said waφed image obtained in means (b) based on the second shift values to obtain an iteratively waφed image in which spatial locations of pixels of said waφed image obtained in means (b) are varied in relation to said second shift values; and
(e) means for subtracting the iteratively waφed image from said first digital image.
42. The system of Claim 41 , further comprising: means for determining third shift values between pixels of said first digital image and pixels of said iteratively waφed image; and means for waφing said iteratively waφed image based on the third shift values to obtain a further iteratively waφed image in which spatial locations of pixels of said iteratively waφed image are varied in relation to said third shift values.
43. The system of Claim 41, wherein said means (a) comprises: means for selecting plural template regions of interest (ROIs) in the first digital image and corresponding search area ROIs in the second image, means for determining shift values between pixels centered in each template ROI and pixels centered in a respective search area ROI and exhibiting highest cross-correlation with respect to said template ROI, means for determining shift vectors for each of said shift values, a first cumulative histogram of said shift vectors for pixels in one lung, and a second cumulative histogram of said shift vectors for pixels in the other lung, means for selecting, based on characteristics of said first and second cumulative histograms of shift vectors, plural of said shift values for derivation of fitted shift values to serve as said first shift values, and means for performing two-dimensional fitting on the shift values selected in the preceding step to derive said fitted shift values which serve as said first shift values; and wherein said waφed image obtained by said means (b) is obtained using said fitted shift values.
44. The system of Claim 43, wherein: said means for selecting plural shift values comprises, means for selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said means for performing two-dimensional fitting comprises, means for performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and means for performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said means (b) comprises, means for using said first set of fitted shift values to waφ said second digital image to obtain a first waφed image, means for using said second set of fitted shift values to waφ said second digital image to obtain a second waφed image, means for producing first and second subtraction images between said first digital image and said first waφed image, and between said first digital image and said second waφed image, respectively, means for using said subtraction images to derive the waφed image for further processing in said means (c) and (d).
45. The system of Claim 44, wherein the means for using said subtraction images to derive the waφed image for further processing in said means (c) and (d) comprises: means for determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, means for selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said means (c), and means (d), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, means for selecting for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, performing further two-dimensional fitting on the further selected fitted shift values, and waφing said second digital image using said further fitted shift values to produce said waφed image for further processing in said means (c) and means (d).
46. The system of Claim 44, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
47. The system of Claim 44, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
48. The system of Claim 44, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
49. The system of Claim 46, wherein said means for selecting plural shift values comprises, means for determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
50. The system of Claim 47, wherein said means for selecting plural shift values comprises, means for determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
51. The system of Claim 53, wherein said means for performing two-dimensional fitting on the shift values comprises: means for using a two-dimensional nth order polynomial function.
52. The system of Claim 53, wherein said means (c) comprises: means for selecting plural template regions of interest (ROIs) in the first image and corresponding search area ROIs in the waφed image obtained in means (b) , means for determining shift values between pixels centered in each template ROI in said first image and pixels centered in a respective search area ROIin said waφed image obtained in means (b) and exhibiting highest cross-correlation with respect to said template ROI, means for determining shift vectors for each of said shift values determined in the preceding means, means for determining a first cumulative histogram of said shift vectors determined in the preceding means for pixels in one lung, and a second cumulative histogram of said shift vectors determined in the preceding means for pixels in the other lung , means for selecting plural of said shift values for derivation of fitted shift values based on characteristics of said first and second cumulative histograms of shift vectors, and means for performing two-dimensional fitting on the shift values selected in the preceding means to derive said fitted shift values; and wherein said waφed image obtained in said means (d) is obtained using said fitted shift values.
53. The system of Claim 52, wherein: said means for selecting plural shift values comprises, means for selecting, based on said first histogram, first and second sets of pixels in said one lung with shift vectors within respective first and second predetermined ranges of angles in said first histogram of shift vectors, and selecting, based on said second histogram, third and fourth sets of pixels in said other lung with shift vectors within respective third and fourth predetermined ranges of angles in said second histogram; said means for performing two-dimensional fitting comprises, means for performing two-dimensional fitting on the shift values of the first and third sets of pixels to derive a first set of fitted shift values, and means for performing two-dimensional fitting on the shift values of the second and fourth sets of pixels to derive a second set of fitted shift values; and said means (d) comprises, means for using said first set of fitted shift values to waφ said waφed image obtained in said means (b) to obtain a first twice-waφed image, using said second set of fitted shift values to waφ said waφed image obtained in said means (b) to obtain a second twice-waφed image, means for producing first and second subtraction images between said first digital image and said first twice-waφed image and between said first digital image and said second twice-waφed image, respectively, means for using said subtraction images to derive the iteratively waφed image for further processing in said means (e).
54. The system of Claim 53, wherein the means for using said subtraction images to derive the iteratively waφed image for further processing in said means (e) comprises: means for determining a histogram of pixel values in lung regions of said first and second subtraction images and determining which of said first and second subtraction images exhibits a narrower histogram of pixel values in each lung region, when one of said subtraction images exhibits narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, means for selecting the waφed image from which the subtraction image having said narrower histograms was produced for further processing in said means (e), and when one of said subtraction images does not exhibit narrower histograms of pixel values in both right and left lung regions in comparison to a histograms of pixel values in the respective right and left lung regions of the other subtraction image, means for selecting for further two-dimensional fitting the fitted shift values of pixels of each respective lung region in the respective first and second waφed images exhibiting the narrower histogram of pixel values, means for performing further two-dimensional fitting on the further selected fitted shift values to derive said second shift values, and means for waφing said waφed image obtained in means (b) using said further fitted shift values serving as said second shift values to produce said iteratively waφed image for further processing in said means (e).
55. The system of Claim 53, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
56. The system of Claim 53, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a ┬▒ 90┬░ range of angles of a peak in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
57. The system of Claim 53, wherein said means for selecting plural shift values comprises, means for selecting, based on said first histogram, a first set of pixels in said one lung with shift vectors within a 0-180┬░ range of angles in said first histogram of shift vectors and selecting said second set of pixels as the remaining pixels in said one lung, and means for selecting, based on said second histogram, a third set of pixels in said other lung with shift vectors within a 0-180┬░ range of angles in said second histogram of shift vectors and selecting said fourth set of pixels as the remaining pixels in said other lung.
58. The system of Claim 55, wherein said means for selecting plural shift values comprises, means for determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
59. The system of Claim 55, wherein said means for selecting plural shift values comprises, means for determining a peak in said first histogram of shift vectors as a shift vector angle having an accumulated shift value greater than twice the average of accumulated shift values of said first histogram.
60. The system of Claim 42, wherein said means for performing two-dimensional fitting on the shift values of said means (c) comprises: means for performing linear inteφolation of shift values.
PCT/US1999/004290 1998-04-02 1999-04-02 Iterative warping for temporal subtraction of radiographs WO1999054705A2 (en)

Priority Applications (4)

Application Number Priority Date Filing Date Title
EP99917300A EP1295244A4 (en) 1998-04-02 1999-04-02 Method, system and computer readable medium for iterative image warping prior to temporal subtraction of chest radiographs in the detection of interval changes
AU35452/99A AU3545299A (en) 1998-04-02 1999-04-02 Method, system and computer readable medium for iterative image warping prior totemporal subtraction of chest radiographs in the detection of interval changes
CA002326776A CA2326776C (en) 1998-04-02 1999-04-02 Method, system and computer readable medium for iterative image warping prior to temporal subtraction of chest radiographs in the detection of interval changes
JP2000545001A JP2003530722A (en) 1998-04-02 1999-04-02 Method, system and computer readable medium for iterative image distortion prior to temporal subtraction of chest radiographs in detecting temporal changes

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US09/053,798 US6067373A (en) 1998-04-02 1998-04-02 Method, system and computer readable medium for iterative image warping prior to temporal subtraction of chest radiographs in the detection of interval changes
US09/053,798 1998-04-02

Publications (2)

Publication Number Publication Date
WO1999054705A2 true WO1999054705A2 (en) 1999-10-28
WO1999054705A3 WO1999054705A3 (en) 2002-12-27

Family

ID=21986613

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US1999/004290 WO1999054705A2 (en) 1998-04-02 1999-04-02 Iterative warping for temporal subtraction of radiographs

Country Status (6)

Country Link
US (1) US6067373A (en)
EP (1) EP1295244A4 (en)
JP (1) JP2003530722A (en)
AU (1) AU3545299A (en)
CA (1) CA2326776C (en)
WO (1) WO1999054705A2 (en)

Families Citing this family (69)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8788020B2 (en) 1998-10-23 2014-07-22 Varian Medical Systems, Inc. Method and system for radiation application
US6937696B1 (en) 1998-10-23 2005-08-30 Varian Medical Systems Technologies, Inc. Method and system for predictive physiological gating
US7158610B2 (en) * 2003-09-05 2007-01-02 Varian Medical Systems Technologies, Inc. Systems and methods for processing x-ray images
US7043066B1 (en) 1998-11-05 2006-05-09 Arch Development Corporation System for computerized processing of chest radiographic images
US6904163B1 (en) * 1999-03-19 2005-06-07 Nippon Telegraph And Telephone Corporation Tomographic image reading method, automatic alignment method, apparatus and computer readable medium
US6941323B1 (en) * 1999-08-09 2005-09-06 Almen Laboratories, Inc. System and method for image comparison and retrieval by enhancing, defining, and parameterizing objects in images
US6594378B1 (en) * 1999-10-21 2003-07-15 Arch Development Corporation Method, system and computer readable medium for computerized processing of contralateral and temporal subtraction images using elastic matching
US6775400B1 (en) * 1999-10-29 2004-08-10 Acuson Corporation Medical diagnostic ultrasonic imaging method and apparatus for suppressing electronic noise
US6940998B2 (en) 2000-02-04 2005-09-06 Cernium, Inc. System for automated screening of security cameras
US6690840B1 (en) * 2000-02-08 2004-02-10 Tektronix, Inc. Image alignment with global translation and linear stretch
US7310437B2 (en) * 2000-03-08 2007-12-18 Fujifilm Corporation Image processing method and system, and storage medium
JP2002324238A (en) * 2001-04-26 2002-11-08 Fuji Photo Film Co Ltd Method and device for positioning image
US7769430B2 (en) * 2001-06-26 2010-08-03 Varian Medical Systems, Inc. Patient visual instruction techniques for synchronizing breathing with a medical procedure
US6795521B2 (en) * 2001-08-17 2004-09-21 Deus Technologies Llc Computer-aided diagnosis system for thoracic computer tomography images
GB0125774D0 (en) * 2001-10-26 2001-12-19 Cableform Ltd Method and apparatus for image matching
US6678399B2 (en) * 2001-11-23 2004-01-13 University Of Chicago Subtraction technique for computerized detection of small lung nodules in computer tomography images
US7127093B2 (en) * 2002-09-17 2006-10-24 Siemens Corporate Research, Inc. Integrated image registration for cardiac magnetic resonance perfusion data
US7260254B2 (en) * 2002-11-25 2007-08-21 Mirada Solutions Limited Comparing images
US20040122787A1 (en) * 2002-12-18 2004-06-24 Avinash Gopal B. Enhanced computer-assisted medical data processing system and method
US20040122709A1 (en) * 2002-12-18 2004-06-24 Avinash Gopal B. Medical procedure prioritization system and method utilizing integrated knowledge base
US20040122719A1 (en) * 2002-12-18 2004-06-24 Sabol John M. Medical resource processing system and method utilizing multiple resource type data
US7187790B2 (en) * 2002-12-18 2007-03-06 Ge Medical Systems Global Technology Company, Llc Data processing and feedback method and system
US20040122708A1 (en) * 2002-12-18 2004-06-24 Avinash Gopal B. Medical data analysis method and apparatus incorporating in vitro test data
US20040122706A1 (en) * 2002-12-18 2004-06-24 Walker Matthew J. Patient data acquisition system and method
US20040122702A1 (en) * 2002-12-18 2004-06-24 Sabol John M. Medical data processing system and method
US20040122707A1 (en) * 2002-12-18 2004-06-24 Sabol John M. Patient-driven medical data processing system and method
US7490085B2 (en) * 2002-12-18 2009-02-10 Ge Medical Systems Global Technology Company, Llc Computer-assisted data processing system and method incorporating automated learning
US20040122705A1 (en) * 2002-12-18 2004-06-24 Sabol John M. Multilevel integrated medical knowledge base system and method
US20040122704A1 (en) * 2002-12-18 2004-06-24 Sabol John M. Integrated medical knowledge base interface system and method
US20040122703A1 (en) * 2002-12-19 2004-06-24 Walker Matthew J. Medical data operating model development system and method
US7373275B2 (en) * 2003-03-14 2008-05-13 Rudolph Technologies, Inc. System and method of point matching measured positions to template positions
US8571639B2 (en) * 2003-09-05 2013-10-29 Varian Medical Systems, Inc. Systems and methods for gating medical procedures
US20050053267A1 (en) * 2003-09-05 2005-03-10 Varian Medical Systems Technologies, Inc. Systems and methods for tracking moving targets and monitoring object positions
JP4439882B2 (en) * 2003-11-14 2010-03-24 キヤノン株式会社 Radiation image processing apparatus and processing method
JP2005149181A (en) * 2003-11-17 2005-06-09 Konica Minolta Medical & Graphic Inc Medical image information management system
US20050105789A1 (en) * 2003-11-17 2005-05-19 Isaacs Hugh S. Method and apparatus for detecting, monitoring, and quantifying changes in a visual image over time
JP2005312776A (en) * 2004-04-30 2005-11-10 Canon Inc Radiograph acquisition device
US8111947B2 (en) * 2004-06-08 2012-02-07 Canon Kabushiki Kaisha Image processing apparatus and method which match two images based on a shift vector
WO2006093523A2 (en) * 2004-07-15 2006-09-08 Kenji Suzuki Computerized scheme for distinction between benign and malignant nodules in thoracic low-dose ct
US20060136259A1 (en) * 2004-12-17 2006-06-22 General Electric Company Multi-dimensional analysis of medical data
US20060136417A1 (en) * 2004-12-17 2006-06-22 General Electric Company Method and system for search, analysis and display of structured data
ATE438161T1 (en) * 2005-01-19 2009-08-15 Koninkl Philips Electronics Nv IMAGE PROCESSING SYSTEM AND METHOD FOR ARRANGEMENT OF IMAGES
US7623731B2 (en) * 2005-06-20 2009-11-24 Honda Motor Co., Ltd. Direct method for modeling non-rigid motion with thin plate spline transformation
US7822224B2 (en) 2005-06-22 2010-10-26 Cernium Corporation Terrain map summary elements
US7889896B2 (en) * 2005-08-18 2011-02-15 Hologic, Inc. Patient worklist management in digital radiography review workstations
US20070078873A1 (en) * 2005-09-30 2007-04-05 Avinash Gopal B Computer assisted domain specific entity mapping method and system
US7769216B2 (en) * 2005-12-29 2010-08-03 Hologic, Inc. Facilitating comparison of medical images
DE102006029718A1 (en) * 2006-06-28 2008-01-10 Siemens Ag Organ system`s e.g. brain, images evaluating method for detecting pathological change in medical clinical picture, involves extending registration for area to extended registration, such that another area is detected
US7697899B2 (en) * 2006-08-31 2010-04-13 Broadcom Corporation RFIC with on-chip acoustic transducer circuit
WO2008141293A2 (en) * 2007-05-11 2008-11-20 The Board Of Regents Of The University Of Oklahoma One Partner's Place Image segmentation system and method
US20090082637A1 (en) * 2007-09-21 2009-03-26 Michael Galperin Multi-modality fusion classifier with integrated non-imaging factors
CN101903885A (en) * 2007-12-18 2010-12-01 皇家飞利浦电子股份有限公司 Consistency metric based image registration
JP5361194B2 (en) * 2008-01-09 2013-12-04 キヤノン株式会社 Image processing apparatus, image processing method, and computer program
US20090196479A1 (en) * 2008-01-31 2009-08-06 Raghav Raman Method and apparatus for computer-aided diagnosis filtered prioritized work item list
JP5016518B2 (en) * 2008-02-21 2012-09-05 富士フイルム株式会社 Alignment apparatus and program thereof
US10667727B2 (en) * 2008-09-05 2020-06-02 Varian Medical Systems, Inc. Systems and methods for determining a state of a patient
US20100061596A1 (en) * 2008-09-05 2010-03-11 Varian Medical Systems Technologies, Inc. Video-Based Breathing Monitoring Without Fiducial Tracking
US20100098305A1 (en) * 2008-10-16 2010-04-22 Riverain Medical Group Llc Multi-image correlation
JP5523791B2 (en) * 2008-10-27 2014-06-18 株式会社東芝 X-ray diagnostic apparatus and image processing apparatus
US8571261B2 (en) * 2009-04-22 2013-10-29 Checkvideo Llc System and method for motion detection in a surveillance video
JP5586917B2 (en) * 2009-10-27 2014-09-10 キヤノン株式会社 Information processing apparatus, information processing method, and program
WO2011135484A1 (en) 2010-04-26 2011-11-03 Koninklijke Philips Electronics N.V. Device and method for providing subtraction images
WO2012106580A2 (en) 2011-02-04 2012-08-09 The Regents Of The University Of California Methods and apparatus for computer-aided radiological detection and imaging
US10540783B2 (en) * 2013-11-01 2020-01-21 Illumina, Inc. Image analysis useful for patterned objects
DE102015200850B4 (en) * 2015-01-20 2019-08-22 Siemens Healthcare Gmbh Method for evaluating medical image data
CN106611411B (en) * 2015-10-19 2020-06-26 上海联影医疗科技有限公司 Method for segmenting ribs in medical image and medical image processing device
US11151726B2 (en) * 2018-01-10 2021-10-19 Canon Medical Systems Corporation Medical image processing apparatus, X-ray diagnostic apparatus, and medical image processing method
US11049606B2 (en) 2018-04-25 2021-06-29 Sota Precision Optics, Inc. Dental imaging system utilizing artificial intelligence
JP2020174861A (en) * 2019-04-17 2020-10-29 キヤノン株式会社 Information processing device, information processing method, and program

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5359513A (en) * 1992-11-25 1994-10-25 Arch Development Corporation Method and system for detection of interval change in temporally sequential chest images
US5623560A (en) * 1992-11-27 1997-04-22 Fuji Photo Film Co., Ltd. Method for adjusting positions of radiation images
US5629988A (en) * 1993-06-04 1997-05-13 David Sarnoff Research Center, Inc. System and method for electronic image stabilization

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5319549A (en) * 1992-11-25 1994-06-07 Arch Development Corporation Method and system for determining geometric pattern features of interstitial infiltrates in chest images
JPH1083442A (en) * 1996-09-09 1998-03-31 Canon Inc Picture synthesis method
US6714689B1 (en) * 1995-09-29 2004-03-30 Canon Kabushiki Kaisha Image synthesizing method

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5359513A (en) * 1992-11-25 1994-10-25 Arch Development Corporation Method and system for detection of interval change in temporally sequential chest images
US5623560A (en) * 1992-11-27 1997-04-22 Fuji Photo Film Co., Ltd. Method for adjusting positions of radiation images
US5629988A (en) * 1993-06-04 1997-05-13 David Sarnoff Research Center, Inc. System and method for electronic image stabilization

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
MOSHFEGHI M., Elastic Matching of Multimodality Medical Images. CVGIP. May 1991. Vol. 53, No. 3, pages 271-282, XP002922620. *
See also references of EP1295244A2 *

Also Published As

Publication number Publication date
EP1295244A2 (en) 2003-03-26
JP2003530722A (en) 2003-10-14
CA2326776C (en) 2008-06-17
CA2326776A1 (en) 1999-10-28
AU3545299A (en) 1999-11-08
WO1999054705A3 (en) 2002-12-27
EP1295244A4 (en) 2009-06-03
US6067373A (en) 2000-05-23

Similar Documents

Publication Publication Date Title
US6067373A (en) Method, system and computer readable medium for iterative image warping prior to temporal subtraction of chest radiographs in the detection of interval changes
US5982915A (en) Method of detecting interval changes in chest radiographs utilizing temporal subtraction combined with automated initial matching of blurred low resolution images
Ishida et al. Iterative image warping technique for temporal subtraction of sequential chest radiographs to detect interval change
US6724925B2 (en) Method and system for the automated delineation of lung regions and costophrenic angles in chest radiographs
US5359513A (en) Method and system for detection of interval change in temporally sequential chest images
US6594378B1 (en) Method, system and computer readable medium for computerized processing of contralateral and temporal subtraction images using elastic matching
Karssemeijer Adaptive noise equalization and recognition of microcalcification clusters in mammograms
van Engeland et al. A comparison of methods for mammogram registration
US6018590A (en) Technique for finding the histogram region of interest based on landmark detection for improved tonescale reproduction of digital radiographic images
US9659390B2 (en) Tomosynthesis reconstruction with rib suppression
WO2002043562A2 (en) Pulmonary nodule detection using multiple-templates matching
US7043066B1 (en) System for computerized processing of chest radiographic images
EP2631873B1 (en) Image alignment of breast images
Yoshida Local contralateral subtraction based on bilateral symmetry of lung for reduction of false positives in computerized detection of pulmonary nodules
US8229190B2 (en) Visual enhancement of interval changes using temporal subtraction and pattern detector
US20080193000A1 (en) Visual Enhancement of Interval Changes Using Temporal Subtraction
Loeckx et al. Non-rigid image registration using a statistical spline deformation model
Good et al. Generalized procrustean image deformation for subtraction of mammograms
Correa et al. Digital chest radiography: comparison of unprocessed and processed images in the detection of solitary pulmonary nodules.
JP4638042B2 (en) Computer processing system for chest X-ray images
Loeckx et al. Temporal subtraction of thorax CR images
Seghers et al. Visual enhancement of interval changes using a temporal subtraction technique
Gauch et al. Segmentation via the Intensity Axis of Symmetry

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A2

Designated state(s): AU CA JP

AL Designated countries for regional patents

Kind code of ref document: A2

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

121 Ep: the epo has been informed by wipo that ep was designated in this application
DFPE Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)
ENP Entry into the national phase

Ref document number: 2326776

Country of ref document: CA

Ref country code: CA

Ref document number: 2326776

Kind code of ref document: A

Format of ref document f/p: F

ENP Entry into the national phase

Ref country code: JP

Ref document number: 2000 545001

Kind code of ref document: A

Format of ref document f/p: F

WWE Wipo information: entry into national phase

Ref document number: 1999917300

Country of ref document: EP

AK Designated states

Kind code of ref document: A3

Designated state(s): AU CA JP

AL Designated countries for regional patents

Kind code of ref document: A3

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

WWP Wipo information: published in national office

Ref document number: 1999917300

Country of ref document: EP