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.
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 for WMHSs is also significant for the following reasons: 1) 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.
With respect to prostate histology, Begelman [49] presented a method for nuclei segmentation on hematoxylin and eosin (H&E) 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.
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 40× magnification) contains approximately 200,000×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×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.
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.
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.
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. Unfortunately, constructing LCPDFs that 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.
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.
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 of potential functions, tailoring LCPDFs to model a specific process devolves 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.
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).
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.
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.
The present invention relates to a method and system for detecting biologicaily relevant structures in a hierarchical fashion, beginning at a low-resolution and proceeding to higher levels of resolution.
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.
It is an object of the present invention to provide probabilistic pairwise Markov models (PPMMs) to classify biologically relevant structures.
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.
It is an object of the present invention to provide an automated system for detecting CaP glands in WMHSs.
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.
It is an object of the present invention to provide a method for incorporating asymmetric costs into MRF estimation.
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.
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.
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 so the cancerous glands are consolidated into regions.
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.
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.
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.
In accordance with the above object and others, the invention is directed to a novel classification approach which weighs the importance of these structures.
a) depicts H&E stained prostate histology section; black ink mark provided by pathologist indicates CaP extent.
b) depicts gland segmentation boundaries.
c) depicts a magnified view of the white box depicted in
d) depicts centroids of cancerous glands before MRF iteration.
e) depicts centroids of cancerous glands after MRF iteration.
a) depicts seed pixel with CB.
b) depicts IB, CB, and CR at step j in the growing procedure.
c) depicts IB, CB, and CR at step j+1.
a)-(d) depict H&E stained WMHSs.
e)-(h) depict initial labels for area=0.25.
i)-(l) depict PPMM GDS results with sensitivity adjusted to 72%.
m)-(p) depict Potts GDS results with sensitivity adjusted to 72%.
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.
b) depicts area under corresponding ROC curves in
a) depicts T2-w MR image with cancerous region as delineated by a radiologist.
b) depicts intensities indicate the probability of cancer for each pixel assuming uniform priors.
c) depicts initial classification after MLE.
d)-(f) depict final classification after WICM for
g)-(i) depict final classification after WMPMMC for
a) depicts luminance channels of two different HER2+ BC histopathology studies.
b) depicts corresponding results for initial region-growing based lymphocyte detection.
c) depicts preliminary Bayesian refinement showing detected BC nuclei and detected lymphocyte nuclei.
d) depicts final lymphocyte detection result after the MRF pruning step.
a) depicts HER2+ breast cancer histopathology image.
b) depicts corresponding Voronoi Diagram constructed using the automatically detected lymphocyte centers as vertices of the graph.
c) depicts corresponding Delaunay Triangulation constructed using the automatically detected lymphocyte centers as vertices of the graph.
d) depicts corresponding Minimum Spanning Tree constructed using the automatically detected lymphocyte centers as vertices of the graph.
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.
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 of modalities (including, but not limited to, ultrasound, radiography, PET, CT).
An embodiment of the present invention is directed to a system specifically designed to operate at low-resolution (0.008 mm 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.
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.
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 MRF 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.
In a preferred embodiment, the present invention incorporates analysis or gland 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.
The invention is described more fully herein. All references cited are hereby incorporated by reference in their entirety herein.
Let the set S={1, 2, . . . , N} reference N sites to be classified. Each site s ε S has two associated random variables: Xs ε Λ indicating its state (class) and Ys εD representing its D-dimensional feature vector. Particular instances of Xs and Ys are denoted by the lowercase variables xs ε Λ and ys εD. Let X=(X1, X2, . . . , XN) and Y=(Y1, Y2, . . . , YN) refer to all random variables Xs and Ys in aggregate. The state spaces of X and Y are the Cartesian products Ω=ΛN and D×N. Instances of X and Y are denoted by the lowercase variables x=(x1, x2, . . . , xN)ε Ω and y=(y1, y2, . . . , yN) εD×N.
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 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} ε E}. If P is a probability measure defined over Ω then the triplet (G, Ω, P) is called a random field.
These concepts are best understood in the context of an example. The graph in FIG. 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 η5={1, 4}. There are six one-element cliques ={{1}, {2}, {3}, {4}, {5}, {6}}, six two-element cliques 2=E, and one three-element clique 3={{1, 4, 5}}. The set is the union of 1, 2, and 3. The set of possible states for each Xs is Λ={b, w}, where b and w represent black and white, respectively. The random variable Ys ε reflects the diameter of site s. For this example we have x=(w, b, b, b, w, b) and y=(0.9, 1.5, 1.6, 1.4, 1.1, 1.3).
Finally, we establish a convention for representing probabilities. Let P(·) indicate the probability of event {·}. For instance, P(Xs=xs) and P(X=x) signify the probabilities of the events {Xs=xs} and {X=x}. Whenever appropriate we will simplify such notations by omitting the random variable, e.g. P(x)≡P(X=x). Let p(·) represent a probability PDF; for example, pg might indicate a Gaussian PDF. The notations P(·) and p(·) are useful in differentiating P(xs) which indicates the probability that {Xs=xs} from pg(xs) which refers to the probability that a Gaussian random variable assumes the value xs.
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 x ε Ω:
The first term in (2.2.1) reflects the influence of the feature vectors. It can be simplified by assuming that all Ys are conditionally independent and identically distributed given their associated Xs. This assumption implies that if the class Xs of site s is known then 1) the classes and features of the remaining sites provide no additional information when estimating Ys and 2) the conditional distribution of Ys is identical for all s ε S. As a result we have
where the use of the single PDF pf indicates that P(ys|xs) is identically distributed across S. 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.
The random field (G, Ω, P) is a Markov Random Field (MRF) if its local conditional probability density functions satisfy the Markov property:
P(xs|x−s)=P(xs|xη
where x−s=(x1, . . . , xs−1, xs+1, . . . , xN), xη
Given the LCPDFs for all s ε 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.
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(x)>0 for all x ε Q satisfies the Markov property if and only if it can be expressed as a Gibbs distribution:
where Z=ΣxεΩ Vc (x) is the normalizing constant and Vc are positive functions, called clique potentials, that depend only on those xs such that s ε c. The following derivation reveals the simplified forms of the LCPDFs:
where represents {c ε : s ε c} and −s represents {c ε : s ∉ c}. This derivation also proves that all Gibbs systems are Markovian.
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.
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.
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. P(xs|xη
where V1 and V2 are the site-invariant potential functions for one- and two-element cliques, respectively. Furthermore, V2 is symmetric in the sense that V2(xs, xr)=V2(xr, xs); this symmetry is needed to ensure the stationarity of P(xs|xη
Before continuing we introduce a new notation. Let p0|m(xs|xη
Consider the following manipulation of p0|m:
In essence, (3.1.2) provides a means for representing p0|m in terms of the single pairwise distribution p0,1, where P1|0 and p0 are a conditional and marginal distribution of p0,1. 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.
Let (G, Ω, P) be a random field such that P(X=x)>0 for all x ε Ω. If p0,1(xs, xr) is a PDF such that p0,1(xs, xr)=p0,1(xr, xs) and
for all s ε S, then X is a MRF.
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:
V
1(xs)=p0(xs)1−|η
and
V
2(xs, xr)=p0,1(xs, xr).
where the symmetry of p0,1 ensures that v2 is also symmetric. Inserting these into (3.1.1) yields
Consider the case of a continuous MRF where Xs ε Λ≡. Let p0,1 be a (necessarily symmetric) Gaussian distribution:
where σ ε + and γ ε[−1,1]. Equation (3.3.1) reflects a bivariate Gaussian that has been rotated such that its principal axes coincide with the lines xs=xr and xs=−xr, i.e. a rotation of π/4. Its requisite marginal and conditional distributions are
Inserting these results into (3.1.3) yields the following:
where ξ=1+γ2(|ηs|−1). The LCPDFs are themselves Gaussian distributions. Thus, by designating p0,1 as a Gaussian distribution, we have created an instance of a Gaussian MRF.
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:
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.
The Potts model can be reformulated as a PPMM as follows:
Note that this probability function is symmetric.
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: ∀x ε Ω, s ε S
P(Xs=xs|Xr=xr: r≠s, r ε AS)=P(Xs=xs|Xr=xr: r ε ηs ∩ A).
This states that a MRF with sites S is a strong-MRF if every A 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).
To understand the connection between PPMMs and strong-MRFs again consider a stationary MRF X with LCPDFs given by p0|m, i.e. P(xs|x−s)=p0|m(xs|xη
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: p0,1(xs, xr)=P(xs, xr). However, this equality does hold if (and only if) X is a strong-MRF. Consequently, only when X is an strong-MRF can p0,1 (x,s, xr) be replaced with P(xs, xr). The overarching principle is that two different joint distributions can share an identical conditional distribution.
Others aspects of strong-MRFs also provide insight into PPMMs. Paget's decomposition of strong-MRFs into pairwise clique probabilities agrees with our representation when X is a strong-MRF. However, our formulation remains valid when X is not a strong-MRF. 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.
The present invention provides a means for incorporating asymmetric costs into MRF 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 MRF optimization routine.
Given an observation of the feature vectors Y, we would like to estimate the states X. Bayesian classification advocates selecting the estimate {circumflex over (x)} ε Ω that minimizes the conditional risk (expected cost)
R(X|{circumflex over (x)},y)=E[C(X,{circumflex over (x)})|y]=ΣxεΩC(x,{circumflex over (x)})P(x|y), (5.2.1)
where E indicates expected value and C(x,{circumflex over (x)}) is the cost of selecting labels {circumflex over (x)} when the true labels are x. In the following subsections we will consider the two most prevalent cost functions used with MRFs.
The most ubiquitous cost function (though this cost is rarely expressed explicitly) is
C
MAP(x,{circumflex over (x)})=ΠsεS[1−δ(xs−{circumflex over (x)})], (5.3.1)
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
Minimizing (5.3.2) over {circumflex over (x)} is equivalent to maximizing the second term. Thus, we have maximum a posteriori (MAP) estimation, which advocates selecting the {circumflex over (x)} that maximizes the a posteriori probability.
When the cardinality of Ω is small the search for the optimal {circumflex over (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(xs|x−s,y)P(x−s|y) (5.3.3)
Increasing the first term in (5.3.3) necessarily increases P(x|y). This suggests a global optimization strategy that sequentially visits each site s ε S and determines the label xs ε Λ that maximizes P(xs|x−s,y). ICM converges to a local, and not global, maximum of P(x|y); however, in many situations the local maximum is more appropriate. The ICM algorithm is as follows:
Typically, the initial labeling x0 is provided by the maximum likelihood estimate of P(y|x). This probability becomes tractable under the typical assumption that all Ys are conditionally independent given their corresponding Xs, i.e. P(y|x)=ΠsεSP(ys|xs). The ICM iteration usually converges in approximately six or seven iterations.
As an alternative to MAP estimation, Marroquin suggested using the following cost function
C
MPM(x,{circumflex over (x)})=ΣsεS[1−δ(xs−{circumflex over (x)}s)]. (5.4.1)
This function counts the number of sites in {circumflex over (x)} that are mislabeled. Inserting (5.4.1) into (5.2.1) yields
The distributions P(xs|y) are called the posterior marginals. Minimizing (5.4.2) is equivalent to independently maximizing each of these posterior marginals with respect to its corresponding {circumflex over (x)}s. Hence, this type of estimation is termed maximum posterior marginals (MPM).
Since the cardinality of Ω is usually large, P({circumflex over (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:
In theory, the final labels xk, which represent a sample from the distribution P(x|y), are independent of the starting conditions x0. However, since this algorithm must eventually terminate, some dependence upon initial conditions does exist. As with ICM, x0 is typically initialized to the MLE. Determining the number of iterations 1 needed for the Markov chain to reach equilibrium will depend upon the distribution P(x|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 maximizes P(xsk|x−sk,ys), the Gibbs sampler selects the label randomly from this distribution.
After reaching equilibrium the amount of time the chain spends in any state x is given by P(x|y). Thus, the posterior marginal P(xs|y) can be estimated as follows:
where m−l is a number of iterations past equilibrium needed to generate an accurate estimate. The value for m, like l, is typically chosen empirically. Given this estimate of P(xs|y), the maximizer can be determined by testing all possible states of xs.
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→0+ define a weighting function. Let {tilde over (X)} be a random field that differs with X only with respect to its probability measure, which is defined as follows:
where the constant Z ensures that ΣxεΩP({tilde over (X)}=x)=1. For convenience we will use following simplified notation: {tilde over (P)}(x)≡P({tilde over (X)}=x). Finally, note that (6.1) implies
To allow certain classification decisions to be weighted more heavily than others CMAP can be generalized as follows:
C
WMAP(x,{circumflex over (x)})=ΠsεSα(xs)[1−δ(xs−{circumflex over (x)}s)]. (6.3)
Consequently, if the estimate {circumflex over (x)} contains any erroneous labels it accrues a cost of ΠsεSα(xs). Inserting (6.3) into (5.2.1) yields
This result indicates that minimizing RWMAP(X|{circumflex over (x)},y) is equivalent to minimizing RMAP({tilde over (X)}|{circumflex over (x)},y); and consequently, the optimal labeling is the MAP estimate of {tilde over (X)}. Since the a posteriori probability of {tilde over (X)} corresponds to the weighted a posteriori probability of X, we refer to this type of estimation as weighted MAP (WMAP) estimation. Note that if α(xs)≡1, WMAP estimation reduces to MAP estimation.
Implementing the WMAP estimation of X is trivial: we simply perform MAP estimation on {tilde over (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:
Modifying ICM to perform WMAP estimation only requires replacing P(xs|x−s,ys) with {tilde over (P)}(xs|x−s,ys). Recognizing that {tilde over (P)}(xs|x−s,ys)∝α(xs)P(xs|x−s,ys) reveals the final form used in step 6.
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).
The cost function CMPM can be generalized as follows:
C
WMPM(x,{circumflex over (x)})=ΠrεSα(xr)ΣsεS[1−δ(xs−{circumflex over (x)}s)]. (7.1)
In this case mislabeling a site whose true label is xs has an associated cost of ΠsεSα(xs). 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,{circumflex over (x)})=ΣsεS[1−δ(xs−{circumflex over (x)}s)], developing a scheme to perform the necessary estimation is problematic.
Inserting (7.1) into (5.2.1) yields
Thus, the WMPM estimate of X is equivalent to the MPM estimate of {tilde over (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:
The only difference between the Gibbs sampler and the weighted Gibbs sampler is the replacement in step 6 of P(xs|x−sk,ys) with
{tilde over (P)}(xs|x−sk,ys)=(α(ω)/Zs)P(xs|x−sk,ys), (7.3)
where Zs ensures ΣωεΛ{tilde over (P)}(ω|x−sk,ys)=1. The remainder of the WMPMMC procedure is identical to that of the MPMMC algorithm.
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.
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.
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 glands. This information is incorporated via a Markov random held (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.
We first provide a basic outline of the algorithm and then discuss each component in detail.
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.
The basic algorithm, illustrated in
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 Ys ε indicates the square root of the area of gland s. The label Xs of gland s is one of two possible classes: Xs ε Λ ≡{ω1, ω2}, where ω1 and ω2 indicate malignancy and benignity. Two glands are neighbors if the distance between their centroids is less than 0.7 mm.
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):
Step 1: Initialize the current region (CR) to the specified seed pixel and establish a 12σg×12σg bounding box centered about it. Initialize the current boundary (CB) to the 8-connected pixels neighboring CR (
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 (
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.
The final segmented regions may overlap; this is resolved by discarding the region with the lower boundary strength.
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 ω1 and benign ω2 glands. Using the equivalent square root of gland area (SRGA), the pdfs pf(ys|ω1) and pf(ys|ω2), which correspond to the PDF pf in (2.2.2), can be accurately modeled with a weighted sum of Gamma distributions:
where y>0 is the SRGA, αε[0, 1] is the mixing parameter, k1, k2>0 are the shape parameters, θ1, θ2>0 are the scale parameters, and Γ 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 ε[0, 1] are labeled malignant; the remainder are classified as benign. The centroids of glands labeled as malignant are shown with blue dots in
PPMMs require the probability p0,1 or equivalently, p1|0 and p0. For the binary classes the general forms of these probabilities are
The required symmetry of p0,1 necessitates that c(1−a)=(1−c)b, yielding c=b/(1+b−a). In the case of the Potts model we have a=eβ/(1+eβ), b=1/(1+eβ), and c=1/2.
Values for a and 6 (or β for the Potts model) can obtained from training data. Though (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.
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.
ICM is based on the following reformulation of (2.2.1):
P(x|y)=P(xs|x−s,y)P(x−s|y)∝P(xs|xη
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 ε S and determines the label xs ε Λ that maximizes the first term of (9.6.3). In the case of binary classes this reduces to following decision:
where Twicm=1/2. Note that P(ω1|xη
An examination of (9.6.4) immediately suggests varying Twicm as a means for favoring one class over the other. As Twicm decreases, (9.6.4) will increasingly prefer ω1. By contrast, increasing Twicm, will favor ω2. Since the value of Twicm implicitly weights the importance of each class, we refer to the modified algorithm as weighted iterated conditional modes (WICM).
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.
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
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 pf. The graph structure connecting the glands is determined from the gland 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 (a and b) and Potts model (β).
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 τarea=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(xs|xη
It is worth mentioning that classification performance could conceivably be adjusted by fixing Twicm=1/2 (i.e. using ICM) and varying τarea. That is, since modifying τarea 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.
The training and test sets were established using a leave-one-out strategy. The performances of the systems were varied by adjusting Twicm (for each algorithm independently) until their respective sensitivities approximately equaled 0.72. That is, when Twicm=0.6 the PPMM GDS detected 72 percent of the malignant glands; to achieve the same sensitivity, the Potts GDS required Twicm=0.6.
The training and testing data—which were identical for both systems—consisted of 20 independent trials using randomized 3-fold cross-validation. By varying Twicm from zero to one, ROC curves were generated for each trial.
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.
The intensities in
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).
An embodiment of the present invention provides a computer-aided diagnosis (CADx) scheme (
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 (
The initial lymphocyte detection is refined by 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 (
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 (
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 (
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%.
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).
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.
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:
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
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
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 ωb 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:
where y>0 is the SRGA, λε[0, 1] is the mixing parameter, k1,k2>0 are the shape parameters, θ1,θ2>0 are the scale parameters, and F is the Gamma function. Note, we use f 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. Those glands whose probabilities exceed the predetermined threshold ρ are labeled malignant; the remainder are classified as benign (
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.
Let S={s1, s2, . . . , sN} represent a set of N unique sites corresponding to the N segmented glands. Let each site sεS have an associated random variable Xsε{ωm,ωb} indicating its state as either malignant or benign. To refer collectively to the states of all glands we have X={Xs:sεS}. Each state Xs is unknown; we only observe an instance of the random variable Ys ε RD representing the D dimensional feature vector associated with gland s. Though our algorithm is extensible to any number of features, currently D=1 with Ys being the SRGA. To collectively refer to the entire scene of feature vectors we have Y={Ys:sεS}.
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:rεS,r≠s,{r,s}εE}. The set of all local neighborhoods establishes a neighborhood structure: η={ηs:sε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
To simplify notation we use Pr{Xr=xr,Xs=xs}≡p(xr,xs) for indicating the probability of a specific event, where xr,xsε{ωm,ωb}. If X is a MRF with respect to the neighborhood structure η, then X satisfies the local Markov property p(xs|x−s)=p(xs|xη
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.
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:
The Potts model disregards all cliques having more or less than two elements, i.e. if |c|≠2 we have Vc(x)=1, where |·| signifies cardinality. Such generic models are unnecessary; assuming that all Xs are i.i.d. and all Xr given Xs are i.i.d. for every rεηs, we can determine the appropriate Vc directly from the data. To our knowledge, the following equations represent the first proposed means for incorporating arbitrary pmfs into the MRF structure:
V
s(xs)=p(xs)1−|η
V
{r,s}(x)=p(xr,xs) for {r,s}ε C. (14)
The functions Vc for higher-order cliques are identically one. The validity of (13) and (14) can be seen by inserting them into (11):
The determination of p(xs) and p(xr,xs) 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,s}(ω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(xs) is the marginal mass function of p(xrx,xs). Since p(xr,xs) is symmetric, both marginals are identical.
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(x|y)∝f(y|x)p(x), where ∝ 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(xs|x72
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 (
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 mm2 per pixel. The approximate image dimensions are 2.1×3.2 cm, i.e. 2100×3200 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.
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 ρ increases from zero to one, yielding the receiver operator characteristic (ROC) curve (solid) in
Qualitative Results: Qualitative results for two WMHSs were previously presented as
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 2100×3200 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.
This application claims priority to U.S. Provisional Patent Application No. 61/093,884 filed Sep. 3, 2008, the disclosure of which is hereby incorporated by reference in its entirety.
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/US09/04979 | 9/3/2009 | WO | 00 | 6/20/2011 |
Number | Date | Country | |
---|---|---|---|
61093884 | Sep 2008 | US |