A method for detecting malignancy-associated changes. A sample of cells is obtained and stained to identify the nuclear DNA material. The sample is imaged with a digital microscope. Objects of interest are identified in the sample of cells based on the intensity of the pixels that comprise the object versus the average intensity of all pixels in the slide image. An exact edge is located for each object and variations in the illumination intensity of the microscope are compensated for. A computer system calculates feature values for each object and, based on the value of the features, a determination is made whether the cell exhibits malignancy-associated changes or not.
|
5. A method of predicting whether a patient will develop cancer, comprising the steps of:
obtaining a sample of apparently normal cells from a patient; determining whether the cells in the sample exhibit malignancy associated changes by: (1) staining the cells in the sample; (2) obtaining an image of the cells with a digital microscope and recording the image in a computer system; (3) analyzing the stored image of the cells to identify the nuclei; and (4) computing a set of feature values for each nucleus found in a subsample of the sample, and from the feature values determining whether the nucleus exhibits a malignancy associated change; determining a composite value based on the nuclei determined to exhibit malignancy-associated changes and the total nuclei in the subsample; comparing the composite value to a predetermined threshold; and predicting whether the patient will develop cancer if the composite value exceeds the predetermined threshold.
1. A method of predicting whether a patient will develop cancer, comprising the steps of:
obtaining a sample of apparently normal cells from the patient; determining whether the cells in the sample exhibit malignancy associated changes by: (1) staining the nuclei of the cells in the sample; (2) obtaining an image of the cells with a digital microscope and recording the image in a computer system; (3) analyzing the stored image of the cells to identify the nuclei; (4) computing a set of feature values for each nucleus found in the sample and from the feature values determining whether the nucleus exhibits a malignancy associated change; and determining a total number of nuclei in the sample that exhibit malignancy-associated changes and from the number predicting whether the patient will develop cancer, wherein the step of determining a total number of nuclei comprises: determining a ratio of nuclei determined to exhibit malignancy-associated changes to the total identified nuclei; comparing the ratio to a predetermined threshold; and predicting that a patient will develop cancer if the ratio exceeds the predetermined threshold. 2. The method of
dilating or contracting the true edge of the object before calculating the features.
3. The method of
6. The method of
|
The present application is a continuation of application Ser. No. 08/644,893, filed May 10, 1996, now U.S. Pat. No. 5,889,881, which is a continuation-in-part of application Ser. No. 08/425,257 filed Apr. 17, 1995, now abandoned which is a continuation of Ser. No. 08/182,453 filed Jan. 10, 1994 now abandoned which is a continuation-in-part of Ser. No. 07/961,596 filed Oct. 14, 1992, now abandoned, the disclosure of which are expressly incorporated herein by reference. The benefit of the priority of the filing dates of the above-identified applications is hereby claimed under 35 U.S.C. §120.
The present invention relates to image cytometry systems in general, and in particular to automated systems for detecting malignancy-associated changes in cell nuclei.
The most common method of diagnosing cancer in patients is by obtaining a sample of the suspect tissue and examining it under a microscope for the presence of obviously malignant cells. While this process is relatively easy when the location of the suspect tissue is known, it is not so easy when there is no readily identifiable tumor or precancerous lesion. For example, to detect the presence of lung cancer from a sputum sample requires one or more relatively rare cancer cells to be present in the sample. Therefore patients having lung cancer may not be diagnosed properly if the sample does not accurately reflect the conditions of the lung.
Malignancy-associated changes (MACs) are subtle changes that are known to take place in the nuclei of apparently normal cells found near cancer tissue. In addition, MACs have been detected in tissue found near precancerous lesions. Because the cells exhibiting MACs are more numerous than the malignant cells, MACs offer an additional way of diagnosing the presence of cancer, especially in cases where no cancerous cells can be located.
Despite the ability of researchers to detect MACs in patients known to have cancer or a precancerous condition, MACs have not yet achieved wide acceptance as a screening tool to determine whether a patient has or will develop cancer. Traditionally, MACs have been detected by carefully selecting a cell sample from a location near a tumor or precancerous lesion and viewing the cells under relatively high magnification. However, it is believed that the malignancy-associated changes that take place in the cells are too subtle to be reliably detected by a human pathologist working with conventional microscopic equipment, especially when the pathologist does not know beforehand if the patient has cancer or not. For example, a malignancy-associated change may be indicated by the distribution of DNA within the nucleus coupled with slight variations in the shape of the nucleus edge. However, nuclei from normal cells may exhibit similar types of changes but not to the degree that would signify a MAC. Because human operators cannot easily quantify such subtle cell changes, it is difficult to determine which cells exhibit MACs. Furthermore, the changes which indicate a MAC may vary between different types of cancer, thereby increasing the difficulty of detecting them.
The present invention is a system for automatically detecting malignancy-associated changes in cell samples. The system includes a digital microscope having a CCD camera that is controlled by and interfaced with a computer system. Images captured by the digital microscope are stored in an image processing board and manipulated by the computer system to detect the presence of malignancy-associated changes (MACs). At the present state of the art, it is believed that any detection of MACs requires images to be captured at a high spatial resolution, a high photometric resolution, that all information coming from the nucleus is in focus, that all information belongs to the nucleus (rather than some background), and that there is an accurate and reproducible segmentation of the nucleus and nuclear material. Each of these steps is described in detail below.
To detect the malignancy-associated changes, a cell sample is obtained and stained to identify the nuclear material of the cells and is imaged by the microscope. The stain is stoichiometric and specific to DNA only. The computer system then analyzes the image to compute a histogram of all pixels comprising the image. First, an intensity threshold is set that divides the background pixels from those comprising the objects in the image. All pixels having an intensity value less than the threshold are identified as possible objects of interest while those having an intensity value greater than the threshold are identified as background and are ignored.
For each object located, the computer system calculates the area, shape and optical density of the object. Those objects that could not possibly be cell nuclei are ignored. Next, the image is decalibrated, i.e., corrected by subtracting an empty frame captured before the scanning of the slide from the current frame and adding back an offset value equal to the average background light level. This process corrects for any shading of the system, uneven illumination, and other imperfections of the image acquisition system. Following decalibration, the images of all remaining objects must be captured in a more precise focus. This is achieved by moving the microscope in the stage z-direction in multiple focal planes around the approximate frame focus. For each surviving object a contrast function (a texture feature) is calculated. The contrast function has a peak value at the exact focus of the object. Only the image at the highest contrast value is retained in the computer memory and any object which did not reach such a peak value is also discarded from further considerations.
Each remaining in-focus object on the image is further compensated for local absorbency of the materials surrounding the object. This is a local decalibration which is similar to that described for the frame decalibration described above, except that only a small subset of pixels having an area equal to the area of a square into which the object will fit is corrected using an equivalent square of the empty frame.
After all images are corrected with the local decalibration procedure, the edge of the object is calculated, i.e., the boundary which determines which pixels in the square belong to the object and which belong to the background. The edge determination is achieved by the edge-relocation algorithm. In this process, the edge of the original mask of the first contoured frame of each surviving object is dilated for several pixels inward and outward. For every pixel in this frame a gradient value is calculated, i.e., the sum and difference between all neighbor pixels touching the pixel in question. Then the lowest gradient value pixel is removed from the rim, subject to the condition that the rim is not ruptured. The process continues until such time as a single pixel rim remains. To ensure that the proper edge of an object is located, this edge may be again dilated as before, and the process repeated until such time as the new edge is identical to the previous edge. In this way the edge is calculated along the highest local gradient.
The computer system then calculates a set of feature values for each object. For some feature calculations the edge along the highest gradient value is corrected by either dilating the edge by one or more pixels or eroding the edge by one or more pixels. This is done such that each feature achieves a greater discriminating power between classes of objects and is thus object specific. These feature values are then analyzed by a classifier that uses the feature values to determine whether the object is an artifact or is a cell nucleus. If the object appears to be a cell nucleus, then the feature values are further analyzed by the classifier to determine whether the nucleus exhibits malignancy-associated changes. Based on the number of objects found in the sample that appear to have malignancy-associated changes and/or an overall malignancy-associated score, a determination can be made whether the patient from whom the cell sample was obtained is healthy or harbors a malignant growth.
The foregoing aspects and many of the attendant advantages of this invention will become more readily appreciated [as the same becomes better understood] by reference to the following detailed description, when taken in conjunction with the accompanying drawings, wherein:
As described above, the present invention is a system for automatically detecting malignancy-associated changes (MACs) in the nuclei of cells obtained from a patient. From the presence or absence of MACs, a determination can be made whether the patient has a malignant cancer.
A block diagram of the MAC detection system according to the present invention is shown in FIG. 1. The system 10 includes a digital microscope 12 that is controlled by and interfaced with a computer system 30. The microscope 12 preferably has a digital CCD camera 14 employing a scientific CCD having square pixels of approximately 0.3 μm by 0.3 μm size. The scientific CCD has a 100% fill factor and at least a 256 gray level resolution. The CCD camera is preferably mounted in the primary image plane of a planar objective lens 22 of the microscope 12.
A cell sample is placed on a motorized stage 20 of the microscope whose position is controlled by the computer system 30. The motorized stage preferably has an automatic slide loader so that the process of analyzing slides can be completely automated.
A stable light source 18, preferably with feedback control, illuminates the cell sample while an image of the slide is being captured by the CCD camera. The lens 22 placed between the sample 16 and the CCD camera 14 is preferably a 20×/0.75 objective that provides a depth of field in the range of 1-2 μm that yields a distortion-free image. In the present embodiment of the invention, the digital CCD camera 14 used is the Microimager™ produced by Xillix Technologies Corp. of Richmond, B.C., Canada.
The images produced by the CCD camera are received by an image processing board 32 that serves as the interface between the digital camera 14 and the computer system 30. The digital images are stored in the image processing board and manipulated to facilitate the detection of MACs. The image processing board creates a set of analog video signals from the digital image and feeds the video signals to an image monitor 36 in order to display an image of the objects viewed by the microscope.
The computer system 30 also includes one or more input devices 38, such as a keyboard and mouse, as well as one or more peripherals 42, such as a mass digital. storage device, a modem or a network card for communicating with a remotely located computer, and a monitor 40.
At step 54, an image of a frame from the slide is captured by the CCD camera and is transferred into the image processor. In this process, the CCD sensor within the camera is cleared and a shutter of the camera is opened for a fixed period that is dependent on the intensity of the light source 18. After the image is optimized according to the steps described below, the stage then moves to a new position on the slide such that another image of the new frame can be captured by the camera and transferred into the computer memory. Because the cell sample on the slide occupies a much greater area than the area viewed by the microscope, a number of slide images are used to determine whether the sample is MAC-positive or negative. The position of each captured image on the slide is recorded in the computer system so that the objects of interest in the image can be found on the slide if desired.
Once an image from the slide is captured by the CCD camera and stored in the image processing board, the computer system determines whether the image produced by the CCD camera is devoid of objects. This is performed by scanning the digital image for dark pixels. If the number of dark pixels, i.e., those pixels having an intensity of the background intensity minus a predetermined offset value, is fewer than a predetermined minimum, the computer system assumes that the image is blank and the microscope stage is moved to a new position at step 60 and a new image is captured at step 54.
If the image is not blank, then the computer system attempts to globally focus the image. In general, when the image is in focus, the objects of interest in the image have a maximum darkness. Therefore, for focus determination the height of the stage is adjusted and a new image is captured. The darkness of the object pixels is determined and the process repeats until the average darkness of the pixels in the image is a maximum. At this point, the computer system assumes that global focus has been obtained.
After performing the rough, global focus at step 62, the computer system computes a histogram of all pixels. As shown in
Returning to
Returning to
The results of step 76 are shown in
Again returning to
If there are objects remaining in the image after the first attempt at removing artifacts at step 76, the computer system then compensates the image for variations in illumination intensity at step 80. To do this, the computer system recalls a calibration image that was obtained by scanning in a blank slide for the same exposure time that was used for the image of the cells under consideration. The computer system then begins a pixel-by-pixel subtraction of the intensity values of the pixels in the calibration image obtained from the blank slide from the corresponding pixels found in the image obtained from the cell sample. The computer system then adds a value equal to the average illumination of the pixels in the calibration image obtained from the blank slide to each pixel of the cell image. The result of the addition illuminates the cell image with a uniform intensity.
Once the variations in illumination intensity have been corrected, the computer system attempts to refine the focus of each object of interest in the image at step 82 (FIG. 2C). The optimum focus is obtained when the object has a minimum size and maximum darkness. The computer system therefore causes the stage to move a predefined amount above the global focus position and then moves in a sequence of descending positions. At each position the CCD camera captures an image of the frame and calculates the area and the intensity of the pixels comprising the remaining objects. Only one image of each object is eventually stored in the computer memory coming from the position in which the pixels comprising the object have the maximum darkness and occupy a minimum area. If the optimum focus is not obtained after a predetermined number of stage positions, then the object is removed from the computer memory and is ignored. Once the optimum focus of the object is determined, the image received from the CCD camera overwrites those pixels that comprise the object under consideration in the computer's memory. The result of the local focusing produces a pseudofocused image in the computer's memory whereby each object of interest is ultimately recorded at its best possible focus.
At a step 84, the computer system determines whether any in-focus objects in the cell image were found. If not, the computer system returns to step 54 shown in
Once an image of the object has been focused, the computer system then compensates for local absorbency of light near the object at a step 85. To do this, the computer system analyzes a number of pixels within a box having an area that is larger than the object by two pixels on all sides. An example of such a box is the box 207 shown in FIG. 6. The computer system then performs a pixel-by-pixel subtraction of the intensity values from a corresponding square in the calibration image obtained from the blank slide. Next the average illumination intensity of the calibration image is added to each pixel in the box surrounding the object. Then the average intensity value for those pixels that are in the box but are not part of the object is determined and this local average value is then subtracted from each pixel in the box that encloses the object.
Once the compensation for absorbency around the object has been made, the computer system then determines a more precise edge of each remaining object in the cell image at step 86. The steps required to compute the edge are discussed in further detail below.
Having compensated for local absorbency and located the precise edge of the object, the computer system calculates a set of features for each remaining object at a step 87. These feature values are used to further separate artifacts from cell nuclei as well as to identify nuclei exhibiting MACs. The details of the feature calculation are described below.
At a step 88, the computer system runs a classifier that compares the feature values calculated for each object and determines whether the object is an artifact and, if not, whether the object is a nucleus that exhibits MACs.
At a step 90, the pseudofocus digital image, the feature calculations and the results of the classifier for each in-focus object are stored in the computer's memory.
Finally, at a step 92, the computer system determines whether further scans of the slide are required. As indicated above, because the size of each cell image is much less than the size of the entire slide, a number of cell images are captured to ensure that the slide has been adequately analyzed. Once a sufficient number of cell images have been analyzed, processing stops at step 94. Alternatively, if further scans are required, the computer system loops back to step 54 and a new image of the cell sample is captured.
As indicated above, before the sample can be imaged by the digital microscope, the sample is stained to identify the nuclear material.
At step 120, the sample is soaked in a bisulfite solution for six minutes followed by a rinse for one minute with distilled water at a step 122. Next, the sample is dehydrated in solutions of 50%, 75% and 100% ethanol for approximately 10 seconds each at a step 124. The sample is then soaked in a final bath of xylene for one minute at a step 126 before a cover slip is applied at a step 128. After the cell sample has been prepared, it is ready to be imaged by the digital microscope and analyzed as described above.
Once the gradients have been calculated for each pixel in the annular ring, the computer system divides the range of gradients into multiple thresholds and begins removing pixels having lower gradient values from the ring. To remove the pixels, the computer scans the object under consideration in a raster fashion. As shown in
Pixels are removed from the annular ring subject to the conditions that no pixel can be removed that would break the chain of pixels around the annular ring. Furthermore, adjacent pixels cannot be removed during the same pass of pixel removal. Once all the pixels are removed having a gradient that is less than or equal to the first gradient threshold, the threshold is increased and the process starts over. As shown in
After locating the precise edge of an object, it is necessary to determine whether those pixels that comprise the edge should be included in the object. To do this, the intensity of each pixel that comprises the newly found edge is compared with its eight neighbors. As shown in
Once the features have been calculated for each in-focus object, the computer system must make a determination whether the object is a cell nucleus that should be analyzed for malignancy-associated changes or is an artifact that should be ignored. As discussed above, the system removes obvious artifacts based on their area, shape (sphericity) and optical density. However, other artifacts may be more difficult for the computer to recognize. To further remove artifacts, the computer system uses a classifier that interprets the values of the features calculated for the object.
As shown in
The particular thresholds used in the binary tree are set by statisticians who compare histograms of feature values calculated on known objects. For example, white blood cells typically have an area less than 50 μm2. Because the present invention treats a red blood cell as an artifact, the binary decision tree can contain a node that compares the area of an object to the 50 μm2 threshold. Objects with an area less than the threshold are ignored while those with an area having a greater area are further analyzed to determine if they are possible MAC cells or artifacts.
In the presently preferred embodiment of the invention, the discriminant functions that separate types of objects are generated by the BMDP program available from BMDP Statistical Software, Inc., of Los Angeles, Calif. Given the discriminant functions and the appropriate thresholds, the construction of the binary tree classifier is considered routine for one of ordinary skill in the art.
Once the binary tree classifier has been developed, it can be supplied with a set of feature values 292 taken from an unknown object and will provide an indication 294 of whether the object associated with the feature data is most likely an artifact or a cell nucleus.
The threshold above which it is likely that a patient exhibiting MACs has or will develop cancer is determined by comparing the MAC scores of a large numbers of patients who did develop cancer and those who did not. As will be appreciated by those skilled in the art, the particular threshold used will depend on the type of cancer to be detected, the equipment used to image the cells, etc.
The MAC detection system of the present invention can also be used to determine the efficacy of cancer treatment. For example, patients who have had a portion of a lung removed as a treatment for lung cancer can be asked to provide a sample of apparently normal cells taken from the remaining lung tissue. If a strong MAC presence is detected, there is a high probability that the cancer will return. Conversely, the inventors have found that the number of MAC cells decreases when a cancer treatment is effective.
As described above, the ability of the present invention to detect malignancy-associated changes depends on the values of the features computed. The following is a list of the features that is currently calculated for each in-focus object.
I.2 Coordinate Systems, Jargon and Notation
Each image is a rectangular array of square pixels that contains within it the image of an (irregularly shaped) object, surrounded by background. Each pixel Pi,j is an integer representing the photometric value (gray scale) of a corresponding small segment of the image, and may range from 0 (completely opaque) to 255 (completely transparent). The image rectangle is larger than the smallest rectangle that can completely contain the object by at least two rows, top and bottom, and two columns left and right, ensuring that background exists all around the object. The rectangular image is a matrix of pixels, Pi,j, spanning i=1, L columns and j=1, M rows and with the upper left-hand pixel as the coordinate system origin, i=j=1.
The region of the image that is the object is denoted by its characteristic function, Ω; this is also sometimes called the "object mask" or, simply, the "mask." For some features, it makes sense to dilate the object mask by one pixel all around the object; this mask is denoted Ω+. Similarly, an eroded mask is denoted Ω-. The object mask is a binary function:
where
and where "(i,j) ε object" means pixels at coordinates: (i, j) are part of the object, and "(i,j) ∉ object" means pixels at coordinates: (i, j) are not part of the object.
II Morphological Features
Morphological features estimate the image area, shape, and boundary variations of the object.
II.1 area
The area, A, is defined as the total number of pixels belonging to the object, as defined by the mask, Ω:
where i, j and Ω are defined in Section I.2 above.
II.2 x_centroid, y_centroid
The x_centroid and y_centroid are the coordinates of the geometrical center of the object, defined with respect to the image origin (upper-left hand corner):
where i and j are the image pixel coordinates and Ω is the object mask, as defined in Section 1.2 above, and A is the object area.
II.3 mean_radius, max_radius
The mean_radius and max_radius features are the mean and maximum values of the length of the object's radial vectors from the object centroid to its 8 connected edge pixels:
where rk is the kth radial vector, and N is the number of 8 connected pixels on the object edge.
II.4 var_radius
The var_radius feature is the variance of length of the object's radius vectors, as defined in Section II.3.
where rk is the kth radius vector, {overscore (r)} is the mean_radius, and N is the number of 8 connected edge pixels.
II.5 sphericity
The sphericity feature is a shape measure, calculated as a ratio of the radii of two circles centered at the object centroid (defined in Section II.2 above). One circle is the largest circle that is fully inscribed inside the object perimeter, corresponding to the absolute minimum length of the object's radial vectors. The other circle is the minimum circle that completely circumscribes the object's perimeter, corresponding to the absolute maximum length of the object's radial vectors. The maximum sphericity value: 1 is given for a circular object:
where rk is the kth radius vector.
II.6 eccentricity
The eccentricity feature is a shape function calculated as the square root of the ratio of maximal and minimal eigenvalues of the second central moment matrix of the object's characteristic function, Ω:
where λ1 and λ2 are the maximal and minimal eigenvalues, respectively, and the characteristic function, Ω, as given by Equation 1. The second central moment matrix is calculated as:
Eccentricity may be interpreted as the ratio of the major axis to minor axis of the "best fit" ellipse which describes the object, and gives the minimal value 1 for circles.
II.7 inertia_shape
The inertia_shape feature is a measure of the "roundness" of an object calculated as the moment of inertia of the object mask, normalized by the area squared, to give the minimal value 1 for circles:
where Ri,j is the distance of the pixel, Pi,j, to the object centroid (defined in Section II.2), and A is the object area, and Ω is the mask defined by Equation 1.
II.8 compactness
The compactness feature is another measure of the object's "roundness." It is calculated as the perimeter squared divided by the object area, giving the minimal value 1 for circles:
where P is the object perimeter and A is the object area. Perimeter is calculated from boundary pixels (which are themselves 8 connected) by considering their 4 connected neighborhood:
where N1 is the number of pixels on the edge with 1 non-object neighbor, N2 is the number of pixels on the edge with 2 non-object neighbors, and N3 is the number of pixels on the edge with 3 non-object neighbors.
II.9 cell_orient
The cell_orient feature represents the object orientation measured as a deflection of the main axis of the object from the y direction:
where Ymoment2 and xycrossmoment2 are the second central moments of the characteristic function Ω defined by Equation 1 above, and λ1 is the maximal eigenvalue of the second central moment matrix of that function (see Section II.6 above). The main axis of the object is defined by the eigenvector corresponding to the maximal eigenvalue. A geometrical interpretation of the cell_orient is that it is the angle (measured in a clockwise sense) between the y axis and the "best fit" ellipse major axis.
For slides of cell suspensions, this feature should be meaningless, as there should not be any a priori preferred cellular orientation. For histological sections, and possibly smears, this feature may have value. In smears, for example, debris may be preferentially elongated along the slide long axis.
II.10 elongation
Features in Sections II.10 to II.13 are calculated by sweeping the radius vector (from the object centroid, as defined in Section II.2, to object perimeter) through 128 discrete equal steps (i.e., an angle of 2π/128 per step), starting at the top leftmost object edge pixel, and sweeping in a clockwise direction. The function is interpolated from an average of the object edge pixel locations at each of the 128 angles.
The elongation feature is another measure of the extent of the object along the principal direction (corresponding to the major axis) versus the direction normal to it. These lengths are estimated using Fourier Transform coefficients of the radial function of the object:
where a2,b2 are Fourier Transform coefficients of the radial function of the object, r(θ), defined by:
II.11 freq_low_fft
The freq_low_fft gives an estimate of coarse boundary variation, measured as the energy of the lower harmonics of the Fourier spectrum of the object's radial function (from 3rd to 11th harmonics):
where an,bn are Fourier Transform coefficients of the radial function, defined in Equation 16.
II.12 freq_high_fft
The freq_high_fft gives an estimate of the fine boundary variation, measured as the energy of the high frequency Fourier spectrum (from 12th to 32nd harmonics) of the object's radial function:
where an,bn are Fourier Transform coefficients of the nth harmonic, defined by Equation 16.
II.13 harmon01_fft, . . . , harmon32 _fft
The harmon01_fft, . . . , harmon32_fft features are estimates of boundary variation, calculated as the magnitude of the Fourier Transform coefficients of the object radial function for each harmonic 1-32:
where an,bn are Fourier Transform coefficients of the nth harmonic, defined by Equation 16.
III Photometric Features
Photometric features give estimations of absolute intensity and optical density levels of the object, as well as their distribution characteristics.
III.1 DNA_Amount
DNA_Amount is the "raw" (unnormalized) measure of the integrated optical density of the object, defined by a once dilated mask, Ω+:
where the once dilated mask, Ω+ is defined in Section I.2 and OD is the optical density, calculated according to [12]:
where IB is the intensity of the local background, and Ii,j is the intensity of the i,j th pixel.
III.2 DNA_Index
DNA_Index is the normalized measure of the integrated optical density of the object:
where iodnorm is the mean value of the DNA amount for a particular object population from the slide (e.g., leukocytes).
III.3 var_intensity, mean_intensity
The var_intensity and mean_intensity features are the variance and mean of the intensity function of the object, I, defined by the mask, Ω:
where A is the object area, Ω is the object mask defined in Equation 1, and {overscore (I)} is given by:
{overscore (I)} is the "raw" (unnormalized) mean intensity.
mean intensity is normalized against iodnorm defined in Section III.2:
III.4 OD_maximum
OD_maximum is the largest value of the optical density of the object, normalized to iodnorm, as defined in Section III.2 above:
III.5 OD_variance
OD_variance is the normalized variance (second moment) of optical density function of the object:
where Ω is the object mask as defined in Section I.2, {overscore (OD)} is the mean value of the optical density of the object:
and A is the object area (total number of pixels). The variance is divided by the square of the mean optical density in order to make the measurement independent of the staining intensity of the cell.
III.6 OD_skewness
The OD_skewness feature is the normalized third moment of the optical density function of the object:
where Ω is the object mask as defined in Section 1.2, {overscore (OD)} is the mean value of the optical density of the object and A is the object area (total number of pixels).
III.7 OD_kurtosis
OD_kurtosis is the normalized fourth moment of the optical density function of the object:
where Ω is the object mask as defined in Section 1.2, {overscore (OD)} is the mean value of the optical density of the object and A is the object area.
IV Discrete Texture Features
The discrete texture features are based on segmentation of the object into regions of low, medium and high optical density. This segmentation of the object into low, medium and high density regions is based on two thresholds: optical density high threshold and optical density medium threshold. These thresholds are scaled to the sample's iodnorm value, based on the DNA amount of a particular subset of objects (e.g., lymphocytes), as described in Section III.2 above.
By default, these thresholds have been selected such that the condensed chromatin in leukocytes is high optical density material. The second threshold is located half way between the high threshold and zero.
The default settings from which these thresholds are calculated are stored in the computer as:
Ahigh is the area of the pixels having an optical density between 0 and 18, Amed. is the area of the pixels having an optical density between 18 and 36 and Alow is the area of the pixels having an optical density greater than 36. Together the areas Ahigh, Amed and Alow sum to the total area of the object. The actual thresholds used are these parameters, divided by 100, and multiplied by the factor iodnorm/100.
In the following discussion, Ωlow, Ωmed, and Ωhigh are masks for low-, medium-, and high-optical density regions of the object, respectively, defined in analogy to Equation 1.
IV.1 lowDNAarea, medDNAarea, hiDNAarea
These discrete texture features represent the ratio of the area of low, medium, and high optical density regions of the object to the total object area:
where Ω is the object mask as defined in Equation 1, and A is the object area.
IV.2 lowDNAamnt, medDNAamnt, hiDNAamnt
These discrete texture features represent the total extinction ratio for low, medium, and high optical density regions of the object, calculated as the value of the integrated optical density of the low-, medium-, and high-density regions, respectively, divided by the total integrated optical density:
where Ω is the object mask as defined in Equation 1, and OD is the optical density as defined by Equation 21.
IV.3 lowDNAcomp, medDNAcomp, hiDNAcomp, mhDNAcomp
These discrete texture features are characteristic of the compactness of low-, medium-, high-, and combined medium- and high-density regions, respectively, treated as single (possibly disconnected) objects. They are calculated as the perimeter squared of each region, divided by 4π (area) of the region.
where P is the perimeter of each of the optical density regions, defined in analogy to Equation 13, and A is the region area, defined in analogy to Equation 2.
IV.4 low_av_dst, med_av_dst, hi_av_dst, mh_av_dst
These discrete texture features represent the average separation between the low-, medium-, high-, and combined medium- and high-density pixels from the center of the object, normalized by the object mean_radius.
where Ri,j is defined in Section II.7 as the distance from pixel Pi,j to the object centroid (defined in Section II.2), and the object mean_radius is defined by Equation 5.
IV.5 lowVSmed_DNA, lowVShigh_DNA, lowVSmh_DNA
These discrete texture features represent the average extinction ratios of the low-density regions, normalized by the medium-, high-, and combined medium- and high-average extinction values, respectively. They are calculated as the mean optical density of the medium-, high-, and combined medium- and high-density clusters divided by the mean optical density of the low density clusters.
where OD is the region optical density defined in analogy to Equation 21, Ω is the region mask, defined in analogy to Equation 1, and A is the region area, defined in analogy to Equation 2.
IV.6 low_den_obj, med_den_obj, high_den_obj
These discrete texture features are the numbers of discrete 8-connected subcomponents of the objects consisting of more than one pixel of low, medium, and high density.
IV.7 low_cntr_mass, med_cntr_mass, high_cntr_mass
These discrete texture features represent the separation between the geometric center of the low, medium, and high optical density clusters (treated as if they were single objects) and the geometric center of the whole object, normalized by its mean radius.
where mean_radius of the object is defined by Equation 5, the object's centroid is defined in Section II.2, Ω is the region mask defined in analogy to Equation 1, and A is the region area defined in analogy to Equation 2.
V Markovian Texture Features
Markovian texture features are defined from the co-occurrence matrix, Δλ,μ of object pixels. Each element of that matrix stands for the conditional probability of the pixel of grey level λ occurring next (via 8-connectedness) to a pixel of grey level μ, where λ, μ are row and column indices of the matrix, respectively. However, the computational algorithms used here for the calculation of Markovian texture features uses so-called sum and difference histograms: Hls and Hmd, where Hls is the probability of neighboring pixels having grey levels which sum to 1, and Hmd is the probability of neighboring pixels having grey level differences of m, where an 8-connected neighborhood is assumed. Values of grey levels, l, m, used in the sum and difference histogram are obtained by quantization of the dynamic range of each individual object into 40 levels.
For completeness, the formulae that follow for Markovian texture features include both the conventional formulae and the computational formulae actually used.
V.1 entropy
The entropy feature represents a measure of "disorder" in object grey level organization: large values correspond to very disorganized distributions, such as a "salt and pepper" random field:
V.2 energy
The energy feature gives large values for an object with a spatially organized grey scale distribution. It is the opposite of entropy, giving large values to an object with large regions of constant grey level:
V.3 contrast
The contrast feature gives large values for an object with frequent large grey scale variations:
V.4 correlation
A large value for correlation indicates an object with large connected subcomponents of constant grey level and with large grey level differences between adjacent components:
where {overscore (Iq)} is the mean intensity of the object calculated for the grey scale quantized to 40 levels.
V.5 homogeneity
The homogeneity feature is large for objects with slight and spatially smooth grey level variations:
V.6 cl_shade
The cl_shade feature gives large absolute values for objects with a few distinct clumps of uniform intensity having large contrast with the rest of the object. Negative values correspond to dark clumps against a light background while positive values indicate light clumps against a dark background:
V.7 cl_prominence
The feature cl_prominence measures the darkness of clusters.
VI Non-Markovian Texture Features
These features describe texture in terms of global estimation of grey level differences of the object.
VI.1 den_lit_spot, den_drk_spot
These are the numbers of local maxima and local minima, respectively, of the object intensity function based on the image averaged by a 3×3 window, and divided by the object area.
where
and where
and I is the object intensity, Ω is the object mask, and A is the object area.
VI.2 range_extreme
This is the intensity difference between the largest local maximum and the smallest local minimum of the object intensity function, normalized against the slide DNA amount, iodnorm, defined in Section III.2. The local maxima, maxi',j' and minima, mini',j', are those in Section VI.1 above.
VI.3 range_average
This is the intensity difference between the average intensity of the local maxima and the average intensity of the local minima, normalized against the slide DNA amount value, iodnorm, defined in Section III.2 above. The local maxima, maxii',j' and minima, mini',j', values used are those from Section VI.1 above.
VI.4 center_of_gravity
The center_of_gravity feature represents the distance from the geometrical center of the object to the "center of mass" of the optical density function, normalized by the mean_radius of the object:
This gives a measure of the nonuniformity of the OD distribution.
VII Fractal Texture Features
The fractal texture features are based on the area of the three-dimensional surface of the object's optical density represented essentially as a three-dimensional bar graph, with the vertical axis representing optical density, and the horizontal axes representing the x and y spatial coordinates. Thus, each pixel is assigned a unit area in the x-y plane plus the area of the sides of the three-dimensional structure proportional to the change in the pixel optical density with respect to its neighbors. The largest values of fractal areas correspond to large objects containing small subcomponents with high optical density variations between them.
The difference between fractal1_area and fractal2_area is that these features are calculated on different scales: the second one is based on an image in which four pixels are averaged into a single pixel, thereby representing a change of scale of fractal1_area. This calculation needs the additional mask transformation: Ωi2,j2 represents the original mask Ω with 4 pixels mapped into one pixel and any square of 4 pixels not completely consisting of object pixels is set to zero. Ωi,j represents Ωi2,j2 expanded by 4 so that each pixel in Ωi2,j2 is 4 pixels in Ωi,j.
VII.1 fractal1_area
where OD*i,j is the optical density function of the image scaled by a factor common to all images such that the possible optical density values span 256 levels.
VII.2 fractal2_area
This is another fractal dimension, but based on an image in which four pixel squares are averaged into single pixels, thereby representing a change of scale of fractal1_area in Section VII.1 above.
where,
with L2, M2 as integers, and OD*i2,j2 is a scaled optical density function of the image, with 4 pixels averaged into one.
VII.3 fractal_dimen
The fractal_dimen feature is calculated as the difference between logarithms of fractal1_area and fractal2_area, divided by log 2. This varies from 2 to 3 and gives a measure of the "fractal behavior" of the image, associated with a rate at which measured surface area increases at finer and finer scales.
VIII Run Length Texture Features
Run length features describe texture in terms of grey level runs, representing sets of consecutive, collinear pixels having the same grey level value. The length of the run is the number of pixels in the run. These features are calculated over the image with intensity function values transformed into 8 levels.
The run length texture features are defined using grey level length matrices,
VIII.1 short0_runs, short45_runs, short90_runs, short135_runs
These give large values for objects in which short runs, oriented at 0°C, 45°C, 90°C, or 135°C, dominate.
VIII.2 long0_runs, long45_runs, long90_runs, long135_runs
These give large values for objects in which long runs, oriented at 0°C, 45°C, 90°C, or 135°C, dominate.
VIII.3 grey0_level, grey45_level, grey90_level, grey135_level
These features estimate grey level nonuniformity, taking on their lowest values when runs are equally distributed throughout the grey levels.
VIII.4 run0_length, run45_length, run90_length, run135_length
These features estimate the nonuniformity of the run lengths, taking on their lowest values when the runs are equally distributed throughout the lengths.
VIII.5 run0_percent, run45_percent, run90_percent, run135_percent
These features are calculated as the ratio of the total number of possible runs to the object's area, having its lowest value for pictures with the most linear structure.
where A is the object's area.
VIII.6 texture_orient
This feature estimates the dominant orientation of the object's linear texture.
where λ'1 is the maximal eigenvalue of the run length pseudosecond moment matrix (calculated in analogy to Section II.9). The run length pseudosecond moments are calculated as follows:
Orientation is defined as it is for cell_orient, Section II.9, as the angle (measured in a clockwise sense) between the y axis and the dominant orientation of the image's linear structure.
VIII.7 size_txt_orient
This feature amplifies the texture orientation for long runs.
where λ'1,λ'2 are the maximal and minimal eigenvalues of the run_length pseudosecond moment matrix, defined in Section VIII.6.
Each of the above features are calculated for each in-focus object located in the image. Certain features are used by the classifier to separate artifacts from cell nuclei and to distinguish cells exhibiting MACs from normal cells. As indicated above, it is not possible to predict which features will be used to distinguish artifacts from cells or MAC cells from non-MAC cells, until the classifier has been completely trained and produces a binary decision tree or linear discriminant function.
In the present embodiment of the invention, it has been determined that thirty (30) of the above-described features appear more significant in separating artifacts from genuine nuclei and identifying cells with MACs. These primarily texture features are as follows:
30 preferred nuclear features | |
1) Area | |
2) mean radius | |
3) OD variance | |
4) OD skewness | |
5) range average | |
6) OD maximum | |
7) density of light spots | |
8) low DNA area | |
9) high DNA area | |
10) low DNA amount | |
11) high DNA amount | |
12) high average distance | |
13) mid/high average distance | |
14) correlation | |
15) homogeneity | |
16) entropy | |
17) fractal dimension | |
18) DNA index | |
19) run 0 percent | |
20) run 45 percent | |
21) run 90 percent | |
22) run 135 percent | |
23) grey level 0 | |
24) grey level 45 | |
25) grey level 90 | |
25) grey level 135 | |
27) run length 0 | |
28) run length 45 | |
29) run length 90 | |
30) run length 135 | |
Although these features have been found to have the best ability to differentiate between types of cells, other object types may be differentiated by the other features described above.
As indicated above, the ability of the system according to the present invention to distinguish cell nuclei from artifacts or cells that exhibit MACs from those that do not depends on the ability of the classifier to make distinctions based on the values of the features computed. For example, to separate cell nuclei from artifacts, the present invention may apply several different discriminant functions each of which is trained to identify particular types of objects. For example, the following discriminant function has been used in the presently preferred embodiment of the invention to separate intermediate cervical cells from small picnotic objects:
cervical cells | picnotic | ||
max_radius | 4.56914 | 3.92899 | |
freq_low_fft | -.03624 | -.04714 | |
harmon03_fft | 1.29958 | 1.80412 | |
harmon04.fft | .85959 | 1.20653 | |
lowVSmed_DNA | 58.83394 | 61.84034 | |
energy | 6566.14355 | 6182.17139 | |
correlation | .56801 | .52911 | |
homogeneity | -920.05017 | -883.31567 | |
cl_shade | -67.37746 | -63.68423 | |
den_drk_spot | 916.69360 | 870.75739 | |
CONSTANT | -292.92908 | -269.42419 | |
Another discriminant function that can separate cells from junk particles is:
cells | junk | ||
eccentricity | 606.67365 | 574.82507 | |
compactness | 988.57196 | 1013.19745 | |
freq_low_fft | -2.57094 | -2.51594 | |
freq_high_fft | -28.93165 | -28.48727 | |
harmon02.fft | -31.30210 | -30.18383 | |
harmon03.fft | 14.40738 | 14.30784 | |
medDNAamnt | 39.28350 | 37.50647 | |
correlation | .27381 | .29397 | |
CONSTANT | -834.57800 | -836.19659 | |
Yet a third discriminant function that can separate folded cells that should be ignored from suitable cells for analysis.
normal interm | rejected objects | ||
sphericity | 709.66357 | 701.85864 | |
eccentricity | 456.09146 | 444.18469 | |
compactness | 1221.73840 | 1232.27441 | |
elongation | -391.76352 | -387.19376 | |
freq_high_fft | -37.89624 | -37.39510 | |
lowDNAamnt | -41.89951 | -39.42714 | |
low_den_obj | 1.40092 | 1.60374 | |
correlation | .26310 | .29536 | |
range_average | .06601 | .06029 | |
CONSTANT | -968.73628 | -971.18219 | |
Obviously, the particular linear discriminant function produced by the classifier will depend on the type of classifier used and the training sets of cells. The above examples are given merely for purposes of illustration.
As can be seen, the present invention is a system that automatically detects malignancy-associated changes in a cell sample. By properly staining and imaging a cell sample, the features of each object found on the slide can be determined and used to provide an indication whether the patient from which the cell sample was obtained is normal or abnormal. In addition, MACs provide an indication of whether cancer treatment given is effective as well as if a cancer is in remission.
While the preferred embodiment of the invention has been illustrated and described, it will be appreciated that various changes can be made therein without departing from the spirit and scope of the invention.
Palcic, Branko, Garner, David M., MacAulay, Calum E., Jaggi, Bruno W., Harrison, S. Alan
Patent | Priority | Assignee | Title |
11280990, | Feb 26 2018 | CALIBER IMAGING & DIAGNOSTICS, INC | System and method for macroscopic and microscopic imaging ex-vivo tissue |
11681133, | Feb 26 2018 | CALIBER IMAGING & DIAGNOSTICS, INC. | System and method for macroscopic and microscopic imaging ex-vivo tissue |
6684092, | Feb 24 1997 | CALIBER IMAGING & DIAGNOSTICS, INC | System for facilitating pathological examination of a lesion in tissue |
6856458, | Feb 17 1999 | CALIBER IMAGING & DIAGNOSTICS, INC | Tissue specimen holder |
6990239, | Jul 16 2002 | The United States of America as represented by the Secretary of the Navy | Feature-based detection and context discriminate classification for known image structures |
6999624, | Jul 12 2002 | The United States of America as represented by the Secretary of the Navy | Context discriminate classification for digital images |
6999625, | Jul 12 2002 | The United States of America as represented by the Secretary of the Navy | Feature-based detection and context discriminate classification for digital images |
7106892, | Sep 19 2002 | Koninklijke Philips Electronics N V | Display of image data information |
7194118, | Nov 10 2000 | CALIBER IMAGING & DIAGNOSTICS, INC | System for optically sectioning and mapping surgically excised tissue |
7379572, | Oct 16 2001 | University of Chicago | Method for computer-aided detection of three-dimensional lesions |
7526116, | Jan 19 2006 | NICHEVISION INCORPORATED | Automated microscopic sperm identification |
7720272, | Jan 19 2006 | NICHEVISION INCORPORATED | Automated microscopic sperm identification |
7864996, | Feb 17 2006 | CALIBER IMAGING & DIAGNOSTICS, INC | System for macroscopic and confocal imaging of tissue |
7907271, | Jun 15 2007 | Novartis AG | Method and system for standardizing microscope instruments |
7978258, | Aug 31 2007 | Novartis AG | Automatic exposure time selection for imaging tissue |
8027030, | Jun 15 2007 | Novartis AG | Method and system for standardizing microscope instruments |
8120768, | Jun 15 2007 | Novartis AG | Method and system for standardizing microscope instruments |
8121365, | Aug 07 2007 | Novartis AG | Method and system for determining an optimal dilution of a reagent |
8149506, | Feb 17 1999 | CALIBER IMAGING & DIAGNOSTICS, INC | Cassette for facilitating optical sectioning of a retained tissue specimen |
8160348, | Aug 06 2007 | Novartis AG | Methods and system for validating sample images for quantitative immunoassays |
8335360, | May 14 2007 | Novartis AG | Compartment segregation by pixel characterization using image data clustering |
8340372, | Dec 17 2002 | Qinetiq Limited | Image analysis |
8417015, | Aug 06 2007 | Novartis AG | Methods and system for validating sample images for quantitative immunoassays |
8572146, | Aug 17 2010 | Fujitsu Limited | Comparing data samples represented by characteristic functions |
8583718, | Aug 17 2010 | Fujitsu Limited | Comparing boolean functions representing sensor data |
8620854, | Sep 23 2011 | Fujitsu Limited | Annotating medical binary decision diagrams with health state information |
8645108, | Aug 17 2010 | Fujitsu Limited | Annotating binary decision diagrams representing sensor data |
8655037, | May 14 2007 | Novartis AG | Compartment segregation by pixel characterization using image data clustering |
8719214, | Sep 23 2011 | Fujitsu Limited | Combining medical binary decision diagrams for analysis optimization |
8781995, | Sep 23 2011 | Fujitsu Limited | Range queries in binary decision diagrams |
8812943, | Sep 23 2011 | Fujitsu Limited | Detecting data corruption in medical binary decision diagrams using hashing techniques |
8838523, | Sep 23 2011 | Fujitsu Limited | Compression threshold analysis of binary decision diagrams |
8874607, | Aug 17 2010 | Fujitsu Limited | Representing sensor data as binary decision diagrams |
8909592, | Sep 23 2011 | Fujitsu Limited | Combining medical binary decision diagrams to determine data correlations |
8930394, | Aug 17 2010 | Fujitsu Limited | Querying sensor data stored as binary decision diagrams |
9002781, | Aug 17 2010 | Fujitsu Limited | Annotating environmental data represented by characteristic functions |
9052523, | Feb 17 1999 | CALIBER IMAGING & DIAGNOSTICS, INC | Cassette for facilitating optical sectioning of a retained tissue specimen |
9075908, | Sep 23 2011 | Fujitsu Limited | Partitioning medical binary decision diagrams for size optimization |
9128024, | Feb 17 1999 | CALIBER IMAGING & DIAGNOSTICS, INC | Tissue specimen holder |
9138143, | Aug 17 2010 | Fujitsu Limited | Annotating medical data represented by characteristic functions |
9176819, | Sep 23 2011 | Fujitsu Limited | Detecting sensor malfunctions using compression analysis of binary decision diagrams |
9177247, | Sep 23 2011 | Fujitsu Limited | Partitioning medical binary decision diagrams for analysis optimization |
9229210, | Feb 26 2012 | CALIBER IMAGING AND DIAGNOSTICS, INC | Tissue specimen stage for an optical sectioning microscope |
9240043, | Sep 16 2008 | Novartis AG | Reproducible quantification of biomarker expression |
9772486, | Feb 17 1999 | CALIBER IMAGING & DIAGNOSTICS, INC. | System for optically sectioning a tissue specimen |
9865086, | Mar 21 2014 | ST. JUDE MEDICAL, CARDIOLOLGY DIVISION, INC. | Methods and systems for generating a multi-dimensional surface model of a geometric structure |
Patent | Priority | Assignee | Title |
4453266, | Apr 21 1980 | Rush-Presbyterian-St. Luke's Medical Center | Method and apparatus for measuring mean cell volume of red blood cells |
5016283, | Nov 04 1986 | Cell Analysis Systems, Inc.; CELL ANALYSIS SYSTEMS, INC , A CORP OF IL | Methods and apparatus for immunoploidy analysis |
5099521, | Mar 24 1988 | FIRST UNION NATIONAL BANK OF FLORIDA ATTN: JOSEPH D SISTARE, III | Cell image processing method and apparatus therefor |
5889881, | Oct 14 1992 | PERCEPTRONIX MEDICAL INC | Method and apparatus for automatically detecting malignancy-associated changes |
6026174, | Oct 14 1992 | PERCEPTRONIX MEDICAL INC | System and method for automatically detecting malignant cells and cells having malignancy-associated changes |
Executed on | Assignor | Assignee | Conveyance | Frame | Reel | Doc |
Mar 26 1999 | Monogen, Inc. | (assignment on the face of the patent) | / | |||
Sep 13 2001 | ONCOMETRICS IMAGING CORP | MONOGEN, INC | ASSIGNMENT OF ASSIGNORS INTEREST SEE DOCUMENT FOR DETAILS | 012282 | /0586 |
Date | Maintenance Fee Events |
Jun 02 2006 | M2551: Payment of Maintenance Fee, 4th Yr, Small Entity. |
Jul 19 2010 | REM: Maintenance Fee Reminder Mailed. |
Dec 10 2010 | EXP: Patent Expired for Failure to Pay Maintenance Fees. |
Date | Maintenance Schedule |
Dec 10 2005 | 4 years fee payment window open |
Jun 10 2006 | 6 months grace period start (w surcharge) |
Dec 10 2006 | patent expiry (for year 4) |
Dec 10 2008 | 2 years to revive unintentionally abandoned end. (for year 4) |
Dec 10 2009 | 8 years fee payment window open |
Jun 10 2010 | 6 months grace period start (w surcharge) |
Dec 10 2010 | patent expiry (for year 8) |
Dec 10 2012 | 2 years to revive unintentionally abandoned end. (for year 8) |
Dec 10 2013 | 12 years fee payment window open |
Jun 10 2014 | 6 months grace period start (w surcharge) |
Dec 10 2014 | patent expiry (for year 12) |
Dec 10 2016 | 2 years to revive unintentionally abandoned end. (for year 12) |