U.S. patent application number 14/347912 was filed with the patent office on 2014-08-21 for systems and methods for automated screening and prognosis of cancer from whole-slide biopsy images.
This patent application is currently assigned to BOARD OF REGENTS OF THE UNIVERSITY OF TEXAS SYSTEM. The applicant listed for this patent is Sos S. Agaian, Ali Almuntashri, Clara M. Mosquera Lopez, Richard Metzler. Invention is credited to Sos S. Agaian, Ali Almuntashri, Clara M. Mosquera Lopez, Richard Metzler.
Application Number | 20140233826 14/347912 |
Document ID | / |
Family ID | 47996715 |
Filed Date | 2014-08-21 |
United States Patent
Application |
20140233826 |
Kind Code |
A1 |
Agaian; Sos S. ; et
al. |
August 21, 2014 |
SYSTEMS AND METHODS FOR AUTOMATED SCREENING AND PROGNOSIS OF CANCER
FROM WHOLE-SLIDE BIOPSY IMAGES
Abstract
The invention provides systems and methods for detection,
grading, scoring and tele-screening of cancerous lesions. A
complete scheme for automated quantitative analysis and assessment
of human and animal tissue images of several types of cancers is
presented. Various aspects of the invention are directed to the
detection, grading, prediction and staging of prostate cancer on
serial sections/slides of prostate core images, or biopsy images.
Accordingly, the invention includes a variety of sub-systems, which
could be used separately or in conjunction to automatically grade
cancerous regions. Each system utilizes a different approach with a
different feature set. For instance, in the quantitative analysis,
textural-based and morphology-based features may be extracted at
image- and (or) object-levels from regions of interest.
Additionally, the invention provides sub-systems and methods for
accurate detection and mapping of disease in whole slide digitized
images by extracting new features through integration of one or
more of the above-mentioned classification systems. The invention
also addresses the modeling, qualitative analysis and assessment of
3-D histopathology images which assist pathologists in
visualization, evaluation and diagnosis of diseased tissue.
Moreover, the invention includes systems and methods for the
development of a tele-screening system in which the proposed
computer-aided diagnosis (CAD) systems. In some embodiments, novel
methods for image analysis (including edge detection, color mapping
characterization and others) are provided for use prior to feature
extraction in the proposed CAD systems.
Inventors: |
Agaian; Sos S.; (San
Antoinio, TX) ; Lopez; Clara M. Mosquera; (San
Antonio, TX) ; Almuntashri; Ali; (San Antonio,
TX) ; Metzler; Richard; (San Antonio, TX) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
Agaian; Sos S.
Lopez; Clara M. Mosquera
Almuntashri; Ali
Metzler; Richard |
San Antoinio
San Antonio
San Antonio
San Antonio |
TX
TX
TX
TX |
US
US
US
US |
|
|
Assignee: |
BOARD OF REGENTS OF THE UNIVERSITY
OF TEXAS SYSTEM
Austin
TX
|
Family ID: |
47996715 |
Appl. No.: |
14/347912 |
Filed: |
September 26, 2012 |
PCT Filed: |
September 26, 2012 |
PCT NO: |
PCT/US2012/057264 |
371 Date: |
March 27, 2014 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
61539827 |
Sep 27, 2011 |
|
|
|
Current U.S.
Class: |
382/133 |
Current CPC
Class: |
G06T 2207/20036
20130101; G06K 9/46 20130101; A61B 5/4381 20130101; G01N 15/1475
20130101; G06K 9/00127 20130101; G06T 2207/10056 20130101; G06T
2207/30024 20130101; G06T 7/181 20170101; G06T 2207/20081 20130101;
G16H 50/70 20180101; G16H 50/30 20180101; G06T 7/0012 20130101;
G06K 9/00147 20130101; A61B 5/7257 20130101; G06K 9/0014 20130101;
G06T 2207/30081 20130101; G06T 2207/10024 20130101; G06T 5/002
20130101; G06T 5/20 20130101; G16H 50/20 20180101; A61B 5/725
20130101; G06T 2207/30096 20130101; A61B 5/7203 20130101; A61B
5/726 20130101; G06T 7/90 20170101; A61B 5/7267 20130101; G06T 7/13
20170101; A61B 10/0241 20130101; G06K 9/6228 20130101; G16H 30/20
20180101 |
Class at
Publication: |
382/133 |
International
Class: |
G06T 7/00 20060101
G06T007/00; G06K 9/00 20060101 G06K009/00 |
Claims
1. A method, performed by a processing unit, for computing pixels
along object edges and producing a deinterlaced image from an
interlaced source, the method comprising: performing image
filtering on a collected image depending on the nature of noise in
the collected image; smoothing the filtered image using
shape-dependent filters; calculating gradient vectors in the image
using different kernels; selecting an edge angle; determine
threshold values within a local dynamic range, generating several
edge maps, and fusing the generated edge maps together.
2. The method of claim 1, further comprising applying an algorithm
for image edge-preserving contrast enhancement which is based on
HVS, Parameterized Logarithm Image Processing operations, and is
integrated with morphological log-ratio approach in order to come
up with an effective edge detection operator that is sensitive to
edges of dark areas in the image.
3. Methods and systems for computation of fractal-like dimension
for modifying input data/image into data describing the image, the
method comprising: performing image filtering depending on the
nature of noise in the image; binarizing the image by using the
method described in claim 1 or 2, or by using an image bit-plane
decomposition method; calculating the fractal dimension of binary
images by using different grid sizes based on a Differential Box
Counting (DBC) algorithm; and fusing resulting fractal
dimensions.
4. A method for computing the fractal dimension of color images
comprising: performing a color model transformation by using
arbitrary functions operating on the original RGB components of the
image; splitting the image in smaller windows and computing for
each block the probability of having m points/pixels in a hypercube
of the size of the window; estimating the color fractal dimension
by using a weighting function and fitting the curve of logarithm of
the window sizes against the logarithm of the total number of boxes
of each window size needed to cover the image.
5. A method for carcinomas color region mapping comprising:
constructing 2-D projections of the RGB color model of a large
variety of biopsy images including normal and cancerous tissue;
computing 2-D histograms of biopsy images; selecting the most
prominent colors by using an algorithm for local maxima location in
3-D surfaces; constructing a color region mapping, in which each
color corresponds to one class or tissue structure.
6. A method for image color quantization comprising: performing an
image color standardization procedure; calculating the distance
between each pixel and the colors considered in the color map of
claim 4 by using a distance metric; minimizing the calculated
distance; and assigning the nearest reference color or class to the
each pixel.
7. A textural-based system and methods for automatically detecting,
classifying, and grading cancerous regions of a histology image
comprising: performing an image color standardization procedure;
forming texture-based feature vectors by using spatial and
transforms domain image information along with fractal analysis;
selecting the group of features that best describe the images;
training a classifier by using the generated feature vectors;
classifying histology images according to the Gleason grading
system; using the result of classification to determine the Gleason
score of the image; assessing the accuracy of the Gleason
grading/scoring system by using cross-validation methods.
8. The method of claim 7, wherein forming texture-based feature
vectors comprises extraction of a set of features from Fourier and
wavelet transforms such us wavelet energy and entropy, phase
information of coefficients, and spatial domain algorithms such as
statistical spatial filtering, wavelet-based fractal dimension
according to the method of claim 3 or others, and gray level
histogram.
9. A textural-based system and methods for automatically detect,
classify, and grade cancerous regions of a histology image
(particularly a prostate biopsy images) comprising: performing an
image color standardization procedure; describing the image data
and forming texture-based feature vectors by using spatial and
transforms domain image information along with color fractal
analysis as claimed in claim 4; selecting the group of features
that best describe the images; training a classifier by using the
generated feature vectors; classifying histology images according
to the Gleason grading system; using the result of classification
to determine the Gleason score of the image; assessing the accuracy
of the Gleason grading/scoring system by using cross-validation
methods.
10. The description step of the system of claim 9, further
comprising the extraction of a set of features from real,
logarithmic, or complex wavelets and multi-wavelets and a new color
fractal dimension algorithm developed in this invention.
Particularly, transform features include joint probability of
detail coefficients, phase information of coefficients, energy
distribution of detail coefficients, and energy of color channels
of the biopsy image. A new color fractal dimension is also
extracted to complete the feature space.
11. A morphology- and architectural-based system and methods for
automatically detect, classify, and grade cancerous regions of a
histology image (particularly a prostate biopsy images) comprising:
performing an image color standardization procedure; segmenting the
image by applying any proper algorithm including the method
proposed in claims 5 and 6; describing the image data and forming
morphology and architectural-based feature vectors by extracting
the shape, size and arrangement of tissue structures; selecting the
group of features that best describe the images; training a
classifier by using the generated feature vectors; classifying
cancerous images according to grading systems, for example Gleason
grading system for prostate carcinomas; using the result of
classification to determine the Gleason score of the image;
assessing the accuracy of the Gleason grading/scoring system by
using cross-validation methods.
12. The description step of the system of claim 11, further
comprising the extraction of features: color information such as
color distribution entropy and 2-D color histograms along with
tissue morphology and architectural features from segmented tissue
structures.
13. Systems and methods for automatically detecting, classifying
and grading cancerous regions of 3-D histological images
(particularly a prostate image) comprising: performing an image
color standardization procedure, segmenting image by applying any
proper algorithm including the method proposed in claims 5 and 6;
mapping of 2-D images or an image into 3-D space by using
pre-developed mapping algorithms, describing the image data and
forming color-, texture-, morphology-, and architectural-based
feature vectors; training a classifier within an framework by using
generated feature vectors; accessing the cancer grade or scale, for
example, according to the Gleason grading scale; accessing the most
prominent, the second most prominent pattern and the third
prominent pattern in said image data and computing a Gleason 3-D
score as a sum of Gleason 3-D grades of said patterns; accessing
the accuracy of Gleason 3-D grading/scoring system; estimating
volume of a tumor region.
14. The systems as claimed in 13, wherein the mapping from 2-D
image to 3-D images includes algorithms for 3-D reconstruction from
a single 2-D image or from several 2-D slides. The processed images
nay be gray scaled or colored images.
15. The systems of claims 7, 9, 11 and 13 wherein automatically
classifying the biopsy images is based on machine learning
algorithms: linear discriminants, Gaussian models, Multivariate
Adaptive Regression Splines (MARS), Classification an Regression
trees (CART), decision trees, k-nearest neighbor, Bayesian, neural
networks, and/or SVM. Boosting algorithms such as AdaBoost, SVM
Boost, etc. can be used to improve the classification
performance.
16. A method to obtain a multidimensional Gleason grading wherein a
revised Gleason value is generated as follows: Revised Gleason
Value=2D Gleason Grade+3rd Dimensional Core Grade
17. A system and method for automatically detecting, classifying
and grading cancerous regions from digitized Whole-Slide
(particularly a prostate image) comprising: performing an image
color standardization procedure; splitting image into smaller
regions; determining and grading a location of a tumor in the
prostate using one or more of the systems of claim 7 or 9 or 11 or
13; accessing a cancer stage; accessing the most prominent, the
second most prominent pattern and third prominent pattern in said
image data and computing a Gleason score as a sum of Gleason grades
of said patterns; generating a cancer map from whole digitized
slides according to Gleason grading/scoring levels associated with
slides processing block-elements; accessing the accuracy of Gleason
grading/scoring system; accessing the estimating size of a tumor
region.
18. The system according to claim 17, wherein said grade estimating
module is configured for estimating a size of a tumor in said at
least one tissue region or whole digitized slides.
19. The system of claim 17, wherein a complete diagnosis is
performed for providing useful information for pathologists about
the disease severity and the localization, size and volume of
cancerous tissue.
20. The system of claim 17 wherein complete diagnosing is
performed, the outputs of the systems include: tumor extent in the
needle biopsy, number and identification of positive needle cores,
percentage of positive cores for cancer, length of the tumor in
positive cores and length of the core, percentage of needle core
tissue affected by carcinoma and area of tissue positive for
cancer, Gleason grade of each selected cancerous region of
interest, percentage of each Gleason pattern in biopsy specimens,
Gleason score of each positive core, including tertiary pattern
information, the most frequent Gleason score in slides, numbers of
identical Gleason scores in the positive cores, revised Gleason
value (claim 16), and localized cancer map including the grade of
each cancerous patch in all needle biopsy slides.
21. The system of claim 17 further comprising methods for surgery
quality control, which allow a surgeon to be immediately aware of
whether the surrounding tissue of a specimen is or not positive for
cancer, prior the patient is closed up.
22. A system for cancer prognosis, wherein the pathology report
including all the outputs of systems claimed in 7, 9, 11 and 13 is
integrated with patient information and PSA analysis results to
construct feature vectors. The resulting feature vectors are used
to aid pathologists in cancer prognosis, risk assessment and
treatment planning.
23. The system of claim 22, further comprising the use of feature
vectors for prostate cancer risk assessment by implementing risk
calculators according to the patient's available information.
24. The system of claim 22 further comprising a method for cancer
prediction based on 2-D and 3-D histopathology images using cancer
map or Gleason graded block and pixel features to predict the state
of the tissue in missing blocks within the scope of the image or to
produce extrapolated values to know the cancer grade in the
neighborhoods outside the core extracted from the prostate gland
during a needle biopsy procedure.
25. The method of claim 24 wherein Kalman filtering or similar
techniques are used for estimation and prediction of cancer grade
(Gleason grade) outside of the scope of observed histopathology
biopsy images.
26. A system and methods for cancer tele-screening and diagnosis
based on biopsy image and an adjudication scheme, wherein expert
pathologists and CAD systems are integrated.
27. The adjudication methods for prostate cancer diagnoses used in
the system of claim 26 consist of replacing one or more
pathologists by CAD systems and evaluate the agreement among them
to produce a faster and less expensive prostate cancer
diagnosis.
28. The systems of claim 26 wherein said adjudication methods is
configure for using morphology- or textural-based 2-D
classification systems for cancer detection and grading, such that
the CADs integrated into the systems perform independent tissue
assessment.
29. The tele-screening system of claim 26 uses communication
networks such as internet for image and patient data acquisition,
pathologists' report acquisition. The final assessment results are
stored in a centralized database available for pathologist
consultation.
30. A method and system constructing an optimal signal detector
using a learning algorithm in tandem with a prediction algorithm
(which may be one in the same) and data compression algorithm.
31. A method and system for detecting carcinomas in a sample image
based on compression rates obtained using a detector trained on
previously classified carcinoma images as specified in claim
30.
32. A method and system for using detection results to compute a
hard, integer Gleason score of a query image through selection of a
best compressing pattern compressors (PCs) using one of the systems
of claims 30 and 31.
33. A method and system for using detection results to compute a
hard, integer Gleason scores of query image pixels and regions
through selection of a best compressing pattern compressors (PCs)
using one of the systems of claims 30 and 31.
34. A method and system for using detection results to compute a
soft, continuous-scale Gleason score of an image (query),
comprising: compressing query image data using specific pattern PCs
and optionally a self PC according to systems in claims 30 and 31;
computing pattern weights determined by: PC compression rates, PC
compression rate statistics, and/or functions of each pattern PC
compression rates or statistics and query image self PC compression
rates or statistics; combining pattern values and weights into a
single Gleason score.
35. A method and system for using detection results to compute a
soft, continuous-scale Gleason scores of query image pixels and
regions, comprising: compressing query image pixels using specific
pattern PCs and optionally a self PC according to systems in claims
30 and 31; computing pattern weights determined by either by
pattern PC compression rate or by comparison of each pattern PC
prediction to a query image's self PC prediction on a pixel by
pixel basis; combining pattern values and weights into a single,
soft Gleason score on a pixel by pixel basis.
36. A system and method for determining the most likely
distribution of Gleason patterns within a query image, comprising:
the systems of claim in 30, 31 and 33 or 35; determining pairwise
uncertainty statistics between pixel scores or grades; forming the
uncertainty statistics into a matrix; computing the eigenvalues of
said matrix; computing the eigenvectors of said matrix; using all
or some of the above features to estimate the most likely current
distribution of Gleason patterns within the sample; using all or
some of the above features to estimate the most likely future
distribution of Gleason patterns within the sample;
37. A system and method for automatically detecting, classifying
and grading cancerous regions from digitized, Whole-Slides
(particularly a prostate image) comprising: segmenting the image
into prominent features; optionally downsizing the image;
converting the image into a 1 dimensional signal; compressing the
signal using the systems of claims 30 and 31; assigning Gleason
scores to each pixel using the system of claims 33 or 35, and 36;
determining, grading, and scoring a location of a tumor in the
prostate using one or more of the systems of claim 33, 35 or 36;
assessing the and grading a location of a tumor in the prostate
using one or more of the systems of claim 33, 35 or 36; estimating
the current carcinoma pattern distribution of a selected region
using the systems of claim 33, 35 or 36; estimating the future
carcinoma pattern distribution of a selected region using the
systems of claim 33, 35 or 36; generating a cancer map from whole,
digitized slides according to the systems of claim 33, 35 or 36;
assessing the accuracy of Gleason grading/scoring system;
38. The system according to claim 37, wherein said grade estimating
module is configured for estimating a size of a tumor in said at
least one tissue region or whole digitized slides.
39. The system of claim 37, wherein a complete diagnosis is
performed for providing useful information for pathologists about
the disease severity and the localization, size and volume of
cancerous tissue.
40. The system of claim 37 wherein complete diagnosing is
performed, the outputs of the systems include: tumor extent in the
needle biopsy, number and identification of positive needle cores,
percentage of positive cores for cancer, length of the tumor in
positive cores and length of the core, percentage of needle core
tissue affected by carcinoma and area of tissue positive for
cancer, Gleason grade of each selected cancerous region of
interest, percentage of each Gleason pattern in biopsy specimens,
Gleason score of each positive core, including tertiary pattern
information, the most frequent Gleason score in slides, numbers of
identical Gleason scores in the positive cores, revised Gleason
value (claim 16), and localized cancer map including the grade of
each cancerous patch in all needle biopsy slides.
41. The system of claim 37 further comprising methods for surgery
quality control, which allow a surgeon to be immediately aware of
whether the surrounding tissue of a specimen is or not positive for
cancer, prior the patient is closed up.
42. A system for cancer prognosis, wherein the pathology report
including all the outputs of systems claimed in 30 through 41 is
integrated with patient information and PSA analysis results to
construct feature vectors. The resulting feature vectors are used
to aid pathologists in cancer prognosis, risk assessment and
treatment planning.
Description
BACKGROUND OF THE INVENTION
[0001] 1. Field of the Invention
[0002] The invention generally relates to diagnostic procedures for
the detection, analysis, assessment, classification, and grading of
cancer from digitized histopathological images. More particularly,
the invention relates to the use of computer analysis of biopsy
samples to detect and classify/grade and score prostate cancer.
[0003] 2. Description of the Relevant Art
[0004] Prostate cancer is the second killer of American men.
Current statistics indicate that 1 in 6 men will suffer from
prostate cancer, 1 in 32 men will die from prostate cancer.
Although prostate cancer can be diagnosed at early stages, 80% of
men diagnosed with prostate carcinoma are above 65 years old.
[0005] The risk of developing prostate cancer is commonly attached
with the patient's age. In general, doctors recommend males of age
50 or higher to screen for prostate cancer. However, screening
could be recommended for males who are at a higher risk of
developing prostate cancer even at younger age, e.g. 40-45. Risk
factors generally increase the chances of developing prostate
cancer. The major risk factor is age, since 63% of prostate cancer
cases occur in males of age 65 or older. Furthermore, family
history is another major indicator. Men who have a family history
of prostate cancer have higher chances of developing prostate
cancer. Considering race, it is reported that African American men
have a higher rate of developing prostate cancer compared to Asian
and Native Americans who have the lowest rate. Considering food
diet, it has been shown that diets that are higher in fats may lead
to a higher chance in developing prostate cancer.
[0006] Biopsies provide the most significant information to
urologists by indicating the stage and the grade of prostatic
carcinoma using the Gleason grading system.
[0007] Tissue samples are sampled by urologists from different
prostate areas called "cores" using a hollow needle. When biopsy
samples are taken from prostate tissues, they are prepared for
observation under microscope by specialized pathologists. Cores are
sampled from suspicious and non suspicious areas in order to have
prognostic information on the areas where cancerous cells may
spread to. Also, sampling different cores would identify cases
where cancer is slowly growing compared to other areas where cancer
is growing more aggressively.
[0008] There are different methodologies and protocols that
laboratories follow in preparing slides of prostate biopsies. In
one method, for a single core examination, three slides are usually
prepared. Each slide may have serial sections of biopsy slices or
multiple slices sections from the same core. The number of
sectioned slices may vary between 2-4 slices. Assuming that four
slices are extracted for each slide, the total number of slices
needed to prepare three slides is 12 slices. For serial slices,
only the 1st and the 3rd slides are stained using Hematoxylin &
Eosin (H&E) or Keratin 903 (K903) prior to pathologists'
investigation. Staining is a significant step that clearly
differentiates various textural components of the biopsy especially
basal cells of the prostate glands. Malignant cells do not have
these cells and hence, can be observed clearly. The next step is
microscopic investigation where the stained slides 1 and 3 are
checked for malignancy. If any suspicious patterns are observed,
the 2nd slide is then stained using (K903) and observed. However,
if the 3rd slide only shows suspicious patterns, further slices are
sectioned from the same core and stained for further investigation.
In case of a negative diagnosis of the current core, slides from
other cores of the prostate gland are prepared for further
diagnosis.
[0009] Pathologists grade cancer tissues by observing biopsy
samples under the microscope to describe their appearance and
architectural features by assigning a grade using Gleason grading
system. Gleason system of grading is a significant histology
grading method due to its popularity and usage in many medical
institutions as well as in literature. The Gleason grading system
was named after Dr. Gleason in 1977, who devised a system of five
grades to describe and differentiate the unique histology features
of each biopsy image excluding the cytological features. The
Gleason grading system assigns grades from 1 to 5 according five
patterns where grade 1 is the most differentiated pattern, and
grade 5 is the least differentiated one. FIG. 1 illustrates the
different architectures of these patterns labeled with their
corresponding grade. Gleason grade 1: cancerous patterns are
composed of generally consistent, well separated, and closely
packed glands. This pattern resembles the normal prostate biopsy
tissue. Gleason grade 2: tumors of this pattern show less uniformed
glands having size and shape variations, with loose arrangements of
the glands and more gaps in between. Gleason grade 3: tumor glands
in this pattern can be recognizable but with higher variations of
shapes (elongated) and mostly smaller in size. Darker cells' nuclei
can be noticed leaving the glands and invading the surrounding
glands. Gleason grade 4: tumors of this pattern can be recognized
by the loss of glands formation, and the inability to separate
glands from each other. It has an architectural jump compared to
previous grades by having fused and poorly defined structures. The
invading of dark nuclei becomes a dominant pattern in this grade.
Gleason grade 5: tumors of this pattern have no glandular
differentiation and consist of sheets of single cells. They have a
total loss of architecture with more dark nuclei pattern spread all
over the prostate tissue.
[0010] Patterns in Gleason grades 1 and 2 are the least common
patterns and seldom occur in general. Gleason pattern 3 and 4 are
the most common patterns, and their existence indicate a worsen
prognosis in general. Gleason pattern 5 is less common than the
previous pattern and seldom seen as diagnosis usually occurs at an
earlier stage of cancer development.
[0011] During a pathologist's examination of biopsy images, they
assign two Gleason grades that represent the most two predominant
patterns in the image. The first grade indicates the most dominant
pattern, and the second grade indicates the second most predominant
pattern. The addition of these two grades yields the Gleason score.
For instance, a Gleason score of 7: 4+3 indicates that prostate
cancer tumor of pattern 4 is the most dominant pattern, combined
with second predominant pattern of grade 3. It is important to note
here that a Gleason score of 4+3 is not the same as 3+4. Although
this yields a total score of 7, the first and second predominance
grades are different at both cases.
[0012] During the microscopic observation of Gleason patterns in
order to assign the two most predominant grades to the biopsy
sample, pathologists may note the existence of a third pattern, or
a "tertiary pattern". Tertiary patterns account for almost 5% of
patterns belonging to grades 4 or 5 patterns. Hence, pathologists
generally report the Gleason score with the tertiary pattern score
even if it accounts for less than 5% due to its role in predicting
advanced tumor stages at some cases.
[0013] Currently, pathologists visually grade prostate biopsy
tissues using Gleason scoring system. This approach is very
subjective and subject to intra and inter-observer variations.
Also, it is a time consuming task and raises difficulties as far as
spatial resolution is considered especially in the subgroups of
grades 3-5 where further classification is needed. There is a
considerable variation of pathologists' grading of prostate biopsy
images. For instance, at Johns Hopkins Hospital, it was reported
that most "2+2 "biopsies sent for evaluation were actually "3+3"
biopsies that had been under graded by a general pathologist.
Pathologists experience is another factor that contributes to
errors in reporting the correct Gleason grades. Most of these
errors result in a repeat needle biopsy procedure for the patient
that could have been avoided if the original biopsy specimens had
been reviewed by more experienced pathologists. In a study
conducted on 3,789 patients collected from 18 studies, the
following were reported: the exact correlation of Gleason grade was
found in 43% of the total cases, the accuracy of Gleason grade plus
or minus one grade was found in 77% of total cases, the under
grading of prostate carcinoma using needle biopsy occurred in 43%,
and the over-grading cases accounted for 15% of the total
cases.
[0014] There is a clear need of improved systems and methods for
screening, quantitative image analysis, and grading of prostate
cancer, providing reliable information on the lesion's extension
and localization. An accurate diagnosis is a crucial step towards
patient treatment selection and management. Several systems have
been proposed for prostate cancer detection and grading from
digitized biopsy slides. The following patents and patent
applications may be considered relevant to the field of the
invention:
[0015] U.S. Pat. No. 7,761,240 B2 to Saidi et al. discloses systems
and methods are provided for automated diagnosis and grading of
tissue images based on morphometric data extracted from the images
by a computer in a 2-level procedure (detection and grading). The
system proceeds as follows: (1) Input image, (2) Background
removal, (3) Histogram-based objects segmentation, (4) Image-level
and object-level feature extraction, (5) Feature selection and/or
classification. The extracted features include: Fractal dimension
of binarized image by several thresholds, fractal code, variance of
symlet4 wavelet coefficients, color channel histograms and tissue
structures classification. The maximum reported accuracy in the
detection stage is 96% and 77.6% in the grading stage. These
systems do not have explored the integration among the color
channels, which could improve the features extraction steps.
[0016] U.S. Pat. No. 8,139,831 B2 to Khamene et al. discloses a
method for unsupervised classification of histological images
obtained from a slide simultaneously co-stained with NIR
fluorescent and H&E stains. The system proceeds as follows: (1)
Prostate gland units segmentation, (2) Feature extraction, (3)
Principal Component Analysis, (4) Bayesian multiclass
classification (PIN, Gleason grades 1 to 5). (5) Gleason score is
also computed. The extracted features include: gland boundary and
region description, structural descriptors of glands by using
Voronoi, Delaunay and Minimum spanning tree graph features, as long
with textural descriptors (Haralick features and Gabor filters). In
this system, only gland centroids are used for graph construction.
Some other important structures such as cell nuclei were not
explores as a nodes of of Voronoi, Delaunay and MST graphs. A
drawback of this method is that additional NIR stain must be used.
The most widespread in clinical practice is H&E stain.
[0017] U.S. Pat. App No. 2010/0098306 A1 to Madabhushi et al.
relates to computer-aided diagnosis using content-based retrieval
of histopathological (prostate/breast/renal) image features based
on predetermine criteria and their analysis for malignancy
determination. The method classifies a region into normal,
cancerous and suspect. The system proceeds as follows: (1) Input
image, (2) Feature extraction, (3) Manifold learning (feature
dimensionality reduction and classification) for classification and
reclassification at higher magnifications The extracted features
include: first and second order statistics of color channels, graph
features (Voronoi, Delaunay and MST from nuclei centers. The
features are standard deviation, average, ratio minimum to maximum,
entropy of edges and areas of polygons and triangles, Statistics of
Co-adjacency matrix constructed from the glands centers), textural
features (Haar wavelet feature and Gabor filter response, Sobel,
derivatives and Kirsch filter response) and morphology information
(nuclear density and gland morphology). The maximum reported
accuracy is 92.8% in Grade 3 vs. Stroma recognition. However, the
classification between Grade 3 vs. 4 is 76.9%.
[0018] U.S. Pat. App No. 2011/0243417 A1 to Madabhushi et al.
discloses a method and a system for detecting biologically relevant
structures in a hierarchical fashion, beginning at a
low-resolution. The relevant structures are classified using
pairwise Markov models (PPMMs). The invention provides a fast CAD
system for detection/diagnosis of medical images. This invention is
used for CaP detection in whole mount WMHS, but can be used for
other organs. The extracted features include: Graph features
(Voronoi diagram, Dealunay triangulation and MST from lymphocyte
centroids) and Gland features and nuclei density. The performance
of the system is estimated by using AUC of ROC curves. The maximum
reported AUC is 0.812.
[0019] U.S. Pat. App No. 2007/0019854 A1 to Gholap et al. discloses
a method and a system that provide automated screening of prostate
needle biopsy specimen in a digital image of prostatectomy
specimen. The system proceeds as follows: (1) Input image, (2)
Contrast modification (if needed), (3) structure segmentation
(lumen, cells, cytoplasm), (4) filtering of isolated cells in
stromal area, (5) cell dilation, (6) gland reconstruction, (7)
feature extraction, (8) classification using neural networks.
2-level process (disease detection and further grading). The
features are: Number of glands, gland size, shape (circularity,
elongation), arrangement and destruction, stroma area, lymphocytes
presence, among others.
[0020] U.S. Pat. App No 2010/0312072 to Breskin et al. discloses a
method of estimating a grade of a prostate cancer from zinc data
associated with the prostate, the zinc data being arranged gridwise
in a plurality of picture-elements representing a zinc map of the
prostate. The method performs cancer grading of at least one tissue
region based on zinc levels associated with zinc map. In this case,
the assessment of cancer severity is not made from histological
quantitative analysis of tissue images.
[0021] It is therefore desirable to have more effective and
accurate procedures for the detection of prostate cancer using
biopsy samples. Unfortunately, all of the above procedures and
apparatus suffer deficiencies in the degree of quantitative
analysis and diagnostic accuracy in the assessment of human or
animal tissue images.
SUMMARY OF THE INVENTION
[0022] The large number of disparate, differently-scaled Gleason
scoring methods can lead to a variety of misinterpretations of
reading actual cancer aggressiveness. We propose generalized
scoring methods which simultaneously covers all previous scoring
methods by relieving the dependency on only two or three of the
most dominant patterns, extending the range of measured patterns to
all possible Gleason patterns, and combining all observed Gleason
pattern grades into one, continuous-scale Gleason score which is
bounded by the traditional Gleason Grade limits, 1 through 5. These
novel scores provide a more accurate, unambiguous measure of
prostate cancer aggressiveness.
[0023] Only few systems have addressed the biopsy whole slide
processing. While some of the systems and methods mentioned above
are able to classify predefined regions of interest, they do not
apply to a biopsy image as a whole. Certain embodiments disclosed
herein map Gleason scores to every, individual pixel within a
sample image or group of image slides. From each pixel score, local
to global measurements about the entire sample are made
available.
[0024] Furthermore, biopsy assessment through integration of both
whole slide 2-D and 3-D systems is described herein. To the best of
our knowledge 3-D classification systems have not been investigated
for Gleason grading. We propose a novel approach in cancer
diagnosis providing a revised Multidimensional Gleason grading when
2-D and 3-D evaluation results are provided. Furthermore, systems
and methods for cancer tele-screening by integrating CAD systems
and pathologist experience are presented in this invention.
[0025] To tackle the limitations of tissue cancer imaging systems,
the present disclosure addresses and overcomes the aforementioned
limitations of the prior art by providing methods and systems of
automated or partially automated quantitative analysis and
assessment of human or animal tissue images such as liver biopsy,
endometrial biopsy, lung biopsy, lymph node biopsy, renal biopsy,
bladder biopsy, rectal biopsy, skin biopsy, and other serial
sections/slides of prostate core images. The methods and systems
are completely general and independent of the process under
consideration.
[0026] In an embodiment, a complete scheme/framework for detecting
and grading of cancerous lesions through automated quantitative
analysis and assessment of human or animal tissue images is
provided. The framework detects, grades, predicts stages of cancers
on serial sections/slides of core or biopsy images. This output
provides information about disease severity, location, and
distribution. Specifically, the output of the invention provides
both an overall Gleason grade and a "map" of cancerous tissue
locations. From this, it is possible to determine the sizes and
amounts of tissue affected by cancer, the primary sample grade,
secondary sample grades, continuous-scale Gleason score, and
higher-order sample grades. Furthermore, it is possible to
determine the most likely distribution of grades within a
sample.
[0027] The framework of the system consists of a diversity of
sub-systems. Each sub-system could be used separately or in a
combined/fused way. Each system use different approaches. For
instance, one sub-system is based on tissue structural and
morphology features and other sub-systems are based on textural
characteristics of the tissue.
[0028] Also provided are systems and methods for accurate detection
and mapping of disease in whole slide images, generating so called
"cancer maps" by extracting new features and integrating one or
more of the classification systems mentioned above.
[0029] In one embodiment, new methods for image analysis and
description of biopsy images are based on edge detection, transform
based features, color and bit-plans based fractal dimensions, and
color mapping are described.
[0030] In another embodiment, new methods for image analysis and
description of biopsy images are based on machine learning,
sequential prediction, and information analysis.
[0031] In one embodiment, a novel "Hard" Gleason score is generated
by (see FIG. 23): Hard Gleason Score=Most Likely Integer Gleason
Grade.
[0032] In another embodiment, a novel "Soft" Gleason score is
generated by (see FIG. 24): Soft Gleason Score=Most Likely
Continuous-Scale Gleason Grade.
[0033] In another embodiment, a novel "Hard" Gleason map is
generated using the most likely Gleason grade of each sample image
pixel (see FIG. 25).
[0034] In another embodiment, a novel "Soft" Gleason map is
generated using the most likely continuous-scale Gleason grade of
each sample image pixel (see FIG. 26).
[0035] In one embodiment, a novel "Gleason distribution" is
generated using the uncertainty in the most likely Gleason grades
of each sample image pixel.
[0036] The scheme/framework based system and sub-systems have many
advantages. They are fast, simple, inexpensive, and provides a
screening test for cancer that does not give a high percentage of
false negative or false positive results. The method facilitates
early detection of the disease and it can be used for screening of
large populations. The invention systems can be used by
pathologists or without pathologists.
[0037] In other embodiments, a novel revised Gleason value is
generated by (see FIG. 16): Revised Gleason Value=2D Gleason
Grade+3rd Dimensional Core Grade.
[0038] In other aspects, systems and methods for the development of
a tele-screening system are provided in which the proposed
computer-aided diagnosis (CAD) systems are integrated as biopsy
readers according to the system configuration, pathologists--expert
systems level of agreement and diagnosis decision rules.
[0039] Embodiments of the computer-aided system can be used for
both: (1) assisting pathologists in several stages. For example:
localization of cancer regions in a large digitized tissue biopsy,
and grading of detected regions, automated quantitative analysis
and assessment, (2) Independent assessment without pathologist
interaction.
[0040] These and other aspects of the present invention will become
evident upon reference to the following detailed description
BRIEF DESCRIPTION OF THE DRAWINGS
[0041] Advantages of the present invention will become apparent to
those skilled in the art with the benefit of the following detailed
description of embodiments and upon reference to the accompanying
drawings in which:
[0042] FIG. 1 is visual representation of the Gleason grading
system;
[0043] FIG. 2 is a flow diagram of the Canny's edge detection
algorithm;
[0044] FIG. 3 is an illustration of the concept for fusion of edge
detected images;
[0045] FIG. 4 is a comparison of the results between one Canny edge
detected image and the output of the proposed edge detection
algorithm;
[0046] FIG. 5 is an example of results of tissue structure
segmentation using color region mapping in prostate cancer
biopsy;
[0047] FIG. 6 shows an exemplary flow diagram of 2-D Gleason
grading systems;
[0048] FIG. 7 shows exemplary morphology- and architectural-based
Gleason grading systems;
[0049] FIG. 8 is a block diagram of the first textural-based 2-D
Gleason grading/scoring system.
[0050] FIG. 9 is a flow diagram of the second textural-based 2-D
Gleason grading/scoring system using the proposed new Color Fractal
Dimension;
[0051] FIG. 10 shows two examples of color model transformation for
computation of new Color Fractal Dimension;
[0052] FIG. 11 is a general flow diagram of the 3-D Gleason grading
system;
[0053] FIG. 12 presents an example of prostate cancer localization
and cancer mapping from prostate cancer biopsy;
[0054] FIG. 13 illustrates a process to reconstruct 3-D images from
2-D slices;
[0055] FIG. 14 shows an example of 3-D reconstruction of segmented
lumens and epithelial nuclei from a prostate cancer biopsy
image;
[0056] FIG. 15 presents a 3-D core grading system;
[0057] FIG. 16 shows the steps for computation of "Revised Gleason
grade";"
[0058] FIG. 17 is a flow diagram of the integration of CADs with
PSA and patient's information for cancer prognosis;
[0059] FIG. 18 shows the structure of the proposed cancer
tele-screening system;
[0060] FIG. 19 is a method and system for diagnosis of cancer by
integrating two pathologists and one CAD;
[0061] FIG. 20 is a method for diagnosis of cancer by integrating
one pathologist and two CADs;
[0062] FIG. 21 presents a fully automated method for diagnosis of
cancer by integrating three CADs;
[0063] FIG. 22 shows an alternative exemplary classification scheme
based on morphology- and architectural-based Gleason feature
selection;
[0064] FIG. 23 shows exemplary steps for computation of a "Hard
Gleason grade";
[0065] FIG. 24 shows exemplary steps for computation of a "Soft
Gleason grade";
[0066] FIG. 25 shows exemplary steps for computation of a "Hard
Gleason map";
[0067] FIG. 26 shows exemplary steps for computation of a "Soft
Gleason map";
[0068] FIG. 27 shows exemplary results for computation of pattern
weights used for the "Soft Gleason Score";
[0069] FIG. 28 shows a flow diagram for the implementation of a
method for cancer prediction based on histopathology images
grading.
[0070] While the invention may be susceptible to various
modifications and alternative forms, specific embodiments thereof
are shown by way of example in the drawings and will herein be
described in detail. The drawings may not be to scale. It should be
understood, however, that the drawings and detailed description
thereto are not intended to limit the invention to the particular
form disclosed, but to the contrary, the intention is to cover all
modifications, equivalents, and alternatives falling within the
spirit and scope of the present invention as defined by the
appended claims.
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0071] It is to be understood the present invention is not limited
to particular devices or methods, which may, of course, vary. It is
also to be understood that the terminology used herein is for the
purpose of describing particular embodiments only, and is not
intended to be limiting. As used in this specification and the
appended claims, the singular forms "a", "an", and "the" include
singular and plural referents unless the content clearly dictates
otherwise. Furthermore, the word "may" is used throughout this
application in a permissive sense (i.e., having the potential to,
being able to), not in a mandatory sense (i.e., must). The term
"include," and derivations thereof, mean "including, but not
limited to." The term "coupled" means directly or indirectly
connected.
[0072] Embodiments herein relate to automated quantitative analysis
and assessment of human or animal tissue images such as liver
biopsy, endometrial biopsy, lung biopsy, lymph node biopsy, renal
biopsy, bladder biopsy, rectal biopsy, skin biopsy, and other
serial sections/slides of prostate core images. More particularly,
but not exclusively, the invention relates to detection, grading,
prediction, and staging of prostate cancer on serial
sections/slides of prostate core images, or part of biopsy images,
which are illustrated as examples.
[0073] It is to be understood that the present invention is not
limited to particular devices or methods, which may, of course,
vary. It is also to be understood that the terminology used herein
is for the purpose of describing particular embodiments only, and
is not intended to be limiting. As used in this detailed
specification and the appended claims, the singular forms "a",
"an", and "the" include singular and plural referents unless the
content clearly dictates otherwise. Furthermore, the word "may" is
used throughout this application in a permissive sense (i.e.,
having the potential to, being able to), not in a mandatory sense
(i.e., must). The term "include," and derivations thereof, mean
"including, but not limited to." The term "coupled" means directly
or indirectly connected. The abbreviation "CaP" refers to prostate
cancer.
[0074] Described herein is the framework for the development of the
cancer grading systems, the functional classification of components
and methods for image analysis and systems integration.
The Edge Detection Method
[0075] Image edge detection is one of the most effective
preprocessing tools that provide essential image edge information
and characteristics. An edge detector can be defined as a
mathematical operator that responds to the spatial change and
discontinuities in a gray-level (luminance) of a pixel set in an
image. This method can be used in wide areas such as image
segmentation, image categorization, image registration, image
visualization, and pattern recognition. These applications may vary
in their outputs but they all share the common need of precise edge
information in order to carry out the needed tasks successfully.
Medical images edge detection is an important step towards object
recognition of the human organs such as brain soft tissues, and
other different organs. It represents an essential preprocessing
work in medical image segmentation. It has also been used for other
medical applications concerned with the reconstruction of Computed
Tomography (CT) and MRI scans. Conventionally, after image
acquisition takes place, various post processing algorithms are
applied to these images in order to extract more information. This
information can be obtained after edge detection, especially of
tissue borders. This process plays a significant rule in
recognition of pathological alternations especially in MRI
images.
[0076] One of the well-known edge detection algorithms is the Canny
edge detection method. Canny method has proven to be superior over
many of the available edge detection algorithms and thus was chosen
mostly for real time implementation and testing. It is considered
as the modern "standard" in the sense that the validity of all
other algorithms is often checked against it. In Canny algorithm,
the Gaussian function is used to smooth the image prior to edge
detection. The filtering or smoothing operation actually services
two purposes. The first one is the reduction of the effect of
Gaussian noise prior to the detection of pixel intensity changes.
The second purpose is setting the resolution or scale at which
intensity changes are to be detected. These factors contribute
effectively towards a precise edge detection method, overcoming the
problem of detecting false edges resulted from noise sitting on
some pixels.
[0077] Canny algorithm has the limitation of being vulnerable to
noise disturbances and impulsive noise. Hence, there are certain
limitations to its application. Also, the traditional Canny
operator has a precise edge positioning but quite sensitive to
noise and may detect false edges, as well as missing some details
of real edges in an image. This issue is of a great importance in
noisy images where many false edges resulting from noise are
detected. Furthermore, the Canny edge detector cannot detect
branching edges to some extent. Because of the existence of
impulsive noise in medical images, traditional Canny algorithm
doesn't perform well since it can't eliminate some noise patterns,
thus, generating many false edges. It is difficult to design a
general edge detection algorithm which performs well in many
contexts and captures the requirements of subsequent processing
stages. Consequently, over the history of digital image processing,
a variety of edge detectors have been devised which differ in their
mathematical and algorithmic properties.
[0078] Canny operator consists of the main processes shown in FIG.
2. The major steps of Canny edge detection algorithm can be
summarized as follows: (1) Smooth an image with Gaussian filter to
have a noise-free image prior to applying the derivative function;
(2) calculate gradient magnitude and gradient direction; (3)
perform non-maximal suppression to ensure the desired edge with one
single pixel width; (4) Determine two threshold values, and then
select possible edge points and trace edges.
[0079] In one limitation with the Canny method, gradient magnitude
and direction are calculated by pixels within a 2-by-2 neighborhood
is sensitive to noise and may detect false edges. Another
limitation is the failure to detect branched edges and curves
especially in neighborhoods where the change in intensity values is
inconsiderable. Furthermore, non-maximal suppression can have a
negative impact on edge detection precision and may miss some edge
connected points. However, it becomes a very essential task to
detect all of these edges especially in medical scans where a
diagnosis decision is going to be made. The developed algorithms
proposes different modifications to the existing Canny algorithm
that succeeds in detecting more edge details as well as image
denoising when compared to the traditional Canny edge algorithm.
Also, one might note that the upper and lower threshold levels are
fixed at Canny algorithm regardless of the nature of the image. The
determination of such threshold values needs to be an image
dependent procedure where the characteristics of different
images/edge information are utilized for an automated selection of
threshold levels.
[0080] There are a few improved algorithms that were developed
based on the traditional Canny algorithm that has a superior
performance to suit our recognition purposes. Such proposed
algorithms may include replacements/modifications of different
processes in Canny method. In one embodiment of the present
invention, a novel algorithm provides methods for multiple edge
detection. This algorithm has a superior performance in localizing
and detecting edges that cannot be detected using the classical
edge detection techniques, including Canny's method. To this end,
several modifications to the original Canny's algorithms are
introduced: [0081] (1) Improved preprocessing: Gaussian smoothing
filter is replaced by median/adaptive median filter depending on
the nature of noise in the image. This enhances the accuracy of
detecting true edges at areas of interest as well as avoiding the
detection of noise edges that could not be removed using the
Gaussian filter. [0082] (2) Improved calculation of the gradient
vector: introducing different kernel sizes at the traditional
9-pixels neighborhood area replacing the existing 3.times.3 kernel.
This adds more freedom in calculating the magnitude at different
directions other than 45 and 135 degrees. Also, other algorithms
may place certain weights of the gradient at the direction of
interest. [0083] (3) Improved selection of threshold values:
algorithms utilizing adaptive approach in automating the selection
of threshold values may replace the existing fixed threshold values
introduced by Canny. This improves edge preservation at different
areas of the image. [0084] (4) Improved non-maximal suppression
(NMS): Canny NMS process compares pixels in the n-pixels
neighborhood area at eight different directions in order to
determine local maxima. An improved algorithm considers dividing
the area into different sub-areas and applies NMS at each area
before a final decision is made about the local maxima.
[0085] In an embodiment, an approach in edge detection is developed
based on dynamic range enhancement algorithm. The introduced
algorithm succeeds in integrating the concepts of logarithmic ratio
of human visual system with contrast enhancement in order to
develop a new edge detection method. The developed method has the
advantage over other traditional methods in terms of enhancing,
detecting, and segmenting different spots and shapes having
significant dark gray levels. The introduced algorithm is tested on
different classes and shows a superior performance compared to
Canny's edge detection operator. Also, it has the advantage of the
precise segmentation and topology/shape preservation of regions of
interest needed for further measurements and analysis purposes. In
the embodiment algorithm, a non-linear image enhancement which is
based on Human Visual System (HVS) dynamic range is integrated with
morphological log-ratio approach in order to come up with an
effective edge detection operator that is sensitive to edges of
dark areas in the image. In applications where edges of interest
are represented by dark gray levels such as cracks texture, a
pre-processing step is needed where spatial filters such
Gaussian/median filters are used to smooth the image prior to
applying any sort of enhancement to avoid noise enhancement.
[0086] As mentioned previously in an embodiment, a Human Visual
System (HVS)-based contrast enhancement and edge detection
algorithm may be used to analyze images. The introduced algorithm
integrates image enhancement, edge detection and logarithmic ratio
filtering algorithms to develop an effective edge detection method.
Also a parameter ((3) is introduced to control the level of
detected edge details and functions as a primary threshold
parameter. The introduced algorithm functions in tracking and
segmenting significant dark gray levels in an image. Simulation
results have shown the effectiveness of the introduced algorithm
compared to other traditional methods such as Canny's algorithm in
preserving object's topology and shape. The developed algorithm
functions at various applications where measurements and
segmentation of dark gray level spots for classification and
tracking purposes are required such as road cracks, lunar surface
images, and remote objects.
[0087] The Human Visual System (HVS) is responsible for
transferring data into information received by the viewer. The
received information can be characterized by attributes such as
brightness, edges, and color shades. Considering brightness, it can
be shown that it is not a simple function of intensity when
simultaneous contrasts are considered since it relates the
perceived brightness of a given region to its background. The
reason behind this is that our HVS perception is sensitive to
luminance contrast rather than the absolute luminance values
themselves. Also, one might consider eyes sensitivity to white
noise added to a uniform background more than in an area with a
high contrast. Eyes are more sensitive to random noise in smooth
areas of an image than in busy areas with many details. Hence, it
is obvious that HVS is not a simple function of intensity. Weber's
law depends on the ratio of intensity values at two points f (x,y)
and f.sub.o (x,y) more than the difference of these values
according to the following Weber's law:
f - f o f o = .gradient. f f o = .omega. ( Weber ` s constant )
##EQU00001##
[0088] In another application of parameterized logarithmic image
processing (PLIP), images are processed as gray tone function
according to the following expression:
g(i,j)=M-f(i,j)
[0089] To begin with, the concept of parameterized logarithmic
image processing (PLIP) was designed to both maintain the pixel
values inside the a bounded rang [0,M) where M is the maximum value
of the range. For instance, M could be 256 for 8 bit gray level
images. Also, PLIP model processes images from human visual system
(HVS) point of view in a more accurate way. The usual arithmetic
operations such as addition and scalar multiplication are argued
not to be suitable for human visual system. The PLIP model
considers processing images as gray tone functions. The
relationship between the traditional gray level function and the
gray tone function is given in as follows:
f(x,y)=M-{tilde over (f)}(x,y)
[0090] Where M is the maximum value of the range. The various
operations associated with PLIP include the PLIP addition, PLIP
subtraction, and PLIP multiplications. The operation of interest in
the developed edge detection algorithm is the PLIP addition
operation. For two gray tone functions a and b, and an arbitrary
linear function .gamma.(M)=AM+B for constants A and B, the PLIP
addition is defined as:
a .sym. b = a + b - a b .gamma. ( M ) ##EQU00002##
[0091] In order to develop PLIP-based edge detection, Canny
algorithm operations are modified to suit PLIP operations. The
gradient magnitude g and direction .theta. calculations are
replaced by PLIP-based formula as follows:
Magnitude g=M(1-exp((- {square root over ( gx.sup.2+
gy.sup.2)})/M))
Direction .theta.=Arc tan( gx, gy)
[0092] As an example of this embodiment, gx and gy may be defined
in a m.times.n window m,n=1, 2, 3 . . . . The values of m and n may
be arbitrary and should be selected according to the application.
For instance, a window of 3.times.3 of 9 gray tone pixels may be
generated as follows:
f 1 f 2 f 3 f 4 f 5 f 6 f 7 f 8 f 9 ##EQU00003## g x = M - M ( f ~
3 f ~ 6 2 f ~ 9 f ~ 1 f ~ 4 2 f ~ 7 ) , g y = M - M ( f ~ 1 f ~ 2 2
f ~ 3 f ~ 7 f ~ 8 2 f ~ 9 ) g _ = .PHI. ( g ) ##EQU00003.2##
[0093] where .phi.(g) is described as isomorphism operator:
.PHI. ( g ) = - M ln ( 1 - g M ) ##EQU00004##
[0094] This section illustrates the mathematical derivation used in
developing the new edge detector. In current classical edge
detection methods, an edge is detected based on a variation of gray
level intensity levels. The levels of detected edges details are
controlled by the selection of threshold levels where a lower
threshold value corresponds to more edge details and vice versa. In
this introduced approach, edges are detected based on the
difference of adjacent pixels' gray levels instead. The gray level
difference is described for the horizontal and vertical directions
as follows:
.DELTA..sub.xI(i,j)=I(i,j)-I(i-1,j);
.DELTA..sub.yI(i,j)=I(i,j)-I(i,j-1)
[0095] Such difference in gray levels is applied also using the
ratio operation of the same surrounding pixels. The maximum ratio
only should be selected when examining rows falling/rising edges,
and column falling/rising edges.
[0096] Rational morphological operations have been investigated in
order to design rational morphological filters:
R = ( MF 1 MF 2 ) .alpha. , R = i = 1 m a i ( MF 1 MF 2 ) i j = 0 n
b j ( MF 1 MF 2 ) j ##EQU00005##
[0097] Where MF1 and MF2 represent two different classes of filters
and a.sup.i and b.sup.i, i=0, 1, . . . m are constants.
[0098] Other forms of the developed visual filters deploy
logarithmic ratio as the following examples:
( log MF i log MF j ) n ; ( log ( MF i ) .THETA.log ( MF j ) log (
MF i ) .sym. log ( MF j ) ) n ##EQU00006##
[0099] The basic steps of the new algorithm are the following: (1)
Smooth an image with shape dependent filter. (2) Calculate the
joint shape-dependent gradient magnitude and direction. (3) Fusion
of edge images. Some examples of this embodiment are the
following:
[0100] Step 1: For example the Smoothing operator could be the
Gaussian kernel, or
[ - 1 - 1 - 1 - 1 .alpha. - 1 - 1 - 1 - 1 ] , [ - 1 - 1 - 2 - 1 - 1
- 2 .alpha. - 2 - 1 - 1 - 2 - 1 - 1 ] , [ - 1 - 2 - 1 - 2 .alpha. -
2 - 1 - 2 - 1 ] , [ - 1 - 1 - 2 .alpha. - 2 - 1 - 1 ]
##EQU00007##
[0101] Where .alpha.=2, 4, or 8.
[0102] Step 2: Joint shape-dependent gradient magnitude and
direction could be described as:
Magnitude ( x , y , .theta. ) = Max ( cos .theta. .differential. G
.differential. x , sin .theta. .differential. G .differential. y )
##EQU00008##
[0103] Where the gradient magnitude and direction are joined and
calculated using pixels within an m.times.n neighborhood, where
m,n=1, 2, 3 . . . .
[0104] This argument assigns the maximum value of the two functions
to the magnitude at a given point, hence, reduces the magnitude
value at angles with low intensity values. Note that the Canny
method gradient magnitude and direction are calculated by pixels
within a 2.times.2 neighborhood is sensitive to noise and may
detect false edges. Also, a set of defined gradient kernel sets for
derivative approximation can be used effectively for edge
detection, line detection, and consequently for feature extraction.
Such kernels are shape dependent and may vary in sizes and
dimensions. A kernel size could be of n.times.m pixels, for example
of 3.times.5, 5.times.7 or 9.times.5 as shown in the next kernel
matrices.
[ 1 2 0 - 2 - 1 2 2 0 - 2 - 2 1 2 0 - 2 - 1 ] [ 1 2 2 0 - 2 - 2 - 1
2 2 2 2 0 - 2 2 - 2 - 2 2 2 2 4 0 - 4 - 2 2 - 2 2 2 2 2 0 - 2 2 - 2
- 2 1 2 2 0 - 2 - 2 - 1 ] [ 1 2 0 - 2 - 1 2 2 0 - 2 - 2 2 2 2 0 - 2
2 - 2 2 2 4 0 - 4 - 2 2 4 4 2 0 - 4 2 - 4 2 2 4 0 - 4 - 2 2 2 2 2 0
- 2 2 - 2 2 2 0 - 2 - 2 1 2 0 - 2 - 1 ] ##EQU00009##
[0105] Step 3: Fusion of edge images. Conventional Canny algorithm
uses a method of non-maximum suppression (NMS) to process the
smoothed image and make sure each edge is of one-pixel width. In
this method, a fusion of two edge detection techniques is proposed.
The first edge-detected image contains the modified kernel set. The
second edge-detected image contains the same kernel set mentioned
previously in addition to using the modified gradient magnitude.
The fusion of the two edges would guarantee the existence of all
edges that one image may miss over the other one. FIG. 3 shows the
proposed fusion concept.
[0106] There are several algorithms that are developed for edges
fusion. Considering medical images, the criteria for selecting the
fusion operator of interest must preserve edge information in both
images as well as noise detection avoidance. Hence, the average
energy fusion method proposed may be utilized. For two edge images
G.sub.L and G.sub.H, the fused edge image G.sub.F can be defined
as:
G.sub.F,i=.omega..sub.LG.sub.L,i+.omega..sub.HG.sub.H,i
[0107] Where G.sub.L,i is the grey-level value at a given pixel
location Pi in the edge image G.sub.L, G.sub.H,i is the grey-level
value at a given pixel location Pi in the edge image G.sub.H. The
average value of the corresponding coefficients in each edge image
.omega. is defined as:
.omega. = E i E L + E H ##EQU00010##
[0108] Where E.sub.L and E.sub.H represent the low frequency and
high frequency energy of edge image respectively. Such an average
value can judge on the contents of the fused image, and can reduce
the unnecessary details that are contained on one edge image over
the other. On the other hand, the energy E.sub.i contained in an
image at a pixel location P.sub.j can be defined as:
E i = j = 1 N P j N ##EQU00011##
[0109] Application to prostate cancer biopsy images: the algorithm
is simulated on a biopsy image having prostate cancer of Gleason
grad 3. FIG. 4 illustrates simulation results of the algorithm and
shows how the proposed edge detection algorithm succeeds in
detecting edges more precisely compared to canny edge detector.
[0110] As mentioned previously, one might consider integrating a
proper enhancement algorithm as a preprocessing stage for a better
visualization and hence, an effective detection of different edge
details.
[0111] This combination succeeds in providing useful information
about objects and structural information as well as an enhanced
visualization when variation of a constant parameter .alpha. is
considered. In this algorithm, Logarithmic "log" functions replace
the morphological filters operation. That is, the deployment of
PLIP model while the maximum ratio of horizontal adjacent gray
level values is investigated. The new rational expression
becomes:
Rx ( i , j ) = max ( log ( I ( i , j ) ) .alpha. ( I ( i , j - 1 )
) .alpha. , log ( I ( i , j ) ) .alpha. ( I ( i , j + 1 ) ) .alpha.
) ##EQU00012##
[0112] Applying PLIP gray tone function on the row ratio R.sub.x
yields:
gx ( i , j ) = max ( log int ( .alpha. f 3 ) int ( .alpha. f 6 2 )
int ( .alpha. f 9 ) int ( .alpha. f 1 ) int ( .alpha. f 4 2 ) int (
.alpha. f 7 ) , log int ( .alpha. f 1 ) int ( .alpha. f 4 2 ) int (
.alpha. f 7 ) int ( .alpha. f 3 ) int ( .alpha. f 6 2 ) int (
.alpha. f 9 ) ) ##EQU00013##
[0113] Where f denotes pixels in gray tone; int is the integer
operation; a is a predefined constant.
[0114] Similarly, the column ratio becomes:
gy ( i , j ) = max ( log int ( .alpha. f 1 ) int ( .alpha. f 2 2 )
int ( .alpha. f 3 ) int ( .alpha. f 7 ) int ( .alpha. f 8 2 ) int (
.alpha. f 9 ) , log int ( .alpha. f 7 ) int ( .alpha. f 8 2 ) int (
.alpha. f 9 ) int ( .alpha. f 1 ) int ( .alpha. f 2 2 ) int (
.alpha. f 3 ) ) ##EQU00014##
[0115] In order to utilize the above expressions to compute the
magnitude in LIP edge detection, the isomorphic transformation is
applied to have the magnitude components gx and gy. An enhanced
magnitude computation considers the maximum horizontal and vertical
ratios as follows:
Magnitude = M ( 1 - exp ( - max { g _ x M , g _ y M } ) )
##EQU00015##
[0116] Different images were tested in order to verify the
functionality of the introduced algorithm. Images range from
medical images to images of other objects where a different
parameter .alpha. is selected for each image. The introduced
algorithm succeeds in detecting the objects of interest only by
eliminating other unwanted objects that belong to different gray
tone levels. Also, it has the ability to detect unforeseen edges
that no other edge detection algorithm can detect due to the very
slight level of intensity variations. Furthermore, it serves as a
tracking application to track the development of certain gray level
in certain areas.
[0117] Most of the traditional edge detection algorithms share one
objective which is to detect accurate and more edge details. In an
embodiment, the previous algorithms are utilized in image filtering
process. Noisy communication channels and sensors may contribute to
the existence of impulsive noise in images. Such images corruption
can be eliminated by implementing image filtering algorithms.
Various image filtering techniques have been introduced but
median-based filters have proven to attract much attention due to
their ability in preserving edge details as well as simplicity in
removing impulsive noise.
[0118] The following steps summarize the process of image filtering
based on edge detection technique: (1) Edge detection processing
where edge detection operator could be any edge detection algorithm
or one of these described previously in this embodiment including
shape dependent edge detection, Noise-resilient edge detection
algorithm for brain MRI images, or logarithmic-ratio based edge
detection. (2) Locate edge pixels coordinates in order to avoid
such edges during image filtering. (3) Map coordinates to the
original/noisy image in order to save gray level intensity values
at these locations. (4) Detect the existence of impulsive noise on
edges. Detection is based on a rapid change of intensity values
between two adjacent pixels. (5) Apply switching filter algorithm
for image filtering. One filter is a filter of different sizes
m.times.n, e.g. 3.times.3, 5.times.5, 7.times.7 among others square
and rectangular kernels, which is applied depending on impulsive
noise density. The second filter operates on noisy edges which is
weighted filter with less weight at the centers of the impulsive
noise locations. Switching occurs only whenever impulsive noise is
detected on edges at locations determined at step (4)
previously.
[0119] Simulation results of applying the previous algorithm on
different images were studied. Prior to implementing the developed
algorithm, impulsive noise with different densities is added to the
image. The simulation results show a superior performance of the
developed algorithm in terms of preserving edge details while
filtering the impulsive noise.
[0120] In conclusion, an edge preservation filtering technique may
be used for image analysis. The idea behind such a concept is to
avoid blurring edges during the filtering of the images, and hence,
avoiding missing some significant branched edges, especially in
medical application, e.g. biopsy tissue images. Simulation results
showed a superior performance in edge preservation compared to
other filtered images where edges appear to be blurry.
Visual Rational Morphological Filters for Image Enhancement and
Edge Detection
[0121] Considering image contrast enhancement, there are various
methods for image enhancement that operate globally on an image
such as histogram equalization, gamma correction, logarithmic
compression, alpha-rooting, and level or curve-type methods. The
drawback of these algorithms is that they are not always adaptive
and in general, they do not take into consideration local
statistical variations within an image. Also, they do not match the
larger dynamic range the human visual system opposes. Hence, the
need to utilize an image enhancement algorithm where pixels are
processed in local neighborhoods is mandatory. An Edge-Preserving
Contrast Enhancement (EPCE) is developed based on dynamic range
compression. The enhancement part of this algorithm is based on
local mean at each pixel as follows:
I ( x , y ) = 2 1 + - 2 .tau. ( x , y ) / .lamda. ( x , y ) - 1
##EQU00016##
[0122] Where I(x,y) is the output image, .tau.(x,y) is either the V
component of the image in HVS color space or the gray level, and
.lamda.(x,y) provides a measure of the local statistic at a given
window for the adjustment of the transfer function to the local
mean and given by:
.lamda. ( x , y ) = C + ( M - C ) ( .mu. ( x , y ) M ( Qx + N ) )
##EQU00017##
[0123] Where C is an enhancement parameter in the range
0.ltoreq.C<255, M is the maximum value in that range, Q and N
are arbitrary constants. This enhancement algorithm functions in an
adaptive manner by enhancing dark area of the image while the light
area is preserved. In order to enhance edges where a significant
gray level discontinuity exists, the ratio of difference of gray
levels is investigated locally to calculate the ratio magnitude R.
An illustrative example is shown for a 3.times.3 neighborhood as
follows:
I ( i - 1 , j - 1 ) I ( i - 1 , j ) I ( i - 1 , j + 1 ) I ( i , j -
1 ) I ( i , j ) I ( i , j + 1 ) I ( i + 1 , j - 1 ) I ( i + 1 , j )
I ( i + 1 , j + 1 ) ##EQU00018## Rx ( i , j ) = max ( I ( i , j ) I
( i , j - 1 ) , I ( i , j ) I ( i , j + 1 ) ) , Ry ( i , j ) = max
( I ( i , j ) I ( i - 1 , j ) , I ( i , j ) I ( i + 1 , j ) )
##EQU00018.2## R ( i , j ) = Rx ( i , j ) 2 + Ry ( i , j ) 2
##EQU00018.3##
[0124] The above expression is used to emboss edges details and to
yield a new and adaptive dynamic range enhancement as follows:
Ien ( i , j ) = .alpha. ( R ( i , j ) ) 1 + - 2 .tau. ( i , j ) /
.lamda. ( i , j ) - 1 ##EQU00019##
[0125] Where len.sub.(i,j) is the enhanced pixel, R(i,j) is the
ratio coefficient, a is an arbitrary function operating according
to the value of R. The proposed formula succeeds in enhancing and
highlighting edges where gray levels transitions are
significant.
[0126] Recall from the log-ratio approach for morphological
filters. The effect of such ratio is to give more weight to edges
of interest that need to be enhanced to the highest dynamic range
of the HVS based on the local statistics as follows:
I ed ( i , j ) = .beta. ( log .tau. ( i , j ) F ( i , j ) ) I en (
i , j ) ##EQU00020##
[0127] Where .beta. is a normalization coefficient controlling the
enhancement saturation level and acting as a primary edge detection
threshold, F(i,j) is local a filtering operation based on local
neighborhood statistics. A further threshold stage is introduced as
an adaptive process to reflect local statistics by deploying
standard deviation in local windows. Several thresholds may be
used. Particularly, we use a threshold T(i,j) for illustrative
purpose which is calculated as follows [16]:
T(i,j)=.alpha.tm(i,j)+C.sigma.(i,2,5)
[0128] Where .alpha.tm(i,j) is the local alpha-trimmed mean, C is a
global constant, and .sigma.(i,j) is the standard deviation. Once
potential edge gray levels are detected, binarization process takes
place. A summary of the developed algorithm is provided in the
table below:
TABLE-US-00001 Process Algorithm Pre-processing
Gaussian/Median/Order-statistics filtering Ratio magnitude
calculation R(i, j) = {square root over (Rx(i, j).sup.2 + Ry(i,
j).sup.2)}{square root over (Rx(i, j).sup.2 + Ry(i, j).sup.2)}
Adaptive enhancement I en ( i , j ) = .alpha. ( R ( i , j ) ) 1 + e
- 2 .tau. ( i , j ) / .lamda. ( i , j ) - 1 ##EQU00021## HVS-based
edge detection I ed ( i , j ) = .beta. ( log .tau. ( i , j ) F ( i
, j ) ) I en ( i , j ) ##EQU00022## Threshold T Adaptive
thresolding Example: T(i, j) = .alpha. tm(i, j) + C .sigma.(i, j)
Binarization I ed ( i , j ) = { 1 .fwdarw. I ed ( i , j ) .gtoreq.
T 0 .fwdarw. elsewhere ##EQU00023##
[0129] The architecture and characteristics of dark nuclei is an
essential feature in deciding the grade of cancer besides other
glandular features. In this case, the proposed algorithm operates
on the image and succeeds to segment dark nuclei. This result can
be used as an input to feed classification systems of interest such
as fractal analysis based classifier for features extraction.
[0130] In this embodiment, we also present an algorithm for
visualizing edges in colored images as well as edge detection using
the previously developed logarithmic ratio approach. The developed
visualization algorithm has a superior performance in highlighting
more unforeseen color edges than the standard color-grayscale
conversion method can present. Also, by integrating the proposed
algorithm with a logarithmic-ratio based edge detector operator,
the developed algorithm outperforms the standard edge detection
operators in gradient colored images where color boundaries
transitions are hard to detect.
[0131] Most of edge detection operators operate on grayscale images
since color image components impose some sort of complexities in
detecting edges. In grayscale images, an edge can be detected by
investigating the discontinuity in gray levels. However, in color
images there are no clear and defined color edges although color
images carry more information and therefore have the potential of
having more hidden and accurate edges. The limitation factor in
detecting all the edges in gray levels images is due to the fact
that different colors have different hue values but they share the
same intensity value. Hence, an effective approach to convert color
components from several color models such as RGB components into
gray level scale is needed.
[0132] Considering color space, there are different representations
of three dimensional color spaces having different characteristics.
One set of color spaces uses Cartesian coordinates for point's
representation in space. Examples of this type of sets are the
three primary illumination colors Red, Green, and Blue (RGB), the
complementary colors Cyan, Magenta, and yellow (CMY), and the
opponent color representation luminance, in-phase and quadrature
(YIQ). Another space representation utilizes polar/cylindrical
coordinates such as Hue, Saturation, and Lightness (HSL) and Hue,
Saturation and Value (HSV). Generally, each color space can be
converted to another color space representation.
[0133] For color images edge detection, previous researches in this
area utilized vector order statistics, distribution of pixel colors
and fusion techniques of edges detected at different color spaces.
Hence, more focus was given to analyzing and extracting color
information rather than investigating new edge detection
operators.
[0134] Considering real world images, one might consider the
precision of edge maps as a result of human visual system (HVS)
perception. One of the well-known image processing frameworks that
incorporate the characteristics of HVS in image enhancement as well
as edge detection is the Parameterized Logarithmic Image processing
(PLIP) model. In PLIP model, a new set of arithmetic operations are
introduced taking into consideration the characteristics of human
visual system.
[0135] Furthermore, the concept of utilizing ratio operations in
morphological filters design that oppose some characteristics of
HVS was introduced. In previous embodiment, the logarithmic image
processing framework and the concept of rational filters design
were integrated to produce a new edge detection algorithm based on
HVS model. The proposed method provided a new approach in edge
detection focusing on the class of images where the variations and
discontinuities in gray levels follow a gradient pattern. In this
embodiment, we introduce an algorithm for visualizing unforeseen
edges represented by any color model. An example is presented
taking the representation of an image in RGB color space.
[0136] RGB color model consists of three independent color planes
representing the three chromaticity components: Red, Green, and
Blue. In order to specify a pixel RGB value, each color component
is specified independently within the range 0 to 255. RGB color
model is used in many applications such as computers, television,
and electronic systems. Conversion of RGB color space to grayscale
is the common preprocessing step prior to applying edge detection.
The concept of converting RGB to grayscale depends on assigning
weight for each color component based on HVS and its perception of
colors. The common weighting vector used is [L, M, N], where
L+M+N=1.
[0137] Some of the common edge detection techniques for color
images utilize fusion technique of edges detected at different
color space using traditional edge detectors. The first fusion
method finds the gradient and the edge map at each of the three RGB
levels, then fuses the resultant edge maps output. The second
method fuses the gradients at each RGB component, and then finds
the edge map of the final fused gradient.
[0138] The developed algorithm focuses on extracting edge details
that cannot be highlighted for detection using the standard
assigned HVS vector weight. To illustrate the operation of the
algorithm, the RGB color model is considered. Given a tri-chromatic
color image, one might consider designing a rational morphological
filter that takes into consideration the assigned HVS weight for
each component in tri-chromatic space. For instance, for a RGB
image let and Fn and Fm represent two different filters as
follows:
Fn=(L*R,M*G,N*B)
Fm=[.beta.R log.sup..alpha.R(L*R),.beta.G
log.sup..alpha.G(M*G),.beta.B log.sup..alpha.B(N*B)];
[.beta.R log.sup..alpha.R(w.sub.1R-w.sub.2G),.beta.R
log.sup..alpha.R(w.sub.1R-w.sub.2B)];
[.beta.G log.sup..alpha.G(w.sub.2G-w.sub.1R),.beta.G
log.sup..alpha.G(w.sub.2G-w.sub.3B)];
[.beta.B log.sup..alpha.B(w.sub.2B-w.sub.1R),.beta.B
log.sup..alpha.B(w.sub.2B-w.sub.3G)];
[0139] Where .beta. and .alpha. are constants assigned to each
color component, w.sub.i, w.sub.2, w.sub.3 represent different
assigned weight for R,G,B color components respectively. The ratio
of different morphological operations may take any of the following
forms:
.beta. 1 ( max ( Fn ) max ( Fm ) ) ; .beta.2 ( min ( Fn ) max ( Fm
) ) ; .beta.3 ( min ( Fn ) min ( Fm ) ) ; .beta.4 ( max ( Fn ) min
( Fm ) ) ; ##EQU00024##
where .beta. is normalization constant
[0140] The above expressions yield successful results in
investigating the intensity of each RGB component at each pixel
taking into consideration the HVS perception characteristics. The
result of carrying out these different operations produces a gray
level image having more edge details than the standard translation
of RGB to grayscale.
[0141] Image processing algorithms presented herein highlight and
visualize different information that can be obtained from variety
of images including biopsy images as well as real world images. The
developed preprocessing algorithms succeed in detecting, segmenting
and visualizing different image details compared to the state of
the art current image processing algorithms.
Method for Carcinomas Color Region Mapping
[0142] Referring now to another embodiment of the invention, a
method for finding carcinomas dominant colors is provided. The goal
of this method is to reduce the number of distinct colors or define
the most populated color in the biopsy image. The resulting
quantized image should be as visually similar as possible to the
original image for the further tissue structure objects
segmentation. For example we may map a biopsy image containing
millions of colors into an image containing a few colors, which are
closely tied to tissue structures in H&E stained needle biopsy
images are also included in the color map. Additional colors of
interest can be added manually to keep information for further
studies.
[0143] Before performing color region mapping, color
standardization in digital histopathology images is performed.
Various algorithms me be used to make the appearance of the image
similar to those in the training database and tackle the problem of
color constancy. Examples of color normalization algorithm include:
Reinhard's method and its modifications [18], color deconvolution,
histogram matching, among others.
[0144] The general algorithm for constructing color region mapping
comprises the following steps: (1) Survey the colors which compose
biopsy images. (2) Determine the `best` set of colors for a biopsy
image by using algorithms for localization of local maxima in 2-D
color histograms. (3) Generate color table. (4) Assign each color
to a class for segmentation purposes. (5) Map images into the
reduced color space.
[0145] For instance, to construct the color region mapping shown in
the following table, instead of adopting an absolute set of color
features for our segmentation algorithm, the 2-D projections of 3-D
RGB color model of a large variety of prostate needle biopsy images
of different grades are computed to construct 2-D histograms and
generate a good approximation of color region distribution in those
images. This is an important step since the variations in colors of
stained images constitute an aspect to take into account before
applying segmentation procedures based on their color composition,
since these variations are always present even if there are defined
laboratory controls and protocols to prepare tissue core samples.
Based on the calculated 2-D histograms, a color region mapping is
proposed taking the colors that appear more frequently in our image
set and that represent any object of interest in the evaluation of
histopatholy images.
[0146] In more detail, still referring to the carcinomas color
region mapping method, the most prominent colors may found by using
an automatic algorithm for local maxima location in 3-D
surfaces.
[0147] The procedure explained above is valid for finding
representative colors of any image dataset. In this embodiment, the
resulting region mapping for prostate cancerous class of test
images is presented in the next table.
TABLE-US-00002 COLOR COMPONENTS STRUCTURE i COLOR R G B (class) 1
##STR00001## 255 255 255 Lumen and non- tissue areas 2 ##STR00002##
205 146 178 Epithelial cytoplasm 3 ##STR00003## 219 133 168 Stroma
4 ##STR00004## 225 77 109 5 ##STR00005## 159 167 203 Blue mucin 6
##STR00006## 90 37 90 Epithelial nuclei
[0148] It is important to note that any color model can be used to
apply the above method. In this case, color model transformation
must be performed before constructing 2-D projections.
[0149] This embodiment also provides a method for tissue structure
segmentation. Segmentation procedure is carried out through pixel
classification. This automatic process consists of minimizing the
distance between each pixel and the colors considered in the color
map, and assigning to each pixel the class corresponding to the
nearest color. To classify a pixel, an appropriate distance metric
may be used. As an example Euclidean distance is:
d.sub.i(m,n)= {square root over
((r.sub.(m,n)-r.sub.i).sup.2+(g.sub.(m,n)-g.sub.i).sup.2+(b.sub.(m,n)-b.s-
ub.i).sup.2)}{square root over
((r.sub.(m,n)-r.sub.i).sup.2+(g.sub.(m,n)-g.sub.i).sup.2+(b.sub.(m,n)-b.s-
ub.i).sup.2)}{square root over
((r.sub.(m,n)-r.sub.i).sup.2+(g.sub.(m,n)-g.sub.i).sup.2+(b.sub.(m,n)-b.s-
ub.i).sup.2)}
[0150] In the equation above d.sub.i(m,n) is the distance in RGB
color model from pixel (m,n) to the i'th color region defined in
the color map table. FIG. 5 shows the structure segmentation
results based on color information of an image that belongs to
prostate cancer, Gleason grade 3 pattern.
[0151] Method for Pattern Detection
[0152] The novel Gleason grading concept exploits the ability of an
optimal predictor/compressor to compress (in an informational
sense) its predicting class optimally. In this embodiment, an
implementer selects an appropriate machine learning technique and
an appropriate prediction technique (which may be one in the same)
to train a set of predictors on images or pre-processed images
containing a plurality specific Gleason pattern grade. By
pre-processed we imply that operations may have been performed on
the image (such as architectural segmentation/classification,
textural segmentation/classification, feature
segmentation/classification, denoising, scaling, etc.) prior to
predictor training. Subsequent use of the word "image" with regard
to predictors refers to an image that has already been
appropriately pre-processed.
[0153] Each predictor should be trained on a specific type of
Gleason pattern. Possible techniques include: maximum likelihood
estimation, Bayesian estimation, Kalman filtering, prediction by
partial matching, context tree weighting, and others. Paired with
an appropriate entropy encoder, each predictor/encoder pair forms
the best compressor for its respective pattern type compared to the
other predictor/encoder pairs. Hereafter, predictor/encoder pairs
shall be referred to as pattern compressors (PCs) for
simplification. Scoring of a biopsy sample image (hereafter
referred to as the query or query image) amounts to compressing
each query pixel with each PC. Smaller compressed sizes per pixel
imply more significant membership to the pattern class of a PC.
Selection, combination, or other analysis of the all PC compression
rates can then be used to estimate the grade of a particular pixel
or the query image on a whole.
[0154] As an exemplifying embodiment, universal
detector/classifiers (UDCs) constructed from a variable order
Markov model which utilizes a prediction by partial matching (PPM)
algorithm for prediction and an arithmetic encoder for compression
can be used as PCs. Each UDC trains on a given Gleason pattern type
such that each is the best UDC data compressor of that pattern type
on average. Grading of a query sample involves selection,
combination, or other analysis of the classifier's relative
compression performance. Optionally, a PC optimized to the query
image can be used for comparison purposes.
[0155] Method for Hard Gleason Scoring & Mapping
[0156] A "Hard Score" may be assigned to a biopsy image to directly
correspond to the single most-likely pattern observed within a
query image. This is analogous to a pathologist's choice of grade,
where a hard decision must be made between possible Gleason grades.
This is achieved by choosing the predictor/encoder pair which best
compresses a query image on average. The hard score of the sample
is the grade of the pattern for which the predictor/encoder pair
was optimized (see FIG. 23).
[0157] The same hard grading method may be used on a pixel by pixel
basis, assigning a hard Gleason grade to each pixel individually
depending on which trained UDC compresses it best. This results in
a map of pixel grades which can distinguish between regions of
different Gleason patterns within a single biopsy image (see FIG.
24). It should be noted that the average hard map grade is not
necessarily equal to the hard grade as determined above. For
instance, if patterns 3 and 5 are both prevalent within a query
image, then the average hard score is likely to be near 4 despite
grade 4 Gleason patterns being absent within the image. The hard
score, however, will be 3 or 4 depending on which pattern is most
prevalent.
[0158] Method for Soft Gleason Scoring & Mapping
[0159] A biopsy image may also be assigned a "Soft Score," here
denoted s. A soft score can better characterize a biopsy image
which is a member to some extent of one or more predictor/encoder
pair classes. One way to generate such a score is to assign
membership values, or weights w.sub.g, to each possible class found
within an image and to compute the weighted mean. Thus, the score
is a weighted combination of all Gleason grades where the grades
G.sub.i range from 1 to 5 and the weights w.sub.i are determined
from a query image:
s=w.sub.1.times.G.sub.1.sym.w.sub.2.times.G.sub.2.sym.w.sub.3.times.G.su-
b.3.times.w.sub.4.times.G.sub.4.sym.w.sub.5.times.G.sub.5
[0160] In the above equation, the x symbol indicates a scaling
operation defining how much of each pattern grade is within a given
sample. .sym. is a generalized combination operator defining how
the weighted grades should by joined into one, global Gleason
score. Other embodiments may include, but are not constrained
to:
s = G 1 w 1 G 2 w 2 G 3 w 3 G 4 w 4 G 5 w 5 w 1 + w 2 + w 3 + w 4 +
w 5 ##EQU00025## s = w 1 + w 2 + w 3 + w 4 + w 5 w 1 G 1 + w 2 G 2
+ w 3 G 3 + w 4 G 4 + w 5 G 5 ##EQU00025.2## s = arg max ( w g , G
g ) , g = 1 , 2 , 3 , 4 , 5 ##EQU00025.3## s = arg min ( w g , G g
) , g = 1 , 2 , 3 , 4 , 5 ##EQU00025.4## s = PLIP ( w 1 .times. G 1
) + PLIP ( w 2 .times. G 2 ) + PLIP ( w 3 .times. G 3 ) + PLIP ( w
4 .times. G 4 ) + PLIP ( w 5 .times. G 5 ) ##EQU00025.5##
[0161] One embodiment is the form the above equation takes by using
standard multiplication as the scaling operation and addition as
the combination operation. The equation then becomes linear
combination of pattern grades:
s = w 1 .times. G 1 + w 2 .times. G 2 + w 3 .times. G 3 .times. w 4
.times. G 4 .times. w 5 .times. G 5 w 1 + w 2 + w 3 + w 4 + w 5
##EQU00026##
[0162] While traditional binary and tertiary scores prove a score
that is the sum of the weighted grades, the above equation
normalizes the score by the sum of the weights to ensure that the
final measure is a score ranging continuously between 1 and 5.
[0163] The relation to the traditional binary Gleason score can be
described as follows: the traditional binary score assigns a weight
of 1 to the two most prevalent pattern grades and a weight of 0 to
the three remaining pattern grades without normalizing the sum by
the sum of the weights. The modified and previously-mentioned
continuous-scale scores can be similarly defined through
appropriate definition of the pattern weights w.sub.i and the
removal of the normalization factor. The benefit of a normalized
system is that it provides a continuous, threat level grade ranging
from 1 to 5 in correspondence to the traditional Gleason grading
system.
[0164] Often, only grades 3, 4, and 5 patterns are considered in
medical practice. Thus, if we have available to us 3
predictor/encoder pairs, each optimal for a grade 3, 4, or 5
Gleason pattern, then predictor/encoder pair's respective pattern
grade is scaled by its respective weight for averaging:
s = w 3 .times. 3 .times. w 4 .times. 4 .times. w 5 .times. 5 w 3 +
w 4 + w 5 ##EQU00027##
[0165] The weights w.sub.3, w.sub.4, and w.sub.5 are the membership
values determined through statistical analysis of the compression
rates Y.sub.PC obtained for each trained PC on a query image
combined with compression rates obtained by a fourth PC, PC.sub.0,
trained on the query image itself. Compression by PC.sub.0
generates a more optimal self-compression rate which is the amount
of information the query image provides about itself.
[0166] For one embodiment utilizing analysis of variance (ANOVA) as
a statistical comparison tool, confidence interval estimates
.sub.PC for the four PC mean compression rates are t-distributed
around their respective measured mean compression rates, Y.sub.PC,
parameterized by the degrees freedom in the sum squared error
(which depends on both sample size and the number factor levels in
question). When the sample size is large enough, the
t-distributions are practically
.sigma. 2 { Y ^ UDC } = MSE n ##EQU00028##
normal, centered at Y.sub.PC with variance where n is the number of
samples (which is just number of query image pixels) and MSE is the
mean square error as determined through 4-level, single factor
ANOVA testing of the four PCs' compression rates.
[0167] FIG. 27 gives an example of the four PC estimate intervals
as overlapping normal distributions centered at their respective
mean compression rates. The total area of the intersection between
the PC.sub.0 distribution, .sub.4, and a particular Gleason
pattern-trained PC distribution, .sub.g, is a measure of how
similar the curves are to one another and therefore provides a
membership value of the query image to that particular PC and can
be used as a weight. For the example in FIG. 27, the amount of
overlap between the compression rate estimates for PCs 3, 4, and 5
with that of PC.sub.0 are 2.55e-10, 3.51e-4, and 1.23e-4,
respectively. Then, the soft score s is computed:
s = 2.55 e - 10 .times. 3 + 3.5 1 e - 4 .times. 4 + 1.23 e - 4
.times. 5 2.55 e - 10 + 3.5 1 e - 4 + 1.23 e - 4 = 4.26
##EQU00029##
[0168] A 3-level, single factor ANOVA test provides information for
whether or not the three Gleason pattern-trained PCs compress the
query image significantly differently. If they do not, then the
soft score is unreliable and the image either belongs to all
classes equally or is not likely a member of any of the classes (at
least as defined by the training set). In this case, the soft grade
returns "CANNOT BE CLASSIFIED" which may prompt a pathologist to
manually inspect an image and possibly add it to a training set for
the proper PC.
[0169] Like hard Gleason mapping, soft Gleason mapping assigns each
pixel an Gleason grade dependent on the compressibility of a query
image pixel I(i, j) by each PC. However, rather than selecting the
best-compressing PC class, soft grading factors in the performance
of all PC classes to generate an average grade. Thus, the soft
grade of each pixel is a weighted mean of each Gleason class. If
only grades 3, 4 and 5 are considered:
s ij = w 3 , ij .times. 3 + w 4 , ij .times. 4 + w 5 , ij .times. 5
w 3 , ij + w 4 , ij + w 5 , ij , ##EQU00030##
or if a PC of each pattern class is available:
s ij = w 1 , ij .times. 1 + w 2 , ij .times. 2 + w 3 , ij .times. 3
+ w 4 , ij .times. 4 + w 5 , ij .times. 5 w 1 , ij + w 2 , ij + w 3
, ij + w 4 , ij + w 5 , ij . ##EQU00031##
[0170] Computation of the weights w.sub.g,ij may differ
substantially from the method used to compute the global weights
w.sub.g because global image compression statistics may neither be
available for nor relevant to the local grade. However, PC
predictions are based on a pixel's local context and are locally
applicable. In fact, a PC estimates the local likelihood function
defining the probability of observing a particular pixel's value
under the assumption that local the pixels are members of the PC's
class. The observed likelihood models produced by each PC after
observation of a pixel can be compared to the optimal likelihood
function as measured by the a query image's self-PC, PC.sub.0, and
measured as an error.
[0171] A convenient measure of error (and the error for which is
used to prove universality) is the Kullback-Leibler divergence. The
KL divergence measures the amount of extra bits produced using one
probabilistic model instead of a more accurate probabilistic model.
In this case, we compare the likelihood function for a subsequent
observation generated by a PC conditioned on a specific Gleason
pattern to the true likelihood function generated by the PC
conditioned on the query image I after observation of pixel I(i,j).
The KL divergence of PC.sub.g to PC.sub.0 after pixel I(i, j) is
given by
ij PC g = a L PC 0 ( a I ( 1 , 1 ) I ( i , j ) ) .times. log 2 ( L
PC 0 ( a I ( 1 , 1 ) I ( i , j ) ) L PC g ( a I ( 1 , 1 ) I ( i , j
) ) ) ##EQU00032##
[0172] where L.sup.PC.sup.g (a|I(1,1) . . . I(i,j)) is the
likelihood of observing the next pixel with value a conditioned on
the previously observed pixel values if those pixels belong to
grade g tissue and L.sup.PC.sup.0(a|I(1,1) . . . I(i,j)) is the
true likelihood of observing the next pixel with value a.
[0173] The smaller the error, the more likely a pixel belongs to a
PC class. It is reasonable, therefore, to weight by the inverse
error:
w g , ij = 1 ij PC g ##EQU00033##
[0174] Once the weights are computed for each PC level on each
pixel, then the soft map can be calculated using.
[0175] Method for Determining the Most Likely Gleason
Distribution
[0176] An new method provides a measure of the uncertainty about
whether a particular observed query image pixel is of one pattern
grade or another. This uncertainty is here called a transition
probability and represents that probability that a pixel is part of
one specified pattern instead of another specified pattern. This
measurement applied to every pixel within a query image then
provides information about the global structure of the
image--specifically, the uncertainty in the distribution (or
relative population) of pixel grades within a map. Analysis of this
uncertainty combined with the observed soft map distribution allows
for estimation of the most likely Gleason distribution within the
image, which is the amount of each grade pattern within the
image.
[0177] Specifically, the percentages of dominant patterns (as
determined from either hard or soft maps) along with their computed
transition probabilities generate a prediction of the most probable
pattern distribution within the sample.
[0178] Secondary likelihoods are given by the weight of one Gleason
pattern when another is dominant. The total transition probability
from the dominant pattern to the secondary pattern is the sum of
the weights of the secondary pattern where the other is dominant
divided by the total sum of all pattern weights where the other is
dominant Using the 3, 4, and 5 patterns only, the transition
probability T.sub.d.fwdarw.s from a dominant pattern d into
secondary pattern s is then
T d .fwdarw. s = i , j c d , ij w s , ij g , i , j c d , ij w g ,
ij , where ##EQU00034## c dij = { 1 w d , ij > w g .noteq. d ,
ij 0 otherwise , for each d = 3 , 4 , 5 ; s = 3 , 4 , 5 ; and g = 3
, 4 , 5 ##EQU00034.2##
[0179] The transition probabilities from each pattern into another
can be collected into a single transition matrix:
T = [ T 3 .fwdarw. 3 T 4 .fwdarw. 3 T 5 .fwdarw. 3 T 3 .fwdarw. 4 T
4 .fwdarw. 4 T 5 .fwdarw. 4 T 3 .fwdarw. 5 T 4 .fwdarw. 4 T 5
.fwdarw. 5 ] ##EQU00035##
[0180] Define a density of states vector M(k) which contains the
observable membership percentages of each pattern within a biopsy
sample at iteration k. Each iteration describes a more likely
distribution than the last based on the transition probabilities.
If there are n.sub.3 hard grade 3 pixels, n.sub.4 hard grade 4
pixels, and n.sub.5 hard grade 5 pixels in a query image, then the
image's initial density of states is
M ( k = 0 ) = [ % of Grade 3 Patterns % of Grade 4 Patterns % of
Grade 5 Patterns ] = [ n 3 n 3 + n 4 + n 5 n 4 n 3 + n 4 + n 5 n 5
n 3 + n 4 + n 5 ] ##EQU00036##
[0181] A prediction for the density of states at time k is given
by
M(k)=T.sup.kM(0)
[0182] (if and only if the transition matrix is stationary with
respect to k). If all transition probabilities are greater than
zero, then T is ergodic and models an aperiodic Markov chain which
has a unique, steady state solution for any initial density of
states with nonzero entries:
M(.infin.)= V.sub..lamda.=1
[0183] where V.sub..lamda.=1 an eigenvector of T for eigenvalue
.lamda.=1 (which in this case is guaranteed to exist) with entries
that sum to 1. Because the eigenvalue of the eigenvector equals 1,
successive operations of the transition matrix T on V.sub..lamda.=1
can only yield V.sub..lamda.=1, which must be the most likely
distribution of Gleason patterns within a sample comparable to the
query image. Under these assumptions, the most likely distribution
of Gleason patterns within query sample are completely
characterized by observed secondary transition probabilities
irrespective of the initial distribution of patterns observed
within the sample.
[0184] System for Prostate Cancer Gleason Grading/Scoring Based on
Morphology and Architecture of Prostatic Tissue (System 1)
[0185] Referring now to another embodiment of the invention, this
system provides several methods including methods for tissue
objects segmentation, gland reconstruction, morphology and
architectural features extraction, and image classification. The
proposed system proceeds into two phases: (1) the description
(Color, tissue structures, gland and architectural features) that
extracts the characteristics from the biopsy image pattern and (2)
the classification that recognize the Gleason grades by using some
characteristics derived from the first phase by implementing
learning algorithms. The flow diagram of the proposed system is
presented in FIG. 7.
[0186] Referring now to this embodiment of the invention in more
details, the image-level features include, color information such
as color distribution entropy and 2-D color histograms. To compute
2-D histograms feature vector a color model transformation may be
perform. This transformation converts each RGB image into its
equivalent form in a suitable color space; for instance, L*a*b or
HSV (Hue, Saturation and Intensity) color models. Next, 2-D
histograms are computed using n bins at each dimension. The number
of bins may be selected by various heuristic methods, for example,
cross-validation.
[0187] In further detail, still referring to the morphology- and
architectural-based grading system, tissue structure are segmented
using any suitable segmentation method, for example the color
region mapping proposed in another embodiment of this invention.
Having segmented tissue structure and gland units the following
features may be extracted from the objects: First order statistics
of morphology of lumen and gland such as size and shape, gland
boundaries features, nuclei density at several scales, nuclei area
and circularity. In addition, architectural features from graphs,
for example Voronoi, Delaunay and Minimum Spanning Tree (MST)
graphs constructed from color primitives (after segmentation
procedures) may be extracted. Such architectural features include
but are not restricted to number of nodes, statistics and entropy
of edge length and polygon area, eccentricity, fractal dimension,
among others. Once the feature extraction process is finished,
algorithms for feature selection and ranking may be implemented to
choose the features that best represent the tissue for
classification purposes. The selected features are processed such
that they are scaled according to their importance.
[0188] Referring now to the image classification step, several
learning algorithms may be used. The classifiers may be decision
trees, k-nearest neighbor, Bayesian, neural networks, and/or SVM.
Boosting algorithms such us AdaBoost, SVM Boost, etc. can be used
to improve the classification performance. Next table shows the
classification accuracy ranges of the multiclass recognition task
for grading prostate carcinomas regions of interest.
TABLE-US-00003 Gleason Correct classification rate pattern for
different classes of images Grade 3 93.33%-99.9% Grade 4
86.67%-99.9% Grade 5 100%
[0189] Textural-Based System for Prostate Cancer Gleason Grading
and Scoring (System 2)
[0190] In this embodiment, transform is considered for transform
(wavelets, cosine, Fourier, and others) features extraction for
classification of prostate cancer biopsy images. The flow diagram
of this system for a wavelet transform is outlined in FIG. 8.
Images are preprocessed using some of the previously discussed
preprocessing techniques in the previous embodiments prior to image
transform and features extraction. Features are classified using
any suitable learning algorithm for example Support Virtual Machine
(SVM). Classification performance showed that Haar wavelet
transform produced an improved result in terms of energy features
extractions when images are first preprocessed using the concept of
morphological filtering and transform coefficients ratio.
[0191] Wavelet Transform: In wavelet transform, the signal is
decomposed into a family of orthogonal .psi..sub.m,n(x) and
.psi.(x) as follows:
.psi..sub.m,n(x)=2.sup.-m/2.psi.(2.sup.-mx-n) for integers m and
n
[0192] For the construction of the mother function .psi.(x), the
determination of a scaling function .PHI.(x) is needed to satisfy
the two-scale difference equation:
[0193] This yields the wavelet kernel .psi.(x) as follows:
.phi. ( x ) = 2 k h ( k ) .phi. ( 2 x - k ) ; ##EQU00037## .psi. (
x ) = 2 k g ( k ) .phi. ( 2 x - k ) ; ##EQU00037.2##
where g(k)=(-1).sup.k h(1-k)
[0194] For wavelet decomposition to a Jth-level, the relation of
the decomposed coefficients to direct previous coefficients can be
given as:
c.sub.j+1,n=.SIGMA..sub.kc.sub.j,kh(k-2n) and
d.sub.j+1,n=.SIGMA..sub.kc.sub.j,kg(k-2n) where
0.ltoreq.j.ltoreq.J
[0195] For 2-D wavelet transform, the basis functions operate along
the horizontal and vertical directions as a pair of low pass and
high pass filters, h and g, respectively. The coefficients at the
HH, HG, GH, and GG correspond to the Approximated, Horizontal,
Vertical, and Detailed submatrices in the 2D transform space.
[0196] To illustrate the use of a transform in this embodiment, the
Haar wavelet transform is shown as an example. The Haar basis
function can be expressed as follows:
H ( 2 ) = ( + 1 + 1 + 1 - 1 ) [ Haar ] 2 n = H ( 2 n ) = ( H ( 2 n
- 1 ) ( + 1 + 1 ) 2 n - 1 I ( 2 n - 1 ) ( + 1 - 1 ) ) , n = 2 , 3 ,
, ##EQU00038##
[0197] Many different features can be extracted from wavelets and
orthogonal transforms such as cosine, sine, Hartley, Fourier, and
others. Some exemplary features may include the normalized
coefficients of energy and entropy features for each decomposed
submatrix at a given level as follows:
energy = i j x ij 2 nxn ##EQU00039## entropy = - ( 1 nxn ) i j ( x
ij 2 norm 2 ) log ( x ij 2 norm 2 ) ##EQU00039.2##
[0198] Where x.sub.ij represents wavelet coefficients; n is the
dimension of the submatrix, and
norm 2 = i j x ij 2 ##EQU00040##
[0199] Due to the diversity of staining colors at each image, the
performance of these features should be evaluated against each
color component of any color model, for example RGB color. Energy
and entropy features are then extracted for each submatrix and
classified using the Support Vector Machine (SVM) classifier as an
example. Features in the training set contain one class label and
several multiple features. The learning objective is to sort the
training data into groups/classes so that the degree of association
is strong among the features of the same class and weak among other
members of other classes. The set of features in our case represent
the wavelet-based feature vectors.
[0200] As discussed previously, there is significant classification
performance diversity among the different grades. In order to
improve the classification performance, the concept of deploying
the Human Visual System (HVS) rational filters is considered in
both, the spatial domain of color channels, and in the wavelet
transform domain.
[0201] In previous embodiment, we developed an algorithm for
visualizing and detecting edges in a colored image using
logarithmic approach. Some of the proposed morphological filters
include the following:
.beta. 1 ( max ( Fn ) max ( Fm ) ) ; ##EQU00041## .beta. 2 ( min (
Fn ) max ( Fm ) ) ; ##EQU00041.2## .beta. 3 ( min ( Fn ) min ( Fm )
) ; ##EQU00041.3## .beta. 4 ( max ( Fn ) min ( Fm ) )
##EQU00041.4##
[0202] Where .beta.n is a constant used for grayscale
normalization.
[0203] For the class of biopsy images, taking as an example the RGB
color model, the filter considers the logarithmic ratio of the
minimum RGB component to the maximum one. This type of filters has
the ability to segregate the dark basal cells from other tissue
characteristics. It operates on biopsy images as follows:
I new = .beta. log 10 ( min ( RGB ) max ( RGB ) ) .alpha.
##EQU00042##
[0204] Where I.sub.new is the output processed/filtered image.
[0205] In another embodiment, the concept of rational filters in
the transform domain is used. The wavelet transform provides a good
localization property in space-frequency domain. Since the most
significant coefficients information exist mainly in the
approximated and detailed coefficient sub matrices, the ratio of
these two matrices is taken, then the energy of the coefficients
ratios is considered in classification.
[0206] A transform approach for image quality measurement using the
properties of HVS has been proposed. For any wavelet transform, the
rational filter coefficients for the low-low "Approximated" and
high-high "Detailed" frequency filter banks in a 1-D signal can be
given as follows:
C ratio 1 D = i Ci A i Ci D ##EQU00043##
[0207] Similarly, the coefficients ratio, C.sub.ratio, for the 2-D
case is calculated as follows:
C ratio 2 D = i j Cij AA i j Cij DD ##EQU00044##
[0208] Another proposed method is to use the phase information of
the Lo and Hi coefficients ratio, which is expressed for 1-D signal
as:
.theta. i = tan - 1 ( i C i Lo i C i Hi ) ##EQU00045##
[0209] Features extracted from the proposed method are fused with
the wavelet energy feature. The following table tabulates the
classification performance for different classes of prostatic
carcinomas by employing the explained approach.
TABLE-US-00004 Gleason Correct classification rate pattern for
different classes of images Grade 3 77.0%-99.9% Grade 4 85.3%-99.9%
Grade 5 96.0-100%
[0210] The classification performance using the proposed method
succeeded in outperforming the previous method where energy feature
is extracted from some color channels of the biopsy images
directly.
[0211] Another extracted feature is related to the fractal
dimension (FD) measure. One embodiment of the invention provides
two developed algorithms in fractal analysis. The developed
algorithms are based on edge detection information as a
preprocessing stage as well as quad-tree decomposition algorithms.
The classification performance resulting from adopting such
pre-processing algorithm showed an improved performance in some
Gleason grades over the classical implementation of the box
counting method.
[0212] Methods for Computation of Fractal Dimension by Improved
Box-Counting Algorithm
[0213] The architecture and geometrical features of biopsy images
of prostate presents a rich domain for fractal analysis through the
concept of fractal dimension (FD). The fractal dimension provides a
geometrical description of an image by measuring the degree of
complexity and how much space a given set occupies near to each of
its points. The growth of cancer through the development of
different architectures associated with different grades shows
distinct features that can be used for classification using fractal
dimension. Malignant cancerous cells demonstrate darker intensity
gray level values compared to the surrounding tissue parts. In this
part, fractal features through Differential Box Counting (DBC) will
be used.
[0214] Given a set S in Euclidean n-space, the self similarity of S
is obtained if Nr distinct copies of S exist by r ratio. The
following equation provides the fractal dimension measure "D" as
follows:
D = log ( N r ) log ( 1 / r ) ##EQU00046##
[0215] To find the box-counting dimension, the image is binarized
by using any proper approach in order to have two levels of
intensities, either 0 or 1, rather than 256 gray levels
intensities. To find a suitable binary image for fractal dimension
computation several algorithms such as the proposed logarithm edge
detection, gray level thresholding, adaptive thresholding, and
bit-plane decomposition may be used. Then, the binary image is
divided into grids of size r and the number of grid boxes Nr
exhibiting self similarities is counted. The next step is to reduce
the grid size, r2 and find the corresponding number of boxes, Nr2.
Grid sizes are reduced in a pattern following 2n. For instance, a
512.times.512 image yields the following grid sizes, 512, 256, 128,
64, 32, 16, 8, 4, and 2 grid sizes. The resulting fractal
dimensions for each grid size are fused together to obtain a unique
fractal dimension.
[0216] Also, in this embodiment we proposed an edge detection-based
fractal dimension measure (EDFDE). It is utilizing the concept of
edge detection using parameterized logarithmic image processing,
which was developed in previous embodiments. Biopsy images go
through PLIP-based edge detection prior to box counting.
[0217] In summary, feature vectors were extracted from wavelet and
fractal analysis domain. Each extracted feature vector represents a
different architectural property used for classification purposes
feature vectors from wavelet transform provided a high performance
for classifying and grading biopsy images. Before classification
step, feature selection is performed. The objective is to select
features sets that best describes the statistical property of the
biopsy images classification features. There are various algorithms
for features selection and validation of features prior to
classification in literature. For example, the minimum-Redundancy
Maximum-Relevance (mRMR) features selection method may be used.
Using the mRMR method, the features space (N-features) is reduced
to M-features (N>M).
[0218] Referring in more details to the image classification step,
several learning algorithms may be used. The classifiers may be
decision trees, k-nearest neighbor, Bayesian, neural networks,
and/or SVM. Boosting algorithms such us AdaBoost, SVM Boost, etc.
can be used to improve the classification performance. The table
below shows an example of classification of images from different
Gleason grades and highlights the final SVM classification
performance using "leave-one-out" cross-validation.
TABLE-US-00005 Gleason Correct classification rate pattern for
different classes of images Grade 3 93.0%-99.9% Grade 4 89.3%-99.9%
Grade 5 96.0-100%
[0219] Textural-Based System for Prostate Cancer Gleason Grading
and Scoring (System 3)
[0220] In this embodiment, the invention provides a system with
several methods for automatic classification of prostate cancer
carcinomas by using wavelet-based and fractals-based textural
features. Image-level features are extracted from wavelet
decomposition and color fractal analysis of images, as well as
object-level fractal dimension features as it is shown in FIG.
9.
[0221] Referring now to the embodiment in more detail, the
following features are extracted from the cancerous images:
[0222] Transform-Based Features:
[0223] In order to compute the feature vectors, real/complex
wavelets or multi-wavelets or other orthogonal transforms (e.g.
Fourier, cosine, sine, Hartley, Hadamard) may be used in general.
Several wavelet transforms are widely used for quantitative image
processing since they provide an effective tool for
multi-resolution analysis. A particular approach may include
parameterized Slant-Haar transform [19], which is described
below:
[0224] Slant-Haar wavelet is the wavelet transform of choice for
this example. The slant-Haar transform of order 2.sup.n for n=3, 4,
5 . . . , is defined by:
X=S.sub.2.sub.nx
[0225] Where x is an input data vector of size n.
S 2 n = S 2 n ( .beta. 4 , , .beta. 2 n ) ##EQU00047## S 2 n = Q 2
n [ A 2 S 2 , 2 n - 1 I 2 S 2 n - 2 , 2 n - 1 ] ##EQU00047.2##
[0226] For -2.sup.2n-2.ltoreq..beta..sub.2.sub.n.ltoreq.2.sup.2n-2
and n.gtoreq.3
[0227] Where S.sub.2,2.sub.n-1 is a matrix of dimension
2.times.2.sup.n-1 comprised of the first two rows of
S.sub.2.sub.n-1 and S.sub.2.sub.n-1.sub.-2,2.sub.n-1 is a matrix of
dimension (2.sup.n-1-2).times.(2.sup.n-1) comprised of the third to
the (2.sup.n-1) th rows S.sub.2.sub.n-1. The symbol denotes the
operator of the Kronecker product. The matrix A.sub.2 is given
by:
A 2 = 1 2 [ 1 1 1 - 1 ] ##EQU00048##
[0228] And the matrix S.sub.2.sub.n-1 is the recursion kernel
defined by:
Q 2 n = [ 1 b 2 n a 2 n a 2 n - b 2 n 1 0 1 ] ##EQU00049## a 2 n =
3 ( 2 2 n - 2 ) 4 ( 2 2 n - 2 ) - .beta. 2 n , b 2 n = ( 2 2 n - 2
) - .beta. 2 n 4 ( 2 2 n - 2 ) - .beta. 2 n ##EQU00049.2##
[0229] For -2.sup.2n-2.ltoreq..beta..sub.2.sub.n.ltoreq.2.sup.2n-2
and n.gtoreq.3
[0230] Another approach may include parameterized PLIP Slant-Haar
transform. In PLIP Slant-Haar wavelet, the arithmetic operations in
wavelet definitions are replaced by PLIP operations [14]:
g 1 + g 2 .fwdarw. g 1 .sym. g 2 = g 1 + g 2 - g 2 g 2 .gamma. ( M
) ##EQU00050## g 1 - g 2 .fwdarw. g 1 .THETA. g 2 = k ( M ) g 1 - g
2 k ( M ) - g 2 + ##EQU00050.2## g 1 g 2 .fwdarw. g 1 * g 2 = .PHI.
- 1 ( .PHI. ( g 1 ) .PHI. ( g 2 ) ) ##EQU00050.3##
[0231] Where k(M), .gamma.(M) are arbitrary linear functions of the
type .gamma.(M)=AM+B for constants A and B. The parameter .epsilon.
is a very small constant included to avoid division by zero in the
case k(M)=g.sub.2.
[0232] Referring now in more detail to the features which can be
extracted from wavelet coefficients, the following features are
considered: [0233] Statistics of wavelet detail coefficients: Using
a wavelet transformation of the three color components of the image
in one or more color models at a specific scale, the first order
statistics of scaled coefficients of each sub-band are computed and
arranged as a feature vector. Those statistics may include average,
median, standard deviation, high order moments, among others.
[0234] Phase information of approximation and detail coefficients:
After transforming the image using any suitable wavelet basis, the
phase information generate by the approximation and detail
coefficients may be computed as follows:
[0234] .theta. i = tan - 1 ( i C i Lo i C i Hi ) ##EQU00051##
[0235] Joint Probability distribution of wavelet detail
coefficients: The color input image is first decomposed into its
RGB components. Before applying 2-D wavelet filters, color model
transformation can be performed. Next, wavelet decomposition is
performed for each channel separately and joint probability of
scaled detail coefficients is computed for each possible pair of
channels. The joint distribution consists of a set of bins which
counts the number of wavelet details coefficients at a given scale
falling in a given area of the plane:
[0235] p Chi , Chj ( x , y ) = k x y AND ( Ch i = x , Ch j = y )
##EQU00052## .A-inverted. i .noteq. j ; ##EQU00052.2## i , j = 1 ,
2 , 3 ##EQU00052.3##
[0236] In equation 3, i,j=1, 2, 3 indicate the three channels of a
color input image. x, y represent the numerical value of the
coefficients which range, in our experiments, from 0 to 255. The
logical AND function will return 1, when both arguments are true,
and 0 otherwise in order to count the number of coefficient falling
in a given bin. [0237] Wavelet Generalized Energy Distribution: In
order to capture how the wavelet energy is distributed within the
images, each input image (gray scale image or R-G-B components) is
divided into non-overlapping s.times.s blocks. A wavelet
transformation, for example Haar transform, is then applied to each
block independently. The wavelet transform converts the data array
into a series of wavelet coefficients, each of which represents the
amplitude of the wavelet function at a particular location within
the array and for a given wavelet scale. The generalized energy of
approximation sub-band can be computed as:
[0237] E i A = j = 1 N A ij p , i = 1 , , K , P = 1 , 2.3 ,
##EQU00053##
[0238] Here, K is the wavelet decomposition level. N is the number
of the coefficients in the approximation sub-band at each
decomposition level. First order statistics, i.e. mean, standard
deviation, median, minimum to maximum ratio, etc, of the energy
distribution are computed to produce feature values.
[0239] According to an embodiment of this invention, the sub-band
coefficients are processed by a window function, wherein the center
coefficient x.sub.ij is replaced by the result of the function
f.sub.w(x.sub.i,j.fwdarw.m.times.n) in a small neighbor w of size
m.times.n, m,n=1, 2, 3, . . . .
[0240] For example, the function for a 3.times.3 neighbor could
be:
f.sub.w(x.sub.i,j.fwdarw.3.times.3)=.SIGMA.[x.sub.i,j(x.sub.i-1,j-1.sup.-
.alpha..sup.1x.sub.i-1,j.sup..alpha..sup.2x.sub.i-1,j+1.sup..alpha..sup.3x-
.sub.i,j+x.sub.i,j.sup..alpha..sup.4x.sub.i,j.sup..alpha..sup.5x.sub.i,j-1-
.sup..alpha..sup.6x.sub.i+1,j-1.sup..alpha..sup.7x.sub.i+1,j.sup..alpha..s-
up.8x.sub.i+1,j+1.sup..alpha..sup.9)].sup.2
[0241] Where .alpha..sub.i are user-defined constants, for example
0, 1, 2, . . . .
[0242] Or,
f w ( x i , j .fwdarw. 1 .times. 1 ) = x i , j 2 , f w ( x i , j
.fwdarw. 1 .times. 1 ) = Log 2 x i , j , f w ( x i , j .fwdarw. 1
.times. 1 ) = x i , j Log x i , j ##EQU00054## f w ( x i , j
.fwdarw. 1 .times. 2 ) = ( x i , j x i , j - 1 ) 2 ##EQU00054.2##
.theta. i = tan - 1 ( i C i log ( Lo ) i C i log ( Hi ) )
##EQU00054.3## .theta. i = tan - 1 ( i C i log ( Lo ( x i , j x i ,
j - 1 ) ) i C i log ( Hi ( x i , j x i , j - 1 ) ) )
##EQU00054.4##
[0243] The features above can be computed using the transformed
wavelet coefficients. In the case complex wavelet are used,
x.sub.ij refers to the magnitude of the coefficients.
[0244] Method for Computation of Color Fractal Dimension
[0245] In accordance with an embodiment of this invention, a new
algorithm for computing fractal dimension of color images is
proposed. The modified algorithm is based on a color space
transformation of the histopathology images. To this end, a new
modification in the (R,G,B) color space is used. The proposed
algorithm extracts a single fractal dimension of a color image or a
fractal dimension vector corresponding to the color fractal
dimension of non-overlapping blocks of size m.times.n. In this
algorithm the color image is considered a set of points in a
5-dimensional space S. Each point s.epsilon.S is represented by its
two spatial dimensions and its three color components in the new
color model s.sub.i=(x,y, f.sub.1(R,G,B,.alpha..sub.1),
f.sub.2(R,G,B,.alpha..sub.2), f.sub.3(R,G,B,.alpha..sub.3)) where
f.sub.i is a function of 3 colors and some set of constants
.alpha..sub.i. The function should be adjusted to enhance the
separation between cancer grades. Various examples of the color
transformation are shown in FIG. 10. Color fractal dimension D of
an image can be approximate by:
N ( L ) = g ( L ) m = 1 N 1 m P ( m , L ) .infin.L - D
##EQU00055##
[0246] Where, N(L) is the total number of boxes (size L) needed to
cover the image, and it is multiplied by a function g(L), such as
g(L)=1, L, L.sup.2, log L, L log L. P(m,L) is the probability
matrix having m points included into a hyper-cube (or box) of size
L centered at an arbitrary point of S. Below, the probability
matrix calculation procedure is outlined:
[0247] 1. Count the number of points s.epsilon.S,
s.sub.i=(x,y,f.sub.1(R,G,B,.alpha..sub.1),f.sub.2(R,G,B,.alpha..sub.2),f-
.sub.3(R,G,B,.alpha..sub.3)) with
[0248] d (F,Fc).ltoreq.L in which d (F,Fc) is a distance between
the center pixel Fc and the rest of pixels F in the box, and
f.sub.i, f.sub.ci, is the value of the i-th dimension of point
s.epsilon.S in the 5-D space s=(x,y,r',g',b'). For example, the
distance d (F,Fc) could be:
[0249] The Minkowski distance
d(F,Fc)=|F-F.sub.c|=max(|f.sub.i-f.sub.ci).A-inverted.i=1,2 . . .
,5
[0250] Or, p-norm distance
d ( F , Fc ) = F - F c = ( i = 1 5 f i - f ci p ) 1 / p
##EQU00056##
[0251] 2. For an arbitrary size L normalize the matrix P(m,L) by
using the following relationship
m = 1 N P ( m , L ) = 1 ##EQU00057##
in which N is the number of pixels included in a box of size L.
[0252] Finally, the fractal dimension D is found by the linear
fitting of the negative of the logarithm of L against the logarithm
of its corresponding. The best function g(L), can be found
using
G = max g { min i , j ( d i , j ) } ##EQU00058##
[0253] Where d.sub.i,j is the distance between pairwise averages
A.nu..sub.i with i=3, 4, 5, of the fractal dimension of the set of
class of images of Gleason grades 3, 4 and 5 calculated for each
function g(L).
[0254] Simplifications of the invented algorithm result in the
color fractal dimension introduced by Ivanovici 20].
[0255] Referring now to the image classification step, several
learning algorithms may be used. The classifiers may be decision
trees, k-nearest neighbor, Bayesian, neural networks, and/or SVM.
Boosting algorithms such us AdaBoost, SVM Boost, etc. can be used
to improve the classification performance. Next table shows the
classification accuracy ranges of the multiclass recognition task
for grading prostate carcinomas regions of interest.
TABLE-US-00006 Gleason Correct classification rate pattern for
different classes of images Grade 3 99.9%-100% Grade 4 86.67%-99.9%
Grade 5 100%
3-D Gleason Grading System (System 4)
[0256] In this embodiment, 3-D image is obtained using a single 2-D
tissue image to assign Gleason grade and score to the third
dimension of the reconstructed image. The building blocks of this
part are shown in FIG. 11 and will consider the following
steps:
[0257] (1) Mapping and pre-processing: feature vectors algorithms
related to wavelet transform, and fractal analysis will be
considered. The initial step considers mapping 2-D images into
slices of layers using pre-developed mapping algorithms in order to
have a 3-D voxel representation of a single 2-D image.
Preprocessing algorithms will consider: Haar-HVS based algorithms,
logarithmic RGB model representation, and edge detection-based
algorithms.
[0258] (2) Feature vectors normalization: feature vectors generated
from above algorithms will be normalized.
[0259] (3) Feature Selection and Classification: feature selection
algorithms may be deployed to select dominant features generated as
described above. Several classifiers can be used for recognition.
For example SVM classifier may be used to classify such features
and generate a Gleason score.
[0260] (4) Integration and validation: the final output of this
subsystem may be verified and integrated with the 2-D Gleason
subsystem.
[0261] Referring now in more details to this embodiment of the
invention, the mapping of a single 2-D image into a 3-D space may
be done by using several mapping algorithms. The used algorithm are
grouped in depth cues such as defocus, linear perspective, shading,
patterned textures, symmetric patterns, occlusion, statistical
patterns, among others.
[0262] System for Whole-Slide Processing and Prostate Cancer
Diagnosis: Detection, Grading, Mapping, Scoring (System 5)
[0263] A computer-aided diagnosis system addresses the detection,
Gleason grading, scoring and mapping of prostate carcinoma in whole
digitized slides. The system follows a multi-resolution approach,
in which specific features are extracted at a specific image
magnification. For example: (1) at low resolution color and
textural features are extracted, (2) at medium resolution,
architectural tissue analysis is performed via graph features, and
(3) At high resolution, the morphological and wavelet features are
used for tissue characterization. The system proceeds as follows:
The biopsy specimen is segmented from the background to save
calculation time, a square or rectangular grid of size m.times.n is
superimposed for patch-based cancer diagnosis (the size of the grid
could change depending on the resolution), spatial and transform
domain (multi-resolution analysis) features are extracted and block
classification is performed by using any of the systems of the
previous embodiments.
[0264] A classification system according to an embodiment of the
invention uses multiple features for classification of Whole-Slide
and several computer algorithms which produce the following final
outputs: tumor extent in the needle biopsy, number and
identification of positive needle cores, percentage of positive
cores for cancer, length of the tumor in positive cores and length
of the core, percentage of needle core tissue affected by carcinoma
and area of tissue positive for cancer, Gleason grade of each
selected cancerous region of interest, percentage of each Gleason
pattern in biopsy specimens, Gleason score of each positive core,
including tertiary pattern information, the most frequent Gleason
score in slides, numbers of identical Gleason scores in the
positive cores, and localized cancer map including the grade of
each cancerous patch in all needle biopsy slides as shown in FIG.
12.
[0265] In accordance to this embodiment, the classification
algorithm uses morphology-, textural-based classification
subsystems or a combination of those systems provided in previous
embodiments of the present invention. As a result, a grade of the
Gleason grade scale is assigned to each cancerous patch in the
whole slide and the Gleason score is computed automatically. In
addition, the system is able to detect tertiary patterns present in
the biopsy image. The percentage of extent of each grade found in
the image is also computed.
[0266] The system of studying biopsy samples may be an automated
procedure to process large number of biopsy slides by integrated a
production line model. The automated system output automatically
classifies 2-D and 3-D Gleason grades by using serial
sections/slices of prostate core biopsy images. The automated
system output automatically recognize de Gleason grade, compute the
Gleason score and predicts the likelihood of prostate cancer
reoccurrence after treatment, as well as more accurate risk
assessment, and the decision on the necessity of further biopsy
samples. The automated system output allows improved and efficient
treatment planning.
[0267] In an embodiment, the method and system allow a surgeon to
be immediately aware of whether the surrounding tissue of a
specimen completely surrounds the tumor, prior to the time that the
patient is closed up.
[0268] 3-D Core Grading/Scoring System (System 6)
[0269] Digital 3-D reconstruction of tissue has many uses in
medical research and it can be also used for cancer diagnosis
because more features can be analyzed from high-resolution 3-D
images. For example, additional information regarding the volume of
tumor areas, the distribution of cells, the branching of blood
vessels in a tumor, etc can be extracted to make better and more
accurate cancer diagnosis and prognosis.
[0270] The building blocks of this embodiment will consider the
following:
[0271] (1) Acquiring images: serial sets of images (M sets from N
patients with 12 slices each) need to be pre-labeled by a multiple
experienced pathologists with Gleason score ranging from 3+3 to 5+5
for learning algorithm training.
[0272] (2) Building 3-D model: building a 3-D image from a set of
slices for preprocessing and featuring vectors generation in the
next step.
[0273] (3) Preprocessing and feature vectors generation:
preprocessing algorithms presented earlier in this proposal are
used for the investigation of the core pattern and a possible
feature vectors generation. The extracted features may include:
shape and structure of cells, graph features, `micro-architecture`
of blood vessels, and tumors' volume, shape, regularity, among
others.
[0274] (4) Classification and validation: in this part, a Gleason
grade and score are given to the 3rd dimension of the core based on
SVM classifier or other suitable classifiers with different
validation methods. The labeled grades are investigated and
verified with a pathologist to insure validation of the grade and
its correlation with the 2-D grading subsystem.
[0275] In accordance to an embodiment of this invention, the image
acquisition process may be performed by using a photon-light, a
stereo-light microscope or an electron microscope. The selection of
the microscope should be done according to the required image
resolution.
[0276] To construct 3-D image from 2-D slices, different software
tools can be used. For instance, ImageJ, Fiji, or developed
algorithms in Maltab, among others. The general flow diagram for
3-D reconstruction is shown in FIG. 13.
[0277] According to another embodiment of the invention, tissue
image preprocessing and structure segmentation may be performed
using 2-D slices. FIG. 14 shows an example of 3-D reconstruction of
segmented lumens and epithelial nuclei from a prostate cancer
biopsy image. Feature extraction may be performed from both 2-D and
3-D images.
[0278] According to another embodiment of the present invention,
Gleason grading and scoring is performed in the 3-D tissue image
from longitudinal and transverse views (see FIG. 15). The total
volume of the tumor is also calculated. In this part, the Gleason
score of the tumor will be computed as the combination of the most
predominant grades in both directions. For example, the most
predominant grade could be the average or mode of the grades which
were found in the longitudinal and transverse tissue evaluation. If
the distribution is bimodal or multimodal, the first component of
the 3-D Gleason grade will be the maximum grade among the modal
values and the second component will be the next grade value.
Ternary patterns may be included in the final Gleason score,
especially if those patterns correspond to a high Gleason
grade.
[0279] System for Whole-Slide Processing and Prostate Cancer
Diagnosis: Detection, Grading, Mapping, Scoring (System 7)
[0280] A computer-aided diagnosis system addresses the detection,
Gleason scoring and mapping of prostate carcinoma in whole
digitized slides. The system follows a prediction/compression (PC)
approach, in which predictors are trained to optimally compress
different Gleason patterns. For example: (1) segment the query
image (color, tissue structures, gland and architectural features)
into pattern-defining characteristics, (2) reduce query to quarter
scale using a wavelet transform approximation, (3) unravel 2-D
query into 8 distinct space-filling curves, (4) train a universal
predictor/compressor (UDC) on each Hilbert curve, (5) train 3 UDCs
on collections of already-graded 3, 4, and 5 patterned, archival
images, (6) compress the query image Hilbert curves with each UDC,
(7) statistically analyze compression results, and (8) grade query
image and pixels with Hard and Soft scores computed from the
statistical results.
[0281] Biopsy image preparation is necessary for formatting image
data into a structure usable by the novel universal
compressor-detectors (UDCs). Images need to be simplified and
reduced to support fast data analysis and also need to be
one-dimensionalized for construction and implementation of a UDC.
Specifically, image data must support two phases of UDC action:
[0282] 1. Training and
[0283] 2. Detection.
[0284] Training builds UDC structures on pathologist-graded,
one-dimensionalized, Gleason pattern images. Detection involves
prediction and compression operations by a UDC on a similarly
prepared query image.
[0285] Thus, image preparation for both training images and query
images consists of three significant steps: (1) structural
segmentation for image simplification, (2) image downscaling for
data reduction, and (3) 1D conversion for UDC support. While
different approaches may be taken to perform these steps, the
methods chosen here are designed to facilitate computational speed
while simultaneously maintaining the robustness of the system.
[0286] First, 512.times.512 pixel, H&E-dyed biopsy images are
segmented into six significant structural elements including:
lumen, epithelial cytoplasm, stroma, blue mucin, and epithelial
nuclei. The method simultaneously extracts the important features
found in the various Gleason patterns and reduces the dynamic range
of pixel values (originally 0-255 in each of the red, green, and
blue channels) to 0-5, effectively reducing the image alphabet size
from the approximately 1.68.times.10.sup.7 color combinations of an
RGB image to 6, which is a very practically sized alphabet for a
UDC to handle.
[0287] Images are also reduced to quarter scale (64.times.64
pixels) in order to reduce the number of pixels from 262,144 per
image to 4,096 per image. This action preserves significant
regional statistics within the sample images while alleviating the
amount of data on which a UDC must train or compress. Scale
reduction is achieved using a 2-level integer Haar wavelet
transform to generate a 64.times.64 approximation of the original
512.times.512 segmented images and to preserve the dynamic range of
the alphabet. Other wavelets could be used for the same task;
however, the Haar is chosen because of its ability to be computed
with extreme speed.
[0288] Lastly, the scaled, segmented image is unwrapped into one
dimension to support use with a UDC. For example, this can be
achieved by using several Hilbert curve raster scans. A Hilbert
curve is a fractally-generated, space-filling curve on integer
powers of 2 sided hyper-squares. For example, for every square
digital image with sides of 2' (where p is an integer) pixels a
Hilbert curve 25] can be constructed which traverses all pixels
continuously through adjacent pixels. In fact, 8 unique Hilbert
curves can be constructed for square 2-D images.
[0289] For UDC training (and compression), a 64.times.64 Haar
approximation of an image is unwrapped into its 8 unique Hilbert
curves on which UDCs can be trained. For query image compression, a
UDC compresses each of the 8 Hilbert sequences separately, and the
average compression rates for each pixel are recorded. Thus, if
UDC.sub.g predicts the probability p.sub.ijk.sup.g for pixel (1,1)
from the context of kth Hilbert curve (k=1, 2 . . . 8), then the
observed compressed size in bits of pixel (i, j) using UDC.sub.g
is
Y ij UDC g = k = 1 8 - log 2 ( p ijk g ) 8 ##EQU00059##
[0290] The distributions of average pixel compression rates for a
query image using each UDC are used to determine the Gleason grade
of the query image. In the following sections, reference to a query
image or its pixels specifically refers to the segmented,
down-scaled versions and not the original query image.
[0291] The novel grading system supplies both hard grades (grades
where the decision must be either grades 3, 4, or 5 and no in
between) and soft grades (where the grade is a fuzzy membership
average of multiple classes, generating a continuous scale Gleason
grade) to a query image as a whole based on UDC compression rates.
Query image pixels are also assigned individual hard and soft
grades depending on the "distance" of a trained UDC's compression
model to an optimal compression model after observation of that
pixel. (Note that only images of grades 3, 4, and 5 are considered
as these represent the significant patterns to cancer grading and
were available for our analysis.)
[0292] The system roughly consists of roughly 8 discrete parts,
illustrated in FIG. 9:
[0293] 1. UDC Training: 3 UDCs are trained on respective classes of
images. 1 UDC is trained on a query image under inspection.
[0294] a. UDC.sub.3 is trained on pathologist-graded pure pattern 3
images.
[0295] b. UDC.sub.4 is trained on pathologist-graded pure pattern 4
images.
[0296] c. UDC.sub.5 is trained on pathologist-graded pure pattern 5
images.
[0297] d. UDC.sub.0 is trained on a single query image.
[0298] 2. Query Image Compression: Each UDC compresses an ungraded
query image, I, and each pixel (i,j)'s compressed size in bits is
stored as an observation Y.sub.ij.sup.UDC.sup.g. Concurrently,
pixel errors .epsilon..sub.ij.sup.UDC.sup.g are also computed using
the Kullback-Leibler divergence of a pattern UDC with respect to
UDC.sub.0, which is the optimal compressor for the query image by
construction.
[0299] a. UDC.sub.3 compresses I with compressed pixel sizes
Y.sub.ij.sup.UDC.sup.3.
[0300] b. UDC.sub.4 compresses I with compressed pixel sizes
Y.sub.ij.sup.UDC.sup.4.
[0301] c. UDC.sub.5 compresses I with compressed pixel sizes
Y.sub.ij.sup.UDC.sup.5.
[0302] d. UDC.sub.0 compresses I with compressed pixel sizes
Y.sub.ij.sup.UDC.sup.0.
[0303] e. KL divergence measures the model errors
.epsilon..sub.ij.sup.UDC.sup.3 of UDC.sub.3 to UDC.sub.0 produced
by observation of pixel I.sub.ij.
[0304] f. KL divergence measures the model errors
.epsilon..sub.ij.sup.UDC.sup.4 of UDC.sub.4 to UDC.sub.0 produced
by observation of pixel I.sub.ij.
[0305] g. KL divergence measures the model errors
.epsilon..sub.ij.sup.UDC.sup.5 of UDC.sub.5 to UDC.sub.0 produced
by observation of pixel I.sub.ij.
[0306] 3. Analysis of Variance Testing: ANOVA is used to determine
whether or not the mean compression rates provided by UDC.sub.3,
UDC.sub.4, and UDC.sub.5 are statistically dissimilar.
[0307] 4. Hard Grading: A hard decision H for the query image
Gleason grade is made corresponding to the best-compressing UDC's
class.
[0308] 5. Hard Mapping: The query image is mapped, pixel by pixel,
with a hard grades H.sub.ij corresponding to the minimally
divergent UDC model at each pixel.
[0309] 6. Soft Grading:
[0310] a. If ANOVA testing indicates the UDC compression rates are
not statistically dissimilar, the query image is given a soft grade
S=UNLCASSIFIABLE, meaning the image is not a distinct member of any
UDC class.
[0311] b. If ANOVA testing indicates that at least 2 UDCs are
statistically dissimilar
[0312] i. Tukey's multiple comparison test is used to test the
independence of each UDC's mean compared to the others' means.
[0313] ii. The query image is given a soft grade S using weighted
means based on the intersections of the t-distributions modeling
the mean UDC compression rates.
[0314] 7. Soft Mapping: The query image is mapped, pixel by pixel,
with a soft grades S.sub.ij computed using weighted means based on
the inverse KL divergence of each pixel.
[0315] 8. Prediction of Most Likely Pattern Composition:
[0316] a. The probability of grade transitions from each grade to
the others are computed from the errors associated with each UDC
model.
[0317] b. The final pattern composition (in percent of each
pattern) is computed through steady state analysis of the
transition model.
[0318] Prognosis System for Integration of Biopsy Classifier
Systems with Patient and PSA Information (System 8)
[0319] In one embodiment of the present invention, a computerized
system is developed by introducing a unique recognition and
prognostic system that provides the following: (1) Recognition and
classification of Gleason grading in 2-D biopsy images, 3-D biopsy
images, and 3-D core patterns; (2) faster and accurate
classification with the introduction of image preprocessing
algorithms; and (3) integration of patient personal data and
prostate-specific antigen (PSA) information for more accurate
recognition and prognostics analysis.
[0320] Described herein is a system for the recognition and
classification of Gleason grades of pathological biopsy images
using three sub-systems for histology grading that are integrated
at a later stage with patient information such as PSA information,
age, race, and family history. The general hierarchy of such
integration is illustrated in FIG. 17.
[0321] The first histology subsystem considers the classification
of 2-D biopsy images; one or more of the previous systems may be
used. The second subsystem, considers the classification of 3-D
biopsy images. In this part, a 2-D image is mapped using certain
mapping algorithms to 3-D space to have a 3rd dimension. The
conventional way of visualizing regular 2-D biopsy images is
replaced by a 3-D view to observe histological details in 3-D
space. The objective of this subsystem is to verify the
classification of the 2-D subsystem in terms of analyzing the
unforeseen information that may exist in the 3-D representation.
The classification performance of the 2-D and 3-D subsystem is
going to be fused for accurate recognition performance. The third
subsystem introduces the concept of 3-D core grading of prostate
biopsy specimen. In this subsystem serial images of a stack of
multiple 2-D tissue biopsy slices are captured from a given
suspicious lesion in the prostate gland. The 3rd dimension of this
stack represents the depth of the tumor or "core". It can be
visualized as a 2-D internal image of the prostate gland when a
vertical cut is made on the stack of the 3-D slices. The objective
of this subsystem is to investigate tertiary patterns, and assigns
a malignancy Gleason grade for the core (depth of the specimen).
The final Gleason score resulted from these three subsystems can be
looked at as a Gleason grade composed from 3 scores, the first two
scores represent the Gleason score from the first two subsystems,
and the third score is the core score as follows:
Gleason score=(predominant pattern grade+second predominant pattern
grade+3-D core grade)
[0322] Regarding the PSA and patient information, information from
these blocks will be integrated as an additional feature for
accurate classification and prognostics purposes.
[0323] Biopsy images of prostate cancer form diversities in terms
of structure, color information, edge information, and
visualization of fine details associated with each distinct grade.
A typical image recognition system for classification purpose is
composed of the following three main blocks/phases: (1) Images
preprocessing; (2) features extraction (feature vectors); and (3)
classification. Considering the first block that is concerned with
image preprocessing, image details are collected from images at
different Gleason grades prior to feature vectors extraction phase.
Preprocessing of images provides significant information in terms
of visualization, segmentation, edge information prior to
categorizing images in different cancer grades through feature
vectors generation and classification.
[0324] In one embodiment of the present invention, methods for 2-D
and 3-D grading systems integration are provided. The steps to
achieve a prostate cancer prognosis system as a follows:
[0325] (1) Acquire PSA information attached with each graded biopsy
image as follows: PSA level, PSA-Velocity, % free-to-total PSA.
[0326] (2) Acquire patient information attached with each labeled
biopsy image as follows: age, race, Body-mass Index (BMI), and
family history of cancer.
[0327] (3) Feature Vector Generation and Classifications of the
above information for features integration with biopsy images
features associated with 2-D, 3-D and Core grading classification.
The systems provides in previous embodiments may be used.
[0328] (4) Verification of the whole system using subjective and/or
objective measures. System output may be compared with biopsies
graded by expert pathologists. Mean Opinion Score (MOS) may be used
to measure the performance of the system. Also, DNA Microarrays
results, if available for each patient, may be analyzed and
compared with our system final output as an objective measure.
[0329] The introduction of a new novel Gleason grading concept
considering different aspects of biopsy images ranging from 2-D to
3-D contributes to an enhanced output of more than just an accurate
classification of cancer grade. The system of this embodiment
integrates PSA and patient information along with results of cancer
automatic grading from the proposed histology system. Hence, the
final "expert system" output is expected to provide the following
functions: [0330] Provides an automated Gleason score to assets
pathologist in accurate classification of cancer grades of prostate
biopsy images. [0331] Aids medical doctors by providing prognostic
information about the stage and level of prostate cancer. [0332]
The generated accurate Gleason score contributes directly to more
accurate risk assessment readings using the common available risk
assessment calculators. Accuracy Contributes to more accurate risk
assessment by integrating a more accurate Gleason grade with other
risk assessment variables. [0333] Biopsy planning by predicting
biopsy locations having suspicious malignant tissues.
[0334] Implementation of the method/or system of embodiments of the
invention can involve one/or combination of pre-developed or new
risk assessment models and 2-D and 3-D Gleason grading/scoring
systems. For example, we may use one of following risk assessment
model: Prostate Cancer Prevention Trial Risk Calculator (PCPTRC)
21], D'Amico model 22], Kaftan nomograms 23], and Cancer of the
Prostate Risk Assessment (CAPRA) 24]). In general, a risk
assessment system uses PSA level (blood test), Gleason grade, and T
stage (size of the tumor on rectal exam and/or ultrasound) to group
men as low, intermediate, or high-risk. In all these methods, the
accurate detection, grading and scoring of biopsy images is a
crucial step towards a better risk assessment.
[0335] System and Methods for Estimation, Prediction, and Cancer
Grading Outside of the Histopathology Biopsy Images
[0336] In another aspect of this embodiment, a method is provided
for predicting the Gleason grade of blocks in histopathology biopsy
images. The basic idea is to predict the state of missing blocks or
to produce extrapolated values to know the cancer grade in a
neighborhood outside the core extracted from the prostate gland
during a needle biopsy procedure. To this end, two models are
fitted. The first model is fitted by taking into account the
distribution of Gleason grades within the processed 2-D or 3-D
core. The second model is a pixel-based model, where the pixel
value and other features discussed in previous embodiments can be
extracted to construct a prediction model. Several algorithms, for
example the Kalman filtering schemes 26] or similar techniques can
be used to make the aforementioned predictions. Once the two
predictions can be performed, they can be fused to get a unique
cancer grade for the new predicted blocks. An exemplary procedure
is shown in FIG. 28.
[0337] Referring in more details to the present embodiment, the
model of the histopathology image is given by:
S(m,n)=AS(m-1,n)+Bw(m,n)
[0338] Where S is the state vector, A is a state transition matrix,
B is an input matrix, w is random system noise. S represents the
noisy observed image features, which is directly related to the
morphology and textural features extracted for classification
purposes. The vector w may be modeled as a Gaussian noise vector
with zero mean, which represent variations in calculated features.
As it was mentioned above, pixel- and block-based model fitting may
be performed.
[0339] The observation equation is given by:
r(m,n)=h.sup.TS(m,n)+.nu.(m,n)
[0340] Where h is an observation matrix and v is additive noise
present in observations.
[0341] The next step is to generate a Kalman filter, for prediction
of the cancerous grade of image blocks by estimating the state
vector. In order to execute the 2-D filtering, we divide the 2-D
problem into equivalent 1-D problems by applying line by line
scanning or keeping constant the variable n or second dimension of
the image or block, and then constructing a global state vector for
all values of n.
[0342] For example, the equivalent 1-D Kalman filtering equations
for a 1-D image scanning process are:
(m)=A(m-1)
[0343] P.sub.b (m)=AP.sub.a(m-1)A.sup.T+BQ.sub.wB.sup.T, Q.sub.w is
the correlation matrix of noise vector w.
[0344] K(m)=P.sub.b(m)h.sup.T hP.sub.b(m)h.sup.T+Q.sub..nu.].sup.-1
Q.sub..nu. is the correlation matrix of noise vector v.
(m)=(m)+K(m)r(m)-h.sub.b(m)]
P.sub.a(m)=I-K(m)h]P.sub.b(m)
[0345] The resulting estimated state vector is used as the input to
a previously trained classifier to obtain the Gleason grade of the
2-D image patch. Extensions and generalizations of the theory of
estimation by Kalman filtering may be used to generate estimated
cancer grade of 2-D slices and then reconstruct 3-D images and
assign a Gleason grade using the classifier for 3-D images
developed in previous embodiments.
Cancer Tele-Screening System and Adjudication Schemes (System
9)
[0346] A tele-screening system according to an embodiment of the
invention is provided for classification and diagnosis from
digitized biopsy images in an ideal adjudication scheme. In this
embodiment, the final report will be produce by integrating
pathologists and automatic expert systems to produce a more
accurate cancer diagnosis based on biopsy assessment. FIG. 18 shows
the proposed tele-screening system. Several variations of the ideal
adjudication scheme are addressed in different aspects of the
present invention.
[0347] In one embodiment, an adjudication paradigm for cancer
diagnosis biopsy by integrating one expert system is proposed. For
instance, the configuration of the system is presented in FIG. 19
and includes the following steps: (1) One expert pathologist and
one of the computerized systems (CAD or expert system) proposed in
various embodiments of this invention perform independent
assessment of all the tissue slides extracted by a biopsy procedure
and produce a diagnostic report and annotated images. (2) An
automatic algorithm analyzes the level of agreement between the two
reports including patch analysis of slides compared to pathologist
annotation. The assessment is considered complete if an acceptable
level of agreement previously defined is reached. (3) If the
minimum acceptable level of agreement is not reached, biopsy
samples are assessed by a second pathologist. (4) The automatic
algorithm of step (2) analyzes the level of agreement among the
three reports in a patch-based basis and force decision, for
example by majority voting rule. (5) The misclassified patches are
added to the training database and the classification model is
updated. (6) The final step uses cross-validation to ensure that
the accuracy and other measurement such as sensitivity and
specificity remain in an acceptable interval. For cross-validation,
5-fold, 10-fold, hold-out or leave-one-out methods may be used.
This system is able to decrease the interpretation time of a biopsy
and reduce the cost associated to a second pathologist.
[0348] In another embodiment, an adjudication paradigm for cancer
diagnosis biopsy by integrating two independent expert systems is
proposed. The flow diagram of the method is shown in FIG. 20. The
configuration of the system includes the following steps: (1) One
expert pathologist and one of the computerized systems proposed in
various embodiments of this invention perform independent
assessment of all the tissue slides extracted by a biopsy procedure
and produce a diagnostic report and annotated images. (2) An
automatic algorithm analyzes the level of agreement between the two
reports including patch analysis of slides compared to pathologist
annotation. The assessment is considered complete if an acceptable
level of agreement previously defined is reached. (3) If the
minimum acceptable level of agreement is not reached, biopsy
samples are assessed by a second expert system, which use different
features from the first CAD. For instance, the first CAD may be
based on morphological features and the second one on textural
features. (4) The automatic algorithm of step (2) analyzes the
level of agreement among the three reports in a patch-based basis
and force decision, for example by majority voting rule. (5) The
misclassified patches are added to the training database and the
classification model is updated. (6) The final step uses a
cross-validation to ensure that the accuracy and other measurement
such as sensitivity and specificity remain in an acceptable
interval. For cross-validation, 5-fold, 10-fold, hold-out or
leave-one-out methods may be used. This system is able to decrease
the interpretation time of a biopsy and further reduce the cost
associated to second and third pathologists.
[0349] In another embodiment, an adjudication paradigm for cancer
diagnosis biopsy by integrating three independent expert systems is
provided. The flow diagram of the method is shown in FIG. 21. The
configuration of the system includes the following steps: (1) Two
CAD systems proposed in various embodiments of this invention
perform independent assessment of all the tissue slides extracted
by a biopsy procedure and produce a diagnostic report and annotated
images. (2) An automatic algorithm analyzes the level of agreement
between the two reports including patch analysis of slides compared
to pathologist annotation. The assessment is considered complete if
an acceptable level of agreement previously defined is reached. (3)
If the minimum acceptable level of agreement is not reached, biopsy
samples are assessed by a second expert system, which use different
features from the first CADs. (4) The automatic algorithm of step
(2) analyzes the level of agreement among the three reports in a
patch-based basis and force decision, for example by majority
voting rule. This system is able to decrease the interpretation
time of a biopsy from minutes to seconds and further reduce the
cost associated to second and third pathologists.
[0350] In broad embodiment, the present invention provides several
systems for 2-D and 3-D classification of histology tissue samples.
Such systems use unique feature vectors representing discriminative
attributes of each cancerous pattern, for example Gleason patterns.
In addition, other systems are proposed for cancer prognosis and
cancer tele-screening by integrating one or more of the recognition
systems. Moreover, novel algorithms for image processing such as
edge detection, color fractal dimension and color region mapping
are provided.
ADDITIONAL EMBODIMENTS
[0351] The present invention which is described hereinbefore with
reference to flowchart and/or block diagram illustrations of
methods, systems, devices, simulations, and computer program
products in accordance with some embodiments of the invention and
has been illustrated in detail by using a computer system. For
instance, the flowchart and/or block diagrams further illustrate
exemplary operations of the computer systems and methods of FIGS. 1
to 28. It is also conceivable that each block of the flowchart
and/or block diagram illustrations, and combinations of blocks in
the flowchart and/or block diagram illustrations, may be
implemented by any computer program instructions and/or hardware.
These computer program instructions may be provided to a processor
of a general purpose computer, a microprocessor, a portable device
such as cell phones, a special purpose computer or device, or other
programmable data processing apparatus to produce a device, such
that the instructions, which execute via the processor of the
computer or other programmable data processing apparatus, create
means for implementing the functions specified in the flowcharts
and/or block diagrams or blocks. The computer program may also be
supplied from a remote source embodied in a carrier medium such as
an electronic signal, including a radio frequency carrier wave or
an optical carrier wave.
[0352] Further modifications and alternative embodiments of various
aspects of the invention will be apparent to those skilled in the
art in view of this description. Accordingly, this description is
to be construed as illustrative only and is for the purpose of
teaching those skilled in the art the general manner of carrying
out the invention. It is to be understood that the forms of the
invention shown and described herein are to be taken as examples of
embodiments. Elements and materials may be substituted for those
illustrated and described herein, parts and processes may be
reversed, and certain features of the invention may be utilized
independently, all as would be apparent to one skilled in the art
after having the benefit of this description of the invention.
Changes may be made in the elements described herein without
departing from the spirit and scope of the invention as described
in the following claims.
* * * * *