Segmentation of Stochastic Images With a Stochastic ... - IEEE Xplore

2424

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 5, MAY 2012

Segmentation of Stochastic Images With a Stochastic Random Walker Method Torben Pätz and Tobias Preusser

Abstract—We present an extension of the random walker segmentation to images with uncertain gray values. Such gray-value uncertainty may result from noise or other imaging artifacts or more general from measurement errors in the image acquisition process. The purpose is to quantify the influence of the gray-value uncertainty onto the result when using random walker segmentation. In random walker segmentation, a weighted graph is built from the image, where the edge weights depend on the image gradient between the pixels. For given seed regions, the probability is evaluated for a random walk on this graph starting at a pixel to end in one of the seed regions. Here, we extend this method to images with uncertain gray values. To this end, we consider the pixel values to be random variables (RVs), thus introducing the notion of stochastic images. We end up with stochastic weights for the graph in random walker segmentation and a stochastic partial differential equation (PDE) that has to be solved. We discretize the RVs and the stochastic PDE by the method of generalized polynomial chaos, combining the recent developments in numerical methods for the discretization of stochastic PDEs and an interactive segmentation algorithm. The resulting algorithm allows for the detection of regions where the segmentation result is highly influenced by the uncertain pixel values. Thus, it gives a reliability estimate for the resulting segmentation, and it furthermore allows determining the probability density function of the segmented object volume. Index Terms—Image segmentation, partial differential equations (PDEs), random variables (RVs), stochastic processes, uncertainty.

I. INTRODUCTION MAGE segmentation in real-world applications is typically performed on noisy images (see Fig. 1). The image noise depends on the image acquisition modality [e.g., digital camera, magnetic resonance, computed tomography (CT), and ultrasound], the acquisition parameters (acquisition time, sound frequency, and magnetic-field strength), and extrinsic parameters (illumination and reflection). Image acquisition itself, in general, results from measurements of physical quantities (photon density, time of flight of waves, spin, and absorption) and, according to good scientific practice, must be thus equipped with information about the measurement error. However, the step of quantifying this error is typically omitted

I

Manuscript received August 22, 2011; revised January 06, 2012; accepted January 24, 2012. Date of publication February 10, 2012; date of current version April 18, 2012. This work was supported by the Deutsche Forschungsgemeinschaft under Grant PR 1038/5-1. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Xin Li. The authors are with the School of Engineering and Science, Jacobs University Bremen, 28759 Bremen, Germany, and also with Fraunhofer Institute for Medical Image Computing MEVIS, 28359 Bremen, Germany (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TIP.2012.2187531

Fig. 1. Noisy images acquired with an (left) ultrasound device and (right) CT.

in image processing, leading to a loss of information about the influence of the input error to the result of the image processing steps. The consideration of error propagation is important in, for example, medical applications, where decisions about patients’ treatments are based on quantitative information extracted from radiological images. Treatment decisions for chemotherapy patients are based on the progression of the volume of the segmented lesions from the noisy images. Thus, it is important to equip the extracted data with a reliability estimate or (this is the aim of the presented paper) to be able to compute the probability density function (PDF) of the extracted data dependent on an estimate or modeling of the input noise. A possibility to get information about the error propagation is to generate input samples via the repeated acquisition of the same scene or via sampling from noise modeling. Then, we apply a segmentation method, e.g., random walker segmentation, on every sample and compute stochastic quantities such as mean and variance from the sample results. This is known as the Monte Carlo method [1]. The approach is time consuming due to the high number of samples needed to get accurate results. The main idea of this paper is to avoid the repeated use of classical segmentation methods by extending them to stochastic images, i.e., images whose pixels are random variables (RVs) describing the noisy acquisition process. Then, we have to segment the stochastic images once to get the stochastic information. Dependent on the numerical methods used for the segmentation of the stochastic image, it can be done as fast as computing the solution of a few Monte Carlo samples. The segmentation of images, i.e., the detection of the shape of objects in images, is a widely investigated topic (see, e.g., [2] for an overview). Many segmentation methods are based on partial differential equations (PDEs), e.g., level-set-based segmentation [3], Mumford–Shah segmentation [4], and the related Ambrosio–Tortorelli approximation [5] or random walker segmentation [6]. Furthermore, there are much more segmentation

1057-7149/$31.00 © 2012 IEEE

PÄTZ AND PREUSSER: SEGMENTATION OF STOCHASTIC IMAGE WITH STOCHASTIC RANDOM WALKER METHOD

methods available in the literature based on energy minimization [2] and Markov random fields (MRFs) [7]. In particular, the MRFs should not be confounded with the stochastic approach proposed here. We focus on error propagation, i.e., we propagate the uncertainty information from the input to the output. In contrast, the MRFs try to avoid erroneous segmentation results but give no information on how the input uncertainty influences the result. Random walker segmentation is, in contrast to other segmentation methods based on PDEs, a supervised segmentation method, meaning that the user input influences the segmentation. The user input consists of defining seeds, i.e., regions where the user specifies whether they belong to the object to be segmented or not. From the unseeded pixels, random walks start, and pixels get a probability to belong to the object. The probability thereby depends on the fraction of random walks reaching the seed region of the object. It is possible to become independent of the user input by incorporating priors [8]. We emphasize that, in the original random walker algorithm, there is no information about included uncertainties, although its name might suggest so. Instead, the name “random walker” is chosen because the direction of the random walk is randomly chosen from probabilities based on the deterministic image gradient. In fact, error propagation is hard or even unfeasible for most of the image processing algorithms including the random walker method. Only few publications are dealing with the computation of error estimates, which come with additional simplifying assumptions. For example, Weber and Malik [9] presumed that the error is Gaussian distributed, and Nestares et al. [10] computed bounds for the error only. Grady [6] made a first attempt to include error propagation in random walker segmentation but was able to prove theorems dealing with independent RVs only. He did not provide an algorithm to compute the stochastic result numerically. Based on stochastic PDEs (SPDEs), Preusser et al. [11] developed a method for error propagation in image processing that was extended in [12]. They identified pixels by RVs and called such images stochastic images. The discretization of stochastic images utilizes the generalized polynomial chaos (gPC) expansion developed in [13]. The gPC expansion provides a basis for the space of square integrable RVs (i.e., RVs with finite variance) built by polynomials in RVs with known distribution. This enables us to deal with deterministic coefficients to represent RVs in the gPC, because the basis functions encode the stochasticity. The resulting SPDEs are solved using the generalized spectral decomposition (GSD) [14], yielding to a tractable and fast algorithm for the high-dimensional stochastic problems. In this paper, we use this approach and combine it with the random walker segmentation. To the best of the authors’ knowledge, the combination of interactive segmentation and error propagation is new. The remaining part is organized as follows: In Section II, we review random walker segmentation and highlight the relation to diffusion equations. In Section III, we review stochastic images and methods needed for their discretization and for the discretization of SPDEs. In Section IV, we introduce the stochastic extension of the random walker segmentation. Results obtained with this segmentation method are presented in Section V. Finally, we draw conclusions in Section VI.

2425

2

Fig. 2. Graph generated from a 3 3 image contains nine nodes and 12 edges. Every edge e has a weight denoted by w .

II. REVIEWING RANDOM WALKER SEGMENTATION We consider image defined on the image domain . Because random walker segmentation is based on pixel values only, we need no additional assumptions about the smoothness of the image. For random walker segmentation, we have to define notation for the graph, which represents image . Graph is a pair containing nodes and edges . Since we use undirected graphs, an edge connecting nodes and is denoted by either or . If every edge has weight that describes the costs for using the edge, we have a weighted graph. The sum of all weights of edges connecting to node , i.e., (1) is called the degree of node . A graph corresponding to a 3 3 image is depicted in Fig. 2. For the random walker segmentation, the graph weights are given by (2) is the normalized difference between the image where intensities at nodes and . Thus, is given by

(3) Parameter

is the only parameter that the user chooses.

A. Relation to the Diffusion Equation Doyle and Snell [15] showed that an infinite number of random walks is equivalent to solving a steady-state diffusion equation with prescribed boundary conditions. Moreover, the steady-state diffusion equation is the necessary condition for minimizers of the Dirichlet integral . For the random walker algorithm, we use a weighted discrete version of the Dirichlet integral. Thus, we introduce the combinatorial Laplacian matrix, i.e., if , if and are adjacent nodes, otherwise

(4)

which allows formulating the discrete version of the weighted Dirichlet integral as (5)

2426


Here, is a vector containing all nodes/pixels of the graph , where is the pixel (resp. the image), i.e., set. For the initialization of the random walker algorithm, the for the object and the user prescribes a set of seed points background. These points act as boundary conditions for the Dirichlet problem with a value of 1 on the object seeds and 0 on the background seeds, because the probability that a random walker starting at a seed points reaches the point is 1. The , are the degrees remaining unseeded points, denoted by of freedom. Reordering the nodes accordingly, the discrete Dirichlet integral (5) can be rewritten as

(6) of matrix describes the dependences between the Part , part describes the dependences between the seed points describe the coupling between unseeded pixels, and and the seeded and unseeded pixels, respectively. of this discrete energy is given by Minimizer (7) This system of linear equations can be solved with an iterative solver, e.g., the method of conjugate gradients [16]. Due to the random walks, this segmentation method might be seen as a stochastic method for segmentation in which randomness and uncertainty are included. However, we emphasize that, by using the equivalence to the Dirichlet problem, the method computes deterministic weights for an elliptic PDE on a graph. As the result, every node/pixel gets values that are interpreted as probabilities for reaching the seed region from this specific node/pixel. These values are used to segment regions in the images. In summary, the original random walker algorithm is a method for the segmentation of deterministic images but does not deal with stochastic images or error propagation. The extension to stochastic images having uncertain gray values is presented in the following sections. III. STOCHASTIC IMAGES In PDE-based image processing, image is typof weakly difically an element of the Sobolev space , we use ferentiable functions. For the discretization of . Thus, a fia multilinear finite-element space nite-element representation of is given by (8) where is the value of the th pixel from the pixel and is its shape function (e.g., tent function, see, set e.g., [2]). In stochastic images, a pixel depends on vector of RVs [12], i.e., , and thus on . Here, and in the following, denotes a random event is a -algebra, and is the probability the event space, measure. Due to the dependence of a pixel on a vector of RVs, the pixel is a RV, i.e., the pixel value depends on the random

. In what follows, we use square integrable RVs event only, i.e., we assume that every stochastic pixel is an element . Thus, a stochastic image is an element of the of space . When we fix a random tensor product space , we get a classical image back from the tensor event . Furtherproduct representation because more, a partial evaluation of the tensor product representation yields a RV for every with respect to the spatial position . Note location inside the image domain that the concept of stochastic images can be also combined with other spatial discretizations, e.g., finite-difference schemes. A. Polynomial Chaos Expansion The polynomial chaos expansion, which was introduced in [17] and extended as the gPC in [13], can be used to represent arbitrary RVs with finite second-order moments. Furthermore, Cameron and Martin [18] proved that every second-order RV , i.e., , can be represented by (9) is a sequence of independent RVs with where are polynomials in forming a basis of known PDF , and . In numerical algorithms, one typically uses a gPC apwith a fixed polynomial deproximation in space gree and a fixed number of independent identically distributed , i.e., RVs (10) This finite-dimensional gPC approximation involves mapof events to , where ping is finite dimensional. Thus, the is mapped to a vector of real numbers. stochastic event This makes operations such as integration over the stochastic is transdimensions much easier, because integral . Furthermore, the independence of formed into into the product of the RVs allows decomposing measure , where the measures of the individual RVs, i.e., is the density of the RV . We write , and we omit denoting the dependence of on for the ease of the presentation. The number of polynomial basis functions depends on the number of RVs and the maximal polynomial degree of the approximation. As usual for polynomial bases, the number of . Thus, rapidly grows basis functions is given by with the number of RVs and the polynomial degree of the approximation. To speed up the convergence of the finite-dimensional ap, which are pairwise proximation, one chooses polynomials orthonormal with respect to the corresponding probability measure of , i.e., where if else.

,

(11)

values are products of 1-D normalFor Gaussian RVs , values are norized Hermite polynomials. For uniform RVs, malized products of Legendre polynomials. The polynomials for more sophisticated distributions are derived in [13]. There


it is shown that can be selected from the Askey scheme for polynomials [19]. In this paper, we use uniformly distributed RVs involving Legendre polynomials. Let us emphasize that the choice of uniformly distributed RVs does not mean that the investigations are restricted to uniform RVs. As shown in [18], we can model all second-order RVs as the sum of polynomials in uniform RVs. The expected value, variance, and analogously higher order stochastic moments of the RV are evaluated as

2427

The variance calculation benefits from the precomputation of , because functions are chosen in advance. We obtain other stochastic moments in a way similar to (16). The construction of stochastic images from sample images is based on the Karhunen–Loève decomposition [20], i.e., it is based on the eigenpairs of the covariance matrix of the samples. The uncorrelated RVs obtained from the Karhunen–Loève decomposition are projected on the polynomial chaos using the Galerkin projection. Details about this stochastic input generation can be found in [12] and [21]. C. SFEM

(12)

B. Polynomial Chaos for Stochastic Images Preusser et al. [11] showed that the representation of an image whose pixels are RVs is obtained from (8) by replacing by RVs ; thus (13)

The stochastic finite-element method (SFEM) [22] for the discretization of SPDE is based on the gPC expansion. Due to the tensor product structure of the ansatz space , the ansatz functions are for and . The resulting system of linear equations tohas a block structure. Grouping the equations for fixed gether, the blocks have the same sparsity structure and extent as a classical FEM matrix. The entries of the stochastic mass and stiffness matrix are given by

by a With the gPC expansion (10), we approximate weighted sum of orthonormal multidimensional polynomials, thus leading to

(17)

(14)

i.e., the blocks are (in this simple setting) scaled FEM matrices.

as the representation of a stochastic image, i.e., an image whose a stogray values are RVs. For fixed , we call coefficient collects the stochastic chastic mode of pixel . Set modes of all pixels for fixed and can be visualized as a classical image, because we have a deterministic value for every pixel. Using this idea, a stochastic image can be visualized as (the number of chaos coefficients per pixel) a sequence of are deterministic images. Let us emphasize that coefficients encode the stochastic inforand that the basis functions mation. Thus, we are dealing with deterministic coefficients but obtain stochastic quantities due to the multiplication with the stochastic basis functions. From the gPC expansion of the stochastic image (14), it is straightforward to compute stochastic moments of the images. With the use of the orthonormal set of basis functions—here the Legendre polynomials—we have

D. Generalized Spectral Decomposition The stochastic equation system arising in the discretization of the random walker segmentation on stochastic images is discretized using the GSD [14], which is based on the SFEM. Compared with the SFEM as the classical solution method for SPDE, the GSD allows for a significant speedup and an enormous reduction of the memory requirements. This is achieved by selecting suitable subspaces and a special basis that captures the dominant stochastic effects. In the GSD, is approximated by (18) where

Consequently, the expected value and variance of a stochastic image are

, are deterministic functions; , , are stochastic functions; and is the number of modes of the decomposition. Thus, in contrast to expansion (13), in the GSD, the solutions are represented by deterministic and stochastic functions, which are not fixed a priori. These flexible basis functions allow finding a representation of , ) but nearly the which has significant less modes (i.e., same approximation quality. Nouy [14] showed how to achieve the modes of an optimal of the problem, i.e., approximation in the energy norm such that

(16)

(19)

i.e., the first coefficient of the gPC expansion is the expected value; thus, its determination requires no further computations.

For details about the GSD, proofs for the optimality of the approximation, and implementation details, we refer to [14].

for if

(15)

,

2428


IV. STOCHASTIC RANDOM WALKER SEGMENTATION The extension of the random walker segmentation to stochastic images is straightforward and follows the methods for processing of stochastic images described in [11] and [12]. To this end, we replace the classical images by stochastic images, as described in Section III. The edge weights for the graph are defined by the difference of two RVs, which describe the gray-value uncertainty for the corresponding pixels. We end up with stochastic edge weights, which are RVs, i.e., (20)

as the discrete version of the stochastic Dirichlet problem. This equation requires the same special ordering of vector and matrix as in the deterministic case. Here, the vector is organized by grouping all coefficients for a gPC polynomial together, i.e., . Matrix is a block matrix, with nonzero entries in the diagonal blocks only, i.e., (26) Reordering vector with respect to seeded or unseeded nodes and using the same stochastic ordering scheme for the new vecand yield tors

Note that we do not calculate the exponential of the gPC expansion explicitly. Instead, we compute an approximation of the exponential in the gPC via the Taylor series, i.e., (21) This method for the computation of the exponential of a gPC quantity, along with methods for other computations such as multiplication, division, computation of square roots, etc., are developed in [23]. It is straightforward to generalize the notion of node degrees and the normalization step for the image gradient to stochastic images as follows:

(27) Note that we have not denoted the dependence of , , , and on and to ease the notation. A stochastic minimizer of the discretized stochastic Dirichlet problem is given by (28) Using the gPC expansion for the discretization of the stochastic quantities and the GSD [14] for the solution of the resulting stochastic equation system, i.e., solving in a dominant subspace of the tensor product space , we get solutions that are random quantities.

(22) Here, for the normalization for the image gradient, we use the maximal expected value in the denominator. From the last equations, it is easy to build the stochastic analog of the Laplacian matrix as if , if are adjacent nodes, otherwise.

(23)

Using the gPC expansions of and , i.e., using the projections onto the gPC, we can represent the entries of the stochastic combinatorial Laplacian matrix in the gPC as well, i.e., (24) To define the linear system of equations for the stochastic random walker problem, we use the same procedure as in the deterministic case. An infinite number of random walks on a graph with stochastic weights is equivalent to solving the disare the crete weighted Dirichlet problem, where weights RVs from (20). Using the gPC discretization of stochastic images from Section III and collecting all pixels in vector , we get (25)

V. RESULTS In the following, we demonstrate the benefits of the stochastic extension of the random walker segmentation by applying the method to stochastic images generated from a collection of classamples with sical images. The first data set consists of a resolution of 100 100 pixels from the artificial “street sequence” [24]. Note that we do not consider the images as a sequence; instead, we treat them as five samples of the noisy and uncertain acquisition of the same scene. The second data set consists of 45 samples with a resolution of 300 300 pixels from an ultrasound device, originally used in [12], and the third data set consists of 25 samples with a resolution of 129 129 pixels of a liver mask in front of a varying background. We comRVs for the ultrasound data pute a stochastic image with set and RVs for the street scene and the liver data set. The number of RVs depends on the decay of the eigenvalues of the covariance matrix of the input samples. For all data sets, we use to compute a stochastic image from a polynomial degree these samples. This polynomial degree for the gPC expansion has proven to be a good balance between the accuracy of gPC expansion and the computational effort. The user defines seed points for the segmentation of the image on the expected value of the stochastic image (see Fig. 4). For the stochastic random walker segmentation on the stochastic image, we vary parameter in the experiments (see Fig. 3). Together with the expected value, we are able to show the variance for every pixel and the realizations of the contour encoded in the stochastic result (see


2429

Fig. 3. Expected value and variance of the probabilities that pixels belong to the object. (Red) Expected value of the object boundary is shown in (first and third row) the expected value pictures. (Red boundaries in the second and fourth row) Furthermore, the variance image together with ten Monte Carlo realizations of the stochastic object boundary is shown. A high variance indicates high influence of the input gray-value uncertainty. Even the small deviations of the object boundary might have an important influence on the volume of the object in 3-D. For comparison, we show the results of a classical random walker segmentation in the rightmost column.

Fig. 4). This gives information on how the gray-value uncertainty in the stochastic input image influences the segmentation results. Thus, the variance and the gPC coefficients of the result reveal how the gray-value uncertainty in the input image propagates through the segmentation process and how it influences the result, respectively. Regions with high variance indicate a high influence of the input gray-value uncertainty on the detection of the object. To highlight the uncertainty of the object-boundary detection, we draw ten Monte Carlo realizations of the stochastic object boundary into the variance images in Fig. 3. The contours are close together and visible through a

fluctuating width of the contour line only. It is obvious from the images that the uncertainty changes from the input to the output. In the input data, the gray-value uncertainty spreads over the image, whereas in the segmentation result, the gray-value uncertainty concentrates at the object boundary. Remark: The classical random walker result can be interpreted as a probability map, i.e., the result of random walker segmentation is a probability for pixels to belong to the object. When we apply the stochastic random walker method, the first interpretation of the result is that we compute “probabilities of probabilities,” because we compute the probability distribution

2430


= 25

Fig. 6. PDF of the segmented area from the street image for (black) and (gray) . From the PDF, it is possible to judge on the reliability of the segmentation; a narrow PDF for the segmented volume indicates that the segmentation process is only marginally influenced by the image noise.

= 50

object. Thus, the RV is Fig. 4. Expected (top row) value and (bottom row) variance of the images to be segmented. (Left) Street image. (Right) Ultrasound image. Color coded are the seed regions for (yellow) the interior and (red) the exterior.

specifying the volume of the object (29)

Evaluating the PDF makes it easy to decide whether the segmented volume is strongly influenced by the image noise or not. For a strong influence on the segmented volume, the PDF is broad; otherwise, it is narrow (see Fig. 6). Obviously, also, the choice of parameter influences the profile of the PDF. A smaller leads to a diffuse object boundary and therefore to a broader PDF. A. Comparison With Monte Carlo and Stochastic Collocation

Fig. 5. Realizations of the stochastic object boundary for the stochastic liver . On the right, we highlight a region of the image, where the image using segmentation result is influenced due to the noise in the input image. Dependent on the noise realization, we observe different boundary positions.

= 10

of pixel values. Let us emphasize that the probability interpretation of the classical result is one possible interpretation only. Mathematically, we compute the result of a diffusion problem, and the stochastic extension is, in this interpretation, a diffusion problem with the stochastic diffusion coefficient. The results can be interpreted as the stochastic solution of the stochastic diffusion problem. The analog to the probability interpretation is that we compute the probabilities for belonging to the object in dependence on the gray-value uncertainty. The stochastic object boundary can be visualized by tracking the deterministic object boundary (corresponding to the value of 0.5 in the result image) for realizations of the RVs, i.e., we perform Monte Carlo sampling on the result. Such a visualization is depicted in Figs. 3 and 5. Furthermore, it is possible to visualize other quantities that can be extracted from the stochastic segmentation result, e.g., the volume of the segmented object, which is a stochastic quantity. The PDF of the segmented volume is obtained from the segmentation result by summing up the RVs at pixels, because these RVs specify the “probability” that the pixels belong to the

To verify the intrusive solution of the resulting SPDE via gPC, stochastic finite elements and the GSD method, we compare the GSD solution with the solutions obtained via Monte Carlo sampling and stochastic collocation. The comparison of the expected value and the variance computed via GSD, stochastic collocation, and Monte Carlo sampling is depicted in Fig. 7. There is a small difference between the variances of the three solutions that might be due to the projection of the Laplacian matrix on the polynomial chaos. However, the great benefit of the GSD method is the significantly better performance. B. Performance Evaluation Due to the availability of several implementation possibilities for the numerical solution of SPDEs, we are able to compare the execution times of the approaches. The comparison of the execution times of the GSD method, the Monte Carlo method, and the stochastic collocation with sparse grid [25] and full grid is shown in Table I. Furthermore, we show the execution time of the classical (deterministic) random walker method on the data sets used. This allows evaluating the computational overhead of the stochastic implementation. The GSD method outperforms the sampling-based approaches. This supports the decision to prefer the GSD method. In fact, the stochastic collocation methods suffer from the “curse of dimension” [26], because the execution times exponentially grow with the number of RVs. The benefit of the Monte Carlo method is that it is independent on the number of RVs. Nevertheless, the 1000 samples used


2431

Fig. 7. Comparison of the discretization methods to verify the intrusive discretization. The small difference between the intrusive discretization via the GSD method and the two sampling-based discretization approaches might be due to the projection of the Laplacian matrix on the polynomial chaos.

TABLE I COMPARISON OF THE EXECUTION TIMES (IN SECONDS) OF THE DISCRETIZATION METHODS APPLIED ON THREE DATA SETS

Fig. 8. (Left) Input “doughnut” without noise and (right) noisy input image treated as an expected value of the stochastic image.

in this comparison are a lower bound for the number of runs needed to get accurate results. Remember that the rate of converfor Monte Carlo simulations. Even with gence is this number of runs, the Monte Carlo method is significantly slower than the GSD method. C. Segmentation and Volumetry of an Object In many applications, the noise of pixels in the image is independent from the noise of the neighboring pixels. It is possible to model this kind of stochastic images with the presented approach as well. We have to use one basic RV for every pixel, i.e., we end up with basic RVs for the gPC. To demonstrate the possibility of modeling such kind of images, we use an artificial test example, i.e., a “doughnut” with an area of 60 pixels in front of a constant background with a resolution of 20 20 pixels. We corrupt the image by uniform noise (see Fig. 8) and treat the noisy image as expected value of the stochastic image. This modeling is close to the situation in real applications. There, the real noise-free image is not available, and thus, the sample at hand is the best available estimate of the

expected value. Due to the high number of basic RVs, we restrict the polynomial chaos to a degree of 1, i.e., we are able to capture the effects expressible in uniform RVs only. Using a polynomial coefdegree of order 1, the polynomial chaos has ficients per pixel. Using a polynomial degree of 2, we end up coefficients. Contemporary personal comwith puters cannot store such a high number of stochastic modes. A solution could be the usage of sparse polynomial chaos expansions introduced in [27], but this has to be further investigated in subsequent work. After the initialization of the expected value with the noisy input image, we have to prescribe values for the remaining chaos coefficients of the input image. As we assume the noise at every pixel to be independent from the noise at other pixels, we have to prescribe a value for the coefficient corresponding to the basic , modeling a RV of the pixel. We set this coefficient to around uniform distributed RV with support the expected value given by the noisy input image.

2432


Fig. 9. Left: (Yellow) The object seed points and (red) the background seed points used as the initialization of the stochastic random walker method. Right: Realizations of the stochastic segmentation result. The object-boundary realizations significantly differ for noise realizations.

The result of the stochastic random walker segmentation on this stochastic image is again a stochastic image with the same number of basic RVs. Because a stochastic diffusion equation has to be solved in the random walker method, stochastic information is transported between the pixels. Thus, a pixel in the result image can depend on all basic RVs of the input image. The visualization of polynomial chaos coefficients of the solution is not intuitive and cumbersome, because we have 401 coefficients per pixel. Thus, we use the same visualization techniques as in the previous section. Realizations of the stochastic object boundary and the seed points for the segmentation are depicted in Fig. 9. In applications, features of the segmented object are of interest, e.g., in medical applications, the volume of the object is of interest to get information about the growth or the shrinkage of the segmented lesion. The volume of the segmented object in a stochastic image is a stochastic quantity as well, because it depends on the particular noise realization. Thus, it is possible to visualize the PDF of the object volume. We investigate two possibilities to compute the volume PDF from the stochastic segmentation result. The first method is introduced in Section V. There, the gPC expansions of all pixels are added, and the PDF of the resulting RV is computed via Monte Carlo sampling from this RV. This method can be compared with methods that are able to consider partial volume effects, because there is no binary decision whether pixels belong to the object or not. In fact, we add all the stochastic possibilities of the pixels to belong to the object. The classical method to compute the random walker result inspires the other possibility to compute the volume PDF from the stochastic result. We generate samples from the stochastic segmentation result via Monte Carlo sampling and estimate the volume of the objects given by pixels with value above 0.5 on every sample. The two approaches for the computation of the object’s volume are compared in Fig. 10. Having in mind that the “real” object volume is 60 pixels, both methods slightly overestimate the object volume, but the real object volume is close to the expected value (60.39 for the summation of the RVs and 60.83 for the object thresholding) of both volume PDFs. VI. CONCLUSION We have discussed an extension of the random walker segmentation to stochastic images whose gray values are RVs. The propagation of information about the gray-value

Fig. 10. PDF for both possibilities of the object’s volume computation. The summation of the polynomial chaos at every pixel is shown in gray; the object thresholding is shown in black. The true object’s volume is 60 pixels.

uncertainty in the input image through the segmentation step yields a stochastic quantity as the segmentation result. From the stochastic segmentation result, the PDF of quantities such as the segmented volume can be analyzed, and regions can be identified, where the object boundary is uncertain due to the noise in the input image. The presented algorithm combines the advantages of a supervised segmentation algorithm with the propagation of information about gray-value uncertainty, which was introduced in [11] and [12]. Such information about the influence of gray-value uncertainty is of interest for applications. In fact, aside from the mentioned medical applications, other engineering disciplines and tasks such as quality control or the extraction of quantitative information in general can benefit from the additional error propagation information. Aside from all these achievements concerning the stochastic random walker method, there are open questions and ideas for further research to improve the presented model. A starting point might be the integration of extensions presented for the classical random walker segmentation method (see, e.g., [8]). As already mentioned, the presented model has a limitation in the number of RVs and the polynomial order, because we need to store the representation in the tensor product space. A possibility to avoid this memory-consuming storage of all chaos coefficients is the sparse polynomial chaos expansion introduced in [27]. The idea is to store the important coefficients at every pixel only. Another important problem that should be addressed in subsequent publications is the generation of input data. The method needs a huge number of input samples to generate an accurate stochastic input image, or we have to use the one sample available as a rough approximation of the expected value. In many applications, the generation of input samples is unfeasible. To use the presented framework in this application as well, we will investigate noise models to generate the stochastic input image. The model used in this paper might not be sufficient for many applications, because the noise is correlated and may depend on the gray value. Furthermore, we need efficient methods for the solution of the resulting SPDE. A huge speedup compared with the stochastic


FEM is the GSD method used in this paper. Further ideas to speed up the computation can benefit from tensor approximation methods recently developed for SPDEs [28]. Aside from these directions for future research, we have worked on the extension of level-set-based segmentation approaches [3] to stochastic images. To use level sets on stochastic images, we have to propagate level sets with an uncertain speed. ACKNOWLEDGMENT The authors would like to thank R. M. Kirby, Scientific Computing and Imaging Institute, Utah, for inspiring discussions on the numerical treatment of SPDEs. REFERENCES [1] N. Metropolis and S. Ulam, “The Monte Carlo method,” J. Amer. Stat. Assoc., vol. 44, no. 247, pp. 335–341, Sep. 1949. [2] G. Aubert and P. Kornprobst, Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations. New York: Springer-Verlag, 2006. [3] J. A. Sethian, Level Set Methods and Fast Marching Methods. Cambridge, U.K.: Cambridge Univ. Press, 1999. [4] D. Mumford and J. Shah, “Optimal approximations by piecewise smooth functions and associated variational problems,” Commun. Pure Appl. Math., vol. 42, no. 5, pp. 577–685, Jul. 1989. [5] L. Ambrosio and M. Tortorelli, “Approximation of functionals depending on jumps by elliptic functionals via gamma convergence,” Commun. Pure Appl. Math., vol. 43, no. 8, pp. 999–1036, 1990. [6] L. Grady, “Random walks for image segmentation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 28, no. 11, pp. 1768–1783, Nov. 2006. [7] C. W. Therrien, “An estimation-theoretic approach to terrain image segmentation,” Comput. Vis. Graph. Image Process., vol. 22, no. 3, pp. 313–326, Jun. 1983. [8] L. Grady, “Multilabel random walker image segmentation using prior models,” in Proc. IEEE CVPR, Jun. 2005, vol. 1, pp. 763–770. [9] J. Weber and J. Malik, “Robust computation of optical flow in a multiscale differential framework,” Int. J. Comput. Vis., vol. 14, no. 1, pp. 67–81, Jan. 1995. [10] O. Nestares, D. J. Fleet, and D. J. Heeger, “Likelihood functions and confidence bounds for total-least-squares problems,” in Proc. IEEE Comput. Vis. Pattern Recogn., 2000, vol. 1, pp. 523–530. [11] T. Preusser, H. Scharr, K. Krajsek, and R. M. Kirby, “Building blocks for computer vision with stochastic partial differential equations,” Int J. Comput. Vis., vol. 80, no. 3, pp. 375–405, Dec. 2008. [12] T. Pätz and T. Preusser, “Ambrosio–Tortorelli segmentation of stochastic images,” in Proc. ECCV, K. Daniilidis, P. Maragos, and N. Paragios, Eds., 2010, vol. 6315, pp. 254–267, Berlin, Germany: Springer-Verlag. [13] D. Xiu and G. E. Karniadakis, “The Wiener–Askey polynomial chaos for stochastic differential equations,” SIAM J. Sci. Comput., vol. 24, no. 2, pp. 619–644, 2002. [14] A. Nouy, “A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations,” Comput. Methods Appl. Mech. Eng., vol. 196, no. 45–48, pp. 4521–4537, Sep. 2007. [15] P. G. Doyle and J. L. Snell, Random Walks and Electric Networks. Washington, DC: Math. Assoc. Amer., 1984. [16] M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” J. Res. Nat. Bur. Stand., vol. 49, no. 6, pp. 409–436, Dec. 1952. [17] N. Wiener, “The homogeneous chaos,” Amer. J. Math., vol. 60, no. 4, pp. 897–936, Oct. 1938.

2433

[18] R. H. Cameron and W. T. Martin, “The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals,” Ann. Math., vol. 48, no. 2, pp. 385–392, Apr. 1947. [19] R. Askey and J. Wilson, “Some basic hypergeometric polynomials that generalize Jacobi polynomials,” Mem. Amer. Math. Soc., vol. 319, pp. 1–55, 1985. [20] M. Loève, Probability Theory, 4th ed. New York: Springer-Verlag, 1977. [21] G. Stefanou, A. Nouy, and A. Clement, “Identification of random shapes from images through polynomial chaos expansion of random level set functions,” Int. J. Numer. Methods Eng., vol. 79, no. 2, pp. 127–155, Jul. 2009. [22] R. G. Ghanem and P. D. Spanos, Stochastic Finite Elements: A Spectral Approach. New York: Springer-Verlag, 1991. [23] B. J. Debusschere, H. N. Najm, P. P. Pébay, O. M. Knio, R. G. Ghanem, and O. P. Le Maıcirc;tre, “Numerical challenges in the use of polynomial chaos representations for stochastic processes,” SIAM J. Sci. Comput., vol. 26, no. 2, pp. 698–719, 2005. [24] B. McCane, K. Novins, D. Crannitch, and B. Galvin, “On benchmarking optical flow,” Comput. Vis. Image Underst., vol. 84, no. 1, pp. 126–143, Oct. 2001. [25] S. Smolyak, “Quadrature and interpolation formulas for tensor products of certain classes of functions,” Sov. Math. Dokl., vol. 4, pp. 240–243, 1963. [26] E. Novak and K. Ritter, “The curse of dimension and a universal method for numerical integration,” in Multivariate Approximation and Splines. Cambridge, MA: Birkaüser, 1998, pp. 177–187. [27] G. Blatman and B. Sudret, “Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach,” Comptes Rendus Mécanique, vol. 336, no. 6, pp. 518–523, Jun. 2008. [28] B. N. Khoromskij and C. Schwab, “Tensor-structured Galerkin approximation of parametric and stochastic elliptic PDEs,” SIAM J. Sci. Comput., vol. 33, no. 1, pp. 364–385, 2011. Torben Pätz received the Diploma degree in mathematics from the University of Bremen, Bremen, Germany, in 2009. He recently finished his Ph.D. dissertation at Jacobs University Bremen, Bremen, on the segmentation of stochastic images with stochastic partial differential equations. He is currently a Research Associate with Jacobs University Bremen and with the Fraunhofer Institute for Medical Image Computing MEVIS, Bremen. His research interests include uncertainty modeling and propagation in image processing, medical image processing, and modeling and simulation of biomedical processes.

Tobias Preusser studied mathematics at the University of Bonn, Bonn, Germany, and at New York University, New York. He received the Ph.D. degree from the University of Duisburg, Duisburg, Germany, with a thesis on anisotropic geometric diffusion in image processing, and the Habilitation degree in mathematics from the University of Bremen, Bremen, Germany, with a thesis on image-based computing. Currently, he is Professor of Mathematical Modeling of Biomedical Processes with Jacobs University Bremen, Bremen, as well as the Head of the Modeling and Simulation Group and Member of the Management Board with the Fraunhofer Institute for Medical Image Computing MEVIS, Bremen. His research is driven by concrete and complex application problems in medicine and includes modeling and simulation of biomedical processes with partial differential equations (PDEs), mathematical image processing, and scientific visualization based on PDEs.