Estimating targets using statistical properties of observations of known targets

[0002] This is a continuation-in-part of U.S. patent application Ser. No. 09/203,108, filed Nov. 30, 1998.

[0003] This invention relates generally to computer vision, and more particularly, to estimating characteristics of scenes represented by images.

[0004] One general problem in computer vision is how to determine the characteristics of a scene from images representing the underlying scene. Following are some specific problems. For motion estimation, the input is usually a temporally ordered sequence of images, e.g., a “video.” The problem is how to estimate the projected velocities of various things—people, cars, balls, background moving in the video. Another problem deals with recovering real-world three-dimensional (3D) structure from a 2D image. For example, how to recover the shape of an object from a line drawing, a photograph, or a stereo pair of photographs. Yet another problem is how to recover high-resolution scene details from a low-resolution image.

[0005] Humans make these types of estimates all the time, frequently subconsciously. There are many applications for machines to be able to do this also. These problems have been studied by many workers with different approaches and varying success for many years. The problem with most known approaches is that they lack machine learning methods that can exploit the power of modern processors within a general framework.

[0006] In the prior art, methods have been developed for interpreting blocks world images. Other prior art work, using hand-labeled scenes, has analyzed local features of aerial images based on vector codes, and has developed rules to propagate scene interpretations. However, these solutions are for specific one-step classifications, and therefore, cannot be used for solving a general class of low-level vision problems. Methods to propagate probabilities have been used, but these methods have not been put in a general framework for solving vision problems.

[0007] Alternatively, optical flow can be estimated from images by using a quad-tree to propagate motion information across scale. There, a brightness constancy assumption is used, and beliefs about the velocity of the optical flow is presented as a gaussian probability distribution.

[0008] The present invention analyzes statistical properties of a labeled visual world in order to estimate a visual scene from corresponding image data. The image data might be single or multiple frames; the scene characteristics to be estimated could be projected object velocities, surface shapes, reflectance patterns, or colors. The invention uses statistical properties gathered from labeled training data to form “best-guess” estimates or optimal interpretations of underlying scenes.

[0009] The invention operates in a learning phase and an inference phase. During the learning phase statistical properties about the training data are modeled as probability density functions, for example, a mixture of Gaussians. A Markov network is built. During the inference phase beliefs derived from a particular images and the density functions are propagated around the network to make estimates about a particular scene corresponding to the particular image.

[0010] Accordingly during the learning phase, training data for typical images and scenes are synthetically generated. A parametric vocabulary for both images and scenes is generated. The probability of image parameters, conditioned on scene parameters (the likelihood function), is modeled, as is the probability of scene parameters, conditioned on neighboring scene parameters. These relationships are modeled with a Markov network where local evidence is propagated to neighboring nodes during an inference phase to determine the maximum a posteriori probability of the scene estimate.

[0011] Humans perform scene interpretations in ways that are largely unknown and certainly mathematically undeterminable. We describe a visual system that interprets a visual scene by determining the probability of every possible scene interpretation for all local images, and by determining the probability of any two local scenes neighboring each other. The first probability allows the visual system to make scene estimates from local image data, and the second probability allows the these local estimates to propagate. One embodiment uses a Bayesian method constrained by Markov assumptions.

[0012] The method according to the invention can be applied to various low-level vision problems, for example, estimating high-resolution scene detail from a low-resolution version of the image, and estimating the shape of an object from a line drawing. In these applications, the spatially local statistical information, without domain knowledge, is sufficient to arrive at a reasonable global scene interpretation.

[0013] Specifically, the invention provides a method for estimating scenes from images. A plurality of scenes are generated and an image is rendered for each scene. These form the training data. The scenes and corresponding images are partitioned into patches. Each patch is quantified as a vector, and the vectors are modeled as a probability densities, for example, a mixture of gaussian distributions. The statistical relationship between patches are modeled as a Markov network. Local probability information is iteratively propagated to neighboring nodes of the network, and the resulting probability densities at each node, a “belief,” is read to estimate the scene.

[0014] In one application of our invention, it is possible to estimate high-resolution details from a blurred, or lower-resolution image. A low-resolution image is the input “image” data, and the “scene” data are the image intensities of the high-resolution details. The invention can also be used to estimate scene motion from a sequence of images. In this application, the image data are the image intensities from two successive images of the sequence, and the scene data are successive velocity maps indicating projected velocities of the visible objects at each pixel position. Another application for our invention is shading and reflectance disambiguation.

[0015] The invention can also be applied to other estimation problems where training and target data can be modeled with probability density functions, for example, in speech recognition, seismological studies, medical diagnostic signals such as EEG and EIKG. In addition, the probability representations can be either discrete or continuous in either the learning or inference phases.

[0016] FIG. 1aflow diagram of a method for estimating targets;

[0017] FIG. 1bis a flow diagram of details of a method for estimating scenes from images according to the invention;

[0018] FIG. 2 is a graph of a network for propagating beliefs of the method;

[0019] FIG. 3 is a graph of a true underlying joint probability relating a scene variable x with an image variable y;

[0020] FIG. 4ais a histogram of scene values observed in training data;

[0021] FIG. 4bis an initial mixture of gaussian distributions fit to the distribution shown in the histogram of FIG. 4a;

[0022] FIG. 4cis a pruned version of the fit of FIG. 4b;

[0023] FIGS. 5a-5cshow a mixture of gaussian fits to conditional probabilities observed in training data;

[0024] FIGS. 6a-6dgraphs of probabilities at various nodes of the network;

[0025] FIG. 6eis the product of the probabilities shown in FIGS. 6a-6d;

[0026] FIG. 6fis a graph of a conditional density;

[0027] FIG. 6gis the probability that is passed in a message;

[0028] FIG. 7 is a graph of the probabilities combined to form the belief of a node;

[0029] FIG. 8 is a graph of the initial probabilities;

[0030] FIGS. 9-10 are graphs of probabilities after one and two iterations; and

[0031] FIG. 11 is a graph of a Markov network having four scene and image nodes.

[0032] Introduction

[0033] We described a method that uses statistical properties of a labeled visual world to estimate characteristics of a scene using either a single image or multiple images. Scene characteristics to be estimated can include projected velocities of objects in the scenes, surface shapes of objects, reflectance patterns, or colors. Our general method can be applied to a number of low-level vision problems. Our method can also be used to model statistical properties of other complex digital signals such as human speech, or seismographs, for example.

[0034] As shown in FIG. 1a, our general method 1 begins with training data 2. The training data are observations of known targets. The training data 2 can randomly be generated. In step 10, we select a representation, discrete 11 or continuous 12, for modeling the training data. In step 20, using either the discrete or continuous representation, we learn statistical relationships about the training data. These relationships can be expressed as either discrete or continuous probability functions (P_d21, or P_c22), for example, vectors and matrices, or Gaussian mixtures.

[0035] After the learning phase, inferences about unknown target can be made. In step 30 using either P_d21 or P_c22, we infer a likely target 31 from observations 32 of the unknown target. Our inference is made by locally propagating of beliefs in a Markov network where the nodes in the network represent the observations, and the beliefs are discrete or continuous statistical representations of the belief. Step 30 can be repeated for other targets that resemble the training data.

[0036] Generate Random Scene and Image for Training Data

[0037] As shown in greater detail in FIG. 1b, during a training phase, our general method 100 generates random scenes x_i(known targets), and corresponding images y_i(observations) as training data 111 in step 110. The random scenes and rendered images can be generated synthetically using computer graphics. The synthetic images should be somewhat typical of unknown images to be processed by the system.

[0038] Partition Scene and into Patches

[0039] In step 120, the scenes and corresponding images are partitioned into local patches 121. The partitioning can be a patchwork of squares that overlay the scene and image. The patches can have multiple sizes, and the patches can be redundantly overlaid. For example, the patches can be formed at multiple of a gaussian pyramid. The pyramid can have, for example, five levels of resolution - from fine to coarse. In addition, the patches can represent image information as viewed through oriented filters of different orientations.

[0040] All patches of a given set of criteria, such as resolution and orientation, but differing spatially, are said to be of the same class, and are assumed to be drawn from the same statistical distribution. The sizes of the patches are small enough to allow modeling, yet large enough to convey meaningful information about the overall scene.

[0041] Quantify Patches as Vectors

[0042] In step 130, principle components analysis (PCA) is used to determine a representation for each patch. Each patch is expressed as a linear combination of basis functions. We represent the patches 121 as a low-dimensional vector 131. For example, each scene patch can be represented as a five-dimensional vector, and each image patch as a seven-dimensional vector. In other words, we represent the patches of the random training data, scenes and images respectively, as points in, say, a five-dimensional and seven-dimensional space.

[0043] Model Probability Densities of Training Data

[0044] Step 140 models the probability densities of all the training data in these low-dimensional spaces with mixtures of gaussian distributions. We use the training data to estimate probabilities of local patches in a most general form as:

[0045] P(scene), P(image\scene) and P(neighboring scene\scene).

[0046] More formally, we model three probability densities 141:

[0047] (1) the prior probability of each scene element x, there is a different prior probability for each class of scene element;

[0048] (2) the conditional probability of a scene element x, given the associated image element y: that is, P(y\x); and

[0049] (3) the conditional probability of a scene element x₁and a neighboring scene element, x₂, that is P(x₁\x₂).

[0050] A neighboring element may be one that is adjacent in spatial position, but also may be one that is nearby in any one of the class attributes, such as scale or orientation.

[0051] It may be useful to modify the training data so that it has a probability distribution that is easier to fit with the mixtures of gaussian distributions. For real images, many distributions of interest have a very sharp spike at the origin. This peak is hard to fit and manipulate with mixtures of gaussians. From our statistical analysis of labeled visual data, we can find the prior probabilities of the scene data. We can than make a second pass through the training data, and randomly omit each training sample with a probability inversely proportional to the prior probability of the scene data. This gives a biased set of data which will have an easier-to-model probability distribution.

[0052] Build Markov Network

[0053] In the last step 150 of the learning phase, the patches and their associated probability densities are organized into a Markov network 200 that represents the statistical relationship of scenes and images. In the Markov network, the nodes represent the low-dimensional vectors, nodes x_ithe scenes, and nodes y_ithe images. The edges connecting the nodes represent the statistical dependencies between the nodes.

[0054] Also, if a gaussian pyramid is used, then nodes at a given resolution level can be connected to spatially neighboring nodes at the same level as well as to nodes at the same spatial position in an adjacent resolution level. Furthermore, connections can be made to scene elements varying in some other dimension, such as oriented filter orientation.

[0055] These connections assist in removing spatial artifacts while estimating scenes. The connected Markov network 200 allows each scene node to update its belief based on accumulated local evidence gathered from other nodes during the following inference phase. The belief is the combined probability densities that will form the final best estimate.

[0056] Iteratively Propagate Beliefs and Read Best Estimate

[0057] During the inference phase, we estimate an unknown target scene 171 from unknown observations or an image 172. Based on rules described below, step 160 iteratively propagates the Bayesian “belief” at each node to neighboring nodes via messages 161. Bayesian or regularization approaches have been used in low-level vision problems in the past. However, to contrast with the prior art, we use labeled training data (scenes and corresponding images), and use strong Markovian assumptions.

[0058] In step 170, the best estimates 171 at each node for the corresponding scene is read given the observed image information. This can be done by examining the probability distribution for the belief at each node, and either taking the mean or the maximum value of that mixture of gaussians distribution. This can tell us what scene value is the best estimate for the true underlying target scene at that position, given the observed image data.

[0059] Example of a 3×3 Markov Network

[0060] FIG. 2 shows a simple 3×3 Markov network 200. For simplicity, we make all data one-dimensional so that the data can be plotted. The “scene data” to be estimated are the 1D x 201 at each node. The 1D image data y 202 that falls on each node will be used to estimate what x is.

[0061] In normal use of the invention, we generate randomly composed computer graphic scenes (known targets) and their corresponding rendered images (observations) in order to create a training set of the images and scenes. Those are used to generate vectors representing training patches of images and scenes from which we gather the desired prior and conditional statistics.

[0062] For this simple example, however, we will make up synthetic data corresponding to the vectors which would represent the training patches of images and scenes. We will make up the underlying joint probability relationships governing images and scenes.

[0063] FIG. 3 shows the joint probability relationship 300 of the variables x and y for the simple example. In FIG. 3, the variable x is along the horizontal axis 301, and the variable y is along the vertical 302. If y is zero, then the variable x can have one of many possible values as shown by the broad distribution of fuzzy horizontal line 303 at the middle of FIG. 3. If the observation y is 2, then x will be somewhere near 3.

[0064] Furthermore, in our simple example, the relationship between the values x of neighboring scene patches is as follows: we always multiply the scene data x by 2 when we go down a “row” 203 of the network 200, and we multiply the scene data x by 1.5 when we go to the right one column 204.

[0065] For our simple example, we make up the image data y which falls on the nodes. Again for simplicity, all the nodes, except node 5 are set with y=0.

[0066] Therefore, all nodes have broad uncertainty as to their values. Node 5 has an observed value y=2. In this case, the observed value of the center node 5 should almost certainly be 3. The Bayesian belief propagation will then take care of propagating that knowledge to all the other nodes in the network 200. The final estimate should be at node 5 where x=3, and other nodes should have their x values scaled up or down by factors of 1.5 or 2 for each right horizontal or downward step, respectively, away from node 5 (and scale by 1/1.5 and 1/2 going in the other directions).

[0067] The example network 200 is a nine node tree where the nodes are numbered sequentially, starting with 1 at the root of the tree. A local scene state of the ith node is x_i, and the image evidence at the ith node is y_i.

[0068] Following the steps of our general method 100 outlined above, we proceed as follows. We collect training data from a computer graphic simulation of the problem. For this example problem, we generate simulated data by drawing from the known joint distributions of y and x, and of x₁and its neighboring node's x_2.

[0069] For the simple 1D problem, it is not necessary to perform principle components analysis (PCA) to reduce the dimensionality of the data collected at each node. We next estimate the desired joint probabilities using a mixture of gaussians probability model, please see Bishop “Neural networks for pattern recognition,” Oxford, 1995.

[0070] FIG. 4ashows histograms of the observed values of x, FIG. 4ba mixture of gaussian distributions fit to that prior probability density, and FIG. 4ca pruned version of that mixture of gaussians. For reasons described below, we prune after each multiplication, or fitting of probabilities.

[0071] FIGS. 5a-cshow mixtures of gaussians fit to some of the needed conditional probabilities 141. We fit a mixture of gaussian to model conditional probabilities P(a\b) by using the joint data for a and b co-occurrences, and weighting the points by 1/P(b), given P(a, b)/P(b) =P(a\b). FIG. 5ashows a mixture of gaussian fit to the probability density y given x, FIG. 5ba mixture of gaussian fit to the probability densities of the neighbor to the right of x, given the value of x, a line with a slope of 1/1.5. FIG. 5cshows a mixture of gaussian fit to the probability densities of the neighbor below x, given the value of x, a line with a slope of 1/2.

[0072] We iteratively calculate the belief at each node, according to the rules described below. The first step is to determine what message to pass from each node to each of its neighbors.

[0073] FIGS. 6a-6dshow graphically each of the probabilities to be multiplied together to generate the message that node 5 passes to the node above it, node 4, on the first iteration. FIG. 6ais the probability from the image, FIG. 6bis from node 2, FIG. 4cfrom node 6, FIG. 6dfrom node 8.

[0074] FIG. 6eis the product of the distributions shown in FIGS. 6a-6d.The dimensionality of the distribution shown in FIG. 6eis then augmented to equal the dimensionality of the distribution shown in FIG. 6fby keeping the distribution constant in the dimensions contained in 6fbut not in 6e.That augmented distribution is then multiplied by the conditional density shown in FIG. 6f,and then marginalized along the dimensions of the distribution contained in FIG. 6e.The result is the message 161 sent from node 5 to node 4, shown in FIG. 6g.

[0075] FIG. 7 shows graphically the probabilities being multiplied together to compute the message that node 5 sends to node 4 in the order of prior probability, the message from the local image data, from the local image data, from neighboring node 4, node 2, node 6, and node 8, and the final belief (estimate) from the image at node 5 at the end of the first iteration.

[0076] FIGS. 8-10 show the “belief” at each node in the network during the first three iterations of our method. As shown in FIG. 8, no information has yet been passed between the nodes, and each node relies only on its local image information, y to estimate its x value. Because all the nodes, except node 5, saw y=0, they have received little information about their x values, and their beliefs about their x values are very broad distributions. Node 5 knows its x value is near 3 because this is implied by y=2. The belief shown at each node is P(y\x) P(x), for the appropriate value of y at each node.

[0077] On the second pass, every node has shared its information with its neighboring nodes as shown in FIG. 9. Nodes 2, 4, 6, and 8 have received informative messages from node 5, the only node that knows what value of x it probably has, and these nodes adjust their beliefs about their values of x accordingly. The distributions shown at each node are P(y\x) P(x) times messages from each of the node's neighbors.

[0078] By the third pass, every node has heard from all neighbors two away, and therefore every node has received the knowledge from node 5. After the third pass, the mean or maximum value of each node's belief is approximately what it should be: node 5 has x has the approximate value of 3, with the other x values scaled up or down by 1.5 for jumps to the right, and 2 for jumps down.

[0079] Mixture Pruning

[0080] A probability mixture of N gaussians, when multiplied by a probability mixture of M gaussians, yields a mixture of NM gaussians. Thus, when mixtures of gaussians are multiplied together, the number of gaussians rises quickly and we must prune the gaussians. One could simply threshold out gaussians with very small weights from the mixture, however, this can give inaccurate mixture fittings. A better pruning process is given in Appendix A. Appendix B details the step of a process to determine the Kulbeck-Leibler distance between the pruned and unpruned mixtures.

[0081] Factorize the Joint Probability

[0082] The details of factorizing the joint probabilities that are used for propagating local evidence to neighboring nodes is described with reference to FIG. 11.

[0083] The network shown in FIG. 11 has four scene and image nodes, respectively:

[0084] x₁, . . . , x₄, and y₁, . . . , y₄.

[0085] We seek a factorization of the joint probability that yields local evidence propagation rules. Our factorization uses three probability manipulation rules repeatedly.

[0086] Rule [1] According to elementary probability: P(a, b) =P(a\b)P(b).

[0087] Rule [2] If node b is between nodes a and c, then P(a, c\b)=P(a\b)P(c\b). This is the statement of the conditional independence of a and c, given b.

[0088] Rule [3] If node b is between nodes a and c, then P(c\a,b)=P(c\b). This is the Markov property that lets us summarize our knowledge of the rest of the chain by our knowledge of its nearest node.

[0089] Note, none of these three rules requires the edges connecting the nodes to be directed. This avoids having to make arbitrary choices about causal relationships in the network 200.

[0090] For the maximum a posteriori (MAP) probability estimate of the parameters x₁, x₂, x₃, x₄, we want to determined argmax_{x1,x2, x3, x4}P(x₁, x₂, x₃, x₄\y₁, y₂, y₃, y₄). This conditional probability differs from the joint probability P(x₁, x₂, x₃, x₄, y₁, y₂, y₃, y₄) by a factor that is constant over the arguments varied. Therefore, we can equivalently choose to find argmax_{x1, x2, x3, x4}P(x₁, x₂, x₃, x₄, y₁, y₂, y₃, y₄), which is simpler to determine.

[0091] An alternate useful estimate of each parameter x_iis the mean of the marginalized distribution, P(x_i\y₁, y₂, y₃, y₄). The mean can be found from the joint distribution, P(x₁, x₂, x₃, x₄, y₁, y₂, y₃, y₄), by marginalizing out (integrating over) all the x parameters except for x_i. The marginalization will yield P(x_i, y₁, y₂, y₃, y₄). This is related by a constant scale factor to the distribution P(x_i\y₁, y₂, y₃, y₄), and thus the mean values of the two distributions will be the same. The factorization steps that follow for the MAP estimate also hold for the mean of the marginalized distribution, with the following change: the operation argmax_xjshould be replaced by integration over the variable x_j(∫x_j). The final argmax operation over the belief at a node is replaced by taking the mean of that belief.

[0092] We factorize the joint probability differently for the computation at each node. Each node j will account for P(x_j) at the very end of its calculation, and never passes that quantity to neighboring nodes. This will lead to a uniform local evidence propagation algorithm, and outputs that are always optimal given the number of nodes which have been reported.

[0093] To proceed by example, we will described the four different cases for each of the four nodes in the network 200. First, we describe the factorization to be performed at each node so that argmaxj at node j gives the same value as:

[0094] argmax_{x1, x2, x3, x4}P(x₁, x₂, x₃, x₄\y₁, y₂, y₃, y₄).

[0095] After the four cases, we present general local evidence propagation rules, which perform the computations of each factorization.

[0096] Calculation at Node 1

[0097] By applying rule 1, then rule 2, we have: $\begin{array}{c}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=P\left(x_2,x_3,x_4,y_1,y_2,y_3,y_4|x_1\right)P\left(x_1\right)\\ =P\left(y_1,x_1\right)P\left(x_2,x_3,x_4,y_2,y_3,y_4|x_1\right)P\left(x_1\right)\end{array}$

[0098] Applying rule 1, then 3, we continue the factorization: $\begin{array}{c}P\left(x_2,x_3,x_4,y_2,y_3,y_4|x_1\right)=P\left(x_3,x_4,y_2,y_3,y_4|x_1,x_2\right)P\left(x_2|x_1\right)\\ =P\left(x_3,x_4,y_2,y_3,y_4|x_2\right)P\left(x_2|x_1\right)\end{array}$

[0099] Applying rule 2 twice:

[0100] P(x₃, x₄, y₂, y₃, y₄\x₂)=P(y₂\x₂) P(x₃, y₃\x₂) P(x₄, y₄\x₂)

[0101] Applying rule 1, then rule 3: $\begin{array}{c}P\left(x_3,y_3|x_2\right)=P\left(y_3|x_2,x_3\right)P\left(x_3|x_2\right)\\ =P\left(y_3|x_3\right)P\left(x_3|x_2\right)\end{array}$

[0102] and $\begin{array}{c}P\left(x_4,y_4|x_2\right)=P\left(y_4|x_2,x_4\right)P\left(x_4|x_2\right)\\ =P\left(y_4|x_4\right)P\left(x_4|x_2\right)\end{array}$

[0103] Applying all those substitutions gives: $\begin{array}{c}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=P\left(x_1\right)P\left(y_1|x_1\right)\\ P\left(x_2|x_1\right)P\left(y_2|x_2\right)\\ P\left(x_3|x_2\right)P\left(y_3|x_3\right)\\ P\left(x_4|x_2\right)P\left(y_4|x_4\right)\end{array}$

[0104] Letting the argmax's slide through variables over which they are constant, we have: $\begin{array}{c}\underset{\mathrm{x1},\mathrm{x2},\mathrm{x3},\mathrm{x4}}{\mathrm{argmax}}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=\underset{\mathrm{x1}}{\mathrm{argmax}}P\left(x_1\right)P\left(y_1|x_1\right)\\ \underset{\mathrm{x2}}{\mathrm{argmax}}P\left(x_2|x_1\right)P\left(y_2|x_2\right)\\ \underset{\mathrm{x3}}{\mathrm{argmax}}P\left(x_3|x_2\right)P\left(y_3|x_3\right)\end{array}$

[0105] The above result is for finding the MAP estimate of the joint posterior probability. As described above, to instead find the mean of the marginalized distribution, we would take the mean over x₁of the distribution:

[0106] P(x₁, y₁, y₂, y₃, y₄)=P(x₁)P(y₁\x₁)∫_x2P(x₂\x₁)P(y₂x₂)∫_x3P(x₃\x₂)P(y₃\x₃)

[0107] Generalization

[0108] We used rule 1 to have P(x_a) appear at node a. By rule 2, each edge leaving node a will contribute a factor of the form P( other variables \x_a). Each of those “other variable” strings are broken down again using rules 1 and 2 to simplify any extra conditioning variables by using rule 3.

[0109] This factorizes the joint probability in a way that reflects the network topology, from the point of view of node a. For a 3-node chain, where nodes b and c fork off from node a, we have:

[0110] P(x_a, x_b, x_c)=P(x_a) P(x_b\x_a)P(x_c\x_a).

[0111] Including the images y forking off from each node, we have:

[0112] P(x_a, x_b, x_c, y_a, y_b, y_c)=P(x_a)P(y_a\x_a)P(x_b\x_a)P(y_b\x_b)P(x_c\x_a)P(y_c\x_c).

[0113] Calculation at Node 2

[0114] We use the three manipulation rules to write a different factorization that is used at node 2. Now the only prior probability over a single variable will be P(x₂). $\begin{array}{c}\underset{\mathrm{x1},\mathrm{x2},\mathrm{x3},\mathrm{x4}}{\mathrm{argmax}}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=\underset{\mathrm{x2}}{\mathrm{argmax}}P\left(x_2\right)P\left(y_2|x_2\right)\\ \underset{\mathrm{x1}}{\mathrm{argmax}}P\left(x_1|x_2\right)P\left(y_1|x_1\right)\\ \underset{\mathrm{x3}}{\mathrm{argmax}}P\left(x_3{|}_{\mathrm{x2}}\right)P\left(y_3|x_3\right)\\ \underset{\mathrm{x4}}{\mathrm{argmax}}P\left(x_4|x_2\right)P\left(y_4|x_4\right).\end{array}$

[0115] Calculation at Node 3

[0116] Factorizing P(x₁, x₂, x₃, x₄, y₁, y₂, y₃, y₄) to leave a factor of: $\begin{array}{c}P\left(x_3\right),\underset{\mathrm{x1},\mathrm{x2},\mathrm{x3},\mathrm{x4}}{\mathrm{argmax}}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=\underset{\mathrm{x3}}{\mathrm{argmax}}P\left(x_3\right)P\left(y_3|x_3\right)\\ \underset{\mathrm{x2}}{\mathrm{argmax}}P\left(x_2|x_3\right)P\left(y_2|x_2\right)\\ \underset{\mathrm{x1}}{\mathrm{argmax}}P\left(\mathrm{x1}|x_2\right)P\left(y_1|x_1\right)\\ \underset{\mathrm{x4}}{\mathrm{argmax}}P\left(x_4|x_2\right)P\left(y_4|x_4\right)\end{array}$

[0117] Calculation at Node 4

[0118] Factorizing P(x₁, x₂, x₃, x₄, y₁, y₂, y₃, y₄) to leave a factor of: $P\left(x_4\right),\underset{\mathrm{x1},\mathrm{x2},\mathrm{x3},\mathrm{x4}}{\mathrm{argmax}}P\left(x_1,x_2,x_3,x_4,y_1,y_2,y_3,y_4\right)=\underset{\mathrm{x4}}{\mathrm{argmax}}P\left(x_4\right)P\left(y_4|x_4\right)\mathrm{argmax}\mathrm{x2P}\left(x_2|x_4\right)P\left(y_2|x_2\right)\underset{\mathrm{x1}}{\mathrm{argmax}}P\left(x_1|x_2\right)P\left(y_1|x_1\right)\underset{\mathrm{x3}}{\mathrm{argmax}}P\left(x_3|x_2\right)P\left(y_3|x_3\right)$

[0119] Local Propagation Rules

[0120] A single set of propagation rules will result in each of the four calculations above arriving at the four different nodes.

[0121] During each iteration, each node xj gathers evidence, then passes an appropriate message to each connecting node x_k. The evidence from node k is the most recent message received from it. The evidence from the image y_iis P(y_j\x_j).

[0122] (1) The message sent to node k from node j starts with the product Q(j;k) of the evidences at node j from the nodes other than node k, the later being the node receiving the message. This includes the evidence P(y_j\x_j) of the local node.

[0123] (2) The message sent to node k is then argmax_xjP(x_j\x_k) Q(j;k). We use a different calculation to read-out the optimum x_jfrom node j.

[0124] (3) To find the x_jwhich maximizes P(x₁, x₂, x₃, x₄, y₁, y₂), take argmax_xjover the product of all the evidences at node j multiplied times P(x_j).

[0125] Propagation Rules, Discrete Case

[0126] Our propagation operations may be easier to articulate for the case of discrete probability representations. In this embodiment we use discrete probability representations during both the learning and inference phases. During training, we measure the co-occurrence histograms H(y_j, x_j) and H(x_j, x_k), for node j that is next to node k. From these histograms, we can estimate P(y_j\x_j) and P(x_j\x_k). If we store the concurrence histogram H(a, b) as a matrix with rows indexed by a and columns indexed by b, then P(a\b) is the row-normalized version of that matrix, after adding a small constant to each count for Poisson arrival statistics. Each row sums to one.

[0127] Node j receives a column-vector message from each node. To send a message from node j to node k, node j:

[0128] (1) term-by-term multiplies each incoming message (except that from node k together, and multiplies in the column-vector P(y_j\x_j), then

[0129] (2) “max-matrix-multiplies” that resulting vector times P(x_j\x_k).

[0130] The resulting column vector is the message to node k.

[0131] The term “max-matrix-multiplication” means term-by-term multiply the column vector times each row of the matrix and set the output for that index of the output column vector equal to the maximum of the multiplied products. For the minimum mean square error (MMSE) estimate, we replace the max-matrix-multiplication step by conventional multiplication of a vector times a matrix.

[0132] To read-out the best estimate of x at node j in the discrete probability representation, we term-by-term multiply the most recent message from each connecting node, multiply in column-vector P(y_j\x_j), and the column-vector P(x_j). The index which maximizes that resulting column vector is the best estimate of x, that is in the target scene.

[0133] Super-resolution Problem

[0134] In one application of our invention, we estimate high-resolution details from a blurred, or lower-resolution image. In this application, the image data are the image intensities of the low-resolution image, and the “scene” data are the image intensities of the high-resolution details.

[0135] Our training images start from randomly shaped blobs covered with random surface markings, rendered by computer graphics techniques. To first obtain a bandpassed version of the image, by applying an oriented, bandpassed filter. We apply a spatially varying multiplicative local gain control factor to this bandpassed image. The gain control factor is computed as the square root of the squared and blurred value of the bandpass filtered image. This contrast gain control normalizes the intensities of image edges, easing the burden on the subsequent modeling steps. The resulting image represents the “observed” information.

[0136] We also apply an oriented, high-pass filter to the rendered image, and then apply the spatially varying local gain control factor that was computed from the bandpassed image. The result represents the corresponding target or “scene” information.

[0137] We generated many such image and scene pairs to build the training data. We broke each image and scene pair into patches in a uniform gridwork at a single spatial scale. We applied PCA separately on the image and scene patches to obtain a lower dimensional representation for each patch.

[0138] We determined the needed conditional and prior probabilities from the training data, and fit mixtures of gaussians to that data. We propagated the local information to obtain an estimated high resolution image. In this embodiment, we use a continuous probability representation during both the learning and inference phase.

[0139] Hybrid Probability Density Representation

[0140] Our method of propagating probability densities described above can be improved upon, in terms of processing speed. The continuous representation of the probability densities allows good fitting of the input image data during the learning phase. A discrete representation allows fast propagation during the inference phase. We will now described a hybrid method that allows both good fitting and fast propagation.

[0141] In the hybrid case, we evaluate the prior and conditional distributions 141 of FIG. 1bonly at a discrete set of scene values, different for each node 201 in the Markov network 200 of FIG. 2. The scenes values are a sampling of those scenes which render to the image at that node. This focuses the computation on the locally feasible scene interpretations. The conditional probability P(x_j\x_k) expresses the ratios of the gaussian mixtures P(x_j, x_k) and P(x_k) evaluated at scene samples at node j and node k, respectively. The conditional probability P(y\x_k) is the probability P(y_k, x_k)/P(x_k) evaluated at the scene samples of node k, for the observed image value, y_k, there.

[0142] Thus, instead of multiplying mixtures of gaussians together in order to propagate information between nodes as described above, we multiply together probability samples at a discrete set of points in the scene domain. The set of scene samples is different at each node of the network 200, and depends on the image information falling on that node. The belief at one node is a set of probability weights on each of the scene samples at that node. The belief propagation from nodej to node k now involves a point-by-point multiplication of the vector of discrete samples from Q(j:k) times each row of the link matrix, P(x_j\x_k) (from propagation rule 1). Following propagation rule 2, the values of the resultant vector is the maximum of the products of each row. For example, if we use between ten and fifty samples, then it is possible to describe the underlying scene adequately while still reducing the processing time. Using a discrete representation during the inference improves processing speed from twenty-four hours to about ten minutes, over two orders of magnitude.

[0143] To select the scene samples, we can condition the mixture P(y, x) on the image element y observed at each node, and sample the scene elements x from the resulting mixture of gaussians. Better results can be obtained by only using those scenes from the training set whose images most closely matched the image observed at that node. This avoids one gaussian mixture modeling step.

[0144] Our scene estimation method in a super sampling application provides high quality zooming of low-resolution images.

[0145] Our invention can also be used to estimate scene motion from a sequence of images. In this application, the image data are the image intensities from two successive images of the sequence, and the scene data are successive velocity maps indicating projected velocities of the visible objects at each pixel position.

[0146] Another application for our invention is shading and reflectance disambiguation. Images can arise from shading effects on a surface, as well as from reflectance changes on the surface itself. For example, an image of a shaded surface could arise from the shaded surface itself, or from a flat surface painted to look like the shaded surface (e.g., a flat picture of it). The image data for that application would be the images themselves. The underlying scene data to be estimated would be the underlying surface shapes and reflection patterns. This method could be used to make a best estimate of the 3D scene and paint pattern depicted by an image.

[0147] Our invention can also be used to make estimations about other complex digital signals. For example, our invention can be used for any signal that can be statistically sampled, and represented as a probability density function, such as speech, seismic data, diagnostic medical data, and so forth.

[0148] In this description of the invention, we used specific terms and examples. It is to be understood that various other adaptations and modifications may be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.