CA2577547A1 - Method and system for discriminating image representations of classes of objects - Google Patents

Abstract

A method of identifying the location of a lesion in an image and evaluating whether the identified lesion is more likely to be cancerous, benign or uncertain is provided where the image includes a plurality of pixels, each pixel having a particular intensity I in the range of 0 <I <2N, where N is an integer > 1. The method includes defining one or more landmark pixels within the lesion; in the, growing a cluster around the landmark pixel(s) which contains the lesion for each intensity value in a set of possible intensity values; and at each of the intensity levels in the set, constructing a region of interest such that the region of interest is a minimal polygon containing the cluster at that intensity level. The method further includes computing a value of a characteristic of the minimal polygon at each intensity level, determining a number that is characteristic of changes in the area or volume of the regions of interest changes over the range of intensity level in the set, and determining whether the lesion is more likely to be cancerous, benign or uncertain in response to the number that is characteristic to the changes in the area or volume of the regions of interest changes, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

Description

TITLE OF THE INVENTION

METHOD AND SYSTEM FOR DISCRIMINATING IMAGE REPRESENTATIONS OF
CLASSES OF OBJECTS

CROSS-REFERENCE TO RELATED APPLICATIONS
[0001] This application claims the benefit of U.S.
Provisional Applications Nos. 60/601,981, filed on August 17, 2004, and 60/647,756, filed on January 31, 2005.

[0002] This work was supported in part by at least one grant R44CA85101 issued by the National Cancer Institute (NCI). The government may have certain rights in this invention.

COMPUTER PROGRAM LISTING APPENDIX SUBMITTED ON A COMPACT DISC
Submitted herewith are two identical, or duplicate, compact discs. The material on these compact discs is incorporated herein by reference. Each of these identical compact discs contains the following files:

File Date Size Description Created LevelSetMapper.java Jan 14, 10195 Java source code. Main entry point 2005 for algorithm. Sets parameters for algorithm, loads image data into memory, and calls methods in LevelSetFinder.
LevelSetFinder.java Jan 4, 4844 Java source code. Initializes 2005 components for the computation of intensity steps. Calls methods in LevelSetHelper.
LevelSetHelper.java Jan 4, 27314 Java source code. Contains main 2005 algorithm and methodology details.
Calls methods in ClusterMapS and BasicSegmenter for connected component analysis and segmentation of cluster.
ClusterMapS.java Jan 4, 43293 Java source code. Contains 2005 supplemental algorithms for connect component analysis.
BasicSegmenter.java Oct 9, 2442 Java source code. Contains 2003 supplemental algorithms for segmenting the image.
Segmentat.ionBuffer.java Oct 9, 21313 Java source code. Work-horse for 2003 BasicSegmentation.
ClusterDispersion.java Apr 19, 1473 Java source code. Computes an 2005 adjustment to intensity step values based on cluster dispersion.
ClusterCondition.java Mar 23, 1207 Java source code. Describes the 2005 cluster attributes and conditions upon which a cluster is removed from consideration.
Lsvuparams.xml Mar 23, 815 Input parameters that control various 2005 aspects of the algorithm.

BACKGROUND OF THE INVENTION

Field of Invention

[0003] The present invention relates to a method, procedure, and set of computer codes for discriminating image representations of classes of objects, such as MRI signals of benign and malignant lesions, where the images have been transformed to meet constraints of the system, e.g., windowing/leveling to display a lesion in the center of a range of displayed intensities.

Prior art

[0004] In many prior applications, data are represented in three dimensions where two of the dimensions (x,y) represent spatial location on a grid of fixed size and the third dimension (w) is a representation of original source data. As an example, in magnetic resonance imaging (MRI) systems, the w dimension is image intensity, which represents signal strength from the imaging system. In the MRI example, the source signal is transformed to an image intensity signal to show a medically appropriate representation of the MRI

image data in a way that satisfies limitations in intensity values imposed by the display system, e.g. 28 =256 discrete intensity levels. Linear transformations on the w-axis, such as windowing and leveling, are typically used to show the medically relevant portions of the image in the central portion of the displayed intensity scale. In many of these applications, the end-user could benefit from being able to distinguish and discriminate objects (e.g., lesions) within the source data on the basis of the original signal strength.
However, the transformations from source signal to image intensity signal, which may vary from case to case, make this comparative analysis difficult and subject to error. In other applications data are represented in 4 dimensions, where 3 of the dimensions represent spatial location (x,y,z) and the fourth dimension (w) is a representation of the source data.
All discussions and descriptions for the invention in 2 spatial dimensions are readily extended to 3 spatial dimensions. While the term pixel is frequently used to refer to 2-dimensions and the term voxel is frequently used to refer to 3-dimensions, in this application we use pixel to refer to 2-dimensions and 3-dimensions.

[0005] In the example of medical MRI, the source signals from gadolinium-enhanced images of malignant lesions are frequently stronger than the source signals from gadolinium enhanced images of benign lesions. However, after the source data have been transformed to image intensities that have been adjusted to optimize medical diagnosis, where this optimization differs from case to case, it is difficult for the radiologist to evaluate the strength of the original magnetic resonance signal on the basis of the images that are presented.

[0006] An object of the present invention is to provide a system, method and computer program product for evaluating whether a lesion in an image is cancerous, benign, or of an uncertain nature based on the intensities of the pixels in the image.

[0007] It is a further object of the present invention to provide a system, method and computer program product for evaluating a lesion which can discriminate classes of objects based on intensity and spatial relationships of pixels, e.g., gradients and shapes within or near the objects.

[0008] It is a further object of the present invention to provide a system, method and computer program product to allow a radiologist to evaluate the strength of the original image signal after the source data have been transformed to image intensities that have been adjusted to optimize medical diagnosis.

[0009] A further object is to simultaneously obtain diagnostic information from a substantial portion of a patient's body, possibly presenting multiple lesions.

[0010] Another object is to obtain highly reliable diagnostic information based on image data obtained from images corresponding to one, or at most two, time frames.

[0011] The scope and content of the present invention is not intended to be limited by or to the above mentioned objects.

[0012] The invention provides a novel method of evaluating images of a body region, wherein at least one of the images shows an abnormality indicative of a possible lesion, comprising the steps of:
determining locations of pixels in each image that show such abnormality;
for each of a set of intensity levels, I, determining a contour around the cluster containing the pixels at the locations determined in the step of determining;
defining a function, F, that discriminates a distinct characteristic of each contour in a nested sequence;
defining a function, G, used to characterize the changes in the value of the function F over the range of contours at each intensity level; and identifying a lesion as being more likely to be benign or more likely to be cancerous based on at least one threshold value for the function G or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

[0013] Briefly, the invention is characterized in particular by the use of a method of evaluating whether a lesion in an image is cancerous, benign or uncertain, where the image includes a plurality of pixels, each pixel having a particular intensity I in the range of 0<_ I< 2N, where N is an integer > 1. One embodiment of the method includes the steps of defining a landmark pixel within the lesion in the image, growing a cluster around the landmark pixel which contains the lesion for each intensity value in a set of possible intensity values; and at each intensity level in the set, constructing a region of interest such that the region of interest is a minimal polygon containing the cluster at that intensity level. The method further includes the steps of computing a value of a characteristic of the minimal polygon at each intensity level in the set, determining a number related to changes in the characteristic values over the range of intensity levels, and determining whether the lesion is more likely to be cancerous, benign or uncertain in response to the number related to changes in the characteristic values over the range of intensity levels.

[0014] An evaluation process according to the invention may typically start with a number of images of different body region slices. A preliminary manual or automated operation may be carried out to exclude from further processing those images that do not show evidence of a lesion.

[0015] For example, all images could be examined visually and only those images that appear to show an abnormality or abnormality could then be evaluated according to the invention.

[0016] Alternatively, all images could be evaluated according to the invention and those images for which the evaluation indicates that no lesion is present could be discarded. For example, images could be discarded if there is no evidence of an abnormality, or if evaluations of each area of interest, or abnormality, in an image generate a number of changes in the size of the image contour that is below a certain threshold value. That value could be selected on the basis of experience, but would be expected to be in the range of 3 to 5.

BRIEF DESCRIPTION OF THE DRAWINGS

[0017] For better understanding of the invention and additional objects and advantages thereof, reference is made to the following detailed description and accompanying drawing of a preferred embodiment, wherein

[0018] Fig. 1 is a diagram of the imaging device or picture archiving system and computer system according to one embodiment of the present invention.

[0019] Fig. 2 is a flowchart of the method to initialize and compute IMax in accordance with one embodiment of the present invention;

[0020] Fig. 3 is a flowchart of the method to compute IMin and number steps, in accordance with one embodiment of the present invention;

[0021] Fig. 4 illustrates the clusters and bounding boxes for a malignant lesion in accordance with one embodiment of the present invention;

[0022] Fig. 5 illustrates the clusters and bounding boxes for a benign lesion in accordance with one embodiment of the present invention;

[0023] Fig. 6 is a graph showing that the invention was better at discriminating benign from malignant breast lesions on MRI on the basis of image intensity than 34 of 35 radiologists tested in a multi-center reader study;

[0024] Fig. 7 is a block diagram illustrating a computer program product in accordance with one embodiment of the present invention; and

[0025] Fig. 8 illustrates one example of a method of using histogram analysis to determine Imax and Imin.

[0026] Fig. 9 is a replica of a presentation illustrating various features and advantages of the invention.

[0027] Each of Figs. 10-13 is a reproduction of a respective one of the panels of Fig. 9.

[0028] Fig. 14 is a screen capture illustrating a diagnostic display produced according to the invention.
DETAILED DESCRIPTION OF THE INVENTION

[0029] ' The current invention is a method, procedure, and set of computer codes for discriminating image representations of classes of objects, such as MRI signals of benign and malignant lesions, where the images have been transformed to meet constraints of the system, e.g., windowing/leveling, to display a lesion in the center of a range of displayed intensities. The methods, procedures, and computer codes of the invention discriminate classes of objects on the basis of intensity and spatial relationships of pixels, e.g., gradients and shapes, within or near the objects.

[0030] The invention will be described in terms of its application in medical -radiology for discriminating benign from malignant lesions on gadolinium-enhanced magnetic resonance images (MRI) on the basis of image intensity values where the image data being analyzed has 256 discrete intensity values and has been subjected to prior windowing and leveling operations according to known techniques. The invention is first described in terms of images corresponding to a 2-dimensional spatial slice. The extension of the invention to a set of 2-dimensional slices that comprise a 3-dimensional data set is described later. It is assumed that windowing and leveling is "reasonably consistent" between cases that are to be discriminated, conforming to standard medical practice, and for each case, a "landmark" showing the approximate location of the lesion, is known. It is understood that the invention may be applied to any imaging system in which the goal is to evaluate the image intensity and spatial relationships of the pixels in the image, within the skill of the ordinary artisan.

[0031] Using standard thresholding and clustering algorithms, a cluster is grown around the landmark for each possible intensity value, which according to one embodiment starts with the highest (e.g., 255) and ending with the lowest (0). The clusters around the landmark form a nested, monotonically increasing (but not necessarily strictly increasing) sequence. At each possible intensity level, a region-of-interest (ROI) is constructed around the cluster in a particular shape such that the ROI is the minimal shape containing the cluster. According to one embodiment, the ROI
is a minimal rectangular box, or rectangular hull formed around the cluster. Other shapes may be used within the skill of the ordinary artisan. The ROIs also form a nested, monotonically increasing (but not necessarily strictly increasing) sequence. According to one embodiment of the present invention, where the ROI is a rectangular box, for each ROI in the sequence, we compute the area of the ROI by multiplying width by height. If the shape for the ROI is not a rectangular box, the area is computed using a different formula, depending on the ROI shape. If the characterization of the ROI being used is not the area, then a different formula may be used. As an example of a possible characterization other than area, in ultrasound, the ratio of width to height is important and this ratio can be used as the chosen characteristic. Further, if the ROI is depicted in 3-dimensions, instead of 2-dimensions, volume of the ROI may be used instead of area.

[0032] A plot of ROI area vs. intensity level is a step function -- the plot of ROI area vs. intensity may remain constant for several intensity levels and then "step" up to a larger size. The number of steps has been found to be highly predictive of whether the lesion is benign or malignant using images from a variety of MRI 'imaging systems and protocols.
Moreover, the number of steps has been found to show a high degree of independence from other discriminatory features and to be useful as a component of a computer-aided-diagnosis or computer-aided-detection system. In the specific example shown here, an image of a lesion is interpreted as being benign if the number of steps is less than or equal to 9 and is interpreted as being cancer if the number of steps is greater than 9. These thresholds may be adjusted as appropriate by an ordinarily skilled artisan. Additionally, other numbers related to the characterization of the ROI may be used.

[0033] While the number of distinct ROIs is a function of shape and gradient of a lesion, it is relatively insensitive to transformations of intensity values, such as windowing and leveling, provided that these transformations are not extreme (e.g., the leveling cannot have reduced the image to a few intensities).

[0034] One embodiment of the present invention can be alternatively described in a more general mathematical context as follows: A cluster is a set of connected pixels. A

contour at level L is constructed by first generating a binary threshold image where a pixel in the threshold image has value 1 if the corresponding pixel in the grey-scale image has value >= L and has value 0 otherwise. A contour at level L is the set of pixels at the boundary of 0's and l's on the binary image. The Outer Contour at level L is the contour at level L
that encloses the landmark and is furthest from the landmark.
The ROI at level L is a geometric object having a specified shape, such as a square or rectangular box, that is of minimal size around a cluster or contour.

[0035] 1. Determine location of pixels in lesion.
When we refer to "pixel" we mean the picture element at a specific location in the coordinate system of the image.

[0036] 2. A landmark within the lesion is selected, either manually or automatically within the lesion. Clusters around the landmark are determined for each level L in a subset of possible intensity levels as determined by a predefined set of rules, and Outer Contours are determined for the cluster at each of the L's. For example, each intensity level within the image may be used, or some subset thereof, e.g., every other or every third intensity level may be sampled and used. In a more general context, other sets of closed paths around the landmark could be defined using other methods that are known, within the skill of the ordinary artisan.

[0037] 3, Define a function, F, from the space of Outer Contours to the space of real numbers. In the specific method described above, for each L we determine the Outer Contour and define the function to be the area of the rectangle, F(C) =Area (B), where B is the ROI defined to be the minimal rectangle around the Outer Contour. In a more general context, the ROI B could be another polygonal shape around the cluster that forms a nested sequence over the range of intensity levels, and F could-be any function that discriminates distinct elements defining characteristics of the ROI in the nested sequence, within the skill of the ordinary artisan.

[0038] 4. Define a function, G, on the range of F
over the set of Outer Contours {C}. In the specific method described above, G({RangeFc})= M, where M is the number of distinct elements in the Range (i.e., the number of times F, the area, changes values). In a more general context, G could be another function used to characterize the function F of step 3, within the skill of the ordinary artisan. Further, it is possible to only consider steps in the Outer Contours in a portion of the range, to consider the density of steps, or other appropriate functions, as will be readily understood by those of ordinary skill in the art.

[0039] 5. Define a feature, i.e., whether the lesion is more likely cancerous, benign, or uncertain, based on the function G. In the specific method described above a single threshold is set at 9 to discriminate benign from malignant lesions. In the more general context, a different threshold could be used or multiple thresholds or another distinguishing characteristic of G could be used to indicate different likelihoods of being benign or malignant, within the skill of the ordinary artisan.

[0040] According to one embodiment, the present invention is implemented on a computer connected to an imaging device or Picture Archiving system, such as a MRI device or other suitable imaging device or hospital PAC system (see Fig.
1). For purposes of this disclosure, we will refer interchangeably to a computer to mean a computer which is separate from the imaging device, or one which is integrated in the imaging device, wherein communication between the user and the computer (i.e., input device and display) is through the imaging device console, such as an MRI console. According to this embodiment, the computer has an input device (e.g., keyboard and mouse), a display device (e.g., monitbr), and a processor. The processor can be a known system, including a storage device, a central processing unit (CPU), and other known components (not shown). The computer can be implemented separately, or as part of the imaging or archiving device. In the latter case, the display and input device of the imaging or archiving device could be used to interface with the computer, rather than separate components.

[0041] Source data consists of pixel intensities of an image derived from the MRI signal captured after use of contrast agent (POST) and pixel intensities of an image derived from the MRI signal captured before use of contrast agent (PRE). Pixel intensities of a subtraction (SUB) image are obtained by subtracting the pixel intensities of PRE from pixel intensities of the POST (Fig. 2, step 1). If there are multiple post contrast images, the post contrast image that corresponds to peak enhancement is used. We indicate below whether SUB or POST is used for each step in the current implementation.

[0042] According to one embodiment of the present invention, parameters are set to: Q=25 mm2, N=4, (Fig. 2, step 2), where Q is a lower bound on the area of the lesion, and N

is determined heuristically to approximate the point at which the cluster effectively grows into background noise. The meaning of the number N is explained as follows: A minimum size of the lesion, E, is obtained by first constructing the Outer Contours at each intensity level, L, starting with the intensity level of the landmark and decrementing, until a level is reached for which the area within the Outer Contour first exceeds Q, the lower bound set by parameter. As intensity level L is further decremented, the area within the Outer Contour increases, ultimately encompassing the entire image, including background tissue outside of the lesion. For each of these Outer Contours, the mean gradient along the corresponding path on the Post image is computed. The level IMax, which corresponds to the maximum mean gradient, is selected and the area within Outer Contour of level IMax is the minimum area of the lesion. As the index L is decremented beyond IMax, the area within the Outer Contours increases.
When the area first exceeds N times the minimum area of the lesion, the Outer Contour is assumed to have extended beyond the lesion and grown into the background tissue.

[0043] The "step feature" is a measurement of a grouping of enhancing pixels on the SUB image, as determined by a landmark, L, defined to be a single pixel in an enhancing group. (Fig. 2, step 3). In general, different landmarks within the same enhancing group will produce different step feature values. The landmark that is used can either be determined by an expert using image and contextual information or determined automatically from image processing and/or histogram analysis. In the implementation according to one embodiment, histogram analysis is used to identify pixels intensities that are likely to be part of a lesion, and cluster analysis is used to identify collections of enhancing pixels that are likely to comprise the lesion. The centroid or other identified region of the cluster of enhancing pixels can be used to identify the landmark. In the implementation according to another embodiment, a radiologist draws a rectangular, or other shaped, ROI around the lesion and the landmark is the center point of the ROI. This ROI is input to the processor by the input device.

[0044] We now describe the step feature algorithmically, and assume for this discussion that there are 256 possible pixel intensity levels on the images, ranging from 255 (highest) to 0(lowest) . Let I(L) denote the pixel intensity at the landmark, each pixel having a particular intensity I in the range of 05 I< 255. According to another embodiment of the invention, each pixel may have a particular intensity I in the range of 0 _ I< 2N1 where N is an integer >
1, which would include image displays with 128, 512 or 1024 intensity levels. Starting with level I=I(L) and decrementing I at each step, we construct the cluster of pixels that are 4-connected to L and have intensity level>=I. A cluster is 4-connected if there is a 4-neighbor path from every pixel in the cluster to every other pixel in the cluster where pixels are 4-neighbor if they are side-by-side or one directly above the other. Other forms of connectivity, including, but not limited to, 8-neighbor in 2-dim and 6-neighbor, 18-neighbor or 26-neighbor in 3-dim can also be used. (See Digital Image Processing, Gonzalez and Waitz, 2nd Edition, Adison & Wesley, 1987.) These clusters form a monotonically increasing set {CN, CN-1, ...} , with Function (CN) <=Function (CN_1) <=.... as the index is decremented. When in 2-dimensions, the Function is the Area of the cluster. When in 3-dimensions, the Function may be the Volume of the cluster. Other alternatives also can be used, within the skill of the ordinary artisan. This process is continued until intensity level equals II, where Function(CI=)>=Q, where the Function is Area when in 2-dim, and Volume when in 3-dim. II is the first level at which the Function of the Outer Contour exceeds the lower bound of the lesion as set by the parameter. (Fig. 2, steps 4-9). Step 5 computes the binary threshold image=used to derive the Outer Contour and Step 6 computes the Function (such as area or volume) within the Outer Contour.

[0045] An Imax and an Imin value can be determined using a histogram analysis. Alternatively, according to one embodiment, starting with intensity level J=II and decrementing by J at each step, the mean gradient on the border of the 4-connected set CJ (MeanGrad(J)) is computed using data from POST. (Fig. 2, step 10). The intensity level at which MeanGrad is maximum defines level I,,, (Fig. 2, steps 11-14). For each pixel on the Outer Contour, the gradient of the corresponding pixel in the Post image is computed using standard image processing techniques. The MeanGrad is defined as the mean of this set of computed gradient values. One example of a method of using histogram analysis to determine Imax and Imin is illustrated in Figure 8. A 64x64 pixel subimage containing a lesion was cropped from a breast MRI
post-contrast image. The graph in Figure 8 shows the histogram of pixel intensities within the cropped subimage, after smoothing. Each pixel in the MRI image covers an area approximately equal to .4x.4=.16 mm2. For each intensity level, the approximate area of pixels having that intensity level is computed by multiplying the number of pixels having that intensity level by .16 mm2 . For each intensity level, the approximate area of pixels having that intensity level or greater is computed by summing the areas for all intensity levels greater than or equal to the given intensity level.
Intensity level 196, shown by the vertical bar furthest to the right, is the first intensity level such that the area of pixels greater than or equal to that level exceeds an area of 25 mm2, corresponding to the parameter Q in the embodiment given above. Intensity level 183, which is used as Imax, shown by the middle vertical bar, is the intensity level at which the histogram reaches its maximum in the right peak of the histogram. The area of pixels having values greater than or equal to Imax is computed as described above. Intensity level 74, which is used as Imin, shown by the left vertical bar, is the highest intensity level such that the area of pixels greater than or equal to that level exceeds the area computed from Imax by a factor of 4, corresponding to the parameter N in the embodiment given above.

[0046] Imin is set as the lowest intensity level for which the Function of Clmin exceeds the Function of Clmax by some pre-determined multiple, i.e., Function (Clmin) > N* Function (CzmaX) - (Fig. 3, Steps 1-7.) Alternative criteria for establishing Imin can be determined from context, cluster size, cluster numbers, or histogram analysis, within the skill of the ordinary artisan.

[0047] Starting with level I=Imax and decrementing through I=Imini the minimum bounding ROI BI around CI is constructed and the Functions representing the characteristics ROIs are computed: B=ma, C B=max-iC ..., with Function (BImax) <- Function (Bimax-1) 5 .... Depending upon the changes in Outer Contours from one intensity level to the next lower intensity level, the minimum bounding ROIs may increase or remain constant. Each time that a decrement in intensity level induces a change in the Function of the minimum bounding ROI, a step counter is incremented. The "step feature" is the final value of the step counter which is output as the total number of steps when B(J) > B (old) where B(old) is the previous minimum bounding ROI. A determination is then made as to whether the lesion is more likely to be benign, cancerous or of an uncertain nature, based on the total number of steps. (Fig. 3, steps 8-12.) It is also contemplated that another number related to the changes in the characteristics of ROI can be used instead of the total number of steps.

[0048] Figures 4 and 5 show the contours and bounding ROIs, in this case, rectangles for a malignant and benign lesion, respectively. Box 1 shows the cluster at intensity level Ima,. Growth of the cluster (pixels that have been added from the previous intensity level) is shown in black.

Bounding boxes that have increased from the previous level are shown with a solid border; bounding boxes that have not increased are shown with a dotted border.

[0049] The two lesions have comparable sizes, shapes, and intensity ranges on the MRI images. However, the malignant lesion shows twelve steps (different boxes); the benign lesion shows three such steps.

[0050] Note that growth of the cluster occurs at many intensity levels --- even for the benign case. In noisy images, growth of the cluster will occur at virtually every intensity level, regardless of the characteristics of the underlying physiological object being studied. The step feature effectively filters out many of the small incidental changes in cluster growth and is relatively insensitive to the placement of the landmark.

Validation of invention for discriminating benign from malignant breast lesions

[0051] Data were obtained from a study conducted under grant R44CA85101 issued by the National Cancer Institute (NCI) to develop and evaluate a computer-aided-diagnosis (CAD) system for breast MRI.

[0052] Breast MRI images for 169 cases were collected from 5 institutions using 6 different imaging systems/protocols. 35 radiologists participated in the study, each interpreting 60 cases. Approximately % of the cases read by each radiologist were cancers and approximately '/ of the cases were benign. All cases were biopsy proven. Each case was interpreted by a minimum of 6 and a maximum of 19 radiologists. The multiple readers resulted in a sample of landmarks, and consequently a sample of step function values, for each case.

[0053] The CAD system calculated the step function and converted values into a binary Intensity feature, using a threshold of 9 steps: steps <_ 9 resulted in Intensity = 0;
steps > 9 resulted in Intensity =1. Performances for both radiologists' interpretations of intensity and the step function were evaluated using accuracy of predicted pathology when compared to the true pathology. The step function produced from the invention was more accurate than the radiologists in 34 of the 35 tests. The average accuracy of the step function was 66%. The average accuracy of the radiologists' interpretations of intensity was 56%. The accuracies are shown in Table 1 below and are also graphically illustrated in Fig. 6.

Radiologist Interpretation Invention Output 1 0.43 0.57 0.43 0.58 0.47 0.58 N 0.52 0.58 0.48 0.6 0.62 0.6 D 0.53 0.63 0.5 0.65 0.55 0.65 I 0.55 0.65 0.55 0.65 0.55 0.65 V 0.57 0.65 0.58 0.65 0.48 0.67 0.5 0.67 0.55 0.67 0.57 0.67 R 0.58 0.67 0.58 0.67 0.6 0.67 A 0.62 0.67 0.55 0.68 0.58 0.68 D 0.58 0.68 0.6 0.68 0.63 0.68 1 0.53 0.7 0.58 0.7 0.58 0.7 0 0.62 0.7 0.65 0.7 0.57 0.72 L 0.57 0.72 0.58 0.72 Average 56% 66%

[0054] The goal of the NCI-supported study was to have an unbiased set of test images. In spite of this goal, the cancer lesions in the test set were larger than benign lesions in the test set by a statistically significant amount. An analysis was conducted to determine if the apparent discrimination derived from the number of steps was a direct consequence of lesion size bias in the test set.

[0055] The lesion size as interpreted by the radiologist was estimated to be (Area = width x height) of a rectangular region-of-interest (ROI) drawn by the radiologist around the lesion. To reduce the effect of very large and very small lesions and force consensus among readers on the number of steps for each case, analysis was restricted to evaluations that met the following criteria:

ROI Area >= 300 pixels;
ROI Area < = 1800 pixels; and Number of steps = consensus (mode) number of steps for readers who interpreted that case.
138 cases (69 benign; 69 cancers) had ROIs that met the criteria and non-ambiguous modes for number of steps among readers.

[0056] For each of case, the mean ROI_area was computed for radiologists that met the evaluation criteria.
Mean ROI areas for cancers and benigns were as follows:
Cancers: 987 pixels; and Benigns: 877 pixels.

[0057] While there remained a difference in size, this difference was not statistically significant using 2-tailed t-test with p=.05.

[0058] Using the consensus (mode) steps for each of the 138 cases, the mean number of steps for cancers and benigns were as follows:
Cancers: 10.9; and Benign: 7.3.

This difference was statistically significant with p-value < .0001.

[0059] It is not fully understood why the step feature is so effective at discriminating benign from malignant breast lesions on MRI while being relatively independent of other known discriminatory features. We conjecture that the step feature combines aspects of gradient, shape, intensity and size that are each discriminatory and that the following factors come into play in the step features.

[0060] Benign lesions are believed, in many cases, to have sharper edges (steeper gradients) than malignant lesions.
This is likely due to angiogenesis of malignant lesions that cause their boundaries to appear fuzzy on MRI. On the other hand, the sharp delineations between benign lesions and surrounding tissue result in a steep gradient. The sharper edges of benign lesions cause them to have larger, but fewer, steps before reaching maximum intensity.

[0061] Malignant lesions are believed, in many cases, to have stronger source signals than benign lesions. The strong source signal may give rise to more opportunity for steps.

[0062] When source signals are transformed to image intensities on the display, the benign lesions may be "stretched out" (windowed) to fill a portion of the displayable intensity range that is comparable to that of a malignant lesion. If a malignant lesion has a strong source signal, it does not have to be stretched out as much as a low-intensity benign lesion for display purposes. Stretching out a low-intensity signal has the effect of adding low order bits to the output signal - if an image spans 64 intensity levels it can be represented by 6 bits; if a stretched image spans 256 levels, it requires 8 bits of-data. The inclusion of the two low order bits may introduce noise that results in the addition of more steps, effectively masking out some of the discriminatory effectiveness of the measured gradient.

Evaluating changes in the rectangular ROI derived from the contour instead of changes in the cluster mitigates some of the effect of the noise induced by the stretching operation.
In particular, noise on pixels that lie on or adjacent to the boundary of the cluster, but lie properly within the rectangular ROI will frequently affect the cluster shape, but not the ROI shape.

[0063] Malignant lesions tend to have more spiculations and irregularities than benign lesions. Growth on these spiculations may induce a larger number of changes in the rectangular hull than the growth in the benign lesions.

[0064] As noted above, imaged malignant lesions may be, in general, larger than imaged benign lesions. This may be due to a variety of reasons, including the sharper gradients of benign lesions that may make them visible at a smaller size and the more rapid growth of malignant lesions between patient visits. The larger malignant lesions have more opportunity for steps. The analysis presented above indicates that this is not the only discriminatory factor, but it still may be a contributing factor to the effectiveness of this feature.

[0065] Referring now to Fig. 7, an article of manufacture or a computer program product 1000 of one embodiment of the invention is illustrated. The computer program product 1000 includes a recording medium 1002, such as, a floppy disk, a high capacity read only memory in the form of an optically read compact disk or CD-ROM, a tape, a transmission type media such as a digital or analog communications link, or a similar computer program product.
Recording medium 1002 stores program means 1004, 1006, 1008, and 1010 on the medium 1002 for carrying out the methods for implementing the image evaluation method according to the embodiments of the invention.

[0066] A sequence of program instructions or a logical assembly of one or more interrelated modules defined by the recorded program means 1004, 1006, 1008, 1010, direct the image evaluation system 100 for implementing the image evaluation of the embodiments of the invention.

[0067] According to another embodiment of the invention, diagnostic information can be produced for a large area of the breast, and even the entire breast, based on the acquired source data. This embodiment will be described below with respect to two-dimensional source data obtained from a cross-sectional "slice" through a breast, i.e. a thin region bounded by two parallel planes. This data is constituted by a matrix, or grid array, of pixels measuring perhaps 256 x 256 pixels. Each pixel has a value representative of image intensity at a corresponding point within the matrix.

[0068] According to the embodiment here under consideration, the procedure described earlier herein can be performed using each pixel in the matrix as a starting point, or landmark. As a result, for each pixel, p, there is produced a numerical value, or step value, equal to the number of steps (number of changes in the bounding box size) that occur when the procedure is carried out over the intensity level range Imax (p) to Imin (p) =

[0069] The result is a set of data containing a step value for each pixel in the original matrix.

[0070] Then, adjacent pixel locations having the same step value are grouped into clusters, the step values are divided into ranges, e.g. >9 and 59, and a map of the "slice"
is produced in which each cluster having step values >9 is displayed in a first color, e.g. red, and each cluster having step value 59 is displayed in a second color, e.g. white. Or alternatively, each cluster having step value s 9 could be set transparent, revealing the pixel intensity shown on the original image.

[0071] Of course, it may prove desirable to divide the step values into a larger number of ranges each associated with a distinct color. Each color would represent a different likelihood that the area in question is malignant, where the greater the number of step values, the greater the chances of malignancy.

[0072] According to an alternative of this embodiment, needed computing power and processing time can be reduced by not determining a step value for each pixel. Instead, a representative step value is computed for a cluster of pixels and that representative step value is assigned to each pixel in the cluster.

[0073] More specifically, at the beginning of the procedure, 4-connected clusters of pixels having the highest intensity value are constructed using data from the SUB image.
Each such cluster is evaluated to determine if it satisfies minimum criteria according to user-defined parameters.

Specifically, a cluster must satisfy the following minimum criteria: (1) when the number of pixels in the cluster is multiplied by a factor that relates pixel size to physical size in the subject that has been imaged (for example, each pixel may correspond to a physical area of 0.8 mm2 on the subject), the product must equal or exceed a minimum value, set by a parameter (for example, 25 mm2) ; and (2) for each pixel in the cluster, the percentage enhancement of the POST

over the PRE is determined by first subtracting the pixel value of the PRE from the pixel value of the POST and then dividing the difference by the value of the PRE, and the average percentage enhancement is computed for all pixels within the cluster, and this average must equal or exceed a minimum value, set by a parameter (for example 80%). If no cluster satisfies the minimum criteria, then clusters of the next lower intensity value are tested; this process is continued until an intensity level is found for which a cluster satisfies the minimum criteria. If no cluster satisfies the minimum criteria for any intensity level, then all pixels in the image are assigned a step value of 0.

[0074] If multiple clusters satisfy the minimum criteria at some intensity level in this iterative process, then the cluster of largest size is used as the starting cluster in the procedure and the intensity level at which this cluster was computed is the "starting intensity level." Once the starting cluster has been determined, the cluster is filled-in to include all pixels that lie within the boundary of the cluster. The average gradient around the perimeter of the cluster is then computed using data from the POST image.
The cluster is then grown by adding to it all pixels having intensity values at the next lower intensity level ("starting intensity level" -1) that are adjacent to the starting cluster. The average gradient of the grown cluster is then computed and compared to the average gradient of the starting cluster. If the average gradient of the grown cluster is less than the average gradient of the starting cluster, then Ima" is set equal to "starting intensity level." If the average gradient of the grown cluster is more than the average gradient of the starting cluster, then the grown cluster is used as a new starting cluster and the next lower intensity level is used a a new starting intensity level, and the procedure is continued at this decremented intensity level.

[0075] Once Imax is set, the procedure continues as described in the original patent application by determining Imi.n and computing the step value as the number of changes in bounding boxes between Imax and Imin= All pixels within the filled-in cluster at ImaX are assigned a step value equal to the step value computed in this way. Based on present studies, it is believed that Imax should correspond to a ROI
corresponding to a physical area of not greater than 400mm2 on the subject.

[0076] The process is then repeated using the same procedure, but may or may not exclude from consideration all pixels that have been assigned a step value. By way of example, the given area of the filled-in cluster at Imax may measure 100 mm2 or more on the source data image. The resulting data is then processed as described above to produce the desired map and the corresponding color coded display.

[0077] The above description defines the discriminating feature in terms of a function of clusters on the 2-dimensional slice. The invention can be extended to a set of 2-dimensional slices by defining one or more statistics on the set of functions on the individual slices. For example, one statistic may be the maximum function value over the set of 2-dimensional slices. This statistic may be used to identify the enhancement in the breast that has the most cancer-like characteristics and may be useful for assessing the likelihood that a breast, rather than a specific lesion, contains cancer.
A second statistic may a measure of the spatial distribution of clusters that have function values exceeding some threshold. Such a measure may be, for example, the standard deviation of the vertical location of the center of enhancements exceeding a step count of 9 located within the 3-dimensional breast. This statistic may be useful for identifying breasts that do not contain cancer but show patterns that appear to be suspicious that are spread throughout the breast.

[0078] The procedures according to the invention can be adapted to produce three-dimensional step value maps and resulting displays using techniques known in the art, such as those employed in the CGI field.

[0079] It has been found that a significant advantage of the invention is that it provides information at least comparable to that produced by methods that are used to reveal "blooming effects". This is a kinetic method described in the paper "Is the 'blooming sign' a promising additional tool to determine malignancy in MR mammography", by Fischer et al, Eur. Radiol (2004) 14:394-401, published online on 27 September 2003, Springer-Verlag. This method is based on data from several images obtained over a period of time, for example at one minute intervals over a period of eight minutes. For this, the patient is required to remain still for the entire period and any movement will reduce the reliability of the information obtained. In contrast, the present invention requires, at most, only two sets of source data, one before and one after injection of an enhancement agent. Indeed, the invention could be carried out on the basis of image data acquired in only one time frame after injection.

[0080] In connection with determination of measurement landmarks, it would be advantageous to eliminate, or filter out, data relating to areas that are not of interest. For example, when the area of study is the breast, the outer surface of the breast can produce image data that would be interpreted as points of enhancement. This data can be filtered out. Similarly, data relating to the region within the chest wall can be filtered out.

[0081] During the procedure for determining the value to be used for Imax, which is based on local gradients, a first maximum gradient value may be identified, and then a second maximum gradient value may be encountered in a region of much lower pixel intensities. If this should occur, the Imax value will be determined on the basis of the first maximum gradient encountered. On the other hand, if a second maximum gradient value is encountered in a region of pixel intensities only slightly less than those associated with the first maximum gradient value, then Imax may be determined on the basis of the second maximum gradient encountered.

[0082] Figure 9 is a chart illustrating various features and advantages of the invention and Figures 10-13 are reproductions of individual panels of Figure 9. Figure 14 is a screen capture of a display of diagnostic results produced according to the invention.

[0083] It is to be understood by one of ordinary skill in the art that the methods set forth herein for evaluating images can be combined with other methods of evaluation to provide even better accuracy of results.

[0084] The foregoing description of the specific embodiments will so fully reveal the general nature of the invention that others can, by applying current knowledge, readily modify and/or adapt for various applications such specific embodiments without undue experimentation and without departing from the generic concept, and, therefore, such adaptations and modifications should and are intended to be comprehended within the meaning and range of equivalents of the disclosed embodiments. It is to be understood that the phraseology or terminology employed herein is for the purpose of description and not of limitation. The means, materials, and steps for carrying out various disclosed functions may take a variety of alternative forms without departing from the invention.

[0085] Thus the expressions "means to..." and "means for...", or any method step language, as may be found in the specification above and/or in the claims below, followed by a functional statement, are intended to define and cover whatever structural, physical, chemical or electrical element or structure, or whatever method step, which may now or in the future exist which carries out the recited function, whether or not precisely equivalent to the embodiment or embodiments disclosed in the specification above, i.e., other means or steps for carrying out the same functions can be used; and it is intended that such expressions be given their broadest interpretation.

Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.Rectangle;
import java.awt.Dimension;
import java.awt.image.RenderedImage;
/**
* Computes a simple segmentation based on the following parameters:
*
* <ul>
*
*<li> The outer boundary of the mass is controlled by the amount it * intersects with the rectangle frame (frameSize, frameOverlay) *<li> The interior boundary of the mass is controlled by the desired * size relative to the mass (minimumMasslnterior) *
* </ul>

public class BasicSegmenter extends Segmenter {
private int frameSize;
private double frameOverlay;
private double minimumMasslnterior;
private SegmentationBuffer segmentationBuffer;

public BasicSegmenter(int size, double overlay, double minInt) {
frameSize = size;
frameOverlay = overlay;
minimumMassInterior = minInt;
segmentationBuffer = new SegmentationBuffer();
segmentationBuffer.setFrameSizeAndOverlap(frameSize, frameOverlay);
segmentationBuffer.setMinimumInterior(minimumMassInterior);
}

public Segmentation segment(RenderedImage img) {
SegmentationBuffer buffer = new SegmentationBuffer();
buffer.setFrameSizeAndOverlap(frameSize, frameOverlay);
buffer.setMinimumInterior(minimumMassInterior);
return new Segmentation(new Dimension(img.getWidth(), img.getHei.ght()), buffer.segment(img));
}

/**
* Segment the region of interest within the given image.
public Segmentation segment(RenderedImage img, Rectangle roi) {

RenderedImage subImage = ImageUtil.crop(img, roi);
byte[] segData = segmentationBuffer.segment(subImage);
Segmentation seg = new Segmentation(new Dimension(img.getWidth(), img.getHeight()));
seg.setRoiData(segData, roi);

seg.setThresholdRange(segmentationBuffer.getThresholdRange());
return seg;
}
public Segmentation fastSegment(RenderedImage img, Rectangle roi) {
RenderedImage subImage = ImageUtil.crop(img, roi);
byte[] segData = segmentationBuffer.fastSegment(subImage);
Segmentation seg = new Segmentation(new Dimension(img.getWidth(), img.getHeight()));
seg.setRoiData(segData, roi);
seg.setThresholdRange(segmentationBuffer.getThresholdRange());
return seg;
}
}
* * * * * * *

Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

/**
* A set of property values to be used to distinguish clusters.
public class ClusterCondition {

protected int filledSize = 0;
protected int enhancement = 0;
protected int ppa = 0;
protected int minStepValue = 0;
public ClusterCondition() {}

public void setFilledSize(int i) { filledSize = i; }
public int getFilledSize() { return filledSize; }
public void setEnhancement(int i) { enhancement = i; }
public int getEnhancement() { return enhancement; }
public void setPPA(int i) { ppa = i; }
public int getPPA() { return ppa; }

public void setMinStepValue(int i) { minStepValue = i; }
public int getMinStepValue() { return minStepValue; }
public String toString() {

StringBuffer buffer = new StringBuffer();
buffer.append("Condition: filledSize=").append(filledSize);
buffer.append(" enhancement=").append(enhancement);
buffer.append(" ppa=").append(ppa);
buffer.append(" minStepValue=").append(minStepValue);

return buffer.toString();
}
}
* * * * * * *
// copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.geom.Point2D;
import java.awt.Rectangle;
import java.util.*;
import com.alanpenn.util.Stat;
~**
* This class computes the dispersion of cluster locations. The * amount of dispersion may correlate with pathology; the higher the * dispersion, the more likely we are dealing with noise.

public class ClusterDispersion {
~**
* Dispersion is computed as a function of the y-coordinate of the * center of the clusters.

public static double computeDispersion(Collection clusterRectangles) {
Stat yStat = new Stat();
for (Iterator rects=clusterRectangles.iterator(); rects.hasNext(); ) {
ClusterRectangle cr = (ClusterRectangle) rects.next();
Point2D center = cr.computeCenter();
yStat.addData(center.getY());
}

return yStat.computeStd();
}

/**
* Compute a step-adjustment value based on the dispersion of the * cluster locations, and the passed coeffificient.

public static int computeAdjustment(double dispersion, double coefficientOfDispersion) {

int adjustment = 0;
if (dispersion > 0) adjustment = (int) Math.ceil(coefficientOfDispersion * dispersion);
else adjustment = 5;
return adjustment;
}
}

* * * * * * *

Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.Dimension;
import java.awt.Point;
import java.awt.Color;
import java.awt.image.DataBufferShort;
import java.awt.image.RenderedImage;
import java.awt.image.Raster;
import java.awt.image.BufferedImage;
import java.awt.image.WritableRaster;
import java.awt.Rectangle;
import java.util.*;
import com.alanpenn.util.Stat;
/**
* Used for connected component analysis of gray-scale images.
*
* @see SegmentationBuffer public class ClusterMap5 {
private static final int _E_ = 0, -NE = 1, _N_-- 2, _NW_ = 3, _W_ = 4, _SW_ = 5, _S_ = 6, SE = 7;

public static final int CLOCKWISE = -1, COUNTER CLOCKWISE = 1;

private static final short MAX ZD = 10000;

// "deactivated" pixels are ignored by all methods public static final Integer DEACTIVATED = new Integer(Integer.MAX VALUE);
// these IDs are reserved for customization by the caller public static final Integer SPECIAL_1 = new Integer(Integer.MAX_VALUE-1);
public static final Integer SPECIAL_2 = new Integer(Integer.MAX_VALUE-2);
public static final Integer SPECIAL_3 = new Integer(Integer.MAX VALUE-3);
// background pixels are eligible for cluster inclusion public static final Integer BACKGROUND = new Integer(Integer.MAX VALUE-10);
private static final Point[] d4 = {new Point(1,0), new Point(0,-1), new Point(-1,0), new Point(0,1)};
private static final byte DA = OxOO;
private static final byte BG = Ox01;
private static final byte FG = 0x02;
private static final byte INNER = 0x04;
private static final byte BDRY = 0x08;

private static final int MAX_IMAGE_SIZE = 512;

private int DEBUG = 0;
private Integer[][] idMap;

// used to operate on a specific cluster private byte [] [] tempMap;

// used to compute interior pixels and statistics private byte [] [] interiorMap;

// a counter that marks the next available cluster ID
private int nextID;

// Integer (ID) -> Cluster private Hashtable clusters;

// the raster used to build the connect components private Raster featureMap;

private Raster intensityMap;
private short[][] enhancementMap;
// The default cluster comparator private ClusterComparator comparator;

// the sub-region of the map under consideration private Rectangle roi;

// the last threshold used to compute the connect components private int threshold;

// the size of the cluster map private Dimension mapSize;
public ClusterMap5() {

idMap = new Integer[MAX_IMAGE_SIZE][MAX_IMAGE_SIZE];
tempMap = new byte[MAX_IMAGE_SIZE+2][MAX_IMAGE_SIZE+2];
interiorMap = new byte[MAX_IMAGE_SIZE+2][MAX_IMAGE_SIZE+2];
comparator = new ClusterComparator();
clusters = new Hashtable();
// massity comparator comparator.setType(ClusterComparator.FOUR);
roi = new Rectangle();
nextID = 0;
threshold = 0;
}

/**
* @deprecated public void initializeMap(Dimension d) {

if (d.width > MAX IMAGE_SIZE 11 d.height > MAX IMAGE_SIZE) throw new IllegalArgumentException("Cannot handle images greater than +
MAX_IMAGE_SIZE + "x" +
MAX IMAGE SIZE);

mapSize = d;
clusters.clear();
roi.setBounds(0,0,d.width, d.height);
nextID = 0;
threshold = 0;
}

public int getThreshold() { return threshold; }
public Dimension getSize() { return mapSize; }
public void setMaps(RenderedImage intensityMap, RenderedImage featureMap, short [] [] enhancementMap) {
Dimension d = new Dimension(intensityMap.getWidth(), intensityMap.getHeight());

if (d.width > MAX_IMAGE_SIZE II d.height > MAX_IMAGE_SIZE) throw new IllegalArgumentException("Cannot handle images greater than +
MA.X_IMAGE_SIZE + "x" +
MAX IMAGE SIZE);

if (featureMap.getWidthO != d.width fI featureMap.getHeightO
d.height) throw new IllegalArgumentException("Feature map must have size "+
d.width+"x"+d.height);
if (enhancementMap.length != d.width 11 enhancementMap[0].length !=
d.height) throw new IllegalArgumentException("Enhancement map must have size " +
d.width+"x"+d.height);
this.intensityMap = intensityMap.getData();
this.featureMap = featureMap.getData();
this.enhancementMap = enhancementMap;
mapSize = d;
clusters.clear();
roi.setBounds(0,0,d.width, d.height);
nextID = 0;
threshold = 0;
}
~**
* @deprecated public void setIntensityMap(RenderedImage bi) {

Dimension d = new Dimension(bi.getWidth(), bi.getHeight());
initializeMap(d);
intensityMap = bi.getData();

}

~**
* @deprecated public void setFeatureMap(RenderedImage bi) {
Dimension d = new Dimension(bi.getWidth(), bi.getHeight());
initializeMap(d);
featureMap = bi.getData();
}

~**
* @deprecated public void setEnhancementMap (short [] [] m) {
enhancementMap = m;
}

public int getNumberOfClusters() { return clusters.size(); }
public Integer getID(int x, int y) {

if (x < roi.x 11 y < roi.y 11 x >= roi.x+roi.width 11 y >= roi.y+roi.height) return DEACTIVATED;
else return idMap [x] [y] ;
}

public Integer getID(Point p) {
return getID(p.x, p.y);
}
private Integer[] getNeighborIDs(Integer[] ids, Point p, Integer[][] map) {
if (ids.length < 4) ids = new Integer[4];
ids[O] = getID(p.x, p.y-1);
ids [1] = getID (p.x+l, p.y) ;
ids[2] = getID(p.x, p.y+l);
ids [3] = getID (p.x-1, p.y) ;
return ids;
}
~**
* Initializes the idMap based on the given mask.

public void initialize(BitMap mask, Integer trueID, Integer falseID) {
if (mask.getWidth() != mapSize.width 11 mask.getHeight() != mapSize.height) throw new IllegalArgumentException("Bitmap does not have expected size: "+

mapSize.width+"x"+mapSize.height);
for (int x=0; x<mapSize.width; x++) {

for (int y=0; y<mapSize.height; y++) {
idMap [x] [y] = mask.getBit(x, y) ? trueID : falseID;
}
}
clusters.clear();
nextID = 0;
roi.setBounds(0,0,mapSize.width,mapSize.height);
}

~**
* Initializes the idMap based on the given mask.

public void setIdsFromBitMap(BitMap mask, Integer trueID) {
if (mask.getWidth() != mapSize.width 11 mask.getHeight() != mapSize.height) throw new IllegalArgumentException("Bitmap does not have expected size: "+

mapSize.width+"x"+mapSize.height);
for (int x=0; x<mapSize.width; x++) {
for (int y=O; y<mapSize.height; y++) {
if (mask.getBit(x, y)) idMap [x] [y] = trueID;
}
}
}

~**
* Sets all IDs to BACKGROUND; clears the list of found clusters.
public void initialize(Integer id) {

for(int x=0; x<mapSize.width; x++) {
for (int y=0; y<mapSize.height; y++) {
idMap [x] [y] = id;
}
}

clusters.clear();
nextlD = 0;
roi.setBounds(o,0,mapSize.width,mapSize.height);
}

~**
* Initializes the map to all BACKGROUND IDs.
public void initialize() {
initialize(BACKGROUND);
}

public void initialize(Set set, Integer id) {

for (Iterator points=set.iterator(); points.hasNext(); ) {
Point p = (Point) points.next();
idMap [p.x] [p.y] = id;

}
}
~**
* Initialize all active pixels in the cluster map the given id.

public void initializeActivatedPixels(Integer id) {
for(int x=0; x<mapSize.width; x++) {
for (int y=0; y<mapSize.height; y++) {
if (idMap[x][y] 1= DEACTIVATED) idMap [x] [y] = id;
}
}

clusters.clear();
nextID = 0;
roi.setBounds(0,0,mapSize.width,mapSize.height);
}

public void setRoi(Rectangle r) {
roi.setBounds(r.x, r.y, r.width, r.height);
}

public Rectangle getRoi() { return roi; }
public Cluster getLargestCluster() {
Cluster maxCluster = null;

if (clusters.size() == 0) return null;

int maxArea = 0;
Integer maxID = null;
Iterator i = clusters.values().iterator();
while ( i . hasNext ( ) ) {
Cluster cluster = (Cluster) i.next();
if (cluster.getFilledSize().intValue() > maxArea) {
maxArea = cluster.getFilledSize().intValue();
maxCluster = cluster;
}
}

return maxCluster;
}

public String toString() {

StringBuffer buffer = new StringBuffer();
Iterator i = clusters.values().iterator();
while ( i . hasNext ( ) ) {
Cluster cluster = (Cluster) i.next();
buffer.append("\n\t" + cluster);
}

return buffer.toString();

}

/**
* Classifies the given pixel into FG, BG, or DA; if FG or BG, * further classifies the pixel as BDRY if it lies on the * d4-margin between the FG and BG.

private byte classify(int x, int y, Integer id) {
// start with a "deactivated" classification byte classification = DA;

if (idMap [x] [y] . equals (id) ) {
if ( (x+1 < mapSize.width) &&
( idMap [x+1] [y] . intValue ( ) >= BACKGROUND. intValue ( ) ) ) classification = FG BDRY;
else if ( (y-1 >= 0) &&
( idMap [x] [y-1] . intValue ( ) >= BACKGROUND. intValue ( ) ) ) classification = FG BDRY;
else if ( (x-l >= 0) &&
( idMap [x-1] [y] . intValue ( ) >= BACKGROUND. intValue ( ) ) ) classification = FG BDRY;
else if ( (y+l < mapSize.height) &&
( idMap [x] [y+1 ] . intValue ( ) >= BACKGROUND. intValue ( ) ) ) classification = FG BDRY;
else classification = FG;

} l se if ( idMap [x] [y] . intValue ( ) < DEACTIVATED. intValue ( ) ) {
if (x+1 < mapSize.width && idMap[x+l][y].equals(id)) classification = BG ~ BDRY;
else if (y-1 >= 0 && idMap[x][y-11.equals(id)) classification = BG BDRY;
else if (x-1 >= 0 && idMap[x-1][y1.equals(id)) classification = BG I BDRY;
else if (y+l < mapSize.height && idMap [x] [y+1] .equals (id) ) classification = BG BDRY;
else classification = BG;
}

return classification;
}

/**
* Must look at p's ID and p's neighbors private byte classify(Point p, Integer id) {
return classify(p.x, p.y, id);
}

private void reclassify(byte[][] map, int x, int y, Integer id) {
// interior mass pixels never change if ( ( (map [x] [y] & FG) > 0) &&
( (map [x] [y] & BDRY) == 0) return;

byte newClass = classify(x, y, id);
if (newClass == map[x][y]) return;
else {
map [x] [y] = newClass;
// if (x+l < mapSize.width) J/ reclassify(map, x+l, y, id);
if (y-1 >= 0) reclassify(map, x, y-1, id);
if (x-1 >= 0) J/ reclassify(map, x-1, y, id);
if (y+l < mapSize.height) J/ reclassify(map, x, y+l, id);
if (x+1 < mapSize.width) map [x+1] [y] = classify (x+1, y, id) ;
if (y-1 >= 0) map [x] [y-1] = classify(x, y-1, id) ;
if (x-1 >= 0) map [x-1] [y] = classify(x-1, y, id) ;
if (y+1 < mapSize.height) map [x] [y+1] = classify(x, y+1, id) ;
}
}
/**
* Computes the connected components of the previously supplied * featureMap. All pixels in the image with pixel values below * <threshold> are ignored.

public Hashtable computeConnectedComponents(int threshold) {
this.threshold = threshold;

if (intensityMap == null) throw new NullPointerException("Intensity map cannot be null.");
if (featureMap == null) throw new NullPointerException("Feature map cannot be null.");
Point p = new Point();
Integer [] ids = new Integer [4] ;
for (int y=0; y<roi.height; y++) {
for (int x=0; x<roi.width; x++) {

p.x = roi.x+x; p.y = roi.y+y;

&&
if ( ( idMap [p . x] [p . y] == BACKGROUND) (featureMap.getSample(p.x, p.y, 0) >= threshold) ) {
ids = getNeighborIDs(ids, p, idMap);
Integer minID = computeMin(ids);

if (minID.intValue() < BACKGROUND.intValue()) {
idMap [p . x] [p . y] = minID;
Cluster cluster = (Cluster) clusters.get(minID);
try {
cluster.addPoint(p);
}

catch (NullPointerException ex) {
System.err.println("Cluster + minID + " does not exist!");
System. exit ( 0 ) ;
}

J/ merge all neighboring clusters with ID less than minID
for (int i=0; i<4; i++) {
if ( (ids [i] . intValue ( ) < BACKGROUND. intValue ( ) ) &&
( ! ids [i] . equals (minlD) ) &&
clusters.containsKey(ids[i]) ) {

// ids[i] is merged with minID
merge(cluster, (Cluster) clusters.get(ids[i]));
deleteCluster (ids [i] ) ;
}
}
}
else {
// create a new cluster Integer id = new Integer(nextID++);
Cluster cluster = new Cluster(id);
clusters.put(id, cluster);
cluster.addPoint(p);
idMap [p=x] [p=y] = id;
}
}
}
}

return clusters;
}

/**
* Initialize the tempory map to FG, BG, & BDRY pixels.
*J
public void initializeCluster(Cluster cluster) {
Integer id = cluster.getID();
for (int x=0; x<mapSize.width; x++) {
for (int y=0; y<mapSize.height; y++) {
tempMap [x] [y] = classify(x, y, id) ;
}
}
}

/**
* Add up the gradient for all pixels in the FG, on the OUTER
* BDRY.

public double computeFrontierGradient(Cluster cluster, Raster raster) {
Rectangle b = cluster.getBounds();
Rectangle bounds = new Rectangle();
bounds.setBounds(Math.max(0,b.x-1), Math.max(0,b.y-1), Math.min(mapSize.width, b.width+2), Math.min(mapSize.height,b.height+2));

double sum = 0.0;
int count = 0;
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {
for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
if ( ( (tempMap [x] [y] & FG) > 0) &&
( (tempMap [x] [y] & BDRY) > 0) &&
( (tempMap [x] [y] & INNER) == 0) ) {
// pixel is OUTER BDRY FG
sum += raster.getSample(x, y, 0);
count++;
}
}
}
return (count > 0 ? sum/count : 0);
}

/**
* Adds 4-connected pixels (with value >= threshold) the given * cluster.

public void addToCluster(Cluster cluster, int threshold) {
this.threshold = threshold;

if (intensityMap == null) throw new NullPointerException("Intensity map cannot be null.");
if (featureMap == null) throw new NullPointerException("Feature map cannot be null.");
Point p = new Point();
Integer id = cluster.getID();
boolean repeat = true;

while (repeat) {

Rectangle b = cluster.getBounds();
Rectangle bounds = new Rectangle();
bounds.setBounds(Math.max(0,b.x-1), Math.max(0,b.y-1), Math.min(mapSize.width, b.width+2), Math.min(mapSize.height,b.height+2));

repeat = false;
for (int y=0; y<bounds.height;. y++) {
for (int x=0; x<bounds.width; x++) {
p.x = bounds.x+x; p.y = bounds.y+y;

if ( ( (tempMap [p.x] [p.y] & BG) == BG) &&
( ( tempMap [p . x] [p . y] & BDRY) == BDRY) &&
(featureMap.getSample(p.x, p.y, 0) >= threshold) ) {
cluster.addPoint(p);
idMap [p . x] [p . y] = id;
reclassify(tempMap, p.x, p.y, id);
repeat = true;

}
}
}
}
}

public Set computePointSet(Cluster cluster) {
Set set = new HashSet();

Rectangle bounds = cluster.getBounds();
Point p = new Point();

for (int y=0; y<bounds.height; y++) {
for (int x=0; x<bounds.width; x++) {
p.x = bounds.x+x;
p.y = bounds.y+y;

if (idMap [p.x] [p.y] .equals (cluster.getID () ) ) set.add(new Point(p));
}
}
return set;
}

public void updateBoundary(List list) {
Iterator i = list.iterator();
while (i.hasNext()) {
Cluster cluster = (Cluster) i.next();
cluster.updateBoundary();
}
}
~**
* Sorts the list of clusters. Note that the idMap is rendered * useless after completion of this action.

public List sortClusters(List list) {
Collections.sort(list, comparator);
return list;
}
~**
* Finds the maximum cluster in the given list according to the * global "comparator".

public Cluster findMax(List list) {
if (list.size() == 0) return null;
int max=0;
for (int i=1; i<list.size(); i++) {
if (comparator.compare(list.get(i), list.get(max)) < 0) max }
return (Cluster) list.get(max);

}

~**
* Return a list of all clusters whose bounding box contains the * given point.

public List getContainerClusters(Point p) {
ArrayList list = new ArrayList();

Iterator i = clusters.values().iterator();
while (i.hasNext () ) {
Cluster cluster = (Cluster) i.next();
if (cluster.getBounds().contains(p)) list.add(cluster);
}
return list;
}

~**
* Return a list of all clusters whose bounding box contains the * given point, and the cluster has the minimum size.
*
* @deprecated public List getContainerClusters(Point p, int minSize) {
ArrayList list = new ArrayList();

Iterator i = clusters.values().iterator();
while (i.hasNext()) {
Cluster cluster = (Cluster) i.next();

if (cluster.getFilledSize().intValue() >= minSize &&
cluster.getBounds().contains(p)) list.add(cluster);
}
return list;
}

~**
* Return the maximum-sized cluster.
public Cluster getMaximumCluster() {

Cluster maxCluster = null;

Iterator i = clusters.values().iterator();
while (i.hasNextO) {
Cluster cluster = (Cluster) i.next();
if (maxCluster == null 11 maxCluster.getFilledSize().intValue() <
cluster.getFilledSize().intValue()) maxCluster = cluster;
}

return maxCluster;

}

~**
* Returns the cluster (or null) that includes the given point, * and that has a filled size that is at least minSize and no * greater than maxiSize.
*
* @deprecated public Cluster getContainerCluster(Point p, int minSize, int maxSize) {
Cluster cluster = null;

Integer id = getID(p.x, p.y);

if ( (0 <= id.intValue () ) && (id.intValue () < BACICGROUND.intValue () ) ) {
cluster = (Cluster) clusters.get(id);

if (cluster != null) {

int clusterSize = cluster.getFilledSize().intValue();
if ( (clusterSize < minSize) 11 (clusterSize > maxSize) cluster = null;
}
}
return cluster;
}

~**
* Returns the cluster (or null) that includes the given point, * and that has a filled size that is at least minSize, no greater * than maxiSize, and has the required minimum enhancement.

public Cluster getContainerCluster(Point p, int minSize, int maxSize, int minEnhancement) {

Cluster cluster = null;
Integer id = getID(p.x, p.y);

if ( (0 <= id.intValue()) && (id.intValue() < BACKGROUND.intValue())) {
cluster = (Cluster) clusters.get(id);

if (cluster != null) {

int clusterSize = cluster.getFilledSize().intValue();

if ((clusterSize < minSize) 11 (clusterSize > maxSize) II
(cluster.computeAverageEnhancement() < minEnhancement) ) {
cluster = null;
}
}
}
return cluster;

}

/**
* Returns the cluster (or null) that includes the given point, * and that has a filled size that is at least minSize, no greater * than maxiSize, and has the required minimum enhancement.

public Cluster getContainerCluster(Point p, ClusterConstraints constraints) {

Cluster cluster = null;
Integer id = getID(p.x, p.y);

if ( ( 0 <= id. intValue ( ) ) && ( id . intValue ( ) < BACKGROUND. intValue ( ) ) ) {
cluster = (Cluster) clusters.get(id);

if (cluster != null) {

int clusterSize = cluster.getFilledSize().intValue();
if ( (clusterSize < constraints.getMinClusterSize()) II
(clusterSize > constraints.getMaxClusterSize()) (cluster.computeAverageEnhancement() <
constraints.getMinClusterEnhancement()) ) {

cluster = null;
}
}
}
return cluster;
}

/**
* Returns the cluster closest to p.

public Cluster getClosestCluster(Point p, int size) {
Cluster closestCluster = null;
double minDistance = Double.MAX_VALUE;
double distance;
Point s = new Point();
Iterator i = clusters.values().iterator();
while ( i . hasNext ( ) ) {
Cluster cluster = (Cluster) i.next();
if (cluster.getFilledSize () . intValue () >= size) {
Rectangle bounds = cluster.getBounds();

// lst point s.x = bounds.x; s.y = bounds.y;
double min = distance(s,p);

// 2nd point s.x = bounds.x; s.y = bounds.y+bounds.height-1;
distance = distance(s,p);
if (distance < min) min = distance;

// 3rd point s.x = bounds.x+bounds.width-1; s.y = bounds.y;
distance = distance(s,p);
if (distance < min) min = distance;
// 4th point s.x = bounds.x+bounds.width-1; s.y = bounds.y+bounds.height-i;
distance = distance(s,p);
if (distance < min) min = distance;
if (min < minDistance) {
minDistance = min;
closestCluster = cluster;
}
}
}
return closestCluster;
}

/**
* Adds the interior points for all clusters in the given list.
public void addInteriorPoints(List list) {

Iterator i = list.iterator();
while (i.hasNext()) {
Cluster cluster = (Cluster) i.next();
addInteriorPoints(cluster);
}
}
/**
* Removes this location from further consideration.

public void deactivatePixel(int x, int y) {
idMap [x] [y] = DEACTIVATED;
}
/**
* Removes this cluster's pixels from further consideration.

public void deactivateCluster(Cluster cluster) {
for (int x=0; x<mapSize.width; x++) {
for (int y=0; y<mapSize.height; y++) {
if (cluster.getID () == idMap [x] [y] ) idMap[x][y] = DEACTIVATED;
}
}
}
/**
* Set's the internal ID of the pixel public void setID(int x, int y, Integer id) {
idMap [x] [y] = id;

}

~**
* Sets the interior statistics for each cluster in the given * list.

public void setInteriorPoints(List list) {
Iterator i = list.iterator();
while (i.hasNext () ) setInteriorPoints((Cluster) i.next());
}

~**
* Sets the cluster's interior statistics public void setInteriorPoints(Cluster cluster) {
Stat stat = computelnteriorStat(cluster);
cluster.setInteriorStat(stat);
}

~**
* Find all clusters whose area (including interior) is at least * <minSize>.
*
* @deprecated public List computeListOfClustersBySize(int minSize) {
ArrayList list = new ArrayList();

Iterator i = clusters.values().iterator();
while (i.hasNext () ) {
Cluster cluster = (Cluster) i.next();
int area = cluster.getFilledSize().intValue();
if (area >= minSize) list.add(cluster);
}
return list;
}

~**
* Find all clusters whose area (including interior) is at least * <minSize>.
*
* @deprecated public List computeSortedListOfClustersBySize(int minSize, int minEnhancement) {
TreeMap map = new TreeMap();

Iterator i = clusters.values().iterator();
while (i.hasNext()) {
Cluster cluster = (Cluster) i.next();
int area = cluster.getFilledSize().intValue();
if ( (area >= minSize) && cluster.computeAverageEnhancement() >=
minEnhancement) map.put(new Integer(area), cluster);
}

return new ArrayList(map.values());
}
/**
* Find all clusters whose area (including interior) is at least * <minSize>.

public List computeSortedListOfClustersBySize(ClusterConstraints constraints) {
TreeMap map = new TreeMap();

Iterator i = clusters.values().iterator();
while (i . hasNext ( ) ) {
Cluster cluster = (Cluster) i.next();
int filledArea = cluster.getFilledSize().intValue();

if ( (filledArea >= constraints.getMinClusterSize()) &&
cluster.computeAverageEnhancement() >=
constraints.getMinClusterEnhancement() &&
cluster.getAverageIntensity().doubleValue() >
constraints.getMinSeedPeakIntensity() ) map.put(new Integer(filledArea), cluster);
}
return new ArrayList(map.values());
}

/**
* Adds the points (maybe none) that fall inside the component * with <code>ID==id</code>.

public void addInteriorPoints(Cluster cluster) {
Point p = new Point();
Integer id = cluster.getIDO;
byte bg = 0;
byte fg = 1;
byte interior = 2;
// copy the map Rectangle bounds = cluster.getBounds();
Rectangle bounds2 = new Rectangle(bounds.x-1, bounds.y-1, bounds.width+2, bounds.height+2);

for (int y=0; y<bounds2.height; y++) {
for (int x=0; x<bounds2.width; x++) {
// form outer-frame of bg pixels if (x==0 11 y==o 11 x== bounds2.width-1 11 y == bounds2.height-1) interiorMap [x] [y] = bg;

else {
// translate to image space p.x = bounds2.x + x;
p.y = bounds2.y + y;
interiorMap [x] [y] = idMap [p . x] [p . y] . equal s ( id) ? fg interior;

}
}
}

The cluster is surrounded by "bg". We want to leak the /J background (bg) to surround the cluster. It will not spread to the cluster's interior (if it exists).
boolean done = false;
while (! done) {
done = true;

for (int y=1; y<bounds2.height-1; y++) {
for (int x=1; x<bounds2.width-1; x++) {

if ( interiorMap [x] [y] == interior ) {

// check the north and west neightbors if ( (interiorMap [x] [y-1] == bg) (interiorMap [x-1] [y]
bg) interiorMap [x] [y] = bg;
done = false;
}
}
}
}

for (int y=bounds2.height-2; y>=l; y--) {
for (int x=bounds2.width-2; x>=l; x--) {
if ( interiorMap[x][y] == interior ) {

// check the south and east neighbors if (interiorMap [x] [y+1] == bg 11 interiorMap [x+1] [y]
bg) {

interiorMap [x] [y] = bg;
done = false;
}
}
}
}
}

At this point, all pixels in interiorMap outside the cluster have an ID == bg; pixels inside have an ID == interior. We update the idMap with the new interior pixels.

for (int y=1; y<bounds2.height-1; y++) {
for (int x=1; x<bounds2.width-1; x++) {
// translate to image space p.x = bounds2.x + x;
p.y = bounds2.y + y;

&&
if ( (interiorMap [x] [y] == interior) ( idMap [p . x] [p . y] . intValue ( ) < DEACTIVATED. intValue ( ) ) ) {
cluster.addInteriorPoint(p);

if ( idMap [p . x] [p = y] . intValue ( ) < BACKGROUND. intValue ( ) ) deleteCluster (idMap [p . x] [p . y] ) ;

idMap [p . x] [p . yl = cluster. getID ( ) ;
}
}
}
}
/**
* Computes the statistics of the points that lie in the interior * of the given cluster.

public Stat computeInteriorStat(Cluster cluster) {
Point p = new Point();
Integer id = cluster.getIDO;
byte bg = 0;
byte fg = 1;
byte interior = 2;
Stat stat = new Stat();

Rectangle b = cluster.getBounds();
Rectangle bounds = new Rectangle();
bounds.setBounds(Math.max(0,b.x-1), Math.max(0,b.y-1), Math.min(mapSize.width, b.width+2), Math.min(mapSize.height,b.height+2));

for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {
tempMap [x] [y] J= INNER;
}
}

The cluster is surrounded by id==bg. We want to spread the background (bg) to fill the exterior of the cluster. It // will not spread to the cluster's interior (if it exists).
boolean done = false;
while (! done) {
done = true;

for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {
if ( ( (tempMap [x] [y] & FG) == 0) &&
( (tempMap [x] [y] & INNER) > 0) ) {
pixel is INNER BG

check the north and west neightbors if ( (y == bounds.y) north neighbor is out of bounds (x == bounds.x) west neighbor is out of bounds ((tempMap[x][y-1] & FG) 0) north neighbor is BG
((tempMap[x-1][y] & FG) 0) ) { // west neighbor is BG

change pixel to OUTER
tempMap [x] [y] '= INNER;
done = false;
}
}
}
}

for (int y=bounds.y+bounds.height-1; y>=bounds.y; y--) {
for (int x=bounds.x+bounds.width-1; x>=bounds.x; x--) {
if ( ( (tempMap [x] [y] & FG) == 0) &&
((interiorMap[xl[y] & INNER) > 0) ) {
pixel is INNER BG

// check the south and east neighbors if ( (y == bounds.y+bounds.height-1) south neighbor is out of bounds (x == bounds.x+bounds.width-1) east neighbor is out of bounds ((tempMap[x][y+1l & FG) 0) south neighbor is BG
((tempMap[x+1][y] & FG) 0) ) { // east neighbor is BG

change pixel to OUTER
tempMap [x] [y] '= INNER;
done = false;
}
}
}
}
}

Add INNER BG pixels to stat.
Change all neighbors of OUTER BG BDRY pixels to OUTER.
for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {

if ( ( (tempMap [x] [y] & FG) == 0) ) {
// pixel is BG

if ((tempMap[x] [y] & INNER) > 0) {

pixel is INNER BG
stat.addData(featureMap.getSample(p.x, p.y, 0));
}
else {

// pixel is OUTER BG

{
if ( ( tempMap [x] [y] & BDRY) > 0) pixel is OUTER BG BDRY
mark all neighbors OUTER
if (x+1 < mapSize.width) tempMap [x+1] [y] -INNER;
if (y-1 >= 0) tempMap [x] [y-1] -INNER;
if (x-1 >= 0) tempMap [x-1] [y] &= -INNER;
if (y+1 < mapSize.height) tempMap [x] [y+1] &= -INNER;
}
}
}
}
}
return stat;
}

/**
* Finds the first point with ID equal to "id".

private Point firstBoundaryPoint(Integer id, Rectangle r) {
for (int y=O; y<r.height; y++) {
for (int x=0; x<r.width; x++) {
if (idMap[r.x+x][r.y+y].equals(id)) return new Point(r.x+x, r.y+y);
}
}
return null;
}

/**
* Determines the 8-connect clockwise boundary of cluster id.
* Note that each pixel b on the 8-connected boundary is such * that d4(b,c)=1, where c is the nearest pixel in the cluster.
* Returns the Vector of links.

public ChainCode findEightBoundary (Integer id, Rectangle bounds) {
if (DEBUG > 0) System.out.println("find8Boundary(" + id + ")");

ChainCode cc = new ChainCode();
Point first = firstBoundaryPoint(id, bounds);
if (first == null) return cc;
cc.setFirst(first);
// Find previous neighbor of first point of chain. We do this so that we know when we've returned to our starting point.
int d;
int lastDir;
d = nearestNeighbor(first, id, N, COUNTER CLOCKWISE);
if (d >= 0) {

lastDir = reflect (d) ;
d = lastDir;
Point p = first;
do {
d = nearestNeighbor(p, id, dir(d+2), CLOCKWISE);
p = ChainCode.neighbor(p, d);

cc.addLink(d);
} while ( (d != lastDir) p.equals(first));
}

if (DEBUG > 0) System.out.println("findBBoundary()");
return cc;
}

* Returns the 180 degree reflexion of d, in the range of 0..7.
private int reflect(int d) {
return dir(d+4);
}

/**
* Returns the direction of "d" in the range of 0..7.
private int dir (int d) {
return (d+8) % 8;
}

/**
/** Returns the index of p's nearest neighbor starting at * direction <d> in the dir" direction.

private int nearestNeighbor (Point p, Integer cc, int d, int dir) {

for (int dd=0; dd<8; dd++) {
int rd = dir(d+(dd*dir));
Point r = ChainCode.neighbor(p, rd);
if (onImage (r) && (idMap [r.x] [r.y] .equals (cc) ) ) return(rd);
}
return(-1);
}

~**
* Returns TRUE of the point lies within confines of the image.
private boolean onImage(Point p) {

return (0 <= p.x && p.x < mapSize.width &&
0 <= p.y && p.y < mapSize.height);
}

public class Cluster {
private Integer id;
private ChainCode boundary;
private Rectangle bounds;
public Stat massStat;
public Stat interiorStat;
public Stat enhancementStat;
public Cluster(Integer id) {
this.id = id;
boundary = null;
bounds = new Rectangle();
massStat = new Stat();
interiorStat = new Stat();
enhancementStat = new Stat();
}

public boolean equals(Object obj) {
if (obj instanceof Cluster) {
Cluster cluster = (Cluster) obj;
return (id.equals(cluster.id));
}
else return false;
}

public Integer getID() {
return id;
}

public Integer getUnfilledSize() {
return new Integer(massStat.count());
}

public Integer getFilledSize() {
return new Integer(massStat.count() + interiorStat.count{));
}

public Integer getInteriorSize() {
if (interiorStat == null) return new Integer(0);
else return new Integer(interiorStat.count());
}

public double getPerimeter() {
if (boundary != null) return boundary.d8Length();
else return 0;
}

public Rectangle getBounds() {
return new Rectangle(bounds);
}

public Stat getMassStat() { return massStat; }

public Stat getInteriorStat() { return interiorStat; }
public Double computeCircularity() {

double perimeter = getPerimeter();
int totalArea = massStat.count() + interiorStat.count();
if (totalArea > 0) return new Double((perimeter*peri.meter) / totalArea);
else return new Double(0);
}

public Integer getMaxIntensity() {
return new Integer((int)massStat.getMax());
}

public Integer getMinIntensity() {
return new Integer((int) massStat.getMin());
}

public Double getAverageIntensity() {
if (massStat.count() > 0) return new Double(massStat.computeMean());
else return new Double(0);
}

public Double getInteriorAvgIntensity{) {
if (interiorStat.count() > 0) return new Double(interiorStat.computeMean());
else return new Double(0);
}

/**
* Returns the average enhancement percentage.

public double computeAverageEnhancement() {
return enhancementStat.computeMean();
}

public void addPoint(Point p) {
addPoint(p.x, p.y);
}
public void addPoint(int x, int y) {
if (massStat.count() == 0) bounds = new Rectangle(x, y, 1,1);
else {
// bounds int xMin = Math.min(x, bounds.x);
int xMax = Math.max(x, bounds.x+bounds.width-1);
int yMin = Math.min(y, bounds.y);
int yMax = Math.max(y, bounds.y+bounds.height-1);
bounds.setBounds(xMin, yMin, xMax-xMin+l, yMax-yMin+l);
}

JJ intensity massStat.addData(intensityMap.getSample(x,y,0));
// boundary (we need to update the boundary) boundary = null;

// enhancement if (enhancementMap != null) {
enhancementStat.addData(enhancementMap[x][y]);
}
}
public void addInteriorPoint(Point p) {
// intensity interiorStat.addData(intensityMap.getSample(p.x,p.y,0));
}

public Double computeMassity() {

return new Double(massStat.count() *
massStat.computeMean() *
computeCompactness().doubleValue());
}

public void clear() {

//id = new Integer(-1);
boundary = null;
massStat.clear();
interiorStat.clear();
bounds.setBounds(0,0,0,0);
enhancementStat.clear();
}

public String toString() {

StringBuffer buffer = new StringBuffer();

buffer.append(" id= ).append(id);
buffer.append(" area=").append(massStat.count());
buffer.append(" interiorArea=").append(interiorStat.count());
buffer.append(" intensity=").append(massStat);
if (boundary != null) buffer.append(', perimeter=").append(boundary.getSize());
else buffer.append(" perimeter=null");
buffer.append(" compactness=").append(computeCompactness());
buffer.append(" enhancement=").append(enhancementStat);
return buffer.toString();
}
public void updateBoundary() {
if (boundary null) boundary = findEightBoundary(id, bounds);
}

public ChainCode getBoundary() { return boundary; }
public Double computeCompactness() {
if (bounds.width * bounds.height > 0) {
int totalArea = massStat.count() + interiorStat.count();
return new Double( (double) totalArea / (bounds.width *
bounds.height));
}
else return new Double(0);
}

~**
* @deprecated public Double aggregate() {
if (massStat.count() > 0) return new Double(massStat.count() * massStat.computeMean() *
massStat.count() *
(bounds.width*bounds.height));
else return new Double(0);
}

public Cluster copy() {
Object rval = null;
try {
rval = clone();

catch (CloneNotSupportedException ex) {
ex.printStackTrace();
}
return (Cluster) rval;
}

public void setInteriorStat(Stat s) {

interiorStat = s;
}
}

public class ClusterComparator implements Comparator {

public static final int AREA = 1;
public static final int MAX_INTENSITY = 2;
public static final int AVG_INTENSITY = 3;
public static final int CIRCULARITY = 4;
public static final int COMPACTNESS = 5;
public static final int FOUR = 100;
private int type = AREA;

public void setType(int type) {
this.type = type;
}
public int compare(Object ol, Object o2) {
Cluster c1 = (Cluster) o1;
Cluster c2 = (Cluster) o2;
{
switch (type) case AREA:
return c2.getUnfilledSize().compareTo(cl.getUnfilledSize());
case MAXINTENSITY:
return c2.getMaxIntensity().compareTo(cl.getMaxIntensity());
case AVG_INTENSITY:
return c2.getAverageIntensity().compareTo(cl.getAverageIntensity());
case CIRCULARITY:
return c2.computeCircularity().compareTo(cl.computeCircularity());
case COMPACTNESS:
return c2.computeCompactness().compareTo(cl.computeCompactness());
case FOUR:
return c2.computeMassity().compareTo(cl.computeMassity());
default:
return 0;
}
}
public boolean equals(Object obj) {
return (obj instanceof ClusterComparator);
}

private int foo(int a) {
if (a < 0) return -1;
else if (0 < a) return 1;
else return 0;
}

}

public void restrictSearch(Rectangle r) {
roi = r;
}

public Rectangle getMaxClusterBounds(List list) {
Rectangle r = new Rectangle(0,0,0,0);
Iterator i = list.iterator();
while (i.hasNext()) {
Cluster cluster = (Cluster) i.next();
Rectangle bounds = cluster.getBounds();
if (bounds.width * bounds.height > r.width*r.height) r.setBounds(bounds.x, bounds.y, bounds.width, bounds.height);
}

return r;
}

private Integer computeMin(Integer[] array) {
if (array.length == 0) throw new NullPointerException("Array cannot be null. );
int ndx = 0;
for (int i=l; i<4; i++) {
if (array [i] . compareTo (array [ndx] ) < 0) ndx = i;
}
return array[ndx];
}

~**
* Merges c2 into cl.

private void merge(Cluster cl, Cluster c2) {
Rectangle bounds = c2.getBounds();
for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {
if (idMap [x] [y] . equals (c2 . getID ( ) ) ) {
idMap [x] [y] = cl. getID ( ) ;
cl.addPoint(x,y);
}
}
}
}
private void deleteCluster(Integer id) {
if (clusters.containsKey(id)) clusters.remove(id);
}

public void printCluster(Cluster cluster) {
StringBuffer buffer = new StringBuffer();
int area = 0;

Rectangle bounds = cluster.getBounds();
buffer.append("\n");
for (int y=bounds.y; y<bounds.y+bounds.height; y++) {
for (int x=bounds.x; x<bounds.x+bounds.width; x++) {
if (idMap [x] [y] .equals (cluster.getID () ) ) {
buffer.append("*");
area++;
}
else buffer.append(" ");
}
buffer.append("\n");
}
buffer.append("Area=").append(area).append("\n");
System.out.println(buffer.toString());
}

public Bufferedlmage createBufferedImage(Cluster cluster, BufferedImage image, int sample) {

return createBufferedImage(cluster.getlD(), image, sample);
}

public BufferedImage createBufferedImage(Integer id, BufferedImage image, int sample) {

if (image == null) throw new NullPointerException("image cannot be null.");
WritableRaster raster = image.getRaster();

if (raster.getWidth() != featureMap.getWidth() raster.getHeight() != featureMap.getHeight()) throw new IllegalArgumentException("Incompatible image ("+
image.getWidth()+"x"+image.getHeight());

for (int x=0; x<featureMap.getWidth(); x++) {
for (int y=0; y<featureMap.getHeight(); y++) {
if (idMap [x] [y] == id) raster.setSample(x, y, 0, sample);
else raster.setSample(x, y, 0, 0);
}
}
return image;
}

private double distance(Point pl, Point p2) {
return Math.sqrt(Math.pow(p1.x-p2.x,2.0)+Math.pow(p1.y-p2.y,2.0));
}

}

* * * * * * *
// copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.Point;
import java.awt.Rectangle;
import java.awt.image.*;
import java.util.List;

import com.alanpenn.image.ClusterMapS.Cluster;
import com.alanpenn.util.Stat;

/**
* Not really "intensity" per se, but rather a statistic from the * connected component analysis. We build the connected components * and start counting from the threshold that produces the maximum * gradient boundary, down to the lowest threshold allowed by the * BasicRoiFinder. As we proceed down from the max, the size of the * cluster bounds increases monotonically. We count the number of * expansions along the way down.
*
* The classification is based on 256x256 image size.
public class LevelSetFinder {

public static final Integer UNKNOWN = new Integer(-1);
public static final Integer BENIGN_LIKE = new Integer(0);
public static final Integer CANCER._LIKE = new Integer(l);
private static final int MAX_STEPS = 7;

private LevelSetHelper levelSetHelper;
private Point landmark = new Point();
private int steps;

public static Integer classify(int s) {
if (s < 0) return UNKNOWN;
else if (s < MAXSTEPS) return BENIGN_LIKE;
else return CANCER_LIKE;
}

public LevelSetFinder() {
levelSetHelper = new LevelSetHelper(new BasicSegmenter(2, 0.15, 0.70));
}

public void setImages(RenderedImage p0, RenderedImage pl, BitMap mask) {
levelSetHelper.setImages(p0, pl);
levelSetHelper.setMask(mask);
levelSetHelper.setGradientImage(pl);
landmark.x = -1; landmark.y = -1;

steps = -1;
}
public void setImages(RenderedImage p0, RenderedImage pl) {
setImages(p0, pl, null);
}

public int computeSteps(Point p) {
return computeSteps(p, null);
}

public int computeSteps(Point p, StringBuffer debugBuffer) {
if (debugBuffer != null) debugBuffer.append("IntensityStepFinder: computeSteps("+p+")\n");
if (! landmark.equals(p)) {
landmark.x = p.x; landmark.y = p.y;
Cluster cluster = levelSetHelper.findFirstCluster(p, debugBuffer);
if (cluster != null) levelSetHelper.findAll(cluster, debugBuffer);
}
steps = levelSetHelper.computeStepsFromMaxGradient(debugBuffer);
return steps;
}
public void saveClusterBoundaries(boolean b) {
levelSetHelper.setSaveChainCodeFlag(b);
}

public Point getLandmark() {
return landmark;
}

public int getSteps() {
return steps;
}

public int computeLevelFromCondition(String condition) {
if (landmark.x < 0 11 landmark.y < 0) throw new RuntimeException("Must compute steps first.");
if (condition.equals("first")) return levelSetHelper.getFirstLevel();
else if (condition.equals("max")) return levelSetHelper.computeMaxGradientLevel();
else if (condition.equals("+50")) return levelSetHelper.computeMaxGradientRoiMultiple(1.5);
else if (condition.equals(11+100")) return levelSetHelper.computeMaxGradientRoiMultiple(2);
else if (condition.equals("-25")) return levelSetHelper.computeMaxGradientOffset(.75);
else if (condition.equals("-5011)) return levelSetHelper.computeMaxGradientOffset(.50);
else if (condition.equals("compactdrop")) return levelSetHelper.computeMaxCompactnessDropoffLevel();
else throw new IllegalArgumentException( Don't know how to handle condition: " +
condition);
}

public LevelSetHelper getLevelSetHelper() {
return levelSetHelper;
}

public void setClusterConstraints(ClusterConstraints cc) {
levelSetHelper.setClusterConstraints(cc);
}

public void setMapListener(MapListener lis) {
levelSetHelper.setMapListener(lis);
}

public List getBoundaryListFromMaxGradient() {
return levelSetHelper.getBoundaryListFromMaxGradient();
}

public List computeChainCodeListFromMaxGradient() {
return levelSetHelper.computeChainCodeListFromMaxGradient();
}

public void startStepMapThread(int mapID, StringBuffer debugBuffer) {
if (debugBuffer != null) debugBuffer.append("IntensityStepFinder:computeStepMap()\n");
levelSetHelper.setDebugBuffer(debugBuffer);
levelSetHelper.setMapID(mapID);

new Thread(levelSetHelper).start();
}

public byte[] computeStepMap(int maplD, StringBuffer debugBuffer) {
if (debugBuffer != null) debugBuffer.append("IntensityStepFinder:computeStepMap()\n");
levelSetHelper.setDebugBuffer(debugBuffer);
levelSetHelper.setMapID(mapID);

return levelSetHelper.computeStepMap();
}
}

* * * * * * *
Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.Color;
import java.awt.Dimension;
import java.awt.Rectangle;
import java.awt.Point;

import java.awt.image.*;
import java.util.*;
import java.text.*;
import javax.swing.*;
import com.alanpenn.util.Stat;
import com.alanpenn.image.ClusterMap5.Cluster;
public class LevelSetHelper implements Runnable {
private static final int DEBUG = 1;

public static final int MAX_INTENSITY_LEVEL = 255;
public static final int MIN INTENSITY LEVEL = 7;
private static final Color[] COLORARRAY =
{Color.black, Color.green, Color.green, Color.green, Color.green, Color.green, Color.yellow, Color.yellow, Color.yellow, Color.yellow, Color.yellow, new Color(0, 0, 50), new Color(0, 0, 100), new Color(0, 0, 150), new Color(0, 0, 200), new Color(0, 0, 250), Color.red, Color.red, Color.red, Color.red, Color.red, Color.cyan, Color.cyan, Color.cyan, Color.cyan, Color.cyan, Color.magenta, Color.magenta, Color.magenta, Color.magenta, Color.magenta, Color.pink, Color.pink, Color.pink, Color.pink, Color.pink, Color.white};
private static final String DELIMITER ="\t";
private static final String STATHEADER =
"case\toffset\tindex\tthresh\tx\ty\twidth\theight\tsize\tmax\tavg\tinter\tintMa x \tintAvg\tcompct\tgradnt\tfrmInt\n";

private NumberFormat FORMAT2 = new DecimalFormat(11.0011);
private LevelRoi[] lrs;
private int firstLevel;
private int lastLevel;
private ClusterMap5 clusterMap;
private Segmenter segmenter;
private RenderedImage gradientImage;
private boolean saveChainCode = false;
private Dimension imageDim;
private ImageDialog imageDialog;
private IndexColorModel colorModel;

// the mask to determine the section of the image to process/ignore private BitMap mask;

// when computing steps, monitor progress to this bar private JProgressBar progressBar;

// append details of the operation to this buffer private StringBuffer debugBuffer;

// upon completion of a step map, announce the map to this listener private MapListener mapListener;

the identification (e.g., slice number) for the newly constructed map private int mapID;

// The constraints that guide the connected component algorithm.
private ClusterConstraints clusterConstraints;

// For expediency, a hash of previously computed gradients.
private GradientHash gradientHash;

public LevelSetHelper(Segmenter segmenter) {
this.segmenter = segmenter;
lrs = new LevelRoi[MAX_INTENSITY_LEVEL+1];
clusterMap = new ClusterMapS();
gradientHash = new GradientHash(segmenter);
}

public void setSaveChainCodeFlag(boolean b) { saveChainCode = b; }
public void setClusterConstraints(ClusterConstraints cc) {
clusterConstraints = cc; }
public ClusterConstraints getClusterConstraints() { return clusterConstraints; }

public void setMapListener(MapListener lis) { mapListener = lis; }
public void setDebugBuffer(StringBuffer sb) { debugBuffer = sb; }
public void setMapID(int id) { mapID = id; }

public void setImages(RenderedImage p0, RenderedImage pl) {
imageDim = new Dimension(p0.getWidth(), p0.getHeight());
RenderedImage subImage = ImageUtil.subtract(pO, pl);
short[][] enhancementMap = ImageUtil.computeEnhancementMap(p0, pl);
pi is the peak intensity; sub will be used as the feature // for connected component construction, and the enhancementMap will be used to constrain the seed clusters.
clusterMap.setMaps(pl, subImage, enhancementMap);
}
public void setGradientImage(RenderedImage image) {
gradientImage = image;
gradientHash.setGradientImage(image);
}

public void setMask(BitMap m) { mask = m; }

public Cluster findFirstCluster(Point p, StringBuffer debugBuffer) {
// initialize for (int i=0; i<=MAX_INTENSITY_LEVEL; i++) lrs [i] = null;
firstLevel = -1;
Cluster cluster = null;

set all points on the map to BACKGROUND

if (mask != null) clusterMap.initialize(mask, clusterMap.BACKGROUND, clusterMap.DEACTIVATED);
else clusterMap.initialize(clusterMap.BACKGROUND);
int level = MAX_INTENSITY_LEVEL;
while (level>=clusterConstraints.getMinStartClusterLevel() && cluster null) {

clusterMap.computeConnectedComponents(level);
cluster = clusterMap.getContainerCluster(p, clusterConstraints);
if (cluster != null) {
clusterMap.initializeCluster(cluster);
clusterMap.addInteriorPoints(cluster);
lrs[level] = makeLR(level, cluster, saveChainCode);
firstLevel = level;
}

if (debugBuffer != null && lrs[level] != null) debugBuffer.append("findFirstCluster:"+lrs[level]+"\n");
level--;
}
return cluster;
}

public Cluster findFirstCluster(Point p) {
return findFirstCluster(p, null);
}

public int getFirstLevel() { return firstLevel; }

public void findAll(Cluster cluster, StringBuffer debugBuffer) {

// continue until we get to minStopClusterLevel or 4x the max grad ROI' if (firstLevel < 0) throw new RuntimeException("First level not found.");

if (firstLevel == clusterConstraints.getMinStopClusterLevel()) {
lastLevel = clusterConstraints.getMinStopClusterLevel();
return;
}

assume the first level we found is at the max gradient int maxLevel = firstLevel;

int level = firstLevel-1;
boolean done = (level < clusterConstraints.getMinStopClusterLevel());
while (! done) {

clusterMap.addToCluster(cluster, level);
clusterMap.addInteriorPoints(cluster);
lrs[level] = makeLR(level, cluster, saveChainCode);

reset the gradient marker if the level is legal, the cluster isn't too big, and the new cluster's gradient // exceeds the current max if ( (level >= clusterConstraints.getMinGradientLevel()) &&
(cluster.getFilledSize().intValue() <=
clusterConstraints.getMaxClusterSize()) &&
(lrs [level] .getGradient () > lrs [maxLevel] .getGradient () ) maxLevel = level;

if (debugBuffer != null) debugBuffer.append("findAll:" + lrs[level]+
" new first="+(maxLevel==level)+", hullRatio="+

FORMAT2.format((double)lrs[level].getArea()/lrs[maxLevel].getArea())+
n\nn) .
, // Set termination based on finding the largest cluster whose area is at most 4 times the size of the first cluster if (lrs'[level].getArea() > (4 * lrs [maxLevel] .getArea() ) ) {
done = true;
lastLevel = level+l;

}lse if (level == clusterConstraints.getMinStopClusterLevel()) {
done = true;
lastLevel = level;
}

level--;
}
}
public int computeMaxGradientLevel() {
int level = firstLevel;

for (int i=firstLevel-1; i>=lastLevel; i--) {
if (lrs [i] .getGradient () > lrs [level] .getGradient () ) level = i;
}
return level;
}

public int computeMaxGradientRoiMultiple(double multiple) {
int max = computeMaxGradientLevel();

int level;
for (level=max; level-l>=lastLevel; level--) {
if (lrs [level-1] .getArea() > (multiple*lrs [max] .getArea() ) ) break;
}
return level;
}

public int computeMaxGradientOffset(double multiple) {

int max = computeMa.xGradientLevel();
int level;
for (level=max; level-l>=lastLevel; level--) {
if (lrs [level-1] .getGradient () < (multiple*lrs [max] getGradient () ) ) break;
}
return level;
}

public int computeStepsFromMaxGradient() {
return computeStepsFromMaxGradient(null);
}

~**
*"steps" is the size of the set of rectangular hulls formed by * the clusters between max-gradient and lastLevel. Note that the * minimum number of steps will be 1, in the case that maxgradient * occurs at lastLevel.
*~ .
public int computeStepsFromMaxGradient(StringBuffer debugBuffer) {
Rectangle lastRoi = null;
int count = 0;

for (int level=computeMaxGradientLevel(); level>=lastLevel; level--) {
if (lrs [level] ! = null) {
Rectangle roi = lrs[level].getRoi();
if (! roi.equals(lastRoi)) {
count++;
lastRoi = roi;
if (debugBuffer != null) debugBuffer.append("Rect hull change at level "+level+", cluster size (pixels)="+
lrs [level] . getArea ( ) +"\n" ) ;
}
}
}
return count;
}

~**
* Form a list of chain codes, starting with the one at maximum * gradient, and including those that correspond to increments in * steps.

public List computeChainCodeListFromMaxGradient() {
Rectangle lastRoi = null;
List list = new ArrayList();

for (int level=computeMaxGradientLevel(); level>=lastLevel; level--) {
if (lrs[level] != null) {

Rectangle roi = lrs[level].getRoi();

if (! roi.equals(lastRoi)) {
list. add (lrs [level] . getChainCode ( ) ) ;
lastRoi = roi;
}
}
}
return list;
}

public List getBoundaryListFromMaxGradient() {
ArrayList list = new ArrayList();

for (int level=computeMaxGradientLevel(); level>=lastLevel; level--) {
list.add(lrs [level] .getChainCode () ) ;
}
return list;
}

public String getClusterStats(String id) {
StringBuffer buffer = new StringBuffer();
Rectangle lastRoi = null;
int count = 0;
NumberFormat format = new DecimalFormat(".00");
double rAvg = 0.0;
int avgCount = 0;

for (int level=firstLevel; level>=1astLevel; level--) {
// compute the average gradient avgCount = 0;
rAvg = 0.0;
for (int i=0; i<5; i++) {
int ndx = level-2+i;
if (ndx <= firstLevel && ndx >= lastLevel) {
avgCount++;
rAvg += lrs [ndx] . getGradient ( ) ;
}
}
Rectangle roi = lrs[level].getRoi();
if (I roi.equals(lastRoi)) {
count++;

Stat massStat = lrs[level].getMassStat();
Stat interiorStat = lrs[level].getInteriorStat();
double gradient = lrs[level].getGradient();
double compactness = lrs[level].getCompactness();
double frameIntersection = lrs[level].getFrameIntersection();
buffer.append(id);
buffer.append(DELIMITER).append(count);
buffer.append(DELIMITER).append(level);
buffer.append(DELIMITER).append(roi.x);
buffer.append(DELIMITER).append(roi.y);

buffer.append(DELIMITER).append(roi.width);
buffer.append(DELIMITER).append(roi.height);
buffer.append(DELIMITER).append(massStat.count());
buffer.append(DELIMITER).append((int)massStat.getMax());
buffer.append(DELIMITER);
buffer.append(format.format(massStat.computeMean()));
buffer.append(DELIMITER).append(interiorStat.count());
buffer.append(DELIMITER).append((int)interiorStat.getMax());
buffer.append(DELIMITER);
buffer.append(format.format(interiorStat.computeMean()));
buffer.append(DELIMITER);
buffer.append(format.format(compactness));
buffer.append(DELIMITER);
buffer.append(format.format(gradient));
buffer.append(DELIMITER);
buffer.append(format.format(frameIntersection));
buffer.append("\n");
}
lastRoi = roi;
}

return buffer.toString();
}

public ChainCode getChainCodeAtLevel(int level) {
return lrs[level].getChainCode();
}

private LevelRoi makeLR(int level, Cluster c, boolean saveChain) {
LevelRoi lr = new LevelRoi();

lr.setLevel(level);
try {
lr.setRoi((Rectangle) c.getBounds().clone());
lr.setMassStat((Stat) c.getMassStat()=.clone());
lr.setlnteriorStat((Stat) c.getInteriorStat().clone());
// Due to costly computation time and memory requirements, // we only recompute and save the cluster boundary if we J/ have requested if (saveChain) {
c.updateBoundary();
1r.setChainCode(c.getBoundary());
}
}
catch (CloneNotSupportedException ex) {
System.out.println("Cannot clone this object.");
}

lr.setCompactness(c.computeCompactness().doubleValue());
if (gradientRaster != null) lr.setGradient(clusterMap.computeFrontierGradient(c, gradientRaster));

compute the gradient strength of the boundary; note we do not compute gradient if the given level is below the minGradientLevel if (gradientImage == null 11 level <
clusterConstraints.getMinGradientLevel()) return lr;
else {

double gradientAvg = gradientHash.computeGradient(c.getBounds());
lr.setGradient(gradientAvg);

double d = gradientHash.getFrameIntersection(c.getBounds());
ir.setFrameIntersection(d);
}
return lr;
}

public int computeMaxCompactnessLevel(int start) {
int maxLevel = start;
double maxCompactness = 0.0;

for (int i=start-1; i>=lastLevel; i--) {
if ( (lrs [i] .getArea () > 100) &&
(lrs[i].getCompactness() > maxCompactness) maxLevel = i;
maxCompactness = lrs[i].getCompactness();
}
}
return maxLevel;
}

public int computeMaxCompactnessDropoffLevel() {
int maxGradNdx = computeMaxGradientLevel();
int maxCompNdx = computeMaxCompactnessLevel(maxGradNdx);
double maxComp = lrs[maxCompNdx].getCompactness();

// search for the max drop-off after the max int maxDrop = maxCompNdx;
double maxDropoff = 0.0;
for (int i=maxCompNdx; i-1>=lastLevel; i--) {

double dropoff = lrs [i] .getCompactness () - lrs [i-1].getCompactness();
if (dropoff > maxDropoff) {
maxDropoff = dropoff;
maxDrop }
}

return maxDrop;
}

public Rectangle computeMaxCompactnessDropoff() {
int ndx = computeMaxCompactnessDropoffLevel();

return lrs [ndx] . getRoi ( ) ;
}
private double computeFrameIntersection(Rectangle roi, int width, Raster raster, int threshold) {
int count = 0;
int total = 0;

// horizontal rasters for (int t=0; t<width; t++) {
for (int x=roi.x; x<roi.x+roi.width; x++) {
total++;
if (raster.getSample(x,roi.y+t,0) >= threshold) count++;
total++;
if (raster.getSample(x,roi.y+roi.height-l-t,0) >= threshold) count++;
}
}

vertical rasters for (int t=0; t<width; t++) {
for (int y=roi.y+width; y<roi.y+roi.height-l-width; y++) {
total++;
if (raster.getSample(roi.x+t,y,0) >= threshold) count++;
total++;
if (raster.getSample(roi.x+roi.width-l-t,y,0) >= threshold) count++;
}
}

return (double) count/total;
}

public String getStatHeader() { return STATHEADER; }

public void setProgressBar(JProgressBar b) { progressBar = b; }
public void run() {
byte[] map = computeStepMap();
if (mapListener != null) mapListener.action(new MapEvent(this, mapID, map));
}

public byte [] computeStepMap () {
BufferedImage image = null;
byte[] map = new byte[imageDim.width * imageDim.height];
if ( (debugBuffer != null) && (imageDialog == null) ) {
imageDialog = new ImageDialog(512, 512);
byte[] r = new byte[COLORARRAY.length];

byte[] g = new byte[COLORARRAY.length];
byte[] b = new byte[COLORARRAY.length];
for (int i=0; i<COLORARRAY.length; i++) {
r[i] = (byte) COLORARRAY [i] .getRed () ;
g[i] = (byte) COLORARRAY [i] . getGreen ( ) ;
b[i] = (byte) COLORAR.RAY [i] . getBlue ( ) ;
}
colorModel = new IndexColorModel(8, COLORARRAY.length, r, g, b);
image = new BufferedImage(imageDim.width, imageDim.height, Bufferedlmage.TYPE_BYTE_INDEXED, colorModel);
}

for (int y=0; y<imageDim.height; y++) {
for (int x=0; x<imageDim.width; x++) {
map[x + (y * imageDim.width)] = (byte) 0;
}
}

initialize for (int i=0; i<=MAX_INTENSITY_LEVEL; i++) lrs [i] = null;
firstLevel = -1;

set all points on the map to BACKGROUND (according to the mask) if (mask != null) clusterMap.initialize(mask, clusterMap.BACKGROUND, clusterMap.DEACTIVATED);
else clusterMap.initialize(clusterMap.BACKGROUND);
if (debugBuffer 1= null) {

imageDialog.setImage(clusterMap.createBufferedImage(clusterMap.DEACTIVATED, image, COLORARRAY.length-1));
System.out.println("Deactivated Pixels (mask)");
imageDialog.show();
}

int minStartingLevel =
Math.max(clusterConstraints.getMinStartClusterLevel(), clusterConstraints.getMinGradientLevel());

// monitor the progress if (progressBar 1= null) {
progressBar.setMinimum(0);
progressBar.setMaximum(MAX_INTENSITY_LEVEL - minStartingLevel + 1);
progressBar.setValue(0);
}
BitMap specialMap = new BitMap(imageDim.width, imageDim.height, false);
int levell = MAX INTENSITY LEVEL;

while (levell >= minStartingLevel) {

eliminate pixels from the cluster map that have already been assigned a step value clusterMap.setIdsFromBitMap(specialMap, clusterMap.SPECIAL_1);
// Compute a list of big clusters clusterMap.computeConnectedComponents(levell);
List clusters =
clusterMap.computeSortedListOfClustersBySize(clusterConstraints);
if (debugBuffer != null) debugBuffer.append( LevelSetHelper.computeStepMap() Leve1="+levell+ \n");

if (clusters.size() > 0) {

if (debugBuffer != null) {
imageDialog.setImage(clusterMap.createBufferedImage(clusterMap.SPECIAL_1, image, COLORARRAY.length-1));
System.out.println("Special 1");
imageDialog.show();
}

Cluster cluster = (Cluster) clusters.get(o);
if (debugBuffer != null) {
imageDialog.setImage(clusterMap.createBufferedImage(cluster, image, COLORARRAY.length-1));
System.out.println("Cluster ID="+cluster.getlD()+
" level-1="+levell);
imageDialog.show();
}

//clusterMap.setInteriorPoints(cluster);
clusterMap.addInteriorPoints(cluster);
clusterMap.initializeCluster(cluster);
if (debugBuffer 1= null) {
imageDialog.setImage(clusterMap.createBufferedImage(cluster, image, COLORARRAY.length-1));
System.out.println("Cluster w/ interior points, ID="+cluster.getIDO);
imageDialog.show();
}

lrs[levell] = makeLR(levell, cluster, saveChainCode);
firstLevel = levell;
int maxLevel = levell;
Set pointSet = clusterMap.computePointSet(cluster);

int level2 = levell-1;
boolean done = (level2 <
clusterConstraints.getMinStopClusterLevel());
while (! done) {

if (debugBuffer != null) debugBuffer.append("LevelSetHelper.findAll Leve12=" +
level2+11\n");

clusterMap.addToCluster(cluster, level2);
clusterMap.addInteriorPoints(cluster);
if (debugBuffer != null) {
System.out.println("\tlevel2=11+level2);
imageDialog.setImage(clusterMap.createBufferedImage(cluster, image, COLORARRAY.length-1));
imageDialog.show();
}

lrs[level2] = makeLR(level2, cluster, saveChainCode);

if ( (level2 >= clusterConstraints.getMinGradientLevel()) &&
(cluster.getFilledSize () . intValue () <=
clusterConstraints.getMaxClusterSize()) &&
(lrs [level2] .getGradient () >
lrs[maxLevel].getGradient()) ){
maxLevel = level2;
pointSet = clusterMap.computePointSet(cluster);
}

if (debugBuffer != null) debugBuffer.append("LevelSetHelper.findAll Level=" +
level2+11, size="+
cluster.getFilledSize()+" pixels, hullRatio="+

FORMAT2 . format ( (double) lrs [level2] . getArea ( ) /lrs [maxLevel] .
getArea ( ) ) +
11\n 1l) ;
if (lrs [level2] .getArea() > (4 * lrs [maxLevel] .getArea() ) ) {
done = true;
lastLevel = level2+1;
}
else if (level2 clusterConstraints.getMinStopClusterLevel()) {
done = true;
lastLevel = level2;
}

level2--;
}

update the map int steps = computeStepsFromMaxGradient(debugBuffer);
if (debugBuffer != null) System.out.println("Steps = "+steps);

for (Iterator points=pointSet.iterator(); points.hasNext(); ) {
Point p = (Point) points.next();
map[p.x + (p.y*imageDim.width)] _ (byte) Math.max(map[p.x + (p.y*imageDim.width)], steps);
speci.alMap.setBit(p.x, p.y, true);
}
if (debugBuffer != null) {
WritableRaster raster = image.getRaster();
raster.setDataElements(0,0,imageDim.width, imageDim.height, map);
imageDialog.setImage(image);
System.out.println("display map");
imageDialog.show();
}

if (debugBuffer != null) {
imageDialog.setImage(clusterMap.createBufferedImage(clusterMap.DEACTIVATED, image, COLORARRAY.length-1));
System.out.println("Deactivated Pixels");
imageDialog.show();
}
}
else levell--;
initialize the cluster map, but preserve the deactivated pixels so that we do not include them in future clusters.
if (mask != null) clusterMap.initialize(mask, clusterMap.BACKGROUND, clusterMap.DEACTIVATED);
else clusterMap.initialize(clusterMap.BACKGROUND);
if (progressBar != null) progressBar.setValue(MAX_INTENSITY_LEVEL - levell);
}

return map;
} =

private int max(int i1, int i2, int i3) {
return Math.max(Math.max(i1,i2), Math.max(i2, i3));
}
}

* * * * * * *
// Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.
package com.alanpenn.image;

import java.awt.Point;

import java.awt.Rectangle;
import java.awt.image.*;
import java.awt.geom.*;
import java.io.*;
import java.util.*;
import javax.media.jai.LookupTableJAI;
import javax.vecmath.*;

import com.alanpenn.breast.*;
import com.alanpenn.display3D.*;
import com.alanpenn.image.g4d.*;
import com.alanpenn.parser.CurtainParser;
import com.alanpenn.parser.ParseException;
import com.alanpenn.parser.ParameterParser;
import com.alanpenn.util.*;

/**
* Creation Date: Jan 14, 2005 *
* Maps :
*
* step0l: First version.
*
* step02: LevelSetHelper was altered to allow unconstrained growth of *"seed" clusters, possibly including pixels that were previously * included in the ontogeny of other clusters. In addition, the step * value at a particular pixel was allowed to increase if it was * included in a second or third seed cluster.
*
* step03: Restrict seed clusters additionally by * minSeedPeakIntensity. Also, do not reset previous map entries with * higher step values.
*
* step04: Same as step02 except sets peakOffset to 1; users (peak+l) * time slice.
*
* step05: Same as step02 except sets peakOffset to 2; uses (peak+2) * time slice.
*
public class LevelSetMapper {

private static final String PAR.AM FILE NAME ="/apa/bin/lsv_params.xml";
private static BreastDatabase breastDB;
private LevelSetFinder levelSetFinder;
private IntegerParameter minClusterIntensity;
private IntegerParameter minGradientIntensity;
private IntegerParameter minClusterSize;
private IntegerParameter maxClusterSize;
private IntegerParameter minClusterEnhancement;
private FloatParameter minSeedPeakIntensity;
private Parameters parameters;
private ClusterConstraints clusterConstraints;
private DataBuffer4D buffer;
private MapDataBuffer mapBuffer;
private SagittalSampleModel ssml, ssm2;

private LookupTableJAI 1ut256;
private BufferedImage bil, bi2;
private Raster rasterl, raster2;
private byte [] [] map;

public LevelSetMapper() throws Exception {

breastDB = BreastUtil.loadDatabase();
clusterConstraints = new ClusterConstraints();
File paramFile = new File(PARAM_FILE_NAME);
parameters = ParameterParser.parse(paramFile);
minClusterIntensity = (IntegerParameter) parameters.get("minClusterIntensity");
minGradientIntensity = (IntegerParameter) parameters.get("minGradientIntensity");
minClusterSize = (IntegerParameter) parameters.get("minClusterSize");
maxClusterSize = (IntegerParameter) parameters.get("maxClusterSize");
minClusterEnhancement = (IntegerParameter) parameters.get("minClusterEnhancement");
minSeedPeakIntensity = (FloatParameter) parameters.get("minSeedPeakIntensity");
// set the min cluster intensity level int ivalue = ((Integer) minClusterlntensity.getValue()).intValue();
clusterConstraints.setMinStopClusterLevel(ivalue);
clusterConstraints.setMinStartClusterLevel(ivalue);
// set the min gradient intensity level ivalue = ((Integer) minGradientIntensity.getValue()).intValue();
clusterConstraints.setMinGradientLevel(ivalue);
// set the min cluster enhancement ivalue = ((Integer) minClusterEnhancement.getValue()).intValue();
clusterConstraints.setMinClusterEnhancement(ivalue);
levelSetFinder = new LevelSetFinder();
levelSetFinder.setClusterConstraints(clusterConstraints);
}

public void map(BreastCase breastCase, String mapID, boolean overwrite, int peakOffset) throws IOException {

ImageStore store = breastCase.getlmageStore();
// get all T1-weighted stacks ImageStack stack = null;

List massList = breastCase.getMassList();
if (massList.size() > 0) {
Mass mass = (Mass) massList.get(0);
// Anatomy anatomy = mass.getAnatomy();
stack = breastCase.getImageStack(anatomy, ImageType.Tl);

List stacks = breastCase.getImageStackList(ImageType.Tl);
for (int i=0; i<stacks.size(); i++) {
stack = (ImageStack) stacks.get(i);
Kinetic3p k3p = stack.getKinetic3p();

File file = store.createMapFile(breastCase.getID(), breastCase.getLocalID(), stack.getAnatomy(), mapID);
if (! file.exists() overwrite) {

/J check to make sure there is the appropriate time J/ slice for continuing...
if (k3p.getPeakIndex()+peakOffset <
stack.getNumberOfTimeSlices()) {
MapDataBuffer mapBuffer = computeMapDataBuffer(breastCase, stack, peakOffset);
mapBuffer.setFile(file);
mapBuffer.writeUncompressed();
}
}
else System.out.println(file+" already exists.");
breastCase.releaseImages();
}
}
public void mapAll(String mapID, boolean overwrite, int peakOffset) {
for (Iterator i=breastDB.iterator(); i.hasNext();
BreastCase bc = (BreastCase) i.next();
System.out.println("Case +bc.getID());
try {
map(bc, mapID, overwrite, peakOffset);
System.out.println("\t...mapped case "+bc.getID());
}
catch (Exception ex) {
System.out.println("\t...Problem with case "+bc.getID()+": "+
ex.getMessage());
}
}
}

private MapDataBuffer computeMapDataBuffer(BreastCase breastCase, ImageStack stack, int peakOffset) throws IOException {

String lid = breastCase.getLocalID();
ImageStore store = breastCase.getImageStore();

buffer = stack.readImageData(store, breastCase.getIDO, lid);
Plane plane = new Plane(Plane.SAGITTAL);

switch (buffer.getDataType()) {

case DataBuffer.TYPEBYTE:
ssml = new SagittalSampleModelByte((DataBuffer4DByte)buffer);
ssm2 = new SagittalSampleModelByte((DataBuffer4DByte)buffer);
break;
case DataBuffer.TYPEUSHORT:
ssml = new SagittalSampleModelShort((DataBuffer4DShort)buffer);
ssm2 = new SagittalSampleModelShort((DataBuffer4DShort)buffer);
break;
}
1ut256 = ImageUtil.createLUT(buffer.getMaxValue());
bil = new BufferedImage(ssml.getWidth(), ssml.getHeight(), BufferedImage.TYPEUSHORTGRAY);
bi2 = new BufferedImage(ssml.getWidth(), ssml.getHeight(), BufferedImage.TYPE USHORT_GRAY);

rasterl = Raster.createRaster(ssml, new DataBufferUShort(0), new Point(0,0));
raster2 = Raster.createRaster(ssm2, new DataBufferUShort(0), new Point(0,0));

Tuple4i bufferDim = buffer.getDimensionO;
// scale the images down to 256x256 float cadScaleFactor = (bufferDim.y > 256 ? 256.Of / bufferDim.y 1.Of);

// compute the minimum cluster size (in pixels) Tuple3d voxelSize = buffer.getVoxelSize();
double pixelSize = voxelSize.y * voxelSize.z / Math.pow(cadScaleFactor, 2.0) ;
int minClusterMM = ((Integer)minClusterSize.getValue()).intValue();
int minClusterPixels = (int) Math.round(minClusterMM / pixelSize);
clusterConstraints.setMinClusterSize(minClusterPixels);
int maxClusterMM = ((Integer)maxClusterSize.getValue()).intValue();
int maxClusterPixels = (int) Math.round(maxClusterMM / pixelSize);
clusterConstraints.setMaxClusterSize(maxClusterPixels);

mapBuffer = new MapDataBuffer();
mapBuffer.setDimension(new Point3i(bufferDim.x, (int)(bufferDim.y * cadScaleFactor), (int)(bufferDim.z *
cadScaleFactor)));
mapBuffer.setPlane(Plane.Sagittal);
HashMap parameterMap = new HashMap();
for (Iterator i=parameters.iterator(); i.hasNext(); ) {
Parameter param = (Parameter) i.next();
if (param instanceof IntegerParameter) parameterMap.put(param.getID(), new Double(((IntegerParameter)param).intValue()));
else if (param instanceof FloatParameter) parameterMap.put(param.getID(), new Double(((FloatParameter)param).floatValue()));
}
mapBuffer.setParameters(parameterMap);
map = new byte [bufferDim.x] [] ;

Stat timeStat = new Stat("Running time");

Curtain curtain = null;
try {
curtain = ChestWallUtil.loadChestwallCurtain(breastCase, stack);
}atch (ParseException ex) {
ex.printStackTrace();
}
catch (IOException ex) {
System.out.println("No chestwall found for "+breastCase.getID()+"
Anatomy="+stack.getAnatomy()+"\n");
curtain = null;
}

for (int sliceNumber=0; sliceNumber<buffer.getNumberOfSagittalSlices();
sliceNumber++) {

long startTime = System.currentTimeMillis();
plane.setSliceNumber(sliceNumber);
try {
int maxSliceValue = buffer.computeMaxValue(plane);
ImageUtil.remapLUT(1ut256, 0, maxSliceValue);
ssml.setTimeSlice(0);
ssml.setSpatialSlice(sliceNumber);
bil.setData(rasterl);
RenderedImage preImage = ImageUtil.windowLevel p(bil, lut256);
if (cadScaleFactor != 1) {
preImage = ImageUtil.scale(preImage, cadScaleFactor, cadScaleFactor, 0, 0);
}
Kinetic3p k3p = stack.getKinetic3p();
ssm2.setTimeSlice(k3p.getPeakIndex()+peakOffset);
ssm2.setSpatialSlice(sliceNumber);
bi2.setData(raster2);
RenderedImage peakImage = ImageUtil.windowLevelOp(bi2, 1ut256);
if (cadScaleFactor != 1) peakImage = ImageUtil.scale(peakImage, cadScaleFactor, cadScaleFactor, 0, 0);

/J Stat imageStat = ImageUtil.computelmageStat(peaklmage);
/J int minIntensityLevel = (int) Math.round(imageStat.computeMean() +

(minSeedPeakIntensity.floatValue()*imageStat.computeStd()));
clusterConstraints.setMinSeedPeakIntensity(0);
// mask out the chest wall BitMap bitMask = new BitMap(preImage.getWidth(), preImage.getHeight(), true);
if (curtain != null) {

Line2D line = (Line2D) curtain.getPlaneIntersection(plane);
for (int y=0; y<bitMask.getHeight(); y++) {

for (int x=0; x<bitMask.getwidth(); x++) {
if (line != null) bitMask.setBit(x, y, x<line.getXl()*cadScaleFactor);
else bitMask.setBit(x, y, true);
}
}
}
levelSetFinder.setImages(preImage, peakImage, bitMask);
StringBuffer debugBuffer = null; //new StringBuffer();
map[sliceNumber] = levelSetFinder.computeStepMap(sliceNumber, debugBuffer);

System.out.println("Case "+breastCase.getIDO+
Anatomy="+stack.getAnatomy()+
Slice "+sliceNumber+" of "+buffer.getNumberOfSagittalSlices());

long endTime = System.currentTimeMillis();
timeStat.addData((endTime-startTime)/1000.0);
System.out.println(timeStat);
}
catch (Exception ex) {
ex.printStackTrace();
}
}
mapBuffer.setData(map);
return mapBuffer;
}
public static void main(String[] args) {
try {

LevelSetMapper mapper = new LevelSetMapper();
boolean overwrite = false;

String maplD = "step02";
int peakOffset = 0;

if (args.length == 0) {
mapper.mapAll(maplD, overwrite, peakOffset);
}
else {

// check if first arg is an integer try {
Integer.parseInt(args[0]);
for (int i=0; i<args.length; i++) {
Integer id = new Integer(Integer.parselnt(args[i]));
BreastCase bc = breastDB.getCase(id);
if (bc != null) mapper.map(bc, mapID, overwrite, peakOffset);

}
}
catch (NumberFormatException ex) {

// first arg might be an image store name ImageStore store = breastDB.getImageStore(args[0]);
if (store == null) return;
System.out.println("Processing cases in image store:
"+store);

for (Iterator cases=breastDB.iterator(); cases.hasNext();
{
BreastCase bc = (BreastCase) cases.next();
try {
if (bc.getImageStore().equals(store)) {
mapper.map(bc, maplD, overwrite, peakOffset);
System.gc();
}
}
catch (FileNotFoundException exc) {
System.out.println(exc.getMessage());
}
}
}
}
}
catch (Exception ex) {
ex.printStackTrace();
}
}
}

* * * * * * *
<parameters id="levelsetviewer">
<parameter id="minClusterIntensity" label="MinCluster" min="1" max="10"
step="1" value="3" units="Lev"/>
<parameter id="minGradientIntensity" label="MinGradient" min="1" max="20"
step="1" value="10" units="Lev"/>
<parameter id="minClusterSize" label="MinSize" min="5" max="50" step="5"
value="25" units="mm"/>
<parameter id="maxClusterSize" label="MaxSize" min="10" max="500" step="10"
value="400" units="mm"/>
<parameter id="minClusterEnhancement" label="MinEnhance" min="O" max="500"
step="10" value="80" units=" s"/>
<parameter id="minSeedPeakIntensity" label="MinSeedPeakIntensity" min=".0"
max="3.0" step=".5" value="1.5" units="std"/>
<typein id="coefficientOfDispersion" label="Coef of Dispersion" width="5"
value="0.000" format="0.000" units=""/>
</parameters>

* * * * * * *
Copyright (c) 2005 Alan Penn & Associates, Inc. All rights reserved.

package com.alanpenn.image;

import java.awt.Dimension;
import java.awt.Rectangle;
import java.awt.Point;
import java.awt.GridLayout;
import java.awt.Color;
import java.awt.Graphics2D;
import java.awt.image.BufferedImage;
import java.awt.image.DataBufferByte;
import java.awt.image.WritableRaster;
import java.awt.image.RenderedImage;
import java.awt.image.Raster;
import javax.swing.ImageIcon;
import javax.swing.JLabel;
import javax.swing.JPanel;
import javax.swing.SwingConstants;
/**
* Segments a gray-scale, single band image. The final segmentation * separates the image into three regions:
*
* <ul>
* <li> outside the mass * <li> on the mass margin * <li> inside the mass * </ul>
*
* The runtime behavior of the algorithm is controlled by the * following parameters:
* <ul>
*
*<li> "Outside" is determined by the desired amount of overlap *(intersection) between the mass and the framing edge around the * image (frameSize, frame0verlap).
*
*<li> "On" the mass margin (and subsequently "inside" as well) is * determined by the desired size of the interior * (minimumMassInterior).
*
* </ul>
*
* Some segmentation algorithms further require that we only consider * every other or every third gray level, for example. The default is * to consider all levels. (thresholdMultiple) *
public class SegmentationBuffer {

// this pixel is in the background public static final byte BACKGROUND = 0;
// this pixel is in the foreground public static final byte FOREGROUND = 1;
private static final byte BG = 0;
private static final byte FGO = Ox01;

86 private static final byte FGl = 0x02;

private static final int MAX_IMAGE_SIZE = 512;
public static final int DEFAULT_FRAME_SIZE = 2;

public static final double DEFAULT_FRAME_OVERLAP = 0.15;
public static final double DEFAULT MINIMUM MASS_INTERIOR = 0.60;
private static final int DEBUG = 0;

private ThresholdRange thresholdRange;
private int imageMin;
private int imageMax;
private byte [] [] map;

private RenderedImage image = null;
private ClusterMap clusterMap;
private int frameSize;
private double frameOverlap;
private double minimumMassInterior;

Some segmentation algorithms further require that we only consider every other or every third gray level. The default is to consider all levels.
private int thresholdMultiple = 1;
public SegmentationBuffer() {
frameSize = DEFAULT_FRAME_SIZE;
frameOverlap = DEFAULT_FRAME_OVERLAP;
minimumMassInterior = DEFAULT_MINIMUM_MASS_INTERIOR;
map = new byte[MAX_IMAGE_SIZE][MAX_IMAGE_SIZE];
}
/**
* Set the size of the bounding frame (thickness), and the amount * of overlap required for the outer boundary of the mass.
*
public void setFrameSizeAndOverlap(int size, double overlap) {
frameSize = size;
frameOverlap = overlap;
}

/**
* Set the minimum size of the mass interior, relative to the * overall mass size.

public void setMinimumInterior(double amount) {
this.minimumMassInterior = amount;
}
/**

87 * Restrict the allowable gray levels.

public void setThresholdMultiple(int mult) {
thresholdMultiple = mult;
}

public void setlmage(RenderedImage img) {
this.image = img;

dataBuffer = new DataBufferByte(img.getWidth() * img.getHeight(), BANKS);
initialize the segmentation; label all elements // BACKGROUND
Point p = new Point();
for (p.y=0; p.y<getHeight(); p.y++) {
for (p.x=0; p.x<getWidth(); p.x++) {
//setElem(0, p, BACKGROUND);
map [p ..ac] [p = y] = BG;
}
}

clusterMap = new ClusterMap(img);

imageMin = ImageUtil.getMinimumValue(image);
imageMax = ImageUtil.getMaximumValue(image);
}

/**
* Perform a fast segmentation (mass only) for the given image.
* @returns The byte array containing the segmentation values.

public byte[] fastSegment(RenderedImage img) {
if (img == null) throw new NullPointerException("Image cannot be null.");
setImage(img);
thresholdRange = new ThresholdRange();
thresholdRange.bottom = findBottomThreshold(guessBottomThreshold(0.20));
thresholdRange.top = imageMax + 1;
return segment(thresholdRange);
}

public int findBottomThreshold(int guess) {

// initialize low and high to the minimum and maximum of the image int low = imageMin;
int high = imageMax;
int maxNdx = -1;

// binary search to find the low threshold int count=0;
int threshold = guess;

88 while (low < high) {

if (count > 0) threshold = (int) (((double) low + (double) high) / 2.0);
count++;

if (DEBUG > 0) {
System.out.println("DEBUG(LOW): SettingBottom: low=" + low +
middle=" + threshold +
high=" + high);
}

clusterMap.reset();
clusterMap.setConnectedComponents(threshold);
//clusterMap.dump("Low" + count + " Threshold=" + threshold);
if (clusterMap.getNumberOfClusters() > 0) {
maxNdx = clusterMap.getLargestCluster();
clusterMap.addInteriorPoints(maxNdx);
//clusterMap.dump("Low" + count + " AddIntPoints=" +
//maxNdx);
//if (DISPLAY_CLUSTER) //addClusterImage(maxNdx, "L");
}

if (computePictureFramelntersection(frameSize, maxNdx) >= frameOverlap) low = threshold + 1;
else high = threshold;
}

return low;
}

public ThresholdRange getThresholdRange() { return thresholdRange; }
/**
* Determines the thresholds which satisfy the given segmentation.
*
* @returns the [low, high] threshold range public ThresholdRange findThresholdRange() {
thresholdRange = new ThresholdRange();

initialize low and high to the minimum and maximum of the image imageMin = ImageUtil.getMinimumValue(image);
imageMax = ImageUtil.getMaximumValue(image);
int low = imageMin;
int high = imageMax;
int maxNdx = -1;

binary search to find the low threshold

89 int count=0;

while (low < high) {
count++;
int threshold =(int) (((double) low + (double) high) / 2.0);
if (DEBUG > 0) {
System.out.println("DEBUG(LOW): SettingBottom: low=" + low +
middle=" + threshold +
high=" + high);
}

clusterMap.reset();
clusterMap.setConnectedComponents(threshold);
//clusterMap.dump("Low" + count + " Threshold=" + threshold);
if (clusterMap.getNumberOfClusters() > 0) {
maxNdx = clusterMap.getLargestCluster();
clusterMap.addinteriorPoints(maxNdx);
//clusterMap.dump("Low" + count + " AddIntPoints=" +
//maxNdx);
//if (DISPLAY_CLUSTER) //addClusterImage(maxNdx, "L");
}

if (computePictureFramelntersection(frameSize, maxNdx) >= frameOverlap) low = threshold + 1;
else high = threshold;
}

thresholdRange.bottom = low;
clusterMap.reset();
clusterMap.setConnectedComponents(thresholdRange.bottom);
int bottomPopulation = 0;

if (clusterMap.getNumberOfClusters() > 0) {
maxNdx = clusterMap.getLargestCluster();
clusterMap.addInteriorPoints(maxNdx);
bottomPopulation = clusterMap.getPopulation(maxNdx);
}
else {
System.err.println("ERROR: Number of clusters +
clusterMap.getNumberOfClusters());
return thresholdRange;
}
if (bottomPopulation < 1) {
System.err.println("ERROR: bottomPopulation=" +
bottomPopulation);
return thresholdRange;
}

low = thresholdRange.bottom;
high = imageMax;
count=0;

double remainingPopulation = 0.0;
while (low < high) {

count++;
int threshold = (int) Math.round( ((double)low+high) / 2.0);
// require "threshold" be a multiple of "thresholdMultiple"
while ( ((threshold - thresholdRange.bottom) % thresholdMultiple) 1=
0) threshold--;
if (DEBUG > 0) {
System.out.println("DEBUG(HIGH): SettingBottom: low=" + low +
middle=" + threshold +
high=" + high);
}

clusterMap.reset();
clusterMap.setConnectedComponents(threshold);
int pop = 0;
if (clusterMap.getNumberOfClusters() > 0) {
maxNdx = clusterMap.getLargestCluster();
clusterMap.addInteriorPoints(maxNdx);
//if (DISPLAY_CLUSTER) //addClusterImage(maxNdx, "H");
pop = clusterMap.getPopulation(maxNdx);
}

remainingPopulation = (double)pop/bottomPopulation;
if (DEBUG > 0) {
System.out.println("\tRemaining Pop=" + remainingPopulation);
}

if ( (remainingPopulation <= minimumMassInterior) (low == threshold) ) high = threshold - 1;
else low = threshold;
}

thresholdRange.top = low;
if (DEBUG > 0) System.out.println("Final range=" + thresholdRange.bottom + "," +
thresholdRange.top);
return thresholdRange;
}

/**
* Perform the segmentation for the given image.

*@returns The byte array containing the segmentation values.

public byte [] segment () {
if (image == null) throw new NullPointerException("Image cannot be null.");
return segment(findThresholdRange());
}
/**
* Perform the segmentation for the given image.
*@returns The byte array containing the segmentation values.
public byte[] segment(RenderedImage img) {

if (img == null) throw new NullPointerException("Image cannot be null.");
setImage(img);
return segment(findThresholdRange());
}

/**
* Perform the segmentation for the given region of interest * within the image.
*@returns The byte array containing the segmentation values.
public byte[] segment(RenderedImage img, Rectangle roi) {

if (img == null) throw new NullPointerException("Image cannot be null.");
RenderedImage crop = ImageUtil.crop(img, roi);
setImage(crop);
return segment(findThresholdRange());
}

/**
* Perform the segmentation for the given image.
*@returns The byte array containing the segmentation values.
public byte[] segment(ThresholdRange r) {

if (image == null) throw new NullPointerException("Image cannot be null.");
clusterMap.reset();
clusterMap.setConnectedComponents(r.bottom);
int maxNdx = -1;
if (clusterMap.getNumberOfClusters() > 0) {
maxNdx = clusterMap.getLargestCluster();
clusterMap.addInteriorPoints(maxNdx);
}

int bottomCounter = 0;

// Set the first data buffer Point p = new Point();

for (p.y=O; p.y<getHeight(); p.y++) {
for (p.x=0; p.x<getWidth(); p.x++) {

if (clusterMap.getID(p) == maxNdx) {
//setElem(0, p, FOREGROUND);
map [p. x] [p. y] = FGO;
bottomCounter++;
}
else {
map [p. x] [p - y] = BG;
//setElem(0, p, BACKGROUND);
}
}
clusterMap.reset();
clusterMap.setConnectedComponents(r.top);
maxNdx = clusterMap.getLargestCluster();
clusterMap.addInteriorPoints(maxNdx);
int topCounter = 0;

// set the second data buffer for (p.y=0; p.y<getHeight(); p.y++) {
for (p.x=0; p.x<getWidth(); p.x++) {

if (clusterMap.getID(p) == maxNdx) {
//setElem(l, p, FOREGROUND);
map [p.x] [p-y] 1= FG1;
topCounter++;
}
else setElem(l, p, BACKGROUND);
}
}
if (DEBUG > 0) {
System.out.println("Bottom=" + bottomCounter +
" Top=" + topCounter +
" Top/Bottom=" +
(double)topCounter/bottomCounter);
return makeDataBuffer().getData();
}
/**
* The area of the "picture frame" is the size of a *<code>k</code>-pixel border that circumscribes the image (which * looks like a,,picture frame).

private int computePictureFrameArea(int k) {
return ((getWidth() * k * 2) +
(getHeight() * k * 2) -(4 * k * k) ) ;
}

/**
* Counts the number of pixels with connected component ID=id that * also intersect a k-wide frame around the border of the image.
* @param k Frame thickness * @param id id of the connected component of interest private int computeRawPictureFrameIntersection(int k, int id) {
int,count = 0;

Point p = new Point();
// horizontal rasters for (int t=0; t<k; t++) {
for (p.x=0; p.x<getWidth(); p.x++) {
P=y = t;
if (clusterMap.getID(p) id) count++;
p.y = getHeight O -1-t;
if (clusterMap.getID(p) id) count++;
}
}

J/ vertical rasters for (int t=0; t<k; t++) {
for (p.y=k; p.y<getHeight()-1-k; p.y++) {
p=x=t;
if (clusterMap.getID(p) id) count++;
p.x = getWidth()-1-t;
if (clusterMap.getID(p) id) count++;
}
}

return count;
}

/**
* Computers the percentage of the pixels with connected component * ID=id that also intersect a k-wide frame around the border of * the image.
* @param k frame thickness *@param id id of the connected component of interest private double computePictureFrameIntersection(int k, int id) {
double rval = 0.0;
int a=0, b=0;
if (k > 0) {
a = computeRawPictureFrameIntersection(k, id);
b = computePictureFrameArea(k);

rval = (double) a / b;
}
else {
System.err.println("ERROR: picture frame thickness +
k) ;
}

if (DEBUG > 0) {
System.out.println("computePictureFramelntersection(" +
k + ", " + id + ") a= " + a +
b=" + b + " rval=" + rval);
}

return rval;
}

private Point FourNeighbor(Point p, int k) {
Point r=null;

switch (k) {
case 0:
r = new Point(p.x+l, p.y); break;
case 1:
r = new Point(p.x, p.y-1); break;
case 2:
r = new Point(p.x-1, p.y); break;
case 3:
r = new Point(p.x, p.y+l); break;
default:
throw new IllegalArgumentException("Bad neighbor index:" + k);
}
return r;
}

private Point EightNeighbor(Point p, int k) {
Point r=null;

switch (k) {
case 0:
r = new Point(p.x+l, p.y); break;
case 1:
r = new Point(p.x+l, p.y-1); break;
case 2:
r = new Point(p.x, p.y-1); break;
case 3:
r = new Point(p.x-1, p.y-1); break;
case 4:
r = new Point(p.x-1, p.y); break;
case 5:
r= new Point(p.x-1, p.y+l); break;
case 6:
r = new Point(p.x, p.y+l); break;
case 7:
r = new Point(p.x+1, p.y+l); break;
default:
throw new IllegalArgumentException("Bad neighbor index:" + k);

}
return r;
}

private boolean on4Boundary(int b, Point p) {
if ( (map [p.x] [p.y] & b) 0) return false;

for (int k=0; k<8; k++) {
Point n = EightNeighbor(p, k);

if (n.x >= 0 && n.x < getWidth() &&
n.y >= 0 && n.y < getHeight()) {
if ( (map [n.x] [n.y] & b) == 0) return true;
}
else {
return true;
}
}
return false;
}

private boolean onBBoundary(byte b, Point p) {
if ( (map [p.x] [p.y] & b) 0) return false;

for (int k=0; k<4; k++) {
Point n = FourNeighbor(p, k);

if (n.x >= 0 && n.x < getWidth() &&
n.y >= 0 && n.y < getHeight()) {
if ( (map[n.x] [n.y] & b) 0) return true;
}
{
else return true;
}
}
return false;
}

~**
* Classify the pixels into the following categories:
* 1. mass interior * 2. mass margin * 3. mass exterior *
* We made a change as of June 25, 2003 to not include the border * of the interior in the margin.
*
* Well, we decided two days later to put it back in. Seems our * statistics were not as good with it left out.

private DataBufferByte makeDataBuffer() {

DataBufferByte buffer = new DataBufferByte(getWidth() * getHeight(), 1);
Point p = new PointO;
for (p.y=0; p.y<getHeight(); p.y++) {
for (p.x=0; p.x<getWidth(); p.x++) {
byte b = Segmentation.ROI;
if ( (map [p.x] [p.y] & FGO) > 0) {
if ( (map[p=x] [p=y] & FG1) 0) b 1= Segmentation.MASSMARGIN;
else if (on8Boundary(FG1, p)) b Segmentation.MASS_MARGIN;
else b Segmentation.MASS_INTERIOR;
}

buffer.setElem(0, (p.y*getWidth() + p.x), b) ;
}
}
return buffer;
}

private int getWidth() {
return image.getWidth();
}

private .int getHeight() {
return image.getHeight();
}

/**
* Computes a guess for the bottom threshold.
private int guessBottomThreshold(double fraction) {

int threshold = imageMax;
while ( (threshold>=imageMin) &&
(computeFrameIntersectionWithImage(threshold) < fraction) threshold--;

return threshold;
}

private double computeFrameIntersectionWithImage(int threshold) {
int count = 0;
int total = 0;

Raster raster = image.getData();
// horizontal rasters int x=0, y=0, t=0;
try {

for (t=0; t<Math.min(image.getHeight(),frameS=ize); t++) {
for (x=0; x<image.getWidth(); x++) {
total++;
if (raster.getSample(x,t,0) >= threshold) count++;
total++;
if (raster.getSample(x,image.getHeight()-1-t,0) >=
threshold) count++;
}
}
// vertical rasters for (t=0; t<Math.min(image.getWidth(),frameSize); t++) {
for (y=frameSize; y<image.getHeight()-1-frameSize; y++) {
total++;
if (raster.getSample(t,y,0) >= threshold) count++;
total++;
if (raster.getSample(image.getWidth()-1-t,y,0) >= threshold) count++;
}
}
}
catch (Exception ex) {
System.out.println("(x,y,t) + x + + y + "," 'f' t +
System.out.println("Image dim=" + image.getWidth() + "x" +
image.getHeight());
ex.printStackTrace();
}

return (double) count/total;
}
}

Claims (75)

What Is Claimed Is:

1. A method of evaluating images of a body region, wherein at least one of the images shows an abnormality indicative of a possible lesion, comprising the steps of:

determining locations of pixels in each image that show such abnormality;
for each of a set of intensity levels, I, determining a contour around the cluster containing the pixels at the locations determined in the step of determining;
defining a function, F, that discriminates a distinct characteristic of each contour in a nested sequence;
defining a function, G, used to characterize the changes in the value of the function F over the range of contours at each intensity level; and identifying a lesion as being more likely to be benign or more likely to be cancerous based on at least one threshold value for the function G or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

2. The method according to claim 1, wherein the step of determining a contour comprises the step of defining the cluster using intensity threshold levels around a landmark.

3. The method according to claim 1, wherein the step of defining a function F comprises the step of defining a function F(C) =Area (B), where B is a geometric contour formed around the cluster.

4. The method according to claim 3, wherein the geometric contour is a rectangular box in 2 or 3 dimensions.

5. The method according to claim 1, wherein the at least one threshold is set at 9 to discriminate benign from malignant lesions.

6. The method according to claim 1, wherein the function, G, reflects the number of changes in the value of the function F over the range of contours at each intensity level.

7. A method of evaluating whether a lesion in an image is cancerous, benign or uncertain, the image comprising a plurality of pixels, each pixel having a particular intensity I in the range of 0 <= I< 2N, where N is an integer > 1, comprising the steps of:

defining a landmark pixel within the lesion in the image;
growing a cluster around the landmark pixel which contains the lesion for each intensity value in a set of possible intensity values in the range 0 to 2N - 1;

at each intensity levels in the set, constructing a region of interest such that the region of interest is a minimal polygon containing the cluster at that intensity level;

computing a value of a characteristic of the minimal polygon at each of the intensity levels in the set;
determining a number relating to changes in the characteristic values over the range of intensity levels; and determining whether the lesion is more likely to be cancerous, benign or uncertain in response to the number related to changes in the characteristic values over the range of intensity levels, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

8. The method according to claim 7, further comprising the step of outputting the changes in the characteristic values, the threshold values and locations of pixels depicting a plurality of lesions within the body region determined in the step of determining, the number of intensity levels within the lesion and an indication of the determination of whether the lesion is more likely to be cancerous, benign or uncertain.

9. The method according to claim 7, wherein the minimal polygon is in 2-dimensions, and the characteristic computed in the step of computing is the area of the minimal polygon.

10. The method according to claim 7, wherein the minimal polygon is in 3-dimensions, and the characteristic computed in the step of computing is the volume of the minimal polygon.

11. The method according to claim 7, wherein 2N = 256.

12. The method according to claim 9, wherein the minimal polygon is a rectangle.

13. The method according to claim10, wherein the minimal polygon is a rectangular box.

14. A method of evaluating images, wherein at least one of the images contains, a lesion, over time comprising the steps of:
gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;
subtracting the pixel intensities of the pre-image from the pixel intensities of the post-image to obtain pixel intensities of a subtraction image;
selecting a region of interest that surrounds the lesion;
setting a landmark within the region of interest, where the landmark may consist of one or more pixels;
starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set {C N, C N-1, ...} , with Function (C N) <=Function (C N-1) <=... , where Function depends in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function (C II) >=Q, where Q is a predetermined value;

determining I min and I max;
starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function (B Imax) <=Function (B Imax-1) <= ...;
incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;
outputting a total number of steps when B(J)>B(previous bounding box); and determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

15. The method according to claim 14, wherein the post-images in one time frame are selected from a plurality of post-images such that the post-images in one time frame correspond to peak enhancement.

16. The method according to claim 14, further comprising setting parameters Q and N such that Q=25 mm2 and N=4.

17. The method according to claim 14, wherein the step of selecting a region of interest comprises inputting a region of interest set by a user through an input device.

18. The method according to claim 14, wherein the step of selecting a region of interest comprises determining the region of interest automatically from an image processing and/or histogram analysis.

19. The method according to claim 10, wherein when 2 dimensions are reflected in the image, M is selected from the group consisting of 4 or 8, and the Function is Area.

20. The method according to claim 10, wherein when 3 dimensions are reflected in the image, M is selected from the group consisting of 6 or 18 or 26, and the Function is Volume.

21. The method according to claim 10, wherein I min and I max are determined as follows:

starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J
(MeanGrad(J)) using data from post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max; and setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function (C Imin) > N*Function (C Imax).

22. The method according to claim 10, wherein I min and I max are determined using a histogram analysis.

23. Software for use with a computer having a memory, an input device for generating images, and a display, the software comprising a computer usable medium having computer readable program code thereon for evaluating images, wherein at least one of the images contains a lesion, comprising:

first program logic for determining locations of pixels in the lesion;

second program logic for, for each of a set of intensity levels, I, determining a contour around the cluster containing the pixels identified in the step of determining;
third program logic for defining a function, F, that discriminates a distinct characteristic of each contour in a nested sequence;

fourth program logic for defining a function, G, used to characterize the changes in the value of the function F over the range of contours at each intensity level; and fifth program logic for identifying a feature responsive to the function G, as being more likely to be benign or more likely to be cancerous based on at least one threshold value for the function G, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

24. The software according to claim 23, wherein the second program logic comprises sixth program logic for defining the cluster using intensity threshold levels around a landmark.

25. The software according to claim 23, wherein the third program logic comprises seventh program logic for defining a function F(C) =Area (B), where B is a geometric contour formed around the cluster.

26. The software according to claim 25, wherein the geometric contour is a rectangular box in 2 or 3 dimensions.

27. The software according to claim 23, wherein the at least one threshold is set at 9 to discriminate benign from malignant lesions.

28. The software according to claim 23, wherein the function, G, reflects the number of changes in the value of the function F over the range of contours at each intensity level.

29. Software for use with a computer having a memory, an input device for generating images, and a display, the software comprising a computer usable medium having computer readable program code thereon for evaluating whether a lesion in an image is cancerous, benign or uncertain, the image comprising a plurality of pixels, each pixel having a particular intensity I in the range of 0 I < 2N, where N is an integer > 1, comprising:

first program logic for defining a landmark pixel within the lesion in the image;
second program logic for growing a cluster around the landmark pixel which contains the lesion for each intensity value in a set of possible intensity values in the range 0 to 2N - 1;
third program logic for, at each intensity level in the set, constructing a region of interest such that the region of interest is a minimal polygon containing the cluster at that intensity level;
fourth program logic for computing a value of a characteristic of the minimal polygon at each of the set of intensity levels in the set;
fifth program logic for determining a number related to changes in the characteristic values over the range of intensity levels; and sixth program logic for determining whether the lesion is more likely to be cancerous, benign or uncertain in response to the number related to changes in the characteristic values over the range of intensity values, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

30 The software according to claim 29, further comprising seventh program logic for outputting the function based on changes in the value of the area of the regions of interest, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region, determined in the step of determining, the number of intensity levels within the lesion and an indication of the determination of whether the lesion is more likely to be cancerous, benign or uncertain.

105 1. The software according to claim 29, wherein the minimal polygon is in 2-dimensions, and the characteristic computed in the step of computing is the area of the minimal polygon.

32. The software according to claim 29, wherein the minimal polygon is in 3-dimensions, and the characteristic computed in the step of computing is the volume of the minimal polygon.

33. The software according to claim 29, wherein 2N = 256.

34. The software according to claim 31, wherein the minimal polygon is a rectangle.

35. The software according to claim 32, wherein the minimal polygon is a rectangular box.

36. Software for use with a computer having a memory, an input device for generating images, and a display, the software comprising a computer usable medium having computer readable program code thereon for evaluating images, wherein at least one of the images contains a lesion, of at least one lesion over time comprising the steps of:
first program logic for gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;
second program logic for subtracting the pixel intensities of the pre-image from the pixel intensities of the post-image to obtain pixel intensities of a subtraction image;

third program logic for selecting a region of interest that surrounds the at least one lesion;
fourth program logic for setting a landmark within the region of interest;
fifth program logic for starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set {C N, C N-1, ...}, with Function (C N) <=Function (C N-1) <=... , where the Function depends at least in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function(C II)>=Q, where Q is a predetermined value;

sixth program logic for determining I max;
seventh program logic for determining I min;

eighth program logic for starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function (B Imax) <=Function (B Imax-1) <=...;

ninth program logic for incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;

tenth program logic for outputting a total number of steps when B(J)>B(previous bounding box); and eleventh program logic for determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

37. The software according to claim 36, wherein the post-images in one time frame are selected from a plurality of post-images such that the post-images in one time frame correspond to peak enhancement.

38. The software according to claim 36, further comprising twelfth program logic for setting parameters Q and N such that Q=25 mm2 and N=4.

39. The software according to claim 36, wherein the first program logic comprises thirteenth program logic for inputting a region of interest set by a user through an input device.

40. The software according to claim 36, wherein the first program logic comprises fourteenth program logic for determining the region of interest automatically from an image processing and/or histogram analysis.

41. The software according to claim 36, wherein when 2 dimensions are reflected in the image, M is selected from the group consisting of 4 or 8, and the Function is Area.

42. The software according to claim 36, wherein when 3 dimensions are reflected in the image, M is selected from the group consisting of 6 or 18 or 26, and the Function is Volume.

43. The software according to claim 36, wherein the sixth and seventh program logic comprise program logic for determining I min and I max as follows:

fifteenth program logic for starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J (MeanGrad(J)) using data from post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max; and sixteenth program logic for setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function(C Imin) >
N*Function(C Imax).

44. The software according to claim 36, wherein the sixth and seventh program logic comprise fifteenth program logic for determining I min and I max using histogram analysis.

45. A computer program product for performing a method of evaluating images, wherein at least one of the images contains a lesion, the computer program product embodied on a computer readable memory and comprising a series of instructions for:

determining locations of pixels in the lesion;
for each of a set of intensity levels, I, determining a contour around the cluster containing the pixels identified in the step of determining;

defining a function, F, that discriminates a distinct characteristic of each contour in a nested sequence;

defining a function, G, used to characterize the changes in the value of the function F over the range of contours at each intensity level; and identifying a feature responsive to the function G, as being more likely to be benign or more likely to be cancerous based on at least one threshold value for the function G, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

46. The computer program product according to claim 45, wherein the instructions for determining a contour comprises instructions for defining the cluster using intensity threshold levels around a landmark.

47. The computer program product according to claim 45, wherein the instructions for defining a function F comprises instructions for defining a function F(C) =Area (B), where B
is a geometric contour formed around the cluster.

48. The computer program product according to claim 45, wherein the geometric contour is a rectangular box in 2 or 3 dimensions.

49. The computer program product according to claim 45, wherein the at least one threshold is set at 9 to discriminate benign from malignant lesions.

50. The computer program product according to claim 45, wherein the function, G, reflects the number of changes in the value of the function F over the range of contours at each intensity level.

51. A computer program product for performing a method of evaluating whether a lesion in an image is cancerous, benign or uncertain, the image comprising a plurality of pixels, each pixel having a particular intensity I in the range of 0 <= I <2N, where N is an integer > 1, the computer program product embodied on a computer readable memory and comprising a series of instructions for:

defining a landmark pixel within the lesion in the image;
growing a cluster around the landmark pixel which contains the lesion for each intensity value in a set of possible intensity values in the range 0 to 2N - 1;
at each intensity level in the set, constructing a region of interest such that the region of interest is a minimal polygon containing the cluster at that intensity level;

computing a value of a characteristic of the minimal polygon at each of the intensity levels in the set;
determining a number related to changes in the characteristic values over the range of intensity levels; and determining whether the lesion is more likely to be cancerous, benign or uncertain in response to the number related to changes in the characteristic values over the range of intensity levels, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

52. The computer program product according to claim 51, further comprising instructions for outputting the function based on changes in the value of the area of the regions of interest determined by the instructions for determining, the number of intensity levels within the lesion and an indication of the determination of whether the lesion is more likely to be cancerous, benign or uncertain.

53. The computer program product according to claim 51, wherein the minimal polygon is in 2-dimensions, and the characteristic computed in the step of computing is the area of the minimal polygon.

54. The computer according to claim 51, wherein the minimal polygon is in 3-dimensions, and the characteristic computed in the step of computing is the volume of the minimal polygon.

55. The computer according to claim 51, wherein 2N = 256.

56. The method according to claim 53, wherein the minimal polygon is a rectangle.

57. The computer according to claim 54, wherein the minimal polygon is a rectangular box.

58. A computer program product for performing a method of evaluating images, wherein at least one of the images contains a lesion, over time, the computer program product embodied on a computer readable memory and comprising a series of instructions for:

gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;

subtracting the pixel intensities of the pre-image from the pixel intensities of the post-images in one time frame to obtain pixel intensities of a subtraction image;
selecting a region of interest that surrounds the lesion;
setting a landmark within the region of interest;
starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set {C N, C N-1, ...} , with Function(C N)<=Function(C N-1) <=... , where the Function depends at least in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function(C II)>=Q, where Q is predetermined value;
determining I max;

determining I min;
starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function(B Imax)<=Function(B Imax-1)<=...;
incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;

outputting a total number of steps when B(J)>B(previous bounding box); and determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

59. The computer program product according to claim 58, wherein the post-images in one time frame are selected from a plurality of post-images such that the post-images in one time frame correspond to peak enhancement.

60. The computer program product according to claim 58, further comprising instructions for setting parameters Q and N
such that Q=25 mm2 and N=4.

61. The computer program product according to claim 58, wherein the instructions for selecting a region of interest comprise instructions for inputting a region of interest set by a user through an input device.

62. The computer program product according to claim 58, wherein the instructions for selecting a region of interest comprise instructions for determining the region of interest automatically from an image processing and/or histogram analysis.

63. The computer program product according to claim 58, wherein when 2 dimensions are reflected in the image, M is selected from the group consisting of 4 or 8, and the Function is Area.

64. The computer program product according to claim 58, wherein when 3 dimensions are reflected in the image, M is selected from the group consisting of 6 or 18 or 26, and the Function is Volume.

65. The computer program product according to claim 58, wherein the steps of determining I max and I min comprise:
starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J
(MeanGrad(J)) using data from the post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max; and setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function(C Imin) > N*Function(C Imax).

66. The computer program product according to claim 58, wherein the steps of determining I max and I min use a histogram analysis.

67. A method of evaluating images, wherein at least one of the images contains a lesion, over time comprising the steps of:

gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;

subtracting the pixel intensities of the pre-image from the pixel intensities of the post-images in one time frame to obtain pixel intensities of a subtraction image;
selecting a region of interest that surrounds the lesion;
setting a landmark within the region of interest;
starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set {C N, C N-1, ...}, with Function(C N)<=Function(C N-1) <=... , where Function depends in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function(C II)>=Q, where Q is a predetermined value;
starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J
(MeanGrad(J)) using data from the post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max;
setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function(C Imin) > N*Function(C Imax);
starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function(B Imax)<=Function(B Imax-1)<=..;

incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;
outputting a total number of steps when B(J)>B(previous bounding box); and determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

68. Software for use with a computer having a memory, an input device for generating images, and a display, the software comprising a computer usable medium having computer readable program code thereon for evaluating images, wherein at least one of the images contains a lesion, over time comprising the steps of:
first program logic for gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;
second program logic for subtracting the pixel intensities of the pre-image from the pixel intensities of the post-images in one time frame to obtain pixel intensities of a subtraction image;
third program logic for selecting a region of interest that surrounds the lesion;
fourth program logic for setting a landmark within the region of interest;
fifth program logic for starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set (C N, C N-1, ...} , with Function(C N)<=Function(C N-1) <=... , where the Function depends at least in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function(C II)>=Q, where Q is a predetermined value;
sixth program logic for starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J (MeanGrad(J)) using data from the post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max;
seventh program logic for setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function(C Imin) >
N*Function(C Imax);

eighth program logic for starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function(B Imax)<=Function(B Imax-1)<=...;
ninth program logic for incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;

tenth program logic for outputting a total number of steps when B(J)>B(previous bounding box); and eleventh program logic for determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

69. A computer program product for performing a method of evaluating images, wherein at least one of the images contains a lesion, over time, the computer program product embodied on a computer readable memory and comprising a series of instructions for:

gathering pixel intensities of post-images in one time frame derived from the signal captured after use of a contrast agent and pixel intensities of a pre-image derived from the signal captured before use of the contrast agent;

subtracting the pixel intensities of the pre-image from the pixel intensities of the post-images in one time frame to obtain pixel intensities of a subtraction image;

selecting a region of interest that surrounds the lesion;
setting a landmark within the region of interest;
starting with level I=I(L) and decrementing I at each step, constructing the cluster of pixels that are M-connected to L, where M is an integer multiple of 2, depending on the number of dimensions reflected in the image, and have intensity level >=I, where I(L) denotes the pixel intensity at the landmark such that the clusters form a monotonically increasing set {C N, C N-1, ...} , with Function (C N)<=Function(C N-1) <=... , where the Function depends at least in part on the number of dimensions reflected in the image, and continuing this process until intensity level equals II, where Function(C II)>=Q, where Q is predetermined value;

starting with level J=II and decrementing by J at each step, computing a mean gradient on the border of C J
(MeanGrad(J)) using data from the post-images in one time frame where the intensity level at which MeanGrad is maximum defines level I max;

setting I min as a highest intensity level for which the area of C Imin exceeds the area of C Imax multiplied by some pre-determined value, Function(C Imin) > N*Function(C Imax);

starting with level I=I max and decrementing through I=I min, constructing a minimum bounding box B I around C I and computing the Functions of the boxes, B Imax ~ B Imax-1~ ..., with Function(B Imax)<=Function(B Imax-1)<=...;
incrementing a step counter each time that a decrement in intensity level induces a change in the Function characteristic of the minimum bounding box;

outputting a total number of steps when B(J)>B(previous bounding box); and determining whether the lesion is more likely to be benign, cancerous or uncertain based on the total number of steps, or based on threshold values and locations of pixels depicting a plurality of lesions within the body region.

70. The method of evaluating images, wherein at least one of the images contains a lesion, as shown and/or described.

71. Software for use with a computer having a memory, an input device for generating images, and a display, the software comprising a computer usable medium having computer readable program code thereon for evaluating images, wherein at least one of the images contains a lesion, as shown and/or described.

72. A computer program product for performing a method of evaluating images, wherein at least one of the images contains a lesion, the computer program product embodied on a computer readable memory and comprising a series of instructions, as shown and/or described.

73. A method for evaluating images containing one or more lesions in a portion of a patient's body, wherein an image consists of a matrix of pixels, each pixel having an intensity indicative of the presence or absence of a lesion, said method comprising the steps of:

identifying a plurality of regions of interest in the image; and for each region of interest, performing the method of claim 1.

74. A method of evaluating whether a plurality of lesions in an image are cancerous, benign or uncertain, the image comprising a plurality of pixels, each pixel having a particular intensity I in the range of 0 <= I< 2N, where N is an integer > 1, said method comprising the steps of:

Identifying a plurality of regions of interest in the image; and for each region of interest, performing the method of claim 8.

75. The software of claim 36 for evaluating images of a plurality of lesions, wherein said third program logic selects a plurality of regions of interest, each surrounding a respective lesion.