Cast shadows and linear subspaces for object recognition

Abstract

The present invention is a method of deriving a reflectance function that analytically approximates the light reflected from an object model in terms of the spherical harmonic components of light. The reflectance function depends upon the intensity of light incident at each point on the model, but excludes light originating from below a local horizon, therefore not contributing to the reflectance because of the cast shadows. This reflectance function is used in the process of machine vision, by allowing a machine to optimize the reflectance function and arrive at an optimal rendered image of the object model, relative to an input image. Therefore, the recognition of an image produced under variable lighting conditions is more robust. The reflectance function of the present invention also has applicability in other fields, such as computer graphics.

Description

BACKGROUND OF THE INVENTION

[0002] 1. Field of Invention

[0003] The present invention relates generally to machine vision, and more particularly, to a method of modeling the shadows cast by various features of an object onto its own surface under certain lighting conditions, to allow a machine to recognize the image of that object when illuminated by a wide variety of lighting conditions.

[0004] 2. Description of Related Art

[0005] It is clear that changes in lighting can have a large effect on how an object, such as a person, looks. Understanding this is essential to building effective object recognition systems.

[0006] U.S. patent application Ser. No. 09/705,507 (Attorney Docket No. NECI 1093, 13919), filed Nov. 3, 2000, entitled “Lambertian Reflectance and Linear Subspaces”, the disclosure of which is hereby incorporated by reference, considered the relationship between the function that describes the lighting intensity, and the reflectance function that describes how much light an object reflects as a function of its surface normal, under a given lighting condition. Representing these functions as spherical harmonics, it was shown that for Lambertian reflectance the mapping from lighting to reflectance is a convolution with a nearly perfect low-pass filter. High-frequency components of the lighting hardly affect the reflectance function. Therefore, nearly all Lambertian reflectance functions could be modeled as some linear combination of nine spherical harmonic components.

[0007] These low-dimensional approximations have two key advantages. First, they are described with few parameters, thus optimization is simpler. Moreover, the images produced under low-frequency lighting are inherently immune to small changes in orientation. A small rotation in the surface normal is equivalent to a small rotation in the light. However, low-frequency light does not appreciably change with a small rotation. The analysis shows this, in that an image produced by an nth order harmonic light component is computed by taking an nth degree polynomial of the components of the surface normal.

[0008] U.S. patent application Ser. No. 09/967,206, incorporated by reference above, expanded the spherical harmonic approach to consider both the effects of broadened specular reflection as well as Lambertian reflection. Thus, the reflectance function that describes how much light an object reflects as a function of its surface normal, under a given lighting condition was made more realistic.

[0009] That work has considered the effects of attached shadows, which arise locally, when some portion of the object is facing away from the light. However, it has not addressed the problems that arise when one part of an object casts a shadow on a different part of the object. This is important, because many objects of interest are characterized by features whose prominence produces cast shadows over a significant portion of the surrounding surface, especially at low angles of illumination (high angles of incidence relative to the surface normal). Indeed, no natural object is convex, and relatively few manufactured objects are, and even fewer when seen in combination. Thus, cast shadows can be expected to play an important role in reflection, and, therefore must be included in harmonic treatments. Brief Summary of the Invention

[0010] The goal of the present invention is to allow a machine, such as a computer, to match an input image of an object, such as a human face, to an object model stored in a database associated with the machine, and to perform this matching with a high degree of accuracy.

[0011] The steps of the recognition method of the present invention include providing a database of three-dimensional object models, providing an input image, positioning each three dimensional model relative to the input image, determining an optimal rendered image of the model that is most similar to the input image, computing a measure of similarity between the input image and the optimal rendered image of each model, and selecting the object whose optimal rendered image is most similar to the input image. The step of determining an optimal rendered image includes deriving a reflectance function for each object that includes the effects of cast shadows, and optimizing that reflectance function to most closely match the input image.

[0012] The method of the present invention includes rendering the images produced by an object model when illuminated only by each of a plurality of spherical harmonic components of light. The summation of these components defines the reflectance function for the object model. The derived reflectance function takes into account the intensity of light components incident on the object model, as well as the portion of incident light which does not reach a given point because that light originates below some local horizon, that point therefore being in a region of cast shadow. The horizon is defined in terms of the polar angle at a given point, and can be a function of the azimuthal direction of the incident light.

[0013] Among the advantages of this approach is the ease of application and computational efficiency over serial methods of dealing with cast shadows, such serial methods being of the type including ray tracing, for example.

DETAILED DESCRIPTION OF THE INVENTION

[0014] Deriving the Reflectance Function

[0015] To derive a reflection function that accounts for some obstruction between the surface point of interest and the light source, consider a small surface region of an object model with surface normal {circumflex over (n)}. Let μ=cosθ, where θ is the polar angle of a ray of light striking the region relative to {circumflex over (n)}. Let φ be the azimuth of the ray, also relative to {circumflex over (n)}. Now if the object is convex, all rays with 0≦θ≦π/2(μ≧0) and 0≦φ≦2π will illuminate the surface. However, in general, only rays whose polar angles satisfy μ(φ)≦μ≦1, 0≦θ≦θ(φ)will contribute to the illumination, and hence to the reflection. Here θ(φ) is just the angle (relative to the normal) that outlines the local horizon as one sweeps through the azimuth. In the usual treatments one takes θ(φ)=π/2, for if θ(φ)˜π/2, so little light intensity per unit area is available that, unless the incoming light is low and highly directional, little error is made. However, for θ(φ) much less than π/2, some compensation must be made for the resulting shadows.

[0016] Three models of cast shadows are considered in the context of the harmonic-decomposition approach. The most general model takes θ(φ) and determines the surface intensity for each Llm lighting component, where L({circumflex over (n)}i)=ΣlmLlmYkn({circumflex over (n)}i) is the illumination on the object. This involves considerable calculation, but will be necessary in many realistic examples.

[0017] The simplest model, the crater model, is to assume that μ(φ)=μ0. In other words, that the region of interest sits in a depression with symmetry about the azimuth. While it is not clear, without more computation, what effective value to choose for μ0 given μ(φ), for μ0 >0 clearly the reduction in incident intensity and shift to higher harmonics can at least be simulated. These features are characteristic of shadows and are impossible to incorporate treating objects as convex. Sketching the derivation of the reflection in this model following the lines and notation of application Ser. No. 09/967,206, incorporated by reference above, one can then consider a somewhat more general model, and finally proceed to the most general case.

[0018] If rCSCX({circumflex over (n)}) is the reflection in the presence of cast-shadows in the crater approximation for a surface normal of {circumflex over (n)}, then

r CSCX ({circumflex over (n)})=u({circumflex over (z)}·{circumflex over (n)})∫d{circumflex over (n)}iL({circumflex over (n)}i)({circumflex over (n)}i·{circumflex over (n)})u({circumflex over (n)}i·{circumflex over (n)}−μ0)

[0019] where L({circumflex over (n)}i) is the incident light intensity in the −{circumflex over (n)}i direction and μ0 defines the effective local horizon of the crater. In the 09/967,206 Application μ0 was taken as zero. Also, u denotes the step function that is one for positive values and zero otherwise. Therefore, the step function including μ as shown in the above equation analytically excludes light originating from any point below the local horizon at the point of interest from contributing to the reflectance.

[0020] Focusing on the Lambertian factor one can find 1$\left({\hat{n}}_i·\hat{n}\right)u\left({\hat{n}}_i·\hat{n}-{\mu }_0\right)=\sum _{\mathrm{lm}}{\left(4\pi /\left(2l+1\right)\right)}^{1/2}{k}_{\mu l}{Y}_{\mathrm{lm}}\left(\hat{n}\right)Y_{\mathrm{lm}}^{*}\left({\hat{n}}_i\right)$

[0021] where Ylm ({circumflex over (n)}) is the lm spherical harmonic function of the direction {circumflex over (n)}, and where 2$k_{\mu l}={\int }_0^{2\pi }d\phi {\int }_{{\mu }_0}^{1}d\mathrm{\mu \mu }{Y}_{10}\left(\theta ,\phi \right)$

[0022] Since

Y l0 (θ,φ)=((2l+1)/4π)1/2Pl(μ)

[0023] 3 $$ k μ ⁢   ⁢ l = ( ( 2 ⁢ l + 1 ) ⁢ π ) 1 / 2 ⁢ ∫ μ 0 1 ⁢ ⅆ μμ ⁢   ⁢ P 1 ⁡ ( μ )

P l (μ)=(2l!)−1dl(μ2−1)l/dμl

[0024] evaluating one can find

k μ0 =π1/2(1−μ02)/2

k μl=π 1/2 (1−μ03)/31/2

[0025] 4 $$ k μ ⁢   ⁢ l = ( ( 2 ⁢ l + 1 ) ⁢ π ) 1 / 2 ⁢ ( ( l - 2 ) ! ⁢ ( - 1 ) ( l - 2 ) / 2 2 l ⁢ ( l 2 - 1 ) ! ⁢ ( l 2 + 1 ) ! + ∑ k = 1 / 2 l ⁢ ( ( 2 ⁢ k ) ! ⁢ ( 2 ⁢ k - ( l - 2 ) - 1 ) ⁢ ( - 1 ) l - k - 1 2 l ⁢ k ! ⁢ ( l - k ) ! ⁢ ( 2 ⁢ k - ( l - 2 ) ) ! ⁢ μ 0 2 ⁢ k - ( l - 2 ) )

[0026] for l≧2, l even, and 5$k_{\mu l}={\left(\left(2l+1\right)\pi \right)}^{1/2}\sum _{k=\left(l+1\right)/2}^l\frac{\left(2k\right)!\left(2k-\left(l-2\right)-1\right){\left(-1\right)}^{l-k-1}}{2^{l}k!\left(2k-\left(l-2\right)\right)!}{\mu }_0^{2k-\left(l-2\right)}$

[0027] for l≧3, l odd. For μ0 =0 the previous results are recovered. For μ0 >0, there is no first-order (in μ0) change in kμl. Second-order changes occur in kμ0 and kμl for 1 even: kμ2 increases as −((5π)1/2/4)μ02 for μ0 small.

Table of kμl
L\μ0	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
0	.886	.887	.851	.806	.744	.665	.567	.452	.319	.168	0
1	1.023	1.022	1.015	.996	.958	.895	.802	.672	.499	.277	0
2	.495	.505	.532	.572	.615	.650	.659	.624	.520	.322	0
3	0	.002	.018	.058	.126	.220	.324	.410	.432	.325	0
4	−.111	−.120	−.143	−.163	−.159	−.109	−.005	.141	.276	.292	0
5	0	−.004	−.026	−.102	−.141	−.189	−.179	−.076	.103	.233	0
6	.050	.051	.074	.069	.020	−.047	−.165	−.221	−.043	.159	0
7	0	.005	.032	.078	.108	.073	−.038	−.153	−.131	.081	0
8	−.029	−.037	−.044	−.016	.054	.114	.083	−.055	−.151	.008	0
9	0	−.006	−.035	−.066	−.047	.041	.116	.052	−.111	−.049	0
10	.019	.026	.024	−.019	−.074	−.054	.057	.109	−.039	−.086	0
11	0	.007	.035	.044	−.014	−.082	−.034	.094	.033	−.101	0
12	−.013	−.019	−.009	.039	.052	−.028	−.084	.030	.082	−.095	0
13	0	−.007	−.031	−.017	.049	.044	−.064	−.041	.093	−.073	0
14	.010	.015	.033	−.041	−.018	.058	.002	−.079	.068	−.049	0

[0028] In the foregoing, the kμl coefficients are shown for the harmonic expansion of the light reflected from a Lambertian surface in the presence of cast shadows in the crater approximation, where the elevation of the local horizon μ(φ)=cos(θ(φ)) is taken to be independent of azimuth, μ(φ)=μ0≧0. In a more general model, the full azimuthal dependence μ(φ) is included, but the Lambertian factor is assumed to have the same form as in the crater model. While this restricts the actual dependence of the reflection on the surface normal, it provides an indication of the role played by higher harmonics.

[0029] Replacing μ0 with μ(φ), the solution is straightforward to obtain from the crater result. Consider kμl, and treat it as a functional of μ(φ) rather than as a function of μ0. Thus 6$k_{\mu l}={\int }_0^{2\pi }d\phi {\int }_{\mu \left(\phi \right)}^{1}d\mathrm{\mu \mu }{Y}_{10}\left(\theta ,\phi \right),$

[0030] where μ=cosθ. For μ(φ)=μ0, kμl, is equal to 2π times a polynomial of degree (l+2) in μ0. For general μ(φ), replace μ0n in the expression for kμl by 7${\mu }_0^n\to {\int }_0^{2\pi }d{\mathrm{\phi \mu }}^n\left(\phi \right)/2\pi$

[0031] If lL is the value of the index l of the largest harmonic of interest, we need only consider n such that 2≦n≦(lL+2):lL=1 for the four-dimensional model (i.e. the 0th and 1st order spherical harmonic components); lL=2 for the nine-dimensional model (i.e., the 0th, 1st and 2nd order spherical harmonic components). Note that the n=0 case is trivial, and that the n=1 case does not enter.

[0032] Even though the table given for kμl assumes μ(φ)=μ0, it is actually useful in the more general case. To see this let 8$\stackrel{\_}{{\mu }^n}={\int }_0^{2\pi }d\phi {\mu }^n\left(\phi \right)/2\pi$

[0033] Then, since {overscore (μn)}≧{overscore (μ)}n for n≧1 note that the actual values of kμl for μ(φ) are indicative of the values given for kμl for μ(φ)=μ0 , where {overscore (μ)}<μ0. While admittedly rough, it is at once apparent that the dependence of μ on φ enhances the harmonic content of the image. Also, since {overscore (μ)}>{overscore (μ2)}>{overscore (μ3)}. . . for most cases of interest, only relatively low values of n (powers of μ) will play a significant role in the kμl.

[0034] To include the full effect of cast shadows in the harmonic representation, one must fully express the expansion of the Lambertian factor. In both the crater model and the more general model an assumption was made to ignore the azimuthal dependence in the unit step function μ( ), the φ dependence of μ(φ). While satisfactory in those cases where μ(φ) is nearly constant, in non-convex regions its variation with azimuth can be significant.

[0035] Though tricky to evaluate, the solution can be written down by analogy with the more general of the two cases evaluated above. If rCSM({circumflex over (n)}) is the reflection in the presence of cast shadows due to a local horizon of μ(φ)=cos(θ(φ)) that is due to an elevation dependent on azimuth φ, then we can write much as above that

r CS ({circumflex over (n)})=u({circumflex over (z)}·{circumflex over (n)})d{circumflex over (n)}iL({circumflex over (n)})({circumflex over (n)}i·{circumflex over (n)})u({circumflex over (n)}i·{circumflex over (n)}−μ(φ))

[0036] 9 $$ ( n ^ i · n ^ ) ⁢ u ⁡ ( n ^ i · n ^ - μ ⁡ ( φ ) ) =   ⁢ ∑ lmm ′ ⁢ k lm ⁢ D m ′ ⁢ m l * ⁡ ( n ^ ) ⁢ Y lm ′ * ⁡ ( n ^ i ) k lm =   ⁢ ∫ 0 2 ⁢ π ⁢ ⅆ φ ⁢ ∫ μ ⁡ ( φ ) 1 ⁢ ⅆ μμ ⁢   ⁢ Y lm ⁡ ( θ , φ ) , μ = cos ⁢   ⁢ θ

[0037] Here Dm′ml(αβγ) is the rotation matrix associated with spherical harmonics, and (β,α) are the polar coordinates of {circumflex over (n)} (γ=0). Angles (αβγ) are the Euler angles of the effective rotation that takes the {circumflex over (z)}-axis from the viewing direction to the direction of the surface normal {circumflex over (n)}.

[0038] The evaluation of klm is complicated by the azimuthal dependence of the elevation μ(φ). One simplification that retains much of the azimuthal or phase coupling introduced by cast shadows is to assume that we can approximate μ(φ) by two values, μ1 for φ1<φ≦φ2 and μ2 for φ2<φ≦φ1+2π. Then the integrals over dφ and dμ separate, and one obtains factors of sin(m(φ2−φ1)/2)/m entering the strength of the contribution to klm from the μ(φ)=φ1 and μ(φ)=φ2 regions. Phases m=pπ/(φ2−φ1), p odd, yield the largest klm. Thus for m=0 10$k_{10}\approx \left({\phi }_2-{\phi }_1\right){\int }_{{\mu }_1}^{{\mu }_2}d\mathrm{\mu \mu }{Y}_{10}\left(\theta ,\phi \right)+2\pi {\int }_{{\mu }_2}^{1}d\mathrm{\mu \mu }{Y}_{10}\left(\theta ,\phi \right),$

[0039] and for m>0 or m<0 11$k_{\mathrm{lm}}\approx e^{im\stackrel{\_}{\phi }}\sin \left(m{\phi }_{21}/2\right)/\left(m/2\right)·{\int }_{{\mu }_1}^{{\mu }_2}d\mathrm{\mu \mu }{\stackrel{\_}{Y}}_{\mathrm{lm}}\left(\theta \right)$

[0040] where {overscore (φ)}=(φ2+φ1)/2,φ21=(φ2−φ2),{overscore (Y)}lm(θ)=e−imφYlm(θ,φ), which is independent of φ.

[0041] Another simplification arises when the local horizon is axially symmetric about an axis {circumflex over (m)} other than the surface normal {circumflex over (n)}; we can call this the megaphone approximation and note its potential utility when only {circumflex over (n)} and {circumflex over (m)}, but not the megaphone aperture μm, vary with surface position. We write

r CSM ({circumflex over (n)}; {circumflex over (m)}, μm)=∫d{circumflex over (n)}iL({circumflex over (n)}i)({circumflex over (n)}i·{circumflex over (n)}ν0)u({circumflex over (n)}i·{circumflex over (m)}−μm)

[0042] (the expression AνB is equal to the greater of (A, B)) and express the three factors in the integrand as 12$L\left({\hat{n}}_i\right)=\sum _{\mathrm{lm}}{L}_{\mathrm{lm}}{Y}_{\mathrm{lm}}\left({\hat{n}}_i\right)$$\left({\hat{n}}_i·\hat{n}\bigvee 0\right)=\sum _{l_1{m}_1}{k}_{l_1}{Y}_{l_1{m}_1}\left(\hat{n}\right)Y_{l_1{m}_1}^{*}\left({\hat{n}}_i\right)$$u\left({\hat{n}}_i·\hat{m}-{\mu }_m\right)=\sum _{l_2{m}_2}{p}_{l_2}{Y}_{l_2{m}_2}\left(\hat{m}\right)Y_{l_2{m}_2}^{*}\left({\hat{n}}_i\right)$

[0043] where 13$k_l={\int }_0^{1}d\mu \mu Y_{\mathrm{l0}}\left(\mu \right)\text{and}{p}_l={\int }_{{\mu }_m}^{1}d\mu Y_{\mathrm{l0}}\left(\mu \right).$

[0044] By this point the analogy with diffuse-specular reflection should be clear. The term Pl has replaced kl(P), but otherwise the solution to that problem can be taken over intact. Thus one finds 14$r_{\mathrm{CSM}}\left(\hat{n};\hat{m},{\mu }_m\right)=\sum _{\mathrm{lm}}{L}_{\mathrm{lm}}\sum _{l_1{m}_1{l}_2{m}_2}I\left(l_1{l}_{2}l;m_1,m_2,m\right)k_{l_1}{p}_{l_2}{Y}_{l_1{m}_1}\left(\hat{m}\right)Y_{l_2{m}_2}\left(\hat{n}\right)$

[0045] where 15$\begin{array}{c}I\left(l_1{l}_{2}l;m_1,m_2,m\right)=\int d{\hat{n}}_i{Y}_{\mathrm{lm}}\left({\hat{n}}_i\right)Y_{l_1{m}_1}^{*}\left({\hat{n}}_i\right)Y_{l_1{m}_2}^{*}\left({\hat{n}}_i\right)\\ =(\left(2l_1+1\right)\left(2l_2=1\right)/\\ {4\pi \left(2l+1\right))}^{1/2}C\left(l_1{l}_{2}l;m_1,m_2,m\right)C\left(l_1{l}_{2}l;000\right)\end{array}$

[0046] where C(l1l2l;m1m2m) are the Clebsch-Gordon(C-G) coefficients. When applicable, this approach is significantly more computable than the general solution above, increasing the speed of computation and/or maximizing the use of processing resources.

[0047] The foregoing expressions for the modification of the reflected intensity in the presence of cast shadows requires that the horizon μ(φ) be determined for each surface position of interest. This requires two steps. First, the surface must be rotated so that the normal {circumflex over (n)} at this position is parallel to the {circumflex over (z)} axis. Then the maximum elevation μ(φ) must be found in each direction of interest φ. (μ(φ)=cosθ(φ) where θ(φ)is the minimum polar angle from the vertical in the direction φ.) The assumption is made that a set of surface positions {r}={(rx,ry,rz)} is given in some known frame of reference.

[0048] The first task is to rotate each surface position r into r′ such that the normal {circumflex over (n)} at the position of interest r0 becomes parallel to {circumflex over (z)} at r′0. This is accomplished by letting r′=Rr, where R is given by 16$R=\hat{z}{\hat{n}}^T+\left(\stackrel{•}{n}\times z\right){\left(\stackrel{•}{n}\times z\right)}^T+\left(\hat{z}\times \left(\stackrel{•}{n}\times z\right)\right){\left(\hat{n}\times \left(\stackrel{•}{n}\times z\right)\right)}^T$

[0049] Note especially that R{circumflex over (n)}={circumflex over (z)}. The three terms in R operate on the component of r parallel to {circumflex over (n)}, to {circumflex over (n)}×{circumflex over (z)}, and to {circumflex over (n)}×({circumflex over (n)}×{circumflex over (z)}), respectively. With this preliminary step out of the way, one can determine the horizon μ(φ)for position r0.

[0050] It is now simplest to pass to a set of variables relative to r′0=Rr0:

r′−r′ 0 =(r′x −r′0x,r′0y,r′z−r′0z)=(dx,dy,h)

d x =(dx2+dy2)1/2 cos φ

d y =(dx2+dy2)1/2 sin φ

[0051] The (dx,dy) give the position of the height h(dx,dy) of the surface above the plane through r0 normal to {circumflex over (n)}. It follows that μ(φ) is given by one of the following:

max(1+dx2/(hcosφ)2)−1/2

dx

max(1+dy2/(hsinφ)2)−1/2

0(zero)

[0052] where only dx or dy where h(dx,dxtanφ) or h(dycotφ,dy), respectively exceed zero enter into the maximum. For h≦0 for all positions along φ,μ=0. Clearly, one will choose to use the first expression for −π/4≦φ≦π/4 and for 3π/4≦φ≦5π/4, and the second for π/4≦φ≦3π/4 and for 5π/4≦φ≦7π/4.

[0053] As an example, see the following illustration computing μ(φ) for a ridge of isosceles-triangular cross section of height h>0 and base b, and unspecified length —somewhat like the nose on a face. At a distance d>b/2 to the left of the centerline of the base of the ridge, taken in the direction φ=±π/2,μ(φ)=(1+(d/hcosφ)2)−1/2 for −π/2≦φ≦π/2 and μ(φ)=0 for π/2≦φ≦3π/2. On the left side of the ridge itself μ(φ)=(1+(b/2h)2)−1/2 for π/2≦φ≦3π/2 and μ(φ)=0 for −π/2≦φ≦π/2. The latter is an example of μ(φ) constant over ranges of φ, whose solution was given above. The former can be approximated over appropriate ranges of φ, or integrated over φ if greater precision is required. Very similar solutions apply for a trapezoidal ridge, with μ(φ)=0 along the top of the ridge.

[0054] The crater model is sufficiently simple to introduce most of the effects due to cast shadows. Odd orders (1=3,5,7, . . . ), which play no role in the case of convex, Lambertian surfaces, become appreciable as the crater deepens. The strength of orders l≧4 oscilate in value as the fit for the various modes waxes and wanes with increasing μ0. From the general solution, note the couplings between the azimuthal frequencies in the horizons and those in the resulting expansions. The actual calculation of the horizon was indicated only briefly, as it has probably been more fully developed in the graphics community.

[0055] The rotation matrices used in the method are discussed by Rose, M. E., Elementary Theory of Angular Momentum. They enable the spherical harmonics in a rotated coordinate system to be expressed in terms of those in the reference system. Thus, 17$D_{m^{\prime }m}^l\left(\mathrm{\alpha \beta \gamma }\right)=\sum _t\frac{{\left(-1\right)}^t{\left(\left(l+m\right)!\left(l-m\right)!\left(l+m^{\prime }\right)!\left(l-m^{\prime }\right)!\right)}^{1/2}}{t!\left(l-t-m^{\prime }\right)!\left(l+m-t\right)!\left(m^{\prime }+t-m\right)!}·e^{-i{m}^{\prime }\alpha }{\left(e^{-im\gamma }\left(\cos \left(\beta /2\right)\right)\right)}^{2l+m-m^{\prime }-2t}{\left(-\sin \left(\beta /2\right)\right)}^{m^{\prime }-m+2t}$

[0056] In the summation over t, t is such that the four factorials in the denominator are well-defined. Thus each of t≧0, l−m′≧t, l+m≧t, t≧m−m′, must be satisfied for each possible t. Also, l≧|m| and l≧|m′|.

[0057] Clearly D000(αβγ)=1.

[0058] Tables for l=1 are commonplace:

	m = 1	m = 0	m = −1
	m′ = 1	e−iaCe−iγ	−e−ia2−1/2s	e−iaSeiγ
	m′ = 0	2−1/2se−iγ	c	−2−1/2seiγ
	m′ = −1	eiaSe−iγ	eia2−1/2s	eiaCeiγ

[0059] where C=(1+cosβ)/2, S=(1−sinβ)/2, c=cos β, s=sinβ. Values for l=2 are also needed, which have been worked out in the following table. The e−im′α and e−mγ factors are omitted the clarity of the table although, of course, these must be included in any solution.

	m = 2	m = 1	m = 0	m = −1	m = −2
m′ = 2	C2	−Cs	(3/8)1/2s2	−Ss	S2
m′ = 1	Cs	C(2c − 1)	−(3/2)1/2sc	S(2c + 1)	−Ss
m′ = 0	(3/8)1/2s2	(3/2)1/2sc	(3c2 − 1)/2	−(3/2)1/2sc	(3/8)1/2s2
m′ = −1	Ss	S(2c + 1)	(3/2)1/2sc	C(2c − 1)	−Cs
m′ = −2	S2	Ss	(3/8)1/2s2	Cs	C2

[0060] Further, a number of indefinite integrals enter in connection with the integration over the elevation (μ=cosθ) of the spherical harmonics. The following table should be self-explanatory:

	∫Ylm(μ)dμ	∫μYlm(μ)dμ
Y00˜1	μ	μ2/2
Y10˜μ	μ2/2	μ3/3
Y1±1˜(1 − μ2)1/2	(μ(1 − μ2)1/2 − cos−1 μ)/2	−(1 − μ2)3/2/3
Y20˜(3μ2 − 1)	−μ(1 − μ2)	μ2(3μ2/2 − 1)/2
Y2±1˜μ(1 − μ2)1/2	−(1 − μ2)3/2/3	((1 − μ2)1/2μ(2μ2 − 1) −
		cos−1 μ)/8
Y2±2˜(1 − μ2)	μ(1 − μ2)/3	μ2(1 − μ2/2)/2

[0061] Application of the Reflectance Function

[0062] The method of deriving the reflectance function described above is used in the object recognition process. To match an input image with an object model in a database, the object model must first be positioned to match the position of the input image. There are standard methods of performing this function known in the art, though a preferred method is disclosed in U.S. patent application Ser. No. 09/538,204 (Attorney Docket No. NECI 1083, 13414), filed Mar. 30, 2000, entitled “Method for Matching a Two Dimensional Image to One of a Plurality of Three Dimensional Candidate Models in a Database”, the disclosure of which is hereby incorporated by reference.

[0063] Once positioned, images of the model are rendered under each of one or more spherical harmonic components of light, based upon the reflectance function described above. These component rendered images can be linearly combined to arrive at an optimal rendered image of the model. An optimal rendered image means a rendered image that most closely approximates the input image. The comparison of the various rendered images to the input image is made by calculating the Euclidian Distance between the images, in a manner that is known in the art. The optimization is performed by a least squares fit technique between the vector that is the reflectance function and the matrix that represents the input image, also known in the art.

[0064] Other steps follow the optimization in the recognition process. One is computing some measure of similarity between the input image and the optimal rendered image for each model. This is accomplished by measuring the Euclidian Distance between each optimal image and the input image. The object model whose optimal rendered image is nearest the input image, assuming some threshold value of similarity must be reached, is a match to the input image.

[0065] While the reflectance function of the present invention is described in application to the field of machine vision, it is not limited to that field. There are numerous parallels between this work and advances in the field of computer graphics, and our invention has applicability there as well. The field of computer graphics is continually striving to produce more realistic rendered images based on object models that exist only in the computer. The reflectance function of the present invention advances that work as well.

[0066] The invention has been described herein with reference to a particular exemplary embodiment. Certain alterations and modifications may be apparent to those skilled in the art, without departing from the scope of the invention. The exemplary embodiment is not meant to be limiting on the scope of the invention, which is defined by the appended claims. claims

Claims

1. A method of recognizing objects under various lighting conditions comprising the steps of: (a) providing a database comprising a plurality of three dimensional models; (b) providing an input image; (c) positioning each three dimensional model relative to the input image; (d) determining, for each three dimensional model, an optimal rendered image which is most similar to the input image, said determining step comprising: (i) deriving a reflectance function that describes an approximation of the set of all possible rendered images that each three dimensional model can produce under all possible lighting conditions, wherein said reflectance function analytically excludes illumination within cast shadows on said model; and (ii) optimizing the reflectance function to determine an optimal rendered image of each model that is most similar to the input image; (e) computing a measure of similarity between the input image and each optimal rendered image; and (f) selecting the three dimensional model corresponding to the optimal rendered image whose measure of similarity is most similar to the input image.

2. The method according to claim 1 wherein the reflectance function includes a local horizon for each point on each said model, wherein no light originating below said local horizon contributes to the reflectance from said point.

3. The method according to claim 2 wherein said local horizon is defined by the polar angle of said local horizion in any direction.

4. The method according to claim 2 wherein said local horizon is variable with the azimuth at each point on said model.

5. A method of deriving a reflectance function of an object model under a variety of lighting conditions, the object model comprising a plurality of surfaces, each of the surfaces defining a normal vector pointing perpendicularly outward, the method comprising the steps of: (a) determining a local horizon for each point on said object model, wherein no light originating from below said local horizon is incident at said point; (b) rendering, for a given orientation of the object model, a plurality of images produced by the object model when illuminated by each of a plurality of spherical harmonic components of incident light, wherein said rendering comprises: (i) calculating the intensity of incident light components upon the object model, relative to the normal at each surface, excluding light originating from below a local horizon at any point on said model; and (ii) calculating the intensity of light components reflected by the model toward the observer; and (c) defining a reflectance function for the object as a linear combination of said plurality of images.

6. The method according to claim 5, wherein said plurality of images comprises the images produced by the 0th, 1st, and 2nd order spherical harmonic components of incident light.

7. The method according to claim 5, wherein the step of determining the local horizon for each point on said object model comprises determining the local horizon as a function of the polar angle of the local horizon at each said point.

8. The method according to claim 7, wherein the function measured by the polar angle is dependent on the azimuth.