System and method for image matting

[0001] This invention relates generally to image editing, and more particularly to matting.

[0002] Matting and compositing are frequently used in image editing, 3D photography, and film production. Matting separates a foreground region from an input image by estimating a color F and an opacity α for each pixel in the image. Compositing blends the extracted foreground into an output image, using the matte, to represent a novel scene.

[0003] The opacity measures a ‘coverage’ of the foreground region, due to either partial spatial coverage or partial temporal coverage, i.e., motion blur. The set of all opacity values is called the alpha matte, the alpha channel, simply a matte.

[0004] Matting is described generally by Smith et al., “Blue screen matting, “Proceedings of the 23rd annual conference on Computer graphics and interactive techniques,” ACM Press, pp. 259-268, and U.S. Pat. No. 4,100,569, “Comprehensive electronic compositing system,” issued to Vlahos on Jul. 11, 1978.

[0005] Conventional matting requires a background with known, constant color, which is referred to as blue screen matting. If a digital camera is used, then a green matte is preferred.

[0006] Blue screen matting is the predominant technique in the film and broadcast industry. For example, broadcast studios use blue matting for presenting weather reports. The background is a blue screen, and the foreground region includes the weatherman standing in front of the blue screen. The foreground is extracted, and then superimposed onto a weather map so that it appears that the weatherman is actually standing in front of the map.

[0007] However, blue screen matting is costly and not readily available to casual users. Even production studios would prefer a lower-cost and less intrusive alternative.

[0008] Rotoscoping permits non-intrusive matting, Fleischer 1917, “Method of producing moving picture cartoons,” U.S. Pat. No. 1,242,674. Rotoscoping involves the manual drawing of a matte boundary on individual frames of a movie.

[0009] Ideally, one would like to extract a high-quality matte from an image or video with an arbitrary, i.e., unknown, background. This process is known as natural image matting.

[0010] Recently, there has been substantial progress in this area, Ruzon et al., “Alpha estimation in natural images,” CVPR, vol. 1, pp. 18-25, 2000, Hillman et al., “Alpha channel estimation in high resolution images and image sequences,” Proceedings of IEEE CVPR 2001, IEEE Computer Society, vol. 1, pp. 1063-1068, 2001, Chuang et al., “A bayesian approach to digital matting,” Proceedings of IEEE CVPR 2001, IEEE Computer Society, vol. 2, pp. 264-271, 2001, Chuang et al., “Video matting of complex scenes,” ACM Trans. on Graphics 21, 3, pp. 243-248, July, 2002, and Sun et al, “Poisson matting,” ACM Trans. on Graphics, August 2004.

[0011] Unfortunately, all of those methods require substantial manual intervention, which becomes prohibitive for long image sequences and for non-professional users.

[0012] The difficulty arises because matting from a single image is fundamentally under-constrained. The matting problem considers the input image as a composite of a foreground layer F and a background layer B, combined using linear blending of radiance values for a pinhole camera:
I_p[x,y]=αF+(1−α)B,   (1)

[0013] where αF is the pre-multiplied image of the foreground regions against a black background, and B is the image of the opaque background in the absence of the foreground.

[0014] Matting is the inverse problem of solving for the unknown values of variables (α, F_r, F_g, F_b, B_r, B_g, B_b) given the composite image pixel values (I_Pr, I_Pg, I_Pb). The ‘P’ subscript denotes that Equation (1) holds only for a pinhole camera, i.e., where the entire scene is in focus. One can approximate a pinhole camera with a very small aperture. Blue screen matting is easier to solve because the background color B is known.

[0015] It desired to perform matting using non-intrusive techniques. That is, the scene does not need to be modified. It is also desired to perform the matting automatically. Furthermore, it is desired to provided matting for ‘rich’ natural image, i.e., images with a lot of fine, detailed structure, such as outdoor scenes.

[0016] Most natural image matting methods require manually defined trimaps to determine the distribution of color in the foreground and background regions. A trimap segments an image into background, foreground and unknown pixels. Using the trimaps, those methods estimate likely values of the foreground and background colors of unknown pixels, and use the colors to solve the matting Equation (1).

[0017] Bayesian matting, and its extension to image sequences, produce the best results in many applications. However, those methods require manually defined trimaps for key frames. This is tedious for a long image sequences.

[0018] It is desired to provide a method that does not require user intervention, and that can operate in real-time as an image sequence is acquired.

[0019] The prior art estimation of the color distributions works only when the foreground and background are sufficiently different in a neighborhood of an unknown pixel.

[0020] It is desired to provide a method that can extract a matte where the foreground and background pixels have substantially similar color distributions.

[0021] The Poisson matting of Sun et al. 2004 solves a Poisson equation for the matte by assuming that the foreground and background are slowly varying. Their method interacts closely with the user by beginning from a manually constructed trimap. They also provide ‘painting’ tools to correct errors in the matte.

[0022] A method that acquires pixel-aligned images has been successfully used in other computer graphics and computer vision applications, such as high-dynamic range (HDR) imaging, Debevec and Malik, “Recovering high dynamic range radiance maps from photographs,” Proceedings of the 24th annual conference on Computer graphics and interactive techniques, ACM Press/Addison-Wesley Publishing Co., pp. 369-378, and Branzoi, “Adaptive dynamic range imaging: Optical control of pixel exposures over space and time,” Proceedings of the International Conference on Computer Vision (ICCV), 2003.

[0023] Another system illuminates a scene with visible light and infrared light. Images of the scene are acquired via a beam splitter. The beam splitter directs the visible to a visible light camera and the infrared light to an infrared camera. That system extracts high-quality mattes from an environment with controlled illumination, Debevec et al., “A lighting reproduction approach to live action compositing,” ACM Trans. on Graphics 21, 3, pp. 547-556, July 2002. Similar systems have been used in film production. However, flooding the background with artificial light is impossible for large natural outdoor scenes illuminated by ambient light.

[0024] An unassisted, natural video matting system is described by Zitnick et al., “High-quality video view interpolation using a layered representation,” ACM Trans. on Graphics 23, 3, pp. 600-608, 2004. They acquire videos with a horizontal row of eight cameras spaced over about two meters. They measure depth discrepancies from stereo disparity using sophisticated region processing, and then construct a trimap from the depth discrepancies. The actual matting is determined by the Bayesian matting of Chuang et al. However, that method has the view dependent problems that are unavoidable with stereo cameras, e.g., reflections, specular highlights, and occlusions. It is desired to avoid view dependent problems.

[0025] Matting is a process for extracting a high-quality alpha matte and foreground from an image or a video sequence.

[0026] Conventional techniques require either a known background, e.g., a blue screen, or extensive manual interaction, e.g., manually specified foreground and background regions.

[0027] Matting is generally under-constrained, because not enough information is obtained when the images are acquired.

[0028] The invention provides a system and method for extracting a matte automatically from images of rich, natural scenes illuminated only by ambient light.

[0029] The invention uses multiple synchronized cameras that are aligned on a single optical axis with a single center of projection. Each camera has the identical view of the scene, but a different depth of field. Alternatively, a single camera can be used to acquire images sequentially at different depths of field.

[0030] A first image or video, has the camera focused on the background, a second image or video has the camera focused on the foreground, and a third image or video is acquired by a pinhole camera so that the entire scene is in focus.

[0031] The images are analyzed according to Fourier image formation equations, which are over-constrained and share a single point of view but differ in their plane of focus. We minimize an error in the Fourier image equations.

[0032] The invention solves the fully dynamic matting problem without manual intervention. Both the foreground and background can have high frequency components and dynamic content. The foreground can resemble the background. The scene can be illuminated only by ambient light.

[0036] System Overview

[0037] FIGS. 1 and 2 shows a system 100 and method 200 according to our invention for automatically extracting a matte 141 from images acquired of a scene 110 including a background region (B) 111 having a background depth of field 131, and a foreground region (F) 112 having a foreground depth of field 132. These can be a natural, real word indoor or outdoor scene illuminated only by ambient light.

[0038] Cameras

[0039] The images are acquired 210 by a background camera 101, a foreground camera 102, and a pinhole camera (P) 103. The three cameras 101-103 are aligned on a single optical axis 160, sharing a single virtual center of projection, using first and second beam splitters 151-152. Therefore, all cameras have an identical point of view of the scene 110. The cameras are synchronized and connected to a processor 140.

[0040] The foreground and background cameras have relatively large apertures, resulting in small, non-overlapping depths of fields 131 and 132. That is, the depths of field are substantially disjoint. The pinhole camera has a very small aperture resulting in a large depth of field 133 with the entire scene in focus.

[0041] The foreground camera produces sharp images for the foreground region within about ½ meter of depth z_F of a foreground image plane 162 and defocuses regions farther away. The background camera produces sharp images for the background region with a background plane 161 at a depth z_B from about four meters to infinity and defocuses the foreground region, see FIG. 2. The pinhole camera is nominally focused on the foreground region. It should be noted that other depth of field setting can be used for the foreground and background cameras, depending on the structure of the scene.

[0042] Alternatively, a single camera can be used to acquire three images sequentially with the different aperture settings. This works for relatively static scenes, or for slowly varying scenes if the frame rate is relatively high or the exposure time is relatively short. In this, the camera is the foreground, background, and pinhole camera as the camera settings are changed in turn.

[0043] Our cameras respond linearly to incident radiance. We connect each camera to the processor 140 with a separate FireWire bus 142. The cameras acquire images at 30 frames per second. We equip each camera with a 50 mm lens 104. The pinhole camera is positioned after the first beam splitter 151. The aperture of the pinhole camera is f=12.

[0044] The pinhole camera 103 is focused on the foreground plane 162, because acquiring a correct matte is more important than correctly reconstructing the background. The foreground and background cameras have f=1.6 apertures and are positioned after the second beam splitter 152. Although, each camera receive only half the light of the pinhole camera, the relative large apertures acquire a relatively large amount of illumination. Therefore, the exposure for these two cameras 101-102 is shorter than the exposure for the pinhole camera. As long as the acquired images are not under-exposed or over-exposed, the color calibration process corrects remaining intensity differences between cameras.

[0045] Calibration

[0046] The cameras are calibrated to within a few pixels. Calibration is maintained by software. The optical axes are aligned to eliminate parallax between cameras. Because the focus is different for the different cameras, the acquired images are of different sizes. We correct for this with an affine transformation. We color correct the images by solving a similar problem in color space. Here, the feature points are the colors of an image of a color chart and the affine transformation is a color matrix. We apply color and position correction in real-time to all image sequences.

[0047] Image Sequences

[0048] For videos, each camera produces a 640×480×30 fps encoded image sequence. The sequences of images are processed by the processor 140 performing a matte extraction method 200 according to our invention.

[0049] Method Overview

[0050] FIG. 2 shows a method 200 for automatically extracting a matte according to the invention. Background, foreground, and pinhole sequence of images (videos) 201, 202, 203, respectively, 303 are acquired 310 of the scene 110 by the cameras 101-103. It should be understood that a single camera can be used as well, acquiring images sequentially at the appropriate different depths of field.

[0051] The pixels in each pinhole image are classified as either background, foreground, or unknown by matching neighborhoods around the pixel with corresponding neighborhood of pixels in the background and foreground images. The classification constructs 220 a trimap 221 for each pinhole image. An optimization process 230 is applied to the unknown pixels. The optimizer minimizes an error in classifying the unknown pixels as either background or foreground pixels. This produces the matte 104.

[0052] Scene Model

[0053] We model the scene 110 as a textured foreground plane 162 with partial coverage, and an opaque textured background plane 161. Because the background depth of field is larger than the foreground depth of field, and because there is no parallax between our cameras, the background region with varying depths can still be approximated as a plane for the purpose of matting.

[0054] We pose matting as an over-constrained optimization problem. For each pixel, there are the seven unknown “scene” values, α, F_{r,g,b}, and B_{r,g,b}, and nine constraint values I_P{r,g,b}, I_F{r,g,b}, and I_B{r,g,b} from the images I acquired by the cameras. The ‘P’ subscript denotes the pinhole images, the ‘F’ subscript the foreground-focused images, and the ‘B’ subscript the background-focused images.

[0055] Optimizer

[0056] We solve Fourier image formation equations by minimizing an error in classifying unknown pixels using the optimizer. To accelerate convergence for our optimizer, we construct 220 the trimaps 221 automatically using depth-from-defocus information, and select initial values that are likely near a true solution for the unknowns of the equations.

[0057] Initial foreground values F₀ for the optimizer are determined by automatically assigning known foreground colors to unknown regions. Initial background values B₀ are determined by reconstructing occluded areas from neighboring images, and then ‘painting’ into always occluded regions. Initial alpha coverage values α₀ are determined by solving a pinhole compositing equation using F₀ and B₀.

[0058] Defocus matting is poorly conditioned when the foreground and background have the same color, when the scene lacks high frequency components, or when the images are under-exposed or over-exposed. To avoid local minima and to stabilize the optimizer in these poorly conditioned areas, we add regularization terms to our optimizer.

[0059] The core of our optimizer 230 is the error function, which is invoked a few hundred times per image. Therefore, the challenge in solving the defocus matting by optimization is selecting an error function that is efficient to evaluate and easy to differentiate. Our error function is a sum-squared pixel value error between the acquired images and composite images rendered from the unknowns.

[0060] Evaluating and differentiating the error function naively make the problem intractable. To move towards a global minimum, the optimizer must find the gradient of the error function, i.e., the partial derivatives with respect to each unknown variable.

[0061] For a 320×240 pixel color image sequence at 30 fps, we need to solve for over 13 million unknowns per second. For instance, numerically evaluating the gradient invokes the error function once for each variable. For our method, this involves rendering three full-resolution images. A very fast ray tracer may be able to render the images in three seconds. That means a single call to the error function also takes three seconds. Therefore, it would take years to optimize a few seconds of video using conventional techniques.

[0062] Therefore, we approach the minimization as a graphics-specific problem. We symbolically manipulate expressions to avoid numerical computations. Thus, we provide a very fast approximation to the image synthesis problem, which enables us to evaluate the error function in milliseconds. We replace numerical evaluation of the error derivative with a symbolic derivative based on our synthesis equations, described below.

[0063] Notation

[0064] We use the following notation to compactly express discrete imaging operations. Monochrome images are 2D matrices that have matching dimensions. Image matrices are multiplied component wise, without a matrix multiplication. A multi-parameter image is sampled across camera parameters, such as, wavelength λ, focus, and time t, as well as pixel location.

[0065] We represent the multi-parameter image with a 3D or larger matrix, e.g., C[x, y, λ, z, t]. This notation and our matting method extend to images with more than three color samples and to other parameters, such as polarization, sub-pixel position, and exposure. Expressions, such as C[λ, z], where some parameters are missing, denote a sub-matrix containing elements corresponding to all possible values of the unspecified parameters, i.e., x, y, and t.

[0066] Generally, our equations have the same form in the x and y dimension, so we frequently omit the parameter y. We also omit the z, λ, and t parameters when these parameters do not for a particular equation.

[0067] A convolution F{circle around (×)}G of an image F and a matrix G has the same size as F. The convolution can be determined by extending edge values of F by half the size of G, so that F is well defined near the edges of F.

[0068] A disk(r)[x, y] is 1/πr² times the partial coverage of the pixel [x, y] by a disk of radius r centered on pixel [0, 0]. If the radius r<½, then the disk becomes a discrete impulse δ[x, y] that is one at [0, 0], and zero elsewhere.

[0069] Convolution with an impulse is the identity operation, and convolution with a disk is a ‘blur’ of the input image.

[0070] A vector ‘hat’ (→) above a variable denotes a multi-parameter image ‘unraveled’ into a column vector along its dimensions in order, e.g.,
{right arrow over (F)}[x+W((y−1)+H(λ−1))]=F[x, y, λ],
for an image with W×H pixels and 1-based indexing. This is equivalent to a raster scan order.

[0071] To distinguish the multi-parameter image vectors from image matrices, elements of the unraveled vectors are referenced by subscripts. Linear algebra operators, such as matrix-vector multiplication, inverse, and transpose operate normally on these vectors.

[0072] Defocus Composites

[0073] Equation 1 is the discrete compositing equation for a pinhole camera. We derive an approximate compositing equation for a camera with a non-zero aperture, which differs from a pinhole because some locations appear defocused. In computer graphics, cameras are traditionally simulated with distributed ray tracing.

[0074] Instead, we instead use Fourier optics, which are well suited to our image-based matting problem. Defocus occurs because the cone of rays from a point in the scene intersects the image plane at a disk called the point spread function (PSF).

[0075] FIG. 3 shows the optical geometry of the situation giving rise to a PSF with pixel radius

[0076] $\begin{array}{cc}r=\frac{f}{2\sigma \#}\frac{z_R\left(z_F-f\right)}{z_F\left(z_R-f\right)}-1,& \left(2\right)\end{array}$

[0077] where the camera is focused at depth z_F, a pixel at a depth z_R, # is the f-stop number, f is the focal length, and σ is the width of a pixel.

[0078] Depths z 300 are positive distances in front of the lens 104. A single plane of points perpendicular to the lens axis with pinhole image αF has a defocused lens image given by the convolution α{circle around (×)}F disk(r). Adding the background to the scene complicates matters because the background is partly occluded near foreground object borders.

[0079] Consider a bundle of rays emanating from a partly occluded background to the lens. The light transport along each ray is modulated by the α value, where the ray intersects the foreground plane. Instead of a cone of light reaching the lens from each background point, a cone cut by the image αF reaches the aperture. Therefore, the PSF varies for each point on the background. The PSF is zero for occluded points, a disk for unoccluded points, and a small cut-out of the a image for partly occluded points. We express the PSF values for the following cases.

[0080] Pinhole

[0081] When fσ is very small, or # is very large, r is less than half a pixel at both planes and Equation 1 holds.

[0082] Focused on Background

[0083] When the background 161 is in focus, the PSF is an impulse, i.e., a zero radius disk with finite integral. Rays in a cone from the background B are still modulated by a disk of (1−α) at the foreground plane, but that disk projects to a single pixel in the final image. Only the average value, and not the shape of the α disk intersected affects the final image. The composition equation is:
I_B=(αF){circle around (×)}disk(r_F)+(1−α{circle around (×)}disk(r_F))B.   (3)

[0084] Focused on Foreground

[0085] When the background is defocused and only the foreground is in focus, the PSF varies along the border of the foreground. Here, the correct image expression is complicated and slow to evaluate, therefore, we use the following approximation:
I_F≈αF+(1−α)(B{circle around (×)}disk(r_B)),   (4)
which blurs the background slightly at foreground borders.

[0086] A matte is a 2D matrix α[x, y], and the foreground and background images are respectively 3D matrices F[x, y, α] and B[x, y, α]. We generalize the two-plane compositing expression with a function of the scene that varies over two discrete spatial parameters, a discrete wavelength (color channel) parameter λ, and a discrete focus parameter zε{1, 2, 3 }: C(α,F,B)[x,y,λ,z]=
(αF[λ]){circle around (×)}h[z]+(1−α{circle around (×)}h[z])(B[λ]{circle around (×)}g[z])|_[x,y],   (5)

[0087] where 3D matrices h and g encode the PSFs:

[0088] $\begin{array}{cc}h\left[x,y,z\right]=\{\begin{array}{cc}\delta \left[x,y\right],& z=1\\ \mathrm{disk}\left(r_F\right)\left[x,y\right],& z=2\\ \delta \left[x,y\right],& z=3\end{array}g\left[x,y,z\right]=\{\begin{array}{cc}\delta \left[x,y\right],& z=1\\ \delta \left[x,y\right],& z=2\\ \mathrm{disk}\left(r_B\right)\left[x,y\right],& z=3\end{array}.& \left(6,7\right)\end{array}$

[0089] Constants r_F and r_B are the PSF radii for the foreground and background planes when the camera is focused on the opposite plane.

[0090] Trimap From Defocus

[0091] The trimap 221 segments the pinhole image into three mutually exclusive and collectively exhaustive regions expressed as sets of pixels. These sets of pixels limit the number of unknown pixels and provide initial estimates for our optimizer 230. In contrast with the prior art, we construct 220 the trimaps 221 automatically as follows.

[0092] Areas in the scene that have high-frequency texture produce high-frequency image content in the pinhole image I_P, and either in the foreground image I_F or the background image I_B, but not both. We use this observation to classify pixels into sets of pixels with high-frequency neighborhoods into three regions based on the z values, which appear ‘sharp’.

[0093] Sets ΩB and ΩF contain pixels that are respectively “definitely background” (α=0) and “definitely foreground” (α=1). Set Ω contains “unknown” pixels that may be either foreground, background, or some blend of foreground and background. This is the set over which we solve for extracting the matte using our optimizer.

[0094] Many surfaces with uniform macro appearance actually have fine structural elements like the pores and hair on human skin, the grain of wood, and the rough surface of brick. This allows us to detect defocus for many foreground objects even in the absence of strong macro texture. We use lower thresholds to detect high frequency components in the background, where only macro texture is visible.

[0095] We determine a first classification of the foreground and background regions by measuring a relative strength of the spatial gradients:
Let D=disk(max(r_F, r_B))
Ω_F1=erode(close((|∇I_F|>|∇I_B|){circle around (×)}D>0.6, D)),D)
Ω_B1=erode(close((|∇I_F|<|∇I_B|){circle around (×)}D>0.4, D)),D)   (8,9)
where erode and close are morphological operators used to improve accuracy. The disk is approximately the size of the PSFs. Then, we classify the ambiguous locations either in both ′Ω_F1 and ′Ω_B1 or in neither:
Ω={tilde over (Ω)}_F1∪{tilde over (Ω)}_B1∪(Ω_F1∩Ω_B1).   (10)

[0096] Finally, we enforce the mutual exclusion property:
Ω_F=Ω_F1∩{tilde over (Ω)}
Ω_B=Ω_B1∩{tilde over (Ω)}.   (11,12)

[0097] Minimization Errors in Classifying Unknown Pixels

[0098] We pose matting as an error minimization problem for each image, and solve the problem independently for each image. Assume we know the approximate depths of the foreground and background planes and all camera parameters. These are reasonable assumptions because digital cameras directly measure their parameters. From the lens to sensor distance we can derive the depths to the planes, if otherwise unknown.

[0099] The foreground and background need not be perfect planes, they just need to lie within the foreground and background depth fields. Because the depth of field is related hyperbolically to depth, the background depth field can stretch to infinity.

[0100] Let u=[{right arrow over (α)}^T{right arrow over (B)}^T{right arrow over (F)}^T]^T be the column vector describing the entire scene, i.e., the unknown pixels in the matting problem, and {right arrow over (C)}(u) be the unraveled composition function from Equation 5.

[0101] The unraveled constraints are {right arrow over (I)}=[{right arrow over (I)}_P^T{right arrow over (I)}_B^T{right arrow over (I)}_F^T]^T. The solution to the matting problem is a scene u* for which the norm of the error vector {right arrow over (E)}(u)={right arrow over (C)}(u)−{right arrow over (I)} is minimized according:

[0102] $\begin{array}{cc}\mathrm{Let}Q\left(u\right)=\sum _k\frac{1}{2}{\overrightarrow{E}}_k^2\left(u\right)u^{*}=\underset{u}{\arg \min }Q\left(u\right).& \left(13,14\right)\end{array}$

[0103] Note that the scalar-valued function Q is quadratic because the function Q contains the terms of the form (α[x]F[i])².

[0104] Iterative solvers appropriate for minimizing such a large system evaluate a given scene u and select a new scene u+Δu as a function of the vector {right arrow over (E)}(u), and a Jacobian matrix J(u). The Jacobian matrix contains the partial derivative of each element of the vector with respect to each element of u,

[0105] $\begin{array}{cc}{J}_{k,n}\left(u\right)=\frac{\partial {\overrightarrow{E}}_k\left(u\right)}{\partial u_n}.& \left(15\right)\end{array}$

[0106] The value k is an index into the unraveled constraints, and the value n is an index into the unraveled unknown array. Henceforth, we write {right arrow over (E)} rather than {right arrow over (E)}(u), and so on, for the other functions of u to simplify the notation in the presence of subscripts.

[0107] A gradient descent solver moves opposite the gradient of Q:

[0108] $\begin{array}{cc}\Delta u=-\nabla Q=-\nabla \sum _k\frac{1}{2}{\overrightarrow{E}}_k^2\mathrm{so}\Delta u_n=-\frac{\partial \sum _k\frac{1}{2}{\overrightarrow{E}}_k^2}{\partial u_n}=-\sum _k\left({\overrightarrow{E}}_k\frac{\partial {\overrightarrow{E}}_k}{\partial u_n}\right)\mathrm{hence}& \left(16,17\right)\\ \Delta u=-{\overrightarrow{E}}^{T}J.& \left(18\right)\end{array}$

[0109] The gradient descent solver has a space advantage over other methods like Gauss-Newton and Levenberg-Marquardt because the gradient function does not need to determine the pseudo-inverse of the Jacobian matrix J. This is important because the vectors and matrices involved are very large.

[0110] Let N be the number of unknown pixels and K be the number of constrained pixels. For 320×240 images, the matrix J has about 6×10⁹ elements.

[0111] We now derive a simple expression for the elements of the Jacobian matrix and determine that the matrix is sparse, so determining Δu is feasible when we do not need the non-sparse inverse of the matrix J. By definition, the elements are:

[0112] $\begin{array}{cc}{J}_{k,n}=\frac{\partial \left({\overrightarrow{C}}_k\left(u\right)-{\overrightarrow{I}}_k\right)}{\partial u_n}=\frac{\partial {\overrightarrow{C}}_k\left(u\right)}{\partial u_n}.& \left(19\right)\end{array}$

[0113] To evaluate Equation 19, we expand the convolution from Equation 5. We change variables from packed 1D vectors indexed by k to images indexed by

[0114] $\begin{array}{cc}C\left[x,z,\lambda \right]\sum _s\alpha \left[s\right]F\left[s,\lambda \right]h\left[x-s,z\right]+\left(1-\sum _s\alpha \left[s\right]h\left[x-s,z\right]\right)\sum _{s}B\left[s,\lambda \right]g\left[x-s,z\right].& \left(20\right)\end{array}$

[0115] An examination of this expansion shows that the matrix J is both sparse and simple. For example, consider the case where unknown pixel u_n corresponds to F[i,λ]. In a full expansion of Equation 20, only one term contains F[i, λ], so the partial derivative contains only one term:

[0116] $\begin{array}{cc}\frac{\partial C\left[x,\lambda ,z\right]}{\partial F\left[i,\lambda \right]}=\alpha \left[i\right]h\left[x-i,z\right].& \left(21\right)\end{array}$

[0117] The expressions for the α and B derivatives are only slightly more complicated, with potentially non-zero elements only at:

[0118] $\begin{array}{cc}\begin{array}{c}\frac{\partial C\left[x,\lambda ,z\right]}{\partial \alpha \left[i\right]}=h\left[x-i,z\right]\left(F\left[i,\lambda \right]-\sum _{s}B\left[\lambda ,s\right]g\left[x-s,z\right]\right)\\ =h\left[x-i,z\right]\left(F\left[i,\lambda \right]-\left(B\left[\lambda \right]\otimes g\left[z\right]\right)\left[x\right]\right)\end{array}\begin{array}{c}\frac{\partial C\left[x,\lambda ,z\right]}{\partial B\left[i,\lambda \right]}=g\left[x-i,z\right]\left(1-\sum _s\alpha \left[s\right]h\left[x-s,z\right]\right)\\ =g\left[x-i,z\right]\left(1-\left(\alpha \otimes h\left[z\right]\right)\left[x\right]\right)\end{array}.& \left(22,23\right)\end{array}$

[0119] The summations in the last two cases are just elements of convolution terms that appear in {right arrow over (E)}, so there is no additional cost for computing these values.

[0120] Trust Region and Weights

[0121] The gradient indicate a direction to change u to reduce the error. We use a so-called dogleg trust region scheme to select the magnitude, see Nocedal and Wright, IEEE PAMI 18, 12, pp. 1186-1198, Springer Verlag. The idea is to take the largest step that decreases the error. We begin with a trust region of radius S=1.

[0122] Let u′=max(0, min(1, u+(SΔu/|Δu|)). If |{right arrow over (E)}(u′)|<{right arrow over (E)}(u), then, we assume we have not overshot the minimum and repeatedly double S until the error increases above the lowest level seen this iteration. If |{right arrow over (E)}(u′)|<{right arrow over (E)}(u), then we assume we have overshot and take the opposite action, repeatedly halving S until we pass the lowest error in this iteration. If S becomes very small, e.g., 10⁻¹⁰ or the error norm decreases by less than 0.1%, then we assume that we are at the local minimum and terminate the optimization process.

[0123] Because our initial estimates are frequently good, we weigh the first N elements of Δu by constant β_α=3 to influence the optimizer to take larger steps in α. This speeds convergence without shifting the global minimum.

[0124] The narrow aperture and long exposure used to acquire the pinhole images produce more noise and motion blur than in the foreground and background images I_F and I_B. This prevents over-fitting the noise. This also reduces the over-representation in {right arrow over (E)} of in-focus pixels that occurs because image F and B are in focus in two of the constraint images and defocused in one each.

[0125] Regularization

[0126] In foreground regions that are low frequency or visually similar to the background, there are many values of u that satisfy the constraints. We bias the optimizer towards likely solutions. This is regularization of the optimization problem, which corresponds to having a different prior probability for a maximum likelihood problem. Regularization also to avoid local minima in the error function and stabilizes the optimizer in regions where the global minimum is in a ‘flat’ region that has many possible solutions.

[0127] We extend the error vector {right arrow over (E)} with p new entries, each entry corresponding to the magnitude of a 7N-component regularization vector. Calling these regularization vectors ε, φ, γ, . . . , the error function Q now has the form:

[0128] $\begin{array}{cc}\begin{array}{c}Q\left(u\right)=\sum _k{\overrightarrow{E}}_k^2\\ =\left[\sum _{k=1}^{9K}{\overrightarrow{E}}_k^2\right]+{\overrightarrow{E}}_{9K+1}^2+{\overrightarrow{E}}_{9K+2}^2+\dots \\ =\sum _{k=1}^{9K}{\overrightarrow{E}}_k^2+{\beta }_1\frac{9K}{7N}\sum _n^{7N}{\varepsilon }_n^2+{\beta }_2\frac{9K}{7N}\sum _n^{7N}{\varphi }_n^2+\dots \end{array}& \left(24\right)\end{array}$

[0129] The regularization vectors are e. Each summation over n appears as a new row in the error vector {right arrow over (E)} and the matrix J for some k>9K:

[0130] $\begin{array}{cc}{\stackrel{->}{E}}_k={\left(\beta \frac{9K}{7N}\sum _n{e}_n^2\right)}^{\frac{1}{2}}{J}_{k,n}=\frac{\partial {\stackrel{->}{E}}_k}{\partial u_n}=\frac{\beta }{{\stackrel{->}{E}}_k}\frac{9K}{7N}\sum _i\left[e_i\frac{\partial e_i}{\partial u_n}\right].& \left(25,26\right)\end{array}$

[0131] The factor

[0132] $\frac{9K}{7N}$
makes the regularization magnitude invariant to the ratio of constraints to unknown pixels, and the scaling factor β allows us to control its significance.

[0133] We select regularization vectors that are both easy to differentiate and efficient to evaluate, i.e., the summations over i generally contain only one non-zero term.

[0134] Regularization influences the optimizer to the most likely of many solutions supported by the image data, but rarely leads to an unsupported solution. We use small weights on the order of β=0.05 for each term to avoid shifting the global minimum.

[0135] Coherence

[0136] The spatial gradients are small,

[0137] $\begin{array}{cc}{e}_n=\frac{\partial u_n}{\partial x};{\left({\overrightarrow{E}}^{T}J\right)}_n=-\frac{{\partial }^2{u}_n}{\partial x^2}.T& \left(27\right)\end{array}$

[0138] We apply separate coherence terms to α, F, and B, for each color channel and for directions x and y. The alpha gradient constraints are relaxed at edges in the image. The F gradient constraints are increased by a factor of ten, where |∇α| is large. These constraints allow sharp foreground edges and prevent noise in the foreground image F where it is ill-defined.

[0139] Discrimination

[0140] The value α is distributed mostly at 0 and 1,
e_n=u_n−u_n²;({right arrow over (E)}^TJ)_n=(u_n−u_n²)(1−2u_n)|1≦n≦N.   (28)

[0141] Background Frequencies Should Appear in B:

[0142] $\begin{array}{cc}\mathrm{Let}G=I_B-I_F\otimes \mathrm{disk}\left(r_F\right)e_n=\frac{\partial u_n}{\partial x}-\frac{\partial {\stackrel{->}{G}}_n}{\partial x};{\left({\stackrel{->}{E}}^{T}J\right)}_n=-\frac{{\partial }^2{u}_n}{\partial x^2}❘4N+1\le 7N.& \left(29\right)\end{array}$

[0143] Other Applications

[0144] Artificial Depth of Field

[0145] We can matte a new foreground onto the reconstructed background, but select the point spread functions and transformations arbitrarily. This enables us to render images with virtual depth of field, and even slight translation and zoom.

[0146] Image Filtering

[0147] Defocus is not the only effect we can apply when recompositing against the original background image. Any filter can be used to process the foreground and background separately using the matte as a selection region, e.g., hue adjustment, painterly rendering, motion blur, or deblur.

[0148] Although the invention has been described by way of examples of preferred embodiments, 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.