Login| Sign Up| Help| Contact|

Patent Searching and Data


Title:
SYSTEM AND METHOD FOR ACCURATE AND RAPID IDENTIFICATION OF DISEASED REGIONS ON BIOLOGICAL IMAGES WITH APPLICATIONS TO DISEASE DIAGNOSIS AND PROGNOSIS
Document Type and Number:
WIPO Patent Application WO/2010/027476
Kind Code:
A1
Abstract:
The present invention relates to a method and system for detecting biologically relevant structures in a hierarchical fashion, beginning at a low-resolution and proceeding to higher levels of resolution. The present invention also provides probabilistic pairwise Markov models (PPMMs) to classify these relevant structures. The invention is directed to a novel classification approach which weighs the importance of these structures. The present invention also provides a fast, efficient computer-aided detection/diagnosis (CAD) system capable of rapidly processing medical images (i.e. high throughput). The computer-aided detection/diagnosis (CAD) system of the present invention allows for rapid analysis of medical images the improving the ability to effectively detect, diagnose, and treat certain diseases.

Inventors:
MADABHUSHI ANANT (US)
MONACO JAMES (US)
TOMASZEWSKI JOHN (US)
FELDMAN MICHAEL (US)
BASAVANHALLY AJAY (US)
Application Number:
PCT/US2009/004979
Publication Date:
March 11, 2010
Filing Date:
September 03, 2009
Export Citation:
Click for automatic bibliography generation   Help
Assignee:
UNIV RUTGERS (US)
MADABHUSHI ANANT (US)
MONACO JAMES (US)
TOMASZEWSKI JOHN (US)
FELDMAN MICHAEL (US)
BASAVANHALLY AJAY (US)
International Classes:
G06V10/25
Foreign References:
US20060064248A12006-03-23
US20060291721A12006-12-28
US20030219760A12003-11-27
US20070031005A12007-02-08
US20050286750A12005-12-29
US20080032293A12008-02-07
Attorney, Agent or Firm:
DAVIDSON, DAVIDSON & KAPPEL, LLC (14th FloorNew York, NY, US)
Download PDF:
Claims:
What is claimed is:

Claim 1 : A method for detecting biologically relevant structures in a hierarchical fashion comprising:

a) performing gland segmentation on the luminance channel to produce gland boundaries;

b) calculating morphological features for each gland;

c) classifying and labeling the morphological features for each gland as either malignant or benign;

d) performing MRP iteration using the labeling in step (c) as a the starting point to produce a final labeling; and

e) consolidating glands labeled as malignant into regions.

Claim 2: The method of claim 1 wherein the biologically relevant structure is selected from the group consisting of the prostate, breast, ovary, bladder, kidney, or brain.

Claim 3: The method of claim 1 wherein the classifying and labeling morphological features in step (c) is performed using probabilistic pairwise Markov models.

Claim 4: The method of claim 1 wherein the performing a MRP iteration in step (d) is performed using weighted MAP (WMAP) and weighted MPM (WMPM) estimation to allow certain decisions to be weighted more heavily than others.

Claim 5: The method of claim 1 wherein the performing a MRP iteration in step (d) is performed using weighted ICM estimation. Claim 6: The method of claim 1 wherein the step of performing a MRF iteration in step (d) is performed using weighted MPM Monte Carlo estimation.

Claim 7: The method of claim 1 wherein the step of performing gland segmentation on the luminance channel to produce gland boundaries in step (a) is performed using CIE Lab color space.

Claim 8: The method of claim 1 wherein the step of performing gland segmentation on the luminance channel to produce gland boundaries in step (a) comprises:

a) convolving a luminance image with a Gaussian kernel at multiple scales σg e {0.2, 0.1, 0.05, 0.025}mm to account for varying gland size;

b) identifying seed pixel peaks (maxima) in the smoothed images as candidate gland centers

c) performing a growing procedure on the pixel peaks of step (b) comprising the following steps:

(i) initializing a current region (CR) to the specified seed pixel and establishing a

12σs x 12σ9 bounding box centered about it and initializing a current boundary (CB) to the

8-connected pixels neighboring CR;

(ii) identity the pixel in CB with the highest intensity and removing the highest intensity pixel from CB and adding it to CR and revising CB to include all 8 connected neighbors of the aggregated pixel which are not in CR;

(iii) defining the internal boundary (IB) as all pixels in CR that are 8-connected with the pixels in CB and computing the current boundary strength which is defined as the mean intensity of the pixels in IB minus the mean intensity of the pixels in CB;

(iv) repeating steps (ii) and (iii) until the algorithm attempts to add a pixel outside the bounding box; (v) identifying a iteration step at which the maximum boundary strength was attained and defining the optimal region as CR at this step.

(vi) resolving any overlapping of the final segmented region by discarding the region with the lower boundary strength.

Claim 9: An automated system for detecting CaP glands in WMHSs using the method of claim 1.

Claim 10: A method for detecting biologically relevant structures in a hierarchical fashion comprising:

a) identifying and segmenting gland lumens to produce segmented boundaries;

b) extracting morphological features from the segmented boundaries;

c) providing a Bayesian classifier to model the glandular area;

d) assigning a value of probability of malignancy to each gland;

e) classifying and labeling the gland as malignant if the probability exceeds the predetermined threshold ταreα or labeling the gland as benign if the probability does not exceed the predetermined threshold rαreα .

f) using a MRF algorithm to refine the labeling of step(c).

Claim 1 1 : The method of claim 10 wherein the classifying and labeling a gland in step (e) is performed using probabilistic pairwise Markov models.

Claim 12: A computer-based system for identifying and assessing the biological significance of abnormal regions in medical images comprising: a) segmenting and consolidating biologically relevant structures using the individual features of these structures, the features and classes of those structures that surround them, and the architectural arrangement of the structures; and b) analyzing the biologically relevant structures in a hierarchical manner to assess the biological significance of any abnormal regions c) assigning the biological significance of step (b) to one of two classes: those regions in which the available regions in which the available information is sufficient to make an assessment of the biological significance or those regions that should be examined at the next hierarchical level.

Claim 13: The system of claim 12 wherein the biologically relevant structure is selected from the group consisting of the prostate, breast or ovary.

Claim 14: The method of Claim 12, whereby the said medical image is a digitized histological section or tissue micro-array of the prostate, breast, ovary, bladder, kidney, or brain.

Claim 15: The method of Claim 12, whereby the said medical image is a radiological image including, but not limited to, an x-ray, magnetic resonance, positron emission tomography, or CAT scan of a prostate, breast, ovary, bladder, kidney, or brain.

Claim 16: The method of Claim 12, whereby the said biologically relevant structures are breast calcifications, prostate glands (acini), or breast lesions.

Claim 17: The method of Claim 12, whereby the said biologically relevant structures are cells or lymphocytes of the prostate, breast, or ovary. Claim 18: The method of Claim 12, whereby the said hierarchy consists of images at different resolutions.

Claim 19: The method of Claim 12, whereby the said individual features describe the morphological, textural, and/or intensity (derived from either color or grey-scale images) of the structures.

Claim 20: The method of Claim 12, whereby the said individual feature is area, width, length, curvature, perimeter, aspect ratio, and compactness.

Claim 21 : The method of Claim 12, whereby the said individual feature is first-order, Haralick, wavelet, or Gabor.

Claim 22: The method of Claim 12, whereby said individual feature is an average, a median, a standard deviation, a difference, a Sobel filter, a Kirsch filter, a horizontal derivative, a vertical derivative, or a diagonal derivative of a pre-selected determinant.

Claim 23: The method of Claim 12, wherein said individual feature is calculated from the boundary of either the structure as the ratio of the structure's area to the area of the smallest circle that circumscribes the structure; the standard deviation, variance, or ratio of the maximum to the average distance between the structure center of mass and the boundary; the ratio of the estimated boundary length (the boundary calculated using a fraction of the boundary pixels) to the actual boundary length; ratio of the boundary length to the area enclosed by the boundary; and the sum of the difference between the distance from the center of the structure to the boundary and the average of the distances from the center to the two adjacent points. Claim 24: The method of Claim 12, whereby the said abnormal regions denote cancer, pre- malignant lesions, prostate intraepithelial neoplasia, or the presence of breast microcalcifications.

Claim 25: The method of Claim 12, whereby the said assessment is Gleason grading or Bloom- Richardson grading.

Claim 26: The method of Claim 12, whereby the said assessment is an indicator of cancer recurrence, disease aggressiveness, diagnosis, prognosis, or theragnosis.

Claim 27: The method of Claim 12, whereby the said architectural feature is the density of the structures or the average direction of the major axis of the structures.

Claim 28: The method of Claim 12, wherein said architectural feature is computed from a Delaunay graph is the standard deviation, average, ratio of minimum to maximum, or disorder of the side lengths or areas of the triangles in the Delaunay graph of said image.

Claim 29: The method of Claim 12, wherein said architectural feature is computed from a minimum spanning tree is the standard deviation, average, ratio of minimum to maximum, or disorder of the edge lengths in the minimum spanning tree of said image.

Claim 30: A computer-aided diagnosis system to automatically detect and grade the extent of LI in digitized HER2+ BC histopathology images comprising:

a) identifying candidate image locations that may represent centers of lymphocytic nuclei using a region-growing algorithm to define sharp, well-defined boundaries; b) refining the initial lymphocyte identification of step (a) using a Maximum a Posteriori estimation to incorporate domain knowledge regarding lymphocyte size, luminance, and spatial proximity; . c) further refining the detection result of step (b) using a Markov Random Field model to ensure that an object is more likely to be labeled as a lymphocyte nucleus if it surrounded by other lymphocyte nuclei

d) reiterating the Markov Random Field model of step (c) until convergence by using an Iterated Conditional Modes algorithm.

Claim 31 : The system of claim 30 wherein the Maximum a Posteriori estimation of step (b) is a Bayesian model used in conjunction with object size and luminance to classify individual objects as either lymphocyte or BC nuclei.

Claim 32: A method for detecting biologically relevant structures in a hierarchical fashion comprising:

a) identifying and segmenting biologically relevant structures;

b) calculating morphological or textural features for each biologically relevant structure;

c) classifying each biologically relevant structure using a MRF;

d) consolidating each biologically relevant structure into regions;

e) drawing biological conclusions from consolidated regions.

Claim 33: The method of claim 32 wherein the biologically relevant structure is selected from the group consisting of the prostate, breast, ovary, bladder, kidney, or brain.

Claim 34: The method of claim 32 wherein the classifying each biologically relevant structure using a MRF in step (c) is performed using probabilistic pairwise Markov models. Claim 35: The method of claim 32 wherein the step of identifying and segmenting biologically relevant structures in step (a) is performed using ClE Lab color space.

Claim 36: The method of claim 32 wherein the step of identifying and segmenting biologically relevant structures in step (a) comprises:

a) convolving a luminance image with a Gaussian kernel at multiple scales σg € {0.2, 0.1, 0.05, 0.025}mm to account for varying gland size;

b) identifying seed pixel peaks (maxima) in the smoothed images as candidate gland centers

c) performing a growing procedure on the pixel peaks of step (b) comprising the following steps:

(i) initializing a current region (CR) to the specified seed pixel and establishing a

12σ5 x 12σff bounding box centered about it and initializing a current boundary (CB) to the

8-connected pixels neighboring CR;

(ii) identity the pixel in CB with the highest intensity and removing the highest intensity pixel from CB and adding it to CR and revising CB to include all 8 connected neighbors of the aggregated pixel which are not in CR;

(iii) defining the internal boundary (IB) as all pixels in CR that are 8-connected with the pixels in CB and computing the current boundary strength which is defined as the mean intensity of the pixels in IB minus the mean intensity of the pixels in CB;

(iv) repeating steps (ii) and (iii) until the algorithm attempts to add a pixel outside the bounding box;

(v) identifying a iteration step at which the maximum boundary strenglb was attained and defining the optimal region as CR at this step.

(vi) resolving any overlapping of the final segmented region by discarding the region with the lower boundary strength, Claim 37: An automated system for detecting CaP glands in WMHSs using, the method of claim 32.

Description:
SYSTEM AND METHOD FOR ACCURATE AND RAPID IDENTIFICATION OF DISEASED REGIONS ON BIOLOGICAL IMAGES WITH APPLICATIONS TO DISEASE

DIAGNOSIS AND PROGNOSIS

Cross Reference to Related Applications

[001] This application claims priority to U.S. Provisional Patent Application No. 61/093,884 filed September 3, 2008, the disclosure of which is hereby incorporated by reference in its entirety.

Field of the Invention

[002] The present invention relates to a method and system for detecting, analyzing and consolidating biologically relevant structures in a hierarchical fashion, beginning at a low- resolution and proceeding to higher levels of resolution. The present invention also provides probabilistic pairwise Markov models (PPMMs) to classify these relevant structures. The invention is directed to a novel classification approach which weighs the importance of these structures. The present invention also provides a fast, efficient computer-aided detection/diagnosis (CAD) system capable of rapidly processing medical images (i.e. high throughput). The computer-aided detection/diagnosis (CAD) system of the present invention allows for rapid analysis of medical images the improving the ability to effectively detect, diagnose, and treat certain diseases.

Background of the Invention

[003] The examination of histological specimens obtained from core needle biopsies remains the definitive test for diagnosing prostate cancer (CaP). If cancer is present in these biopsy specimens, a surgeon may perform a radical prostatectomy (excise the prostate). Following prostatectomy, the prostate is sliced into whole-mount histological sections (WMHS). Staging and grading these WMHSs can help determine patient prognosis and treatment. Additionally, the spatial extent of CaP, as established by the analysis of WMHSs, can be registered to other modalities (eg. magnetic resonance imaging), providing a "ground truth" for the evaluation of computer-aided diagnosis (CAD) systems that operate on these modalities. The development of CAD algorithms tor WMHSs is also sigmticant tor the following reasons: i; CAD offers a viable means for analyzing the vast amount of the data present in WMHSs, a time- consuming task currently performed by pathologists, 2) the extraction of reproducible, quantifiable features can help refine our own understanding of prostate histopathology, thereby helping doctors improve performance and reduce variability in grading and detection, and 3) the data mining of quantified morphometric features may provide means for biomarker discovery, enabling, for example, the identification of patients whose CaP has a high likelihood of metastasis or post-treatment recurrence.

[004] With respect to prostate histology, Begelman [49] presented a method for nuclei segmentation on hematoxylin and eosin (RScE) stained prostate tissue samples using a semi- supervised, fuzzy-logic engine. Adiga [50] presented a sophisticated system for the fast segmentation of cell nuclei in multispectral histological images. Doyle et al. used image texture features [51] and features derived from segmented nuclei and glands [52] to determine Gleason grade in core biopsy samples. Exploiting both domain-specific knowledge and low-level textural features, Naik et al. (References [53] and [54] )developed an automated system for segmenting glandular structures and discriminating between intermediate Gleason grades in core samples. To aid in manual cancer diagnosis, Gao [55] applied histogram thresholding to enhance the appearance of cytoplasm and nuclei. Hafiane [56] performed gland structure segmentation using a spatially constrained adaptation of fuzzy c-means and a multiphase level-set algorithm.

[005] The most significant impediment in the development of automated systems for the analysis of MHSs is data volume. A typical whole-mount histological section digitized at 0.25 μm per pixel (equivalent to 4Ox magnification) contains approximately 200, 000 x 60, 000 pixels. This is 800 times the number pixels in a digitized mammogram. Consequently, nearly all previous systems (see References [84] to [95]) restrict their analysis to specific regions of interest (ROIs). These ROIs are typically less than one-thousandth the size of a WMHS. The only previous attempt — aside from our own — to identify CaP in relatively large HSs was presented by Diamond. Diamond divided a single WMHS into small 100 x 100 patches, classifying each patch individually using a single Haralick feature. However, the algorithm requires manual segmentation and classification of the glands. Additionally, execution time was 5.5 hours. [006] Contextual information can be invaluable in estimation tasks. For example, when determining a missing word in a document, clues can be ascertained from adjacent words, words in the same sentence, and even words in other sentences. In image processing tasks such as denoising, texture synthesis, and smoothing, estimating the intensity at any single pixel is facilitated by knowledge of the remaining pixel intensities. Unfortunately, modeling contextual information may be exceedingly difficult, especially if extensive dependencies exist among all sites (a generic term representing entities such as pixels or words). However, in many systems the preponderance of contextual information regarding an individual site is contained in those sites that neighbor it. These neighborhoods are often functions of spatial or temporal proximity. Since the contextual information in such systems can be restricted to local neighborhoods, modeling the systems becomes tractable.

[007] In a Bayesian framework, the restriction of contextual information to local neighborhoods is called the Markov property, and a system of sites that obeys this property is termed a Markov random field (MRF). The merits of MRFs have been demonstrated for a variety of computer vision tasks such as clustering, and texture synthesis. MRFs have demonstrated particular proficiency in image segmentation and object detection tasks. Zhang [33], for example, used MRFs and expectation-maximization to segment magnetic resonance (MR) brain images. Farag [34] used a similar approach to segment multimodal brain and lung images. Li [35]applied MRFs to tumor detection in digital mammography. Employing an adaptive Markov model, Awate [36] presented an unsupervised method for classifying MR brain images. In general, MRF segmentation techniques have evolved into sophisticated algorithms that employ multiresolution analysis and complex boundary models.

[008] MRFs are established through the construction of local conditional probability density functions (LCPDFs). These LCPDFs — one centered about each site — model the local inter- site dependencies of a random process. In combination, these LCPDFs can establish a joint probability density function (JPDF) relating all sites. However, only LCPDFs of certain functional forms will reconstitute a valid JPDF. Specifically, the Gibbs-Markov equivalence theorem states that the JPDF will be valid if (and only if) it, and transitively the LCPDFs, can be represented as Gibbs distributions. Untortunately, constructing L( ^ f υts mat simultaneously meet the following three conditions:

1) satisfy the Markov property;

2) combine to yield a valid JPDF; and

3) sufficiently model an underlying process

is nontrivial. Consequently, current MRF models are generic and/or heuristic, and thus, ignore crucial information. Currently, the computer vision literature advocates two disparate methods for constructing LCPDFs that meet these three conditions: parametric and nonparametric modeling.

[009] Nonparametric modeling primarily focuses on condition (3). It is assumed, though the justification seems mostly empirical, that conditions (1) and (2) will be realized during the relaxation procedure (e.g. stochastic relaxation or iterated conditional modes). For example, in the case of texture reconstruction (i.e. constructing textural images from prescribed textural patterns), the LCPDFs — whose functional forms are initially unconstrained — are estimated from training images; the reconstruction algorithm attempts to rectify these unconstrained LCPDFs to satisfy conditions 1 and 2 during the actual synthesis (reconstruction). That is, the synthesizing procedure implicitly modifies the LDCPFs during texture generation, producing (hopefully) a random texture whose LCPDFs form a valid JPDF and whose appearance resembles the original texture. Regardless of the success of this process, the rectified LCPDFs are themselves not directly accessible (though they could perhaps be estimated by sampling the synthesized texture). Additionally, since we have no insight as to the possible functional forms of the LCPDFs, we can not anticipate whether or not the rectified LCDPFs will be able to describe our specific system. This "black box" type approach removes the insight typically provided by Bayesian modeling.

[010] Parametric methods directly address conditions (1) and (2) by exploiting the Gibbs- Markov equivalence. That is, representing the JPDF as a Gibbs distribution guarantees that the attendant LCPDFs satisfy the Markov property and form a viable JPDF. However, satisfying condition 3 within the Gibbsian framework is difficult. Since Gibbs distributions are expressed as a product ot potential functions, tailoring L(SfUtS to model a speciiic process αevoives into the intelligent selection of these functions. Unfortunately, potential functions are mathematical abstractions, lacking intuition. Consequently, constructing LCPDFs through their selection becomes an ad hoc procedure, usually resulting in generic and/or heuristic models.

[Oil] Random fields, which are simply collections of random variables, provide a means for performing complex classification tasks within a Bayesian framework. Markov random fields — a particular type of random field — have proven useful in a variety of applications such as segmentation, texture synthesis. However, since random fields can assume a large number of states (classes), techniques for classifying them cannot rely on exhaustive searches and must employ more sophisticated techniques. Such techniques include simulated annealing, iterated conditional modes (ICM), loopy belief propagation, and graph cuts. However, these techniques, and nearly all other such methods, perform either maximum posterior marginals (MPM) or maximum a posteriori (MAP) estimation (or a variation of the two).

[012] The weakness of MPM and MAP estimation lies in their inability to weight certain classification decisions more heavily than others — they implicitly weight all decisions equally. This ability is crucial for two reasons. First, some tasks naturally give rise to outcomes that are asymmetrical. For example, when identifying cancerous tissue, mislabeling a malignant lesion is far more egregious than misclassifying a benign structure. Second, many classification systems require the capability of adjusting their performance. For example, the performance of commercial systems for detecting mammographic abnormalities is typically adjusted to the highest detection sensitivity that incurs no more than two false positive per image. Unfortunately, since MRFs are restricted to MPM and MAP estimation, they are ill suited for such tasks.

[013] Though others have addressed these problems [See References 1-10], they have used approaches that leverage low-level features (e.g. Haralick, wavelet, fractal). These features require extensive processing time, precluding the ability to analyze large images such as those found in histological analysis.

Brief Summary of Invention [014) The present invention relates to a method and system tor detecting Dioiogicaiiy relevant structures in a hierarchical fashion, beginning at a low-resolution and proceeding to higher levels of resolution.

[015] It is an object of the present invention to provide a system and method for detecting potentially cancerous regions in prostate histology at low-resolution.

[016] It is an object of the present invention to provide probabilistic pairwise Markov models (PPMMs) to classify biologically relevant structures.

[017] It is an object of the present invention to provide probabilistic pairwise Markov models (PPMMs). PPMMs formulate the LCPDFs in terms of pairwise probability density functions (PDFs), each of which models an interaction between two sites.

[018] It is an object of the present invention to provide an automated system for detecting CaP glands in WMHSs.

[019] In certain embodiments of the present invention, gland size alone is a sufficient feature for yielding an accurate and efficient algorithm for detecting cancerous glands.

[020] It is an object of the present invention to provide a method for incorporating asymmetric costs into MRF estimation.

[021] In accordance of this object, the present invention provides weighted MAP (WMAP) and weighted MPM (WMPM) estimation, generalizations of MAP and MPM estimation that allow certain decisions to be weighted more heavily than others.

[022] In a further embodiment of the present invention, we introduce weighted ICM (WICM), a novel adaptation of ICM capable of WMAP estimation, and weighted MPM Monte Carlo (WMPMMC), an extension of MPM Monte Carlo capable of WMPM estimation.

[023] In accordance with the above objects and others, the invention is directed in part to a system and method for detecting biologically relevant structures wherein: 1) gland segmentation is performed on the luminance channel of a color H & E stained WMHS, producing gland boundaries, 2) the system then calculates morphological features for each gland, 3) the features are classified, labeling the glands as either malignant or benign, 4) the labels serve as the starting point tor the MKt- iteration whicn men produces me tinai labeling, ana DU me cancerous gianαs are consolidated into regions.

[024] An embodiment of the present invention provides a fast, efficient computer-aided detection/diagnosis (CAD) system capable of rapidly processing medical images (i.e. high throughput). The computer-aided detection/diagnosis (CAD) system of the present invention allows for rapid analysis of medical images the improving the ability to effectively detect, diagnose, and treat certain diseases.

[025] The present invention is further directed to the construction of a computer-based system for identifying and assessing the biological significance of suspicious (abnormal) regions in medical images comprising segmenting and consolidating biologically relevant structures using the individual features of these structures, the features and classes of those structures that surround them, and the architectural arrangement of the structures.

[026] In certain embodiments, the invention is directed to a method of detection/grading of malignant glands in a prostate histological section, the detection of malignant regions in prostate magnetic resonance image (MRI), or the detection of micro-calcifications in a mammogram.

[027] In accordance with the above object and others, the invention is directed to a novel classification approach which weighs the importance of these structures.

Brief Description of Drawings

[028] Figure 1 depicts a graphical representation showing six sites and binary states indicating that each site has an associated feature ^ e R.

[029] Figure 2 depicts an overview of a classification algorithm for detecting CaP glands from digitized WMHSs.

[030] Figure 3(a) depicts H&E stained prostate histology section; black ink mark provided by pathologist indicates CaP extent.

[031] Figure 3(b) depicts gland segmentation boundaries.

[032] Figure 3(c) depicts a magnified view of the white box depicted in Figure 3(b). [033J higure J(Cl) depicts centroids ot cancerous glands betore MKt iteration.

[034] Figure 3(e) depicts centroids of cancerous glands after MRF iteration.

[035] Figure 4 depicts region growing examples.

[036] Figure 4(a) depicts seed pixel with CB.

[037] Figure 4(b) depicts IB, CB, and CR at step j in the growing procedure.

[038] Figure 4(c) depicts IB, CB, and CR at step j+1.

[039] Figure 5 depicts results from gland detection algorithm.

[040] Figures 5(a)-(d) depict H&E stained WMHSs.

[041] Figures 5(e)-(h) depict initial labels for _area = 0.25.

[042] Figures 5(i)-(l) depict PPMM GDS results with sensitivity adjusted to 72%.

[043] Figures 5(m)-(p) depict Potts GDS results with sensitivity adjusted to 72%.

[044] Figure 6(a) depicts receiver operator characteristic (ROC) curves for PPMM GDS (solid) and Potts GDS (dashed) algorithms averaged over 20 trials using randomized 3 -fold cross-validation. Dotted ROC curve indicates performance of initial classification based only on gland area.

[045] Figure 6(b) depicts area under corresponding ROC curves in Figure 6(a) for PPMM GDS and Potts GDS.

[046] Figure 7(a) depicts T2-w MR image with cancerous region as delineated by a radiologist. [047] Figure 7(b) depicts intensities indicate the probability of cancer for each pixel assuming uniform priors.

[048] Figure 7(c) depicts initial classification after MLE.

[049] Figures 7 (d)-(f) depict final classification after WICM for T equal to 0.9, 0.55, and 0.05, respectively.

[050] Figures 7 (g)-(i) depict final classification after WMPMMC forτ equal to 0.9, 0.55, and 0.05, respectively.

[051] Figure 8 depicts ROC curve indicating system performance over 15 studies.

[052] Figure 9 depicts a flowchart illustrating the 3 main steps in the CADx scheme for LI- based stratification of HER2+ BC histopathology.

[053] Figure 10(a) depicts luminance channels of two different HER2+ BC histopathology studies.

[054] Figure 10(b) depicts corresponding results for initial region-growing based lymphocyte detection.

[055] Figure 10(c) depicts preliminary Bayesian refinement showing detected BC nuclei and detected lymphocyte nuclei.

[056] Figure 10(d) depicts final lymphocyte detection result after the MRF pruning step.

[057] Figure 1 l(a) depicts HER2+ breast cancer histopathology image.

[058] Figure 1 l(b) depicts corresponding Voronoi Diagram constructed using the automatically detected lymphocyte centers as vertices of the graph. [059] Figure 1 1 (c) depicts corresponding Delaunay Triangulation constructed using the automatically detected lymphocyte centers as vertices of the graph.

[060] Figure 1 l(d) depicts corresponding Minimum Spanning Tree constructed using the automatically detected lymphocyte centers as vertices of the graph.

[061] Figure 12 depicts an example of current segmented region (CR), internal boundary (IB), and current boundary (CB) during a step of the region growing algorithm.

[062] Figure 13 depicts: (a), (g) H&E stained prostate histology sections with black ink mark indicating RT. (b), (h) Gland segmentation boundaries, (c), (i) Magnified views of white boxes from (b), (h). Centroids of cancerous glands before (d), (j) and after (e), (k) MRF iterations, (f), (i) Estimated cancerous regions (green) with HFT (yellow).

[063] Figure 14 depicts a solid ROC curve that describes the performance (obtained over four studies) of the initial gland classification (without MRF) as the probability threshold p varies from zero to one.

Detailed Description of Invention

[064] In accordance with the above objects and others, the invention is directed in part to a system and method for detecting biologically relevant structures wherein: 1) biologically relevant structures are identified and segmented, 2) morphological (and possibly other) features are extracted from the segmentations, 3) the segmented structures are classified using an MRF, and 4) similarly classified structures are consolidated into regions.

[065] In accordance with the above objects and others, the invention is directed in part to a system and method for detecting biologically relevant structures wherein: 1) gland segmentation is performed on the luminance channel of a color H & E stained WMHS, producing gland boundaries, 2) the system then calculates morphological features for each gland, 3) the features are classified, labeling the glands as either malignant or benign, 4) the labels serve as the starting point for the MRF iteration which then produces the final labeling, and 50 the cancerous glands are consolidated into regions. The system and method of the present invention can be applied to the analysis of variety of organs (including, but not limited to, ovarian, bladder, cervical, rectal, breast) imaged using a variety ol modalities (including, but not limited to, uitrasounα, radiography, PET, CT).

[066] An embodiment of the present invention is directed to a system specifically designed to operate at low-resolution (0.008mm per pixel) which will eventually constitute the initial stage of a hierarchical analysis algorithm, identifying areas for which a higher-resolution examination is necessary. A hierarchical methodology provides an effective means for analyzing high density data. Even at low resolutions, gland size and morphology are noticeably different in cancerous and benign regions.

[067] An embodiment of the present invention is directed to probabilistic pairwise Markov models (PPMMs). PPMMs formulate the LCPDFs in terms of pairwise probability density functions (PDFs), each of which models an interaction between two sites. They allow the creation of relatively sophisticated LCPDFs, increasing our ability to model complex processes. Furthermore, since the pairwise probabilities have meaningful interpretations, the functional forms, modeling capabilities, and expected behavior of the available LCPDFs become more apparent. Thus, PPMMs obviate the need for ad hoc representations and return the construction of the LCPDFs to the familiar Bayesian framework.

[068] An embodiment of the present invention is directed to a method of detection/grading of malignant wherein the initial stage of a hierarchical analysis algorithm indicates those areas for which a higher-resolution examination is necessary. At the proposed resolution the predominately visible structures are the glands. Gland size and morphology are noticeably different in cancerous and benign regions, even at low-resolution. In cancerous areas the glands tend to be smaller than those found in benign regions. Also, cancerous glands tend to be proximate to other cancerous glands. This information is included using Markov random fields (MRFs). However, unlike the preponderance of MRP strategies which rely on heuristic formulations, we introduce a straight-forward, mathematically rigorous formulation that allows the MRF to be modeled directly from the data. The cancerous glands are then formed into regions. The methodology of the present invention can be applied to a variety of gland based features and proximity measures. Accordingly, the system and method of the present invention can be used to assist pathologists in examination of histological sections and/or as an independent means of assessment. [069] In a preferred embodiment, the present invention incorporates analysis or gianα structure on low resolution images as the initial stage of a hierarchical approach to detection and grading of prostate cancer, rapid segmentation of gland boundaries and morphological features, and integration of probability distributions into the Markov Random Fields framework, allowing the interactions between neighboring glands to be modeled probabilistically.

[070] The invention is described more fully herein. All references cited are hereby incorporated by reference in their entirety herein.

Review of Markov Random Fields

Random Field Nomenclature

[071] Let the set 5 = {1, 2, . . . , N) reference N sites to be classified. Each site se S has two associated random variables: X 3 e A indicating its state (class) and Y 3 e R D representing its D- dimensional feature vector. Particular instances of X s and Y s are denoted by the lowercase variables χ s e Λ and y s €R D . Let X = (X 1 , X 2 , . . . , XN) and Y = (Y 1 , Y 2 , . . . , YN) refer to all random variables X s and Y s in aggregate. The state spaces of X and Y are the Cartesian products Ω, = A N and R°* N . Instances of X and Y are denoted by the lowercase variables χ = (X 1 , X 2 , . ■ • , XN) e Ω and y = (yi , y 2 , ■ ■ ■ , VN) e R 0 *^.

[072] Let G = [S, E) establish an undirected graph structure on the sites, where S and E are the vertices (sites) and edges, respectively. A clique c is any subset of S which constitutes a fully connected subgraph of G, i.e. each site in c shares an edge with every other site. The set e contains all possible cliques. A neighborhood η s is the set containing all sites that share an edge with s, i.e. η s = {r : r € S, r φ s, {r, s} G E). If P is a probability measure defined over Ω then the triplet (G, Ω, P) is called a random field.

[073] These concepts are best understood in the context of an example. The graph in Figure 1 has sites S = {1, 2, 3, 4, 5, 6} and edges E = {{1, 2} , {1, 4} , {1, 5} , {2, 3} , {2, 6} , {4, 5}}. The neighborhood of site 5, for example, is 77 5 = {1, 4}. There are six one-element cliques CI = {{1} , {2} , {3} , {4} , {5} , {6}}, six two-element cliques Q 2 = E, and one three-element clique 63 = {{1, 4, 5}}. The set e is the union of G 1 , Q 2 , and e 3 . The set of possible states for each X 3 is Λ = {6, w), where b and w represent black and white, respectively. The random variable Y s eK reflects the diameter of site s. For this example we have x = {w, b, b, b, w, b) and Y = (O- Q 1 LS 1 Le 1 IA l-I 1 I-S).

[074] Finally, we establish a convention for representing probabilities. Let P{ ) indicate the probability of event {•}. For instance, P(X s = χ s ) and P(X = χ ) signify the probabilities of the events {X s = χ s } and {X = χ}. Whenever appropriate we will simplify such notations by omitting the random variable, e.g. P(χ) ≡P(X=χ). Let p( ) represent a probability PDF; for example, p g might indicate a Gaussian PDF. The notations P(-) and p(-) are useful in differentiating P( χ s ) which indicates the probability that {X s = χ s } from p g {x s ) which refers to the probability that a Gaussian random variable assumes the value x s .

Maximum a Posteriori Estimation

[075] Given an observation of the feature vectors Y we would like to estimate the states X. The preferred method is maximum a posteriori (MAP) estimation which entails maximizing the following quantity over all xG Ω: oc P(y|x) P(x) . (2.2.1)

[076] The first term in (2.2.1) reflects the influence of the feature vectors. It can be simplified by assuming that all Y 3 are conditionally independent and identically distributed given their associated X s . This assumption implies that if the class X 3 of site s is known then 1) the classes and features of the remaining sites provide no additional information when estimating Y s and 2) the conditional distribution of Y s is identical for all s e S. As a result we have

^ (y|χ) = ^l** ) ses = H p f (y a \x a ) , (2.2.2) sζS

[077] where the use of the single PDF p/ indicates that P{y s \x s ) is identically distributed across 5. The second term in (2.2.1) reflects the influence of the class labels. In general, modeling this high-dimensional PDF is intractable. However, if the Markov property is assumed its formulation simplifies. This is the topic of the following subsections. Random Fields and the Markov Froperty

[0781 The random field (G 1 Ω, P) is a Markov Random Field (MRP) if its local conditional probability density functions satisfy the Markov property:

P{x s \x- s ) = P{x s \x η3 ) , (2.3.1)

[079] where χ - β = (xi , . . . , x 8 -i , x a+ ι , . . . , XN), *η s = (^(l) , • • ■ > x η.(\η M \)), Vs{ι) & S is the i th element of the set η s , and \η s \ is the cardinality of η s . Thus, the Markov property simplifies the forms of the LCPDFs.

[080] Given the LCPDFs for all s G S we can uniquely recover the corresponding joint probability density function P(x) using Brook's Lemma. However, it is unlikely that a group of arbitrarily selected LCPDFs will be consistent in the sense that they define a valid JPDF. That is, consistent LCPDFs are limited to certain functional forms. The Hammersley-Clifford theorem establishes the functional forms of LCPDFs that are both consistent and satisfy the Markov property.

Gibbs Markov Equivalence

[081] The connection between the Markov property and the JPDF of X is revealed by the Hammersley-Clifford (Gibbs-Markov equivalence) theorem. This theorem states that a random field (G, Ω, P) with P(χ) > 0 for all x G Ω satisfies the Markov property if and only if it can be expressed as a Gibbs distribution:

[082] where Z = ∑ xeΩ FLe e ^c ( x ) * s the normalizing constant and V c are positive functions, called clique potentials, that depend only on those χ s such that s£c. The following derivation reveals the simplified forms of the LCPDFs:

Ilcee., ^(x) ΣλεΛ ricee. ^c(Zi , ■ • ■ , x β -i , λ, x β+ i , . , XN ) rue. f cM = P(a; β |χ,J , (2.4.2)

Σλ Λ ricee, v c(xi , - - - , x s -i , λ, χ s +ι , . . . , XN) [083] where G s represents {c€ e : s € c) and C. s represents {c£ £ : s f c). i nis derivation also proves that all Gibbs systems are Markovian.

Probabilistic Pairwise Markov Models

[084] In this section we introduce a means for constructing consistent LCPDFs using probability density functions instead of potentials. We refer to the resulting MRFs as probabilistic pairwise Markov models (PPMMs). We then demonstrate how both Gaussian MRFs and MRFs based on the Potts model can be interpreted as PPMMs. Finally, we compare PPMMs to a subset of MRFs called strong-MRFs.

[085] A PPMM whose pairwise PDFs are modeled as Gaussian distributions will produce a Gaussian MRF (i.e. a MRF whose LCPDFs and JPDFs are Gaussian distributions). These PDFs need not be Gaussian; they can assume any bivariate distribution (e.g. Gamma, exponential, Poisson). This allows the creation of relatively sophisticated LCPDFs, increasing our ability to model complex processes. Furthermore, since the pairwise probabilities have meaningful interpretations, the functional forms, modeling capabilities, and expected behavior of the available LCPDFs become more apparent. Thus, the PPMMs of the present invention obviate the need for ad hoc representations and return the construction of the LCPDFs to the familiar Bayesian framework.

Assumptions and Derivations

[086] We begin by invoking two prevalent assumptions. First, we assume that only the pairwise dependencies among sites are significant. Second, we assume that the LCPDFs are stationary, i.e. is identical for each site s£ S having m = |?7 s | neighbors. With these assumptions (2.4.2) simplifies to the following:

r[x* l χϊ? J - U AeA l / i (λ)l l reil5 l / 2 (λ,i r ) > V A Λ )

[087] where V\ and Vi are the site-invariant potential functions for one- and two-element cliques, respectively. Furthermore, Vi is symmetric in the sense that V 2 (x s , χ r ) = Vi (x r , x s ); this symmetry is needed to ensure the stationarity of P{x s \x. Vs ) across all possible lattice structures. (For regular lattices, such as rectangular grids, Vi need not be symmetric.) [088 J Betore continuing we introduce a new notation. Let represent me stationary distribution m = |77 s |. The naming convention motivating Po\ m results from considering s as stationary site 0 and its neighbors η s as stationary sites 1 . . . m. This same convention will be used to represent other distributions involving these m+1 sites. For example, V $ \Q reflects the conditional probability of the state of stationary site 3 given the state of stationary site 0.

[089] Consider the following manipulation of Po\ m '-

[090] In essence, (3.1.2) provides a means for representing Po\ m in terms of the single pairwise distribution po,i, where Pi|o and Po are a conditional and marginal distribution of Po 11 . Furthermore, the similarities between (3.1.2) and (3.1.1), both in their forms and the assumptions used to derive them, suggest the following theorem.

[091] Let (G, Ω, P) be a random field such that P(X = x) > 0 for all x G Ω. If p o ,i{χ s , χ r ) is a PDF such that po,i( χ s , χ r) = po,i{ χ r, χ s ) and

P ( x Ix ϊ - P^ πreq. PiloH s. )

[092] for all s G S, then X is a MRF.

[093] Proof: It suffices to show that (3.1.3) can be expressed in the form given by (3.1.1). Consider the following one- and two-element potential functions:

and

V 2 (X 3 , Xr) = Pθ,l ( x s, x τ) ■ [094] where the symmetry ol po,i ensures that V2 is also symmetric, inserting these into (3.1.1) yields

p, , x _ Po(X 5 ) 1"1 ^ 1 n reη , po,i( χ ,, χ τ)

Z-Ae A Po(X,) 1 -"" 1 li r e,. Po,i (λ, xr)

Po(^ a ) Ilrei?. Pl |0 (^ I ^) Σλ €Λ Pθ(λ) U r s Pl|θ(»r |λ) '

Gaussian MRF as a PPMM

[095] Consider the case of a continuous MRF where X s € Λ ≡ R. Let po,i be a (necessarily symmetric) Gaussian distribution:

Po ,_ (z β , -c r ) = ^§^ exp {-^7 {x s 2 + 2Υx s x r + xl)} , (3.3.1)

[096] where σ € R + and 7 € [- 1 , 1]. Equation (3.3.1) reflects a bivariate Gaussian that has been rotated such that its principal axes coincide with the lines x s = x r and x s = —x r , i.e. a rotation of π/4. Its requisite marginal and conditional distributions are

} and

[097] Inserting these results into (3.1.3)yields the following:

[098] where ξ = 1 + 7 2 (|r/ β | - 1). The LCPDFs are themselves Gaussian distributions. Thus, by designating po,i as a Gaussian distribution, we have created an instance of a Gaussian MRF. Potts Model as a rrMM

[099] In this subsection we demonstrate how the Potts model, perhaps the most prevalent MRF formulation, can be interpreted as a PPMM. The potential functions of the Potts model are as follows:

v ^ - J ^ if |c| = 2 and x r = x S) c= {r, s} n * λ \

V c ^ - \ l otherwise, ^ ;

[0100] where β > 0 (except in rare instances). Since e β > 1, neighboring pixels with identical states will contribute more to the JPDF and their respective LCPDFs than neighboring pixels with differing states. The degree of contribution is a function of β, with greater values of β producing "smoother" solutions. To see this, consider the MAP estimation in (2.2.1). Increasing β increases the weight of the second term, further encouraging neighboring sites to share the same label.

[0101] The Potts model can be reformulated as a PPMM as follows:

if X s — X τ

( ,x s , x r λ ) = I |Λ|(e"+|Λ|-l) s ~ τ ,., < -

Po,i (a:*, Sr) = \ i ' . , . (3.5.2) otherwise.

( |Λ|(e0+|Λ|-l)

[0102] Note that this probability function is symmetric.

Connection to Strong-MRFs

[0103] A better understanding of PPMMs can be garnered by comparing them to strong- MRFs. Strong-MRFs are MRFs that, in addition to the Markov property, impose the following constraint: Vx e Ω, s G S

P {X s =x s \X r = x r : r ≠ s, r G A C S) = P{X s = x s \X r = x τ : r G η s D A) .

[0104] This states that a MRF with sites 5 is a strong-MRF if every A c S is itself a MRF with the previous graph structure over S persisting among the remaining sites in A. In general this is not true for MRFs since removing sites (i.e. summing them out of the JPDF) usually forms new interactions (edges). [0105] To understand the connection between FPMMs and strong-MKPs again consider a stationary MRP X with LCPDFs given by Bayes Law implies that

[0106] However, it is not necessarily true that the numerators (or denominators) are equal. This further implies that the following equality will not hold in general: p o ,i ( χ s , χ r ) = P(xs, Xr). However, this equality does hold if (and only if) X is a strong-MRP. Consequently, only when X is an strong-MRP can po,i (χ s , Xr) be replaced with P{x s , x r ). The overarching principle is that two different joint distributions can share an identical conditional distribution.

[0107] Others aspects of strong-MRPs also provide insight into PPMMs. Paget' s decomposition of strong-MRPs into pairwise clique probabilities agrees with our representation when X is a strong-MRP. However, our formulation remains valid when X is not a strong-MRP. Furthermore, work by Paget and Bishop has implied that the aforementioned decomposition is only possible for one- and two-element cliques. This suggests that an analogous extension of our PPMMs to higher-order cliques may not be possible.

Incorporation of asymmetric costs into MRF estimation

[0108] The present invention provides a means for incorporating asymmetric costs into MRP estimation. In an embodiment of the present invention, the present invention provides weighted MAP (WMAP) and weighted MPM (WMPM) estimation; they are generalizations of MAP and MPM estimation that allow certain decisions to be weighted more heavily than others. In another embodiment of the present invention, the present invention provides weighted ICM (WICM), a novel adaptation of ICM capable of WMAP estimation, and weighted MPM Monte Carlo (WMPMMC), an extension of MPM Monte Carlo capable of WMPM estimation. Thus, the present invention can be used for any MRP optimization routine. Review of Bayesian Risk

[0109] Given an observation of the feature vectors Y, we would like to estimate the states X. Bayesian classification advocates selecting the estimate xe Ω that minimizes the conditional risk (expected cost)

R(X|x, y) = £ [C (X, x) |y] = ∑ xeΩ C (x, x) P(x|y) , (5.2.1)

[0110] where E indicates expected value and C (x, x) is the cost of selecting labels x when the true labels are x. In the following subsections we will consider the two most prevalent cost functions used with MRPs.

Maximum a Posteriori Estimation [0111] The most ubiquitous cost function (though this cost is rarely expressed explicitly) is

CAMP (x, X) = ILe S I 1 - 5 , (5.3.1)

[0112] where δ is the Kronecker delta. Thus, a cost of 1 is incurred if one or more sites are labeled incorrectly. Inserting (5.3.1) into (5.2.1) yields

R MAP

= 1 - P(x|y) . (5.3.2)

[0113] Minimizing (5.3.2) over x is equivalent to maximizing the second term. Thus, we have maximum a posteriori (MAP) estimation, which advocates selecting the x that maximizes the a posteriori probability.

[0114] When the cardinality of Ω is small the search for the optimal x can be determined by testing all possible states. Since this is rarely the case, other means become necessary. Simulated annealing is an effective stochastic relaxation method. Unfortunately, its high computation complexity diminishes its usefulness. A far less computationally intensive and more popular technique is iterated conditional modes (ICM), a deterministic relaxation scheme. The ICM algorithm is predicated on the following reformulation of the a posteriori probability:

P(x|y) = P (x s |x. β> y) P (x_ s |y) (5.3.3) [0115] Increasing the lirst term in (5.3.3) necessarily increases -^(χ|yj. mis suggests a giooai optimization strategy that sequentially visits each site 565 and determines the label χ s e Λ that maximizes P(x s |χ- S , y). ICM converges to a local, and not global, maximum of P( χ |y); however, in many situations the local maximum is more appropriate. The ICM algorithm is as follows:

Iterated Conditional Modes

Input: Initial labeling x°

Output: Final labeling x fc after iteration k

1. k = 0

2. do

3. k = k + 1

4. x fc = x fc-1

5. for V s G S do , y s )

8. while x fc φ x fc-1

[0116] Typically, the initial labeling χ° is provided by the maximum likelihood estimate of P(y|χ). This probability becomes tractable under the typical assumption that all Y s are conditionally independent given their corresponding X s , i.e. P(y|x) = Il s es The ICM iteration usually converges in approximately six or seven iterations.

Maximum Posterior Marginals

[0117] As an alternative to MAP estimation, Marroquin suggested using the following cost function

CMPM ( X , X ) = ∑, 6S [1 -« ( Ϊ . -Ϊ.)] - (5-4.1)

[0118] This function counts the number of sites in x that are mislabeled. Inserting (5.4.1)into (5.2.1)yields

iWM (X| χ , y) = W5> - «S (χ β β )] ) p(χ|y)

∑ ∑ P(x|y)-∑ ∑ δ (x s -x s ) P(x|y) seS xζΩ seS xeQ

\S\ -∑ P(x s \y) . (5.4.2) sζS [0119] 1 he distributions are caiieα me posterior marginals. Minimizing ip.t.z; is equivalent to independently maximizing each of these posterior marginals with respect to its corresponding χ s . Hence, this type of estimation is termed maximum posterior marginals (MPM).

[0120] Since the cardinality of Ω is usually large, P(x s \y) can not be evaluated directly; and consequently, Marroquin proposed using a Monte Carlo method — MPM Monte Carlo (MPMMC). The basic rationale for his approach is as follows: Stochastic relaxation schemes such as the Metropolis algorithm Metropolis 1953 or the Gibbs Sampler form Markov chains whose equilibrium distribution is -P(x|y). Consider the following algorithm for the Gibbs Sampler:

Gibbs Sampler

Input: Initial labeling x°

Output: Final labeling x fc after iteration k

1. k = 0

2. do

3. k = k + l

4. x fc = x fc~1

5. for V s G 5 do

6. sample x£ from P(X s = ω|x^ s , y s )

7. end for

8. while k < oo

[0121] In theory, the final labels χ k , which represent a sample from the distribution P(χ|y), are independent of the starting conditions x°. However, since this algorithm must eventually terminate, some dependence upon initial conditions does exist. As with ICM, χ° is typically initialized to the MLE. Determining the number of iterations I needed for the Markov chain to reach equilibrium will depend upon the distribution P(χ|y). This is usually chosen in an ad hoc fashion. It is worth mentioning that the only difference between ICM and the Gibbs sampler (other than the termination criterion) occurs in step 6: ICM chooses the label that maximizesP(χ s |χ^ s , y s ), the Gibbs sampler selects the label randomly from this distribution.

[0122] After reaching equilibrium the amount of time the chain spends in any state x is given by P( χ |y). Thus, the posterior marginal can be estimated as follows: P(X 3 =ω\y) « τb ∑÷i'Li ό (x « -ω) , (5.4.3)

[0123] where m-l is a number of iterations past equilibrium needed to generate an accurate estimate. The value for m, like /, is typically chosen empirically. Given this estimate of P(x s \y), the maximizer can be determined by testing all possible states of x s .

Weighted Maximum a Posteriori Estimation

[0124] MAP estimation weights each classification decision equally. In this section we generalize MAP estimation to allow for asymmetric costs. However, before continuing we first introduce several definitions that will prove useful in the subsequent derivations. Let a : Λ → R^ define a weighting function. Let X be a random field that differs with X only with respect to its probability measure, which is defined as follows:

P(X=x) = ^P (X=X) FUs ^ (O . (6.1)

[0125] where the constant Z ensures that ∑ xe π P(X =x) = 1. For convenience we will use following simplified notation: P(χ) = P(X = x). Finally, note that (6.1) implies

[0126] To allow certain classification decisions to be weighted more heavily than others C M AP can be generalized as follows:

Cw M AP (x, x) = Pises a ( χ s) [1 - δ (Xs -Xs)] ■ (6.3)

[0127] Consequently, if the estimate x contains any erroneous labels it accrues a cost of u s es a ( χ s)- Inserting (6.3) into (5.2.1) yields

Rw MA p W*, y)

[0128] This result indicates that minimizing y) is equivalent to minimizing R MAP (X\9., y); and consequently, the optimal labeling is the MAP estimate of X. Since the α posteriori probability of X corresponds to the weighted α posteriori probability of X, we refer to this type of estimation as weighted MAP (WMAP) estimation. Note that ≡ 1, WMAP estimation reduces to MAP estimation.

[0129] Implementing the WMAP estimation of X is trivial: we simply perform MAP estimation on X. Consequently, we are free to use any MAP estimation scheme available for MRFs. For example, consider this following explicit implementation of WMAP using weighted ICM (WICM), an adaptation of ICM:

Weighted Iterated Conditional Modes Input: Initial labeling x°, weights α(ω) Output: Final labeling x fc after iteration k

1. k = 0

2. do

3. k = k + l

4. x k = X^ "1

5. for V s e S do

7. end for

8. while x fc φ x fc~x

[0130] Modifying ICM to perform WMAP estimation only requires replacing P(x s |χ- S , y s ) with reveals the final form used in step 6. Weighted Maximum Posterior Marginals

[0131] In this section we introduce a means for modifying MPM to allow certain classes to be weighted more heavily than others. We call the modified estimation scheme weighted MPM (WMPM). We begin by deriving WMPM. We then demonstrate how WMPM can be implemented using an adaption of the Monte Carlo method introduced by Marroquin. We refer to this adaptation as weighted MPM Monte Carlo (WMPMMC).

[0132] The cost function C M PM can be generalized as follows:

CH / MPM ( X , X ) = lives a M ∑s€s l 1 ~ δ C 7 - 1 )

[0133] In this case mislabeling a site whose true label is x s has an associated cost of u s es a ( χ s)- Notice that the penalty for mislabeling site s varies depending upon the labels of the remaining sites. Though it is possible to remove these dependencies by redefining the cost function as CWMPM2(X, X) = developing a scheme to perform the necessary estimation is problematic.

[0134] Inserting (7.1) into (5.2.1) yields

RWMAP ( X | χ , y ) p ( x l y )

[0135] Thus, the WMPM estimate of X is equivalent to the MPM estimate of X. This implies that any approach used for performing MPM estimation can be modified to perform WMPM estimation. For example, consider the WMPMMC estimation procedure. Altering this to perform WMPM estimation requires a change to the Gibbs sampler as follows: Weighted Liibbs Sampler

Input: Initial labeling x°, weights a(ω)

Output: Final labeling x fc after iteration k

1. k = 0

2. do

3. k = k + l

4. x fc = x*- 1

5. for V s G 5 do

6. sample z£ from (a(ω)/Z s )P(X s = ω\^ s , y s )

7. end for

8. while k < oo

[0136] The only difference between the Gibbs sampler and the weighted Gibbs sampler is the replacement in step 6 of P(x s \ χ -t s , y s ) with

[0137] where Z 5 ensures ∑ ω£Λ ^( ω l x - s> Vs) — 1- The remainder of the WMPMMC procedure is identical to that of the MPMMC algorithm.

Examples of Preferred Embodiments

[0138] Examples of the system and method of embodiments of the present invention are provided below. The examples show the use of PPMM and weighted MRF estimation.

[0139] In Embodiments I and II of this application, we evaluate the algorithms in the context of two MRF-based classification systems: 1) for detecting prostate cancer (CaP) in whole-mount histological sections (WMHSs) and 2) for detecting CaP in T2-weighted 3 Tesla in vivo magnetic resonance imaging (MRI). Specifically, we illustrate how WICM and WMPMMC can be used to vary classification performance, enabling the construction of receiver operator characteristic (ROC) curves.

Embodiment I: Application of PPMMs to CaP Detection on Prostate Histological Sections

[0140] In the present embodiment, we demonstrate our novel classification system for detecting cancerous glands on digitized whole-mount histological sections of the prostate. In the present embodiment of the present invention, we utilize the fact that cancerous glands tend to be proximate to other cancerous gianαs. i nis iniormauon is incorporated via a Markov ranαom neiα (MRF), which we formulate using both the Potts model and the PPMM of the present invention. To demonstrate the advantages of PPMMs of the present invention over typical Markov models (specifically the Potts model), we evaluate both in the context of an algorithm to detect prostate cancer (CaP) on digitized whole-mount histological sections (WMHSs) of radical prostatectomy specimens.

[0141] We first provide a basic outline of the algorithm and then discuss each component in detail.

Algorithm Overview

[0142] Figure 3(a) illustrates a prostate histological (tissue) section. The pinkish hue results from the hematoxylin and eosin (H&E) staining procedure. The black circle delimits the spatial extent of CaP as delineated by a pathologist; since it exists on the physical slide, it remains an indelible part of the digitized image. The numerous white regions are the lumens — holes in the prostate through which fluid flows — of secretory glands. Our system identifies CaP by leveraging two biological properties: 1) cancerous glands (and hence their lumens) tend to be smaller in cancerous compared to benign regions and 2) malignant/benign glands tend to be proximate to other malignant/benign glands.

Method:

[0143] The basic algorithm, illustrated in Figure 2, proceeds as follows: 1) The glands (or, more precisely, the gland lumens) are identified and segmented. Figure 3(b) illustrates the segmented gland boundaries in green. Figure 3(c) shows a magnified view of the white box in Figure 3(b). 2) Morphological features are extracted from the segmented boundaries. Currently, we only consider one feature: glandular area. 3) A Bayesian classifier — using the distribution Pf in (2.2.2) which models the glandular area — assigns a probability of malignancy to each gland. If the probability exceeds the predetermined threshold τ αreα , the gland is labeled as malignant; otherwise it is labeled benign. The blue dots in Figure 3(d) indicate the centroids of those glands classified as malignant. 4) The MRF algorithm refines these labels, thereby enforcing the precept that nearby glands tend to share the same class. The malignant labels resulting from the MKF iteration are shown in Mgure J^e;; me Diue αots indicate xne ceniroiαs oi those glands labeled as cancerous.

Specific Notation for CaP Detection Algorithm

[0144] We recapitulate the gland classification problem in terms of the MRF nomenclature. Let the set S- {1, 2, . . . , N} reference the N glands in a WMHS. The random variable Y s eR indicates the square root of the area of gland s. The label X s of gland s is one of two possible classes: X s eA ≡ {u>i, u>2}, where W 1 and w 2 indicate malignancy and benignity. Two glands are neighbors if the distance between their centroids is less than 0.7 mm.

Gland Segmentation

[0145] Since color information is not needed to identify gland lumens on digitized WMHSs, all segmentation is performed using the luminance channel in CIE Lab color space. The CIE Lab color space is known to be more perceptually uniform than the RGB color space. In the luminance images glands appear as regions of contiguous, high intensity pixels circumscribed by sharp, pronounced boundaries. Our procedure to segment these glands proceeds as follows: The luminance image is convolved with a Gaussian kernel at multiple scales σ g € {0.2, 0.1, 0.05, 0.025}mm to account for varying gland size. Peaks (maxima) in the smoothed images are considered candidate gland centers. These single pixel peaks serve as seeds for the following region growing procedure (which operates on the original image):

[0146] Step 1 : Initialize the current region (CR) to the specified seed pixel and establish a

12σ s x 12σ ff bounding box centered about it. Initialize the current boundary (CB) to the 8- connected pixels neighboring CR (Figure 4(a)).

Step 2: Identity the pixel in CB with the highest intensity. Remove this pixel from CB and add it to CR. Revise CB to include all 8-connected neighbors of the aggregated pixel which are not in

CR (Figure 4(a)-(c)).

Step 3: Define the internal boundary (IB) as all pixels in CR that are 8-connected with the pixels in CB. Compute the current boundary strength which is defined as the mean intensity of the pixels in IB minus the mean intensity of the pixels in CB.

Step 4: Repeat steps 2 and 3 until the algorithm attempts to add a pixel outside the bounding box.

Step 5: Identify the iteration step at which the maximum boundary strength was attained. Define the optimal region as CR at this step. [0147] The final segmented regions may overlap; this is resolved by discarding the region with the lower boundary strength.

Feature Extraction, Modeling, and Bayesian Classification

[0148] Gland area is a feature known to discriminate benign from malignant glands. Since we employ a Bayesian framework, we require estimates of the conditional probability density functions of gland area for both malignant ω\ and benign u>2 glands. Using the equivalent square root of gland area (SRGA), the pdfs and which correspond to the PDF p/ in (2.2.2), can be accurately modeled with a weighted sum of Gamma distributions:

Λi β , k. »> - «^j + (l-α) »«j , (9.4.1)

[0149] where y > 0 is the SRGA, a e [0, 1] is the mixing parameter, Zc 1 , Zc 2 > 0 are the shape parameters, θ\ , 0 2 > 0 are the scale parameters, and F is the Gamma function. A Bayesian classifier uses these PDFs to calculate the probability of malignancy for each gland. Those glands whose probabilities exceed the predetermined threshold τ area e [0, 1] are labeled malignant; the remainder are classified as benign. The centroids of glands labeled as malignant are shown with blue dots in Figure 2(d).

Modeling the Homogeneity of Neighboring Glands as a MRF

Probabilistic Pairwise Markov Models for Binary Classes

[0150] PPMMs require the probability po,i > or equivalently, P 1(O and po. For the binary classes the general forms of these probabilities are

and

Po(Zs) = [Po(W 1 ) P 0 (W 2 ) ] = [c 1 -c] . (9.6.2)

[0151] The required symmetry of Po 11 necessitates that c (1 -α) = (1 -c) b, yielding c = b/ (1+b-a). In the case of the Potts model we have α = e β / (l + e β ), b = 1/ (l + e β ), and c = 1/2. [0152] Values for α and b (or β lor the Potts model) can obtained trom training data. I nougn (9.6.1) is a nonparametric distribution in the sense that there is no assumed model, the limited degrees of freedom allow the use of parametric estimation techniques. Since maximum likelihood estimation is numerically untenable for MRFs, maximum pseudo-likelihood estimation (MPLE) is the preferred alternative. MPLE maximizes the product of the LCPDFs over all samples, and unlike MLE, does not require computing the intractable normalizing factor Z.

Weighted Iterated Conditional Modes for MAP Estimation

[0153] We defined the optimal states x as those that maximize the a posteriori probability. The literature provides several means for performing this maximization. We have selected iterated conditional modes (ICM), a deterministic (as opposed to stochastic) relaxation technique. However, since our system requires the ability to favor certain classification results (i.e. misclassifying a malignant lesion is more serious than misclassifying a benign one), we must slightly adapt the ICM algorithm to incorporate this capability. Before discussing this adaptation we review the basic ICM algorithm.

[0154] ICM is based on the following reformulation of (2.2.1):

P(x|y) = P(z s |x- s , y) P(x- s |y) ex P(x s ηs , y s ) P(x- s \y) . (9.6.3)

[0155] Increasing the first term of (9.6.3) necessarily increases P(x|y). Since both terms depend only on s and its neighborhood, they can be easily evaluated. This suggests a global optimization strategy that sequentially visits each site s e S and determines the label χ s e Λ that maximizes the first term of (9.6.3). In the case of binary classes this reduces to following decision:

χ = ( U 1 if P(ωι \x Vs , y s ) > T wιcm ( 9 6 4 ) s \ W 2 otherwise. v • • /

[0156] where T wιcm ICM, after several iterations, converges to a local maximum of (2.2.1). Since the maximization problem is non- convex, no realizable optimization strategy is guaranteed to identify the global maximum. However, in many situations, the local maximum is actually more appropriate. [0157] An examination ol (9.6.4) immediately suggests varying T W ιcm as a means lor favoring one class over the other. As T wιcm decreases, (9.6.4) will increasingly prefer ωγ. By contrast, increasing T wιcm will favor ω 2 . Since the value of T wιcm implicitly weights the importance of each class, we refer to the modified algorithm as weighted iterated conditional modes (WICM). Figure 3(e) illustrates the result of applying WICM to the gland labels in Figure 3(d); the MRF prior was modeled using our PPMM. It is worth noting that the WICM algorithm can be derived in a more general and mathematically rigorous fashion using Bayesian cost functions.

Experimental Results

[0158] In this section we qualitatively and quantitatively compare the effectiveness of the PPMM and the Potts model by alternately incorporating both into the gland detection algorithm. We begin by discussing the dataset over which the CaP detection performance will be evaluated. Next we address the specific procedures used to train and test the system. Finally, we present the results and discuss their significance.

Data

[0159] The dataset consists of 20 prostate histological sections from 19 patients at two separate institutions (University of Pennsylvania and Queens University in Canada). In 13 cases the cancerous extent was delineated by a pathologist on the physical slide, and consequently, remains in the digital image. An example of this was provided in Figure 3(a). In seven cases the CaP extent was determined after digitization. An example of this is shown in Figure 3(c); the black line demarcating CaP extent is overlaid on the image and was not present during processing.

Experimental Design

Training

[0160] The training procedure begins by segmenting the glands in each training image and then calculating their areas and centroids. Glands whose centroids fall within the pathologist provided truth are labeled as malignant; otherwise they are labeled benign. A MLE procedure uses these labeled samples to estimate the parameters for the mixtures of Gamma distributions used to model p/. I he graph structure connecting tne gianαs is αetermineu lrυm me giαiiu centroids: two glands share an edge if the distance between their.centroids is less than 0.7 mm. A MPLE procedure uses the graph structure and gland labels to estimate parameters for both the PPMM (α and b) and Potts model (/3).

Testing

[0161] For each test image the algorithm segments the glands and then extracts their areas and centroids. The system uses these areas to determine the probability of malignancy for each gland. Those glands whose probability of malignancy is greater than τ αreα = 0.25 (chosen empirically) are classified as malignant; otherwise they are classified as benign. These classification results are passed to the MRF iteration — weighted iterated conditional modes. A graph structure over the glands is established using the methodology described in the training procedure. In addition to this graph structure, the WICM algorithm requires the conditional distribution P(x 3 ηa ); this distribution is provided by either the PPMM or Potts model. Consequently, we have two gland detection systems which we will refer to as PPMM GDS and Potts GDS. The final gland labels from both systems will depend upon the threshold T wιcm used in the WICM alogirithm. Varying T wιcm from zero to one will allow us to adjust classifier performance, yielding receiver operator characteristic (ROC) curves.

[0162] It is worth mentioning that classification performance could conceivably be adjusted by fixing T wιcm = l/2 (i.e. using ICM) and varying τ area . That is, since modifying τ αreα alters the initial conditions, it can cause ICM to converge to a different local maximum of the MAP estimate. In fact, such an approach has been used elsewhere. However, the success of this methodology is predicated on the assumption that the individual modes (local maxima) of the a posteriori distribution correspond to meaningful classification results. We are aware of no justification for this argument. By contrast, the rationale governing WICM can be firmly established using Bayesian risk.

Results

Qualitative Results

[0163] The training and test sets were established using a leave-one-out strategy. The performances of the systems were varied by adjusting T wιcm (for each algorithm independently) until their respective sensitivities approximately equaled 0.72. l hat is, when ' ±w ' icm =U.8 the PPMM GDS detected 72 percent of the malignant glands; to achieve the same sensitivity, the Potts GDS required T wιcm = 0.6. Figures 5(i)-(l) show results for four histological slices using PPMM GDS. Figures 5(m)-(p) provide results for the same slices using the Potts GDS. The blue dots indicate the centroids of those glands labeled as malignant. The black lines enclose the spatial extent of the CaP. Both systems demonstrate remarkable proficiency, especially considering that they only examine glandular area while disregarding factors such as nuclei density, color information from the H&E staining, morphology, and glandular arrangement. Notice that every cancerous delineation (black mark) of appreciable area contains multiple detections (blue dots).

Quantitative Results

[0164] The training and testing data — which were identical for both systems — consisted of 20 independent trials using randomized 3 -fold cross-validation. By varying T wιcm from zero to one, ROC curves were generated for each trial. Figure 6(a) illustrates the average of these ROC curves for both systems (r α7 - = 0.25). The true positive rate is the ratio of malignant glands correctly classified to the total number of malignant glands. The false positive rate is the ratio of benign glands incorrectly classified to the total number of benign glands. For reference, this figure also includes the ROC curve (dotted line) indicating the classification performance if no MRF were used. That is, prior to the MRF iteration a Bayesian classifier assigns an initial probability of malignancy to each gland based solely on its area; the performance capability of this classifier varies as τ area ranges from zero to one, yielding the dotted ROC curve in Figure 6(a). Figure 6(b) provides a bar plot of the average area under the ROC curves for the PPMM GDS (0.87) and Potts GDS (0.83) algorithms. The error bars indicate the standard error of this statistic. The difference in performance is statistically significant.

Embodiment II: CaP Detection in Multi-protocol MRI

[0165] An embodiment of the present invention provides a classification system for detecting CaP in ex vivo MRI. This system uses T2-weighted (T2-w) 3 Tesla MRI from 15 2D slices. Figure 7(a) illustrates a T2-w MR image; the green overlay indicates the cancerous extent as specified by a radiologist. Let the set S = {1, 2, . . . , N} reference the N pixels in the T2-w MR image that reside within the prostate. The random vector Y s € K represents the 21 features associated with pixel s. Each multi-dimensional distribution / {y s \x s ) (one tor each class) was modeled using 21 one-dimensional histograms. The neighborhood η s of a pixel s is the typical 8- connected region. All results were produced using leave-one-out cross-validation.

Results

[0166] The intensities in Figure 7(b) indicate the probability of malignancy at each pixel based solely assuming uniform priors (i.e. MLE). Those pixels labeled whose probability exceeds 1/2 are indicated by the red overlay in Figure 7(c). Figures 7(d)-(f) illustrate the malignant labels after WICM for T equal to 0.9, 0.55, and 0.05, respectively. Figures 7(g)-(i) depict the results when, instead of WICM, we use WMPMMC. The black ROC curve in Figure 8 indicates system performance over all 15 studies as T varies from 0 to 1. The AUCs for WICM and WMPMMC are 0.800 and 0.812, respectively.

Embodiment III: Detecting Lymphocytic Infiltration in Breast Histology

[0167] Breast cancer (BC) is the second leading cause of cancer-related deaths in women, with more than 182,000 new cases of invasive BC predicted in the United States for 2008 alone. Researchers have shown that the presence of lymphocytic infiltration (LI) in histopathology is a viable prognostic indicator for various cancers, including breast cancers that exhibit amplification of the HER2 gene (HER2+ BC).

[0168] An embodiment of the present invention provides a computer-aided diagnosis (CADx) scheme (Figure 9) to automatically detect and grade the extent of LI in digitized HER2+ BC histopathology images.

Automated Lymphocyte Detection Scheme

[0169] We first attempt to identify candidate image locations that could represent centers of lymphocytic nuclei. A region-growing algorithm exploits the fact that lymphocyte nuclei in the luminance channel are identified as continuous, circular regions of low intensity circumscribed by sharp, well-defined boundaries (Figure 10(b)). However, a number of BC nuclei are also detected (i.e. false positives) along with the lymphocyte nuclei. [0170] The initial lymphocyte detection is retineα Dy utilizing Maximum a Posteriori estimation to incorporate domain knowledge regarding lymphocyte size, luminance, and spatial proximity. Specifically, a Bayesian model is used in conjunction with object size and luminance to classify individual objects as either lymphocyte or BC nuclei (Figure 10(c)).

[0171] Since LI is characterized by a high density of lymphocytes, the detection result is further refined by utilizing a Markov Random Field (MRF) model, which ensures that an object is more likely to be labeled as a lymphocyte nucleus if it surrounded by other lymphocyte nuclei. The MRF model is iterated until convergence by using an Iterated Conditional Modes algorithm (Figure 10(d)).

Graph-Based Feature Extraction

[0172] After discarding the BC nuclei in a histopathology image, the remaining lymphocyte nuclei are used as vertices in the construction of 3 graphs, a Voronoi Diagram (Figure 1 1 (b)), Delaunay Triangulation (Figure 11 (c)), and Minimum Spanning Tree (Figure 1 1 (d)). A total of 25 graph-based features are extracted for each image. An additional 25 nuclear features are extracted directly from the lymphocyte nuclei by calculating statistics that describe the number and density of nuclei in the image. The graph-based and nuclear features are combined to create an architectural feature set of 50 features for each histopathology image.

Support Vector Machine Based Classification

[0173] For a total of 41 HER2+ BC histopathology samples, a SVM classifier is used to calculate the ability of the architectural feature set to discriminate between images with high and low levels of LI. Over 100 trials of randomized 3-fold cross-validation, the architectural features produced a mean classification accuracy of 89.71% ± 2.84%. For comparison, the feature extraction and classification were repeated for manually detected lymphocyte nuclei, producing a classification accuracy of 99.59% ± 0.92%.

[0174] The aforementioned embodiments which encompass the construction of a computer- based system capable of identifying and assessing the biological significance of suspicious (abnormal) regions/objects in medical images by segmenting and classifying them using PPMMs and weighted estimation can be extended to the analysis of variety of organs (including, but not limited to, ovarian, bladder, cervical, rectal, breast; imaged using a variety or modalities (including, but not limited to, ultrasound, radiography, PET, CT).

Embodiment IV: Detection of Prostate Cancer from Whole-Mount Histology Images Using

Markov Random Fields

[0175] In a certain embodiment of the present invention, our CAD algorithm proceeds as follows: 1) gland segmentation is performed on the luminance channel of a color H&E stained WMHS, producing gland boundaries, 2) the system then calculates morphological features for each gland, 3) the features are classified, labeling the glands as either malignant or benign, 4) the labels serve as the starting point for the MRF iteration which then produces the final labeling, and 5) the cancerous glands are consolidated into regions.

Methodology:

Gland Segmentation

[0176] In the luminance channel of histological images glands appear as regions of contiguous, high intensity pixels circumscribed by sharp, pronounced boundaries. To segment these regions we adopt a routine first used for segmenting breast microcalcifications[81]. We briefly outline this algorithm. First define the following:

[0177] 1) current region (CR) is the set of pixels representing the segmented region in the current step of the algorithm, 2) current boundary (CB) is the set of pixels that neighbor CR in an 8-connected sense, but are not in CR, and 3) internal boundary (IB) is the subset of pixels in CR that neighbor CB. These definitions are illustrated in Figure 12. The growing procedure begins by initializing CR to a seed pixel assumed to lie within the gland. At each iteration CR expands by aggregating the pixel in CB with the greatest intensity. CR and CB are updated, and the process continues. The algorithm terminates when the L 00 norm from the seed to the next aggregated pixel exceeds a predetermined threshold. That is, the L norm establishes a square bounding box about the seed; the growing procedure terminates when the algorithm attempts to add a pixel outside this box. During each iteration the algorithm measures the boundary strength which is defined as the average intensity of the pixels in IB minus the average intensity of the pixels in CB. After the growing procedure terminates, the region with the greatest boundary strength is selected. Seed pixels are established by rinding peaKS in me image aπer uaussian smoothing. Since gland sizes can vary greatly, we smooth at multiple scales, each of which is defined by the sigma σ g e {0.2,0.1,0.05,0.025} mm of a Gaussian kernel. The growing procedure operates on the original image. The length / of each side of the bounding box used for terminating the segmentation step is tied to the scale: / = 12σ g . The final segmented regions may overlap. In this event the region with the highest boundary measure is retained.

[0178] Figures 13(a) and 13(g) are H&E stained WMHSs with black ink marks providing a rough truth (RT) of CaP extent. Gland segmentation results are shown in Figures 13(b) and 13(h). Figures 13(c) and 13(i) provide magnified views of the regions of interest in Figures 13(b) and 13(h). The centroids of glands whose probability of malignancy exceeds p = 0.15 are marked with green dots in Figures 13(d) and 13(J). This labeling is refined by the MRF iteration, producing the gland centroids shown in Figures 13(e) and 13(k). Figures 13(f) and 13(1) show the aggregation of cancerous glands into regions (green) along with a high-fidelity truth (HFT) of CaP extent (yellow).

Feature Extraction, Modeling, and Bayesian Classification

[0179] Gland area is used to discriminate benign from malignant glands. Since we employ a Bayesian framework, we require estimates of the conditional probability density functions (pdfs) of gland area for both malignant ω m and benign ω h glands. Using the equivalent square root of gland area (SRGA), the pdfs f(y \ ω m ) and f(y \ ω b ) can be accurately modeled with a weighted sum of gamma distributions:

[0180] where y>0 is the SRGA, /lφ,l] is the mixing parameter, A:, ,A: 2 >0 are the shape parameters, #, ,# 2 >0 are the scale parameters, and F is the Gamma function. Note, we use / to indicate a continuous pdf and p to denote a discrete probability mass function (pmf). A Bayesian classifier uses these pdfs to calculate the probability of malignancy for each gland. l hose glands whose probabilities exceed the predetermined threshold p are labeled malignant; the remainder are classified as benign (Figures 13(d) and 13(J)).

Improved Classification Using Markov Random Fields

[0181] In addition to glandular features such as area, a highly indicative trait of cancerous glands is their proximity to other cancerous glands. This can be modeled using MRFs.

Formulation of Gland Proximity as a MRF

[0182] Let . , s N } represent a set of N unique sites corresponding to the N segmented glands. Let each site seS have an associated random variable X s e{ω m b } indicating its state as either malignant or benign. To refer collectively to the states of all glands we have X={X, : seS}. Each state X s is unknown; we only observe an instance of the random variable

F r e R D representing the D dimensional feature vector associated with gland s . Though our algorithm is extensible to any number of features, currently D=I with K 5 being the SRGA. To collectively refer to the entire scene of feature vectors we have Y={Y S : seS}.

[0183] Consider the undirected graph G= {S, E} , where the set S of sites represents the vertices and the set E contains the edges. A local neighborhood η s is defined as follows: η s = {r : reS,r ≠ s, {r,s)eE} . The set of all local neighborhoods establishes a neighborhood structure: η = {η s : s e S} . A clique is a set of the vertices of any fully connected subgraph of G . The set C contains all possible cliques. These concepts are best understood in the context of an example. The graph in Figure 1 has sites 5 = {1,2,3,4,5,6} and edges

E = {{l ,2}, {l ,4}, {l ,5}, {2,3}, {2,6}, {4,5}} . The neighborhood of site 5 , for example, is η 5 = {l ,4} . There are six one-element cliques C 1 = {{l }, {2 }, {3 }, {4}, {5 }, {ό}} , six two-element cliques C 2 =E , and one three-element clique C 3 ={{l,4,5}} . The set C is the union of these three sets. The state X 5 of each site is either black or white, i.e. Λ={b,w}. Our specific neighborhood structure is determined by the distance between gland centroids. If m s denotes the centroid of gland s , then r e η s if . [0184] To simplify notation we use PrjΛf, = x r ,X s = x s ] = p[x r ,x s ) for indicating the probability of a specific event, where x r ,x s e{ω m b }. If X is a MRP with respect to the neighborhood structure η , then X satisfies the local Markov property p[x s \ x_ s )=p\x s \ x η j, where x_ s indicates x without x s and x η ={x ς : s e 77J. Additionally, X is a MRF with respect to η if and only if p(x) is a Gibbs distribution [82]: V c are nonzero functions that depend only on those x s for which sec . The local conditional probabilities follow directly:

J x ,, )- M-*J A*) *"' ' ' Υ ωeλMX) ~ < IaeΛM*- ) c πstc v c( χ ) cτ seicκ( χ ) ΓKM

(H) c stc ωeAc sec ωeAc sec

[0185] This distribution has the same form as the Gibbs distribution for p(x), but now the product is only over those cliques c that contain s .

Integration of Data Derived PMFs into the MRF

[0186] Since it is difficult to derive Gibbs distributions that model a set of training data, generic models are usually assumed. The most prevalent formulation is the Potts model which is defined as follows for two-element cliques:

V ( r r \- \ e'P if *r = *, and {r, j } e C n r >

[e p if x r ≠ x s and Jr 1 JJe C.

[0187] The Potts model disregards all cliques having more or less than two elements, i.e. if \c\ ≠ 2 we have V c {x) = 1 , where | | signifies cardinality. Such generic models are unnecessary; assuming that all X s are i.i.d. and all X r given X s are i.i.d. for every reη s , we can determine the appropriate V c directly from the data. To our knowledge, the following equations represent the first proposed means for incorporating arbitrary pmfs into the MRF structure: y*) = Λ.*J for {r,5}€C. (14)

[0188] The functions V c for higher-order cliques are identically one. The validity of (13) and (14) can be seen by inserting them into (1 1) :

pMYl p ( χ r \ χ s)

[0189] The determination of p(x s ) and p(x r ,x s ) from training data is straight-forward. For example, consider the randomly selected two-element clique {r,s} where both sites are malignant. The probability V^ r ^{ω m m )=p(ω m m ) can be found by examining all permutations of two neighboring glands and determining the percentage in which both are cancerous. The pmf p(x s ) is the marginal mass function of p(x r ,x ^ ). Since p(x r ,x ^ ) is symmetric, both marginals are identical.

Label Estimation and Aggregation

[0190] The goal is to estimate the hidden states X given the observations Y using maximum a posteriori (MAP) estimation, i.e. maximizing p(x \ y) over x . Bayes laws yields p{ χ I yfcfiy I χ )p{ χ ) > where oc signifies proportionality. The Iterated Conditional Modes (ICM) [?] algorithm indicates that the maximization of p(x | y) need not occur at all sites simultaneously; we can perform MAP estimation on each site individually by maximizing p[x s I x η ,yj∞ f{y s | x s )p[x s ,x η J. After estimating each individual state X s , the entire scene of states X is updated. The process iterates until convergence, usually requiring only five or six iterations (Figures 13(e) and 13(k)).

[0191] Our ultimate objective is to delineate the spatial extent of the cancerous regions. Following the neighborhood structure defined above, each gland centroid can be considered the center of a disk of diameter d . If the disks of two centroids overlap, they are considered neighbors. This leads to a simple formulation for cancerous regions: the union of all disks of diameter d centered at the centroids of the malignant glands (Figures 13(f) and 13(1)).

Results:

[0192] Data and CAD Training: The data consists of four H&E stained prostate WMHSs obtained from different patients. An initial pathologist used a black marker to delineate a very rough truth (RT) of CaP extent. An second pathologist performed a more detailed annotation of the digitized slices, producing a high fidelity truth (HFT). The digital images have a resolution of 0.01 mm 2 per pixel. The approximate image dimensions are 2.1x3.2 cm, i.e. 2100x3200 pixels. The training step involves estimating the parameters for the SRGA pdfs f(y \ ω m ) and f(y \ ω b ) using (10) and determining the MRF pmfs in (13) and (14) . The training/testing procedure uses a leave-one-out strategy.

[0193] Quantitative Results: We first assess the ability of the CAD system to discriminate malignant and benign glands. The quality of the gland segmentation is implicit in this performance measure. A gland is considered cancerous if its centroid falls within the HFT. The performance of the classification step above varies as the threshold p increases from zero to one, yielding the receiver operator characteristic (ROC) curve (solid) in Figure 14. Since the resulting classification serves as the initial condition for the MRF iteration, the MRF performance also varies as a function of p , producing the dashed ROC curve in Figure 14. The operating points for the qualitative results in Figure 13 are shown by the ring and dot (p = 0.15 ). We next measure the accuracy of the final CAD regions produced at this operating point by comparing them with the HFT. The sensitivity (the ratio of the cancerous area correctly marked to the total cancerous area) and specificity (the ratio of benign area correctly marked to the total benign area) are 0.8670 and 0.9524 , respectively.

[0194] Figure 14 depicts a solid ROC curve that describes the performance (obtained over four studies) of the initial gland classification (without MRF) as the probability threshold p varies from zero to one. The classification results at each threshold p are passed to the MRF iteration, yielding the dashed ROC curve. The operating points for the qualitative results in Figure 13 are shown by the ring and dot ( p = 0.15 ).

[0195] Qualitative Results: Qualitative results for two WMHSs were previously presented as Figure 13.

Conclusion

[0196] Over a cohort of four studies we have demonstrated in the present embodiment a simple, powerful, and rapid —requiring only four to five minutes to analyze a 2100x3200 image on 2.4 Ghz Intel Core2 Duo Processors— method for the detection of CaP from low-resolution whole-mount histology specimens. Relying only on gland size and proximity, the CAD algorithm highlights the effectiveness of embedding physically meaningful features in a probabilistic framework that avoids heuristics. The above embodiment introduces a novel method for formulating data derived pmfs as Gibbs distributions, obviating the need for generic MRF models.

References

[1] Scott Doyle, Carlos Rodriguez, Anant Madabhushi, John Tomaszeweski, and Michael Feldman, "Detecting prostatic adenocarcinoma from digitized histology using a multi-scale hierarchical classification approach.," ConfProc IEEE Eng Med Biol Soc, vol. 1, pp. 4759^762, 2006.

[2] S. Doyle, A. Madabhushi, M. Feldman, and J. Tomaszeweski, "A boosting cascade for automated detection of prostate cancer from digitized histology,"

MICCAI, pp. 504-51 1, 2006.

[3] M.A. Roula, A. Bouridane, F. Kurugollu, and A. Amira, "A quadratic classifier based on multispectral texture features for prostate cancer diagnosis," in Signal Processing and Its Applications, 2003. Proceedings. Seventh International Symposium on, July 2003, vol. 2, pp. 37-40 vol.2.

[4] James Diamond, Neil H Anderson, Peter H Bartels, Rodolfo Montironi, and Peter W Hamilton, "The use of morphological characteristics and texture analysis in the identification of tissue composition in prostatic neoplasia.," Hum Pathol, vol. 35, no. 9, pp. 1 121-1 131, Sep 2004.

[5] Reza Farjam, Hamid Soltanian-Zadeh, Kourosh Jafari-Khouzani, and Reza A Zoroofi, "An image analysis approach for automatic malignancy determination of prostate pathological images.," Cytometry B Clin Cytom, vol. 72, no. 4, pp. 227-240, JuI 2007.

[6] S Naik, S Doyle, A Madabhushi, J Tomaszewski, and M Feldman, "Gland segmentation and gleason grading of prostate histology by integrating low-, high-level and domain specific information," Workshop on Microscopic Image Analysis with Applications in Biology, 2007. [7] A. Tabesh, M. Teverovskiy, Ho-Yuen Fang, V. F. KΛimar, D. Verbel, A. Rotsianti, ana U. Saidi, "Multifeature prostate cancer diagnosis and gleason grading of histological images," /EEE Transactions on Medical Imaging, vol. 26, no. 10, pp. 1366-1378, Oct. 2007.

[8] S. Doyle, S. Agner, A. Madabhushi, M. Feldman, and J. Tomaszewski, "Automated grading of breast cancer histopathology using spectral clustering with textural and architectural image features," in Proc. 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro ISBI 2008,

14-17 May 2008, pp. 496-499.

[9] M. Roula, J. Diamond, A. Bouridane, P. Miller, and A. Amira, "A multispectral computer vision system for automatic grading of prostatic neoplasia," in Proc. IEEE International Symposium on Biomedical Imaging, 7-10 July 2002, pp. 193-196.

[10] Po-Whei Huang and Cheng-Hsiung Lee, "Automatic classification for pathological prostate images based on fractal analysis.," /EEE Trans Med Imaging, vol. 28, no. 7, pp. 1037-1050, JuI 2009.

[1 1] T. N. Pappas, "An adaptive clustering algorithm for image segmentation," /EEE Transaction on Signal Processing, vol. 40, no. 4, pp. 901-914, April 1992.

[12] C. A. Bouman and M. Shapiro, "A multiscale random field model for bayesian image segmentation," /EEE Transactions on Image Processing, vol. 3, no. 2, pp. 162-]77, March 1994.

[13] J. Besag, "On the statistical analysis of dirty pictures," Journal of the Royal Statistical Society. Series B (Methodological), vol. 48, no. 3, pp. 259-302, 1986.

[14] M. A. T. Figueiredo and J. M. N. Leitao, "Unsupervised image restoration and edge location using compound gauss-markov random fields and the mdl principle;' /EEE Transactions on Image Processing, vol. 6, no. 8, pp. 1089-1 102, Aug. 1997.

[15] R. Paget and 1. D. Longstaff, "Texture synthesis via a noncausal nonparametric multiscale markov random field," /EEE Transactions on Image Processing, vol. 7, no. 6, pp. 925-931, June 1998. [16J AJexey Zalesny and Luc J. Van (Jool, "A compact moαei tor viewpoint dependent texture synthesis," In SMILE Workshop, London, UK, 2001, pp. 124-143, Springer- Verlag.

[17] James Monaco, John E. Tomaszewski, Michael D. Feldman, Mehdi Moradi, Parvin Mousavi, Alexander Boag, Chris Davidson, Purang Abolmaesumi, and Anant Madabhushi, "Probabilistic pairwise markov models: application to prostate cancer detection," Medical Imaging 2009: Image Processing, vol. 7259, no. I, pp. 725903, 2009.

[18] S. Geman and D. Geman, "Stochastic relaxation, gibbs distribution, and the bayesian restoration of images," IEEE Transactions on Pattern Recognition and Machine Intelligence, vol. 6, pp. 721-741, 1984.

[19] J. Yedidia, W. Freeman, and Y. Weiss, "Generalized belief propagation," Advances in Neural Information Processing Systems, pp. 689-695, 2000.

[20] Y. Boykov, O. Veksler, and R. Zabih, "Fast approximate energy minimization via graph cuts," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 23, no. II, pp. 1222- 1239, Nov. 2001.

[21] R.C. Dubes, A.K. Jain, S.G. Nadabar, and CC. Chen, "Mrf model-based algorithms for image segmentation," in Proceedings ofthelOth International Conference on Pattern Recognition, 1990, vol. i, pp. 808-814 voU.

[22] R. Szeliskl, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov, A. Agarwala, M. Tappen, and C. Rother, "A comparative study of energy minimization methods for markov random fields with smoothness-based priors," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 6, pp. 1068-1080, June 2008.

[23] J. Marroquin, S. Miller, and T. Poggio, "Probabilistic solution of ill-posed problems in computational vision," Tech. Rep. AIM-897, MIT Artificial Intelligence Laboratory, Mar. 6 1987.

[24] J. L. Marroquin, B. C. Vemuri, S. Botello, E. Calderon, and A. Femandez-Bouzas, "An accurate and efficient bayesian method for automatic segmentation of brain mri," IEEE Transactions on Medical Imaging, vol. 21, no. 8, pp. 934-945, Aug. 2002.

[25] Jonathan Chappel ow, Satish Viswanath, James Monaco, Mark Rosen, John Tomaszewski, Michael Feldman, and Anant Madabhushi, "Improving supervised classification accuracy using non-rigid multimodal image registration: detecting prostate cancer, zuu», vol. ovo, p. ovnυv; SPIE.

[26] R. Duda, P. Hart, and D. Stork, Pattern Classification, John Wiley & Sons, 2001.

[27] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, "Equation of state calculations by fast computing machines," 1. Chem.

Phys., vol. 21, pp. 1087-1092, June 1953.

[28] T. N. Pappas, "An adaptive clustering algorithm for image segmentation," /EEE Transaction on Signal Processing, vol. 40, no. 4, pp. 901-914, April 1992.

[29] J. Besag, "On the statistical analysis of dirty pictures," Journal of the Royal Statistical Society. Series B (Methodological), vol. 48, no. 3, pp. 259-302, 1986.

[30] M. A T. Figueiredo and 1. M. N. Leitao, "Unsupervised image restoration and edge location using compound gauss-markov random fields and the mdl principle," /EEE Transactions on Image Processing, vol. 6, no. 8, pp. 1089-1 102, Aug. 1997.

[3 I] R. Paget and J. D. Longstaff, "Texture synthesis via a noncausal nonparametric multiscale markov random field," /EEE Transactions on Image Processing, vol. 7, no. 6, pp. 925-931, June 1998.

[32] Alexey Zalesny and Luc J. Van Gool, "A compact model for viewpoint dependent texture synthesis," in SMILE Workshop, London, UK, 2001. pp. 124-143, Springer-Verlag.

[33] Y. Zhang, M. Brady, and S. Smith, "Segmentation of brain rnr images through a hidden markov random field model and the expectation-maximization algorithm," /EEE Transactions on Medical Imaging, vol. 20, no. I, pp. 45-57, Jan 2001.

[34] AA Farag, AS. ΕI-Baz, and G. Gimel'farb, "Precise segmentation of multimodal images," /EEE Transactions on Image Processing, vol. 15, no. 4, pp. 952-968, April 2006.

[35] H.D. Li, M. Kallergi, L.P. Clarke, V.K. Jain, and R.A. Clark, "Markov random field for tumor detection in digital mammography," Medical Imaging, IEEE Transactions on, vol. 14, no. 3, pp. 565-576, Sep 1995.

[36] S. Awate, T. Tasdizen, and R. Whitaker, "Unsupervised texture segmentation with nonparametric neighborhood statistics," in Computer Vision ECCV, 2006, pp. 494-507. [37] C. A Bouman and M. Shapiro, "A multiscale random Held model lor Dayesian image segmentation," IEEE Transactions on Image Processing, vol. 3, no. 2, pp. 162-177, March 1994.

[38] J. Wu and A C. S. Chung, "A segmentation model using compound markov random fields based on a boundary model," IEEE Transactions on Image Processing, vol. 16, no. I, pp. 241- 252, Jan. 2007. [12] J. Besag, "Spatial interaction and the statistical analysis of lattice systems," Journal of the Royal Statistical Society. Series B (Methodological), vol. 36, no. 2, pp. 192-236, 1974.

[39] S. Geman and D. Geman, "Stochastic relaxation, gibbs distribution, and the bayesian restoration of images," IEEE Transactions on Pattern Recognition and Machine Intelligence, vol. 6, pp. 721-741, 1984.

[40] A. A. Efros and T.K. Leung, "Texture synthesis by non-parametric sampling," The Proceedings of the Seventh IEEE International Conference on Computer Vision, vol. 2, pp. 1033-1038, 1999.

[4I] M. Hammersley and P. Clifford, "Markov field on finite graphs and lattices," Unpublished, 1971.

[42] A. Levada, N. Mascarenhas, and A. Tannus, "Pseudolikelihood equations for potts mrf model parameter estimation on higher order neighborhood systems," IEEE Geoscience and Remote Sensing Letters, vol. 5, no. 3, pp. 522-526, July 2008.

[43] J. Chappelow, N. Bloch, N. Rofsky, E. Genega, R. Lenkinski, W. DeWoIf, S. Viswanath, and A. Madabhushi, "Collinarus: Collection of image-derived non-linear attributes for registration using splines," SPIE Medical Imaging, vol. 7260, 2009.

[44] S. Viswanath, B. Bloch, M. Rosen, J. Chappelow, R. Toth, R. Lenkinski N. Rofsky, E. Genega, A. Kalyanpur, and A. Madabhushi, "Integrating structural and functional imaging for computer assisted detection of prostate cancer on multi-protocol in vivo 3 tesla mri," SPIE Medical Imaging, vol. 7260, 2009.

[45] S. Doyle, M. Hwang, S. Naik, A. Madabhushi, M. Feldman, and T. Tomaszewski, "Using manifold learning for content-based image retrieval of prostate histopathology," in MICCAI 2007 Workshop: Content-based Image Retrieval for Biomedical Image Archives: achievements, problems and prospects, 2007. [46] M. J. Donovan, S. Hamann, M. Clayton, F. M Khan, M. Sapir, V. Bayer-Zubek, U. Fernandez, R. Mesa-Tejada, M. Teverovskiy, V. Reuter, P. Scardino, and C. Cordon-Cardo, "Systems pathology approach for the prediction of prostate cancer progression after radical prostatectomy.," J CHn Oncol, vol. 26, no. 24, pp. 3923-3929, Aug 2008.

[47] C. Cordon-Cardo, A. Kotsianti, D. Verbel, M. Teverovskiy, P. Capodieci, S. Hamann, Y. Jeffers, M. Clayton, F. Elkhettabl, F. Khan, M. Sapir, V. BayerZubek, Y. Vengrenyuk, S. Fogarsi, O. Saidi, V. Reuter, H. Scher, M. Kattan, F. Bianco, T. Wheeler, G. Ayala, P. Scardino, and M. Donovan, "Improved prediction of prostate cancer recurrence through systems pathology.," J CHn Invest, vol. 1 17, no. 7, pp. 1876-1883, JuI 2007.

[48] G. Begelman, E. Gur, E. Rivlin, M. Rudzsky, and Z. Zalevsky, "Cell nuclei segmentation using fuzzy logic engine," in ICIP, 2004, pp. V: 2937-2940. [23] U. Adiga, R. Malladi, R. Fernandez-Gonzalez, and CO. de Solorzano, "High-throughput analysis of multispectral images of breast cancer tissue," Image Processing, IEEE Transactions on, vol. 15, no. 8, pp. 2259-2268, Aug. 2006.

[49] S. Doyle, A. Madabhushi, M. Feldman, and J. Tomaszeweski, "A boosting cascade for automated detection of prostate cancer from digitized histology," MlCCAI, pp. 504-51 1, 2006.

[50] S. Doyle, M. Hwang, K. Shah, A. Madabhushi, M. Feldman, and J. Tomaszeweski, '-utomated grading of prostate cancer using architectural and textural image features," in ISBI, 12-15 April 2007, pp. 1284-1287.

[5 I] S Naik, S Doyle, A Madabhushi, J Tomaszewski, and M Feldman, "Gland segmentation and gleason grading of prostate histology by integrating low-, high-level and domain specific information," Workshop on Microscopic Image Analysis with Applications in Biology, 2007.

[52] S. Naik, S. Doyle, S. Agner, A. Madabhushi, M. Feldman, and J.h Tomaszewski, automated gland and nuclei segmentation for grading of prostate and breast cancer histopathology," ISBI, pp. 284-287, May 2008.

[53] M. Gao, P. Bridgman, and. Kumar, "Computer-aided prostate cancer diagnosis using image enhancement andjpeg2000," 2003, vol. 5203, pp. 323-334, SPIE.

[54] A. Hafiane, F. Bunyak, and K. Palaniappan, "Fuzzy clustering and active contours for histopathology image segmentation and nuclei detection," in Advanced Concepts for Intelligent Vision Systems, 2008, pp. 903-yi4.

[55] V. Kumar, A. Abbas, and N. Fausto, Robbins and Cotran Pathologic Basis of Disease, Saunders, 2004. [3 I] R. Duda, P. Hart, and D. Stork, Pattern Classification, John Wiley & Sons, 2001.

[56] S. Banerjee, B. Carlin, and A. Gelfand, Hierarchical Modeling and Analysis for Spatial Data, CRC Press, 2004.

[57] G. Grimmett, "A theorem about random fields," Bulletin of the London Mathematical Society, vol. 5, no. 13, pp. 81-84, 1973. [34] J. Moussouris, "Gibbs and markov random systems with constraints," Journal of Statistical Physics, vol. 10, pp. 1 1-33, Jan. 1974.

[58] D. Geman, "Random fields and inverse problems in imaging," in Lecture Notes in Mathematics. 1991, vol. 1427, pp. 1 13-193, Springer- Verlag. [36] R. B. Potts, "Some generalised order-disorder transformations," Proceedings of the Cambridge Philosophical Society, vol. 48, pp. 106-109, 1952.

[59] R. Paget, "Strong markov random field model," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol 26, no. 3, pp. 408-413, March 2004.

[60] Y. Bishop, S. Fienberg, and P. Holland, Discrete Multivariate Analysis: Theory and Practice, MIT Press, 1975. [39] A. K. Jain, Fundamentals of Digital Image Processing, Prentice Hall, 1989.

[61] S. A. Hojjatoleslami and J. Kittler, "Region growing: a new approach," IEEE Trans, on Image Processing, vol. 7, no. 7, pp. 1079-1084, July 1998. [41] A. Papoulis, Probability, Random Variables, and Stochastic Processes, McGraw-Hili, Inc., New York, NY, 1965.

[62] S. Geman and C. Graffigne, "Markov random field image models and their applications to computer vision," in Proceedings of the International Congress of Mathematicians, 1986, pp. 1496- 1517.

[63] R. Szeliski, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov, A. Agarwala, M. Tappen, and C. Rother, .- comparative study of energy minimization methods for markov random fields," in ECCV, 2006, pp. 16-29.

[64] R.C. Dubes, A.K. Jain, S.G. Nadabar, and CC. Chen, "Mrfmodel-based algorithms for image segmentation," in Proceedings oftheIOth International Conference on Pattern Recognition, 1990, vol. i, pp. 808-814 voU.

[65] R. Maitra, S. R. Roys, and R. P. Gullapalli, "Test-retest reliability estimation of functional mri data," Magnetic Resonance in Medicine, vol. 48, no. 1, pp. 62-70, 2002.

[66] E. Salli, H. J. Aronen, S. Savolainen, A. Korvenoja, and A. Visa, "Contextual clustering for analysis of functional mri data," IEEE Transactions on Medical Imaging, vol. 20, no. 5, pp. 403- 414, May 2001.

[67] R. Baldick, Applied Optimization, Cambridge University Press, 2006.

[68] Begelman, G., Gur, E., Rivlin, E., Rudzsky, M., Zalevsky, Z.: Cell nuclei segmentation using fuzzy logic engine. In: ICIP. (2004) V: 2937-2940

[69] Naik, S., Doyle, S., Agner, S., Madabhushi, A., Feldman, M., Tomaszewski, J.:

Automated gland and nuclei segmentation for grading of prostate and breast cancer histopathology. ISB! (May 2008) 284-287

[70] Gao, M., Bridgman, P., Kumar, S.: Computer-aided prostrate cancer diagnosis using image enhancement and jpeg2000. Volume 5203., SPIE (2003) 323-334

[71] Doyle, S., Madabhushi, A., Feldman, M., Tomaszeweski, J.: A boosting cascade for automated detection of prostate cancer from digitized histology. MICCAI (2006) 504-51 1

[72] Kumar, V., Abbas, A., Fausto, N.: Robbins and Cotran Pathologic Basis of Disease. Saunders (2004)

[73] Hojjatoleslami, S., Kittler, J.: Region growing: a new approach. IEEE Trans, on Image Processing 7(7) (July 1998) 1079-1084

[74] Geman, S., Geman, D.: Stochastic relaxation, gibbs distribution, and the bayesian restoration of images. IEEE Trans, on PAMI 6 (1984) 721-741

[75] Besag, J.: On the statistical analysis of dirty pictures. Journal of the Royal Statistical Society. Series B (Methodological) 48(3) (1986) 259-302

[76] Begelman, G., Gur, E., Rivlin, E., Rudzsky, M., Zalevsky, Z.: Cell nuclei segmentation using fuzzy logic engine. In: ICIP. (2004) V: 2937-2940

[77] Naik, S., Doyle, S., Agner, S., Madabhushi, A., Feldman, M., Tomaszewski, J.: Automated gland and nuclei segmentation for grading of prostate and breast cancer histopathology. ISBI (May 2008) 284-287

[78] Gao, M., Bridgman, P., Kumar, S.: Computer-aided prostrate cancer diagnosis using image enhancement and jpeg2000. Volume 5203., SPIE (2003) 323-334

[79] Doyle, S., Madabhushi, A., Feldman, M., Tomaszeweski, J.: A boosting cascade for automated detection of prostate cancer from digitized histology. MICCAI (2006)

504-51 1

[80] Kumar, V., Abbas, A., Fausto, N.: Robbins and Cotran Pathologic Basis of Disease.

Saunders (2004)

[81] Hojjatoleslami, S., Kittler, J.: Region growing: a new approach. IEEE Trans, on

Image Processing 7(7) (July 1998) 1079-1084

[82] Geman, S., Geman, D.: Stochastic relaxation, gibbs distribution, and the bayesian restoration of images. IEEE Trans, on PAMI 6 (1984) 721-741

[83] Besag, J.: On the statistical analysis of dirty pictures. Journal of the Royal Statistical

Society. Series B (Methodological) 48(3) (1986) 259-302

[84] M. A. Roula, A. Bouridane, F. Kurugollu, and A. Amira, "A quadratic classifier based on multispectral texture features for prostate cancer diagnosis," in Signal Processing and Its

Applications, 2003. Proceedings. Seventh International Symposium on, July 2003, vol. 2, pp. 37-

40 vol.2.

[85] Reza Farjam, Hamid Soltanian-Zadeh, Kourosh Jafari-Khouzani, and Reza A Zoroofi, "An image analysis approach for automatic malignancy determination of prostate pathological images.," Cytometry B Clin Cytom, vol. 72, no. 4, pp. 227-240, JuI 2007.

[86] A. Tabesh, M. Teverovskiy, Ho- Yuen Pang, V. P. Kumar, D. Verbel, A. Kotsianti, and O. Saidi, "Multifeature prostate cancer diagnosis and gleason grading of histological images," /EEE Transactions on Medical Imaging, vol. 26, no. 10, pp. 1366-1378, Oct. 2007.

[87] Po-Whei Huang and Cheng-Hsiung Lee, "Automatic classification for pathological prostate images based on fractal analysis.," /EEE Trans Med Imaging, vol. 28, no. 7, pp. 1037-1050, JuI 2009.

[88] S. F. Patten, J. S. Lee, D. C. Wilbur, T. A. Bonfiglio, T. J. Colgan, R. M. Richart, H. Cramer, and S. Moinuddin, "The autopap 300 qc system multicenter clinical trials for use in quality control rescreening of cervical smears: Ii. prospective and archival sensitivity studies.,"

Cancer, vol. 81, no. 6, pp. 343-347, Dec 1997.

[89] Muhammad Atif Tahir and Ahmed Bouridane, "Novel round-robin tabu search algorithm for prostate cancer classification and diagnosis using multispectral imagery," /EEE Trans Inf

Technol Biomed, vol. 10, no. 4, pp. 782-793, Oct 2006.

[90] K. Jafari-Khouzani and H. Soltanian-Zadeh, "Multiwavelet grading of pathological images of prostate," /EEE Transactions on Biomedical Engineering, vol. 50, no. 6, pp. 697-704, June

2003.

[91] S. Naik, S. Doyle, S. Agner, A. Madabhushi, M. Feldman, and J. Tomaszewski, "Automated gland and nuclei segmentation for grading of prostate and breast cancer histopathology," in Proc.

5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro ISBI 2008,

14-17 May 2008, pp. 284-287.

[92] Yoav Smith, Gershom Zajicek, Michael Werman, Galina Pizov, and Yoav Sherman,

"Similarity measurement method for the classification of architecturally differentiated images,"

Comput. Biomed. Res., vol. 32, no. 1, pp. 1-12, 1999.

[93] G. Begelman, M. Pechuk, Ε. Rivlin, and Ε. Sabo, "System for computer-aided multiresolution microscopic pathology diagnostics," in Proc. IEEE International Conference on

Computer Vision Systems ICVS '06, 04-07 Jan. 2006, pp. 16-16.

[94] R. Stotzka, R. Mnner, P. H. Bartels, and D. Thompson, "A hybrid neural and statistical classifier system for histopathologic grading of prostatic lesions.," Anal Quant Cytol Histol, vol.

17, no. 3, pp. 204-218, Jun 1995.

[95] James Diamond, Neil H Anderson, Peter H Bartels, Rodolfo Montironi, and Peter W

Hamilton, "The use of morphological characteristics and texture analysis in the identification of tissue composition in prostatic neoplasia.," Hum Pathol, vol. 35, no. 9, pp. 1 121-1 131, Sep

2004.