Lecture 10 Notes: Characteristic Strip Expansion, Shape from Shading, Iterative Solutions

6.801/6.866: Machine Vision, Lecture 10 1 Lecture 10: Characteristic Strip Expansion, Shape from Shading, Iterative Solutions In this lecture, we will continue studying our shape from shading problem, this time generalizing it beyond Hapke surfaces, making this framework amenable for a range of machine vision tasks. We will introduce a set of ordinary differential equations, and illustrate how we can construct a linear dynamical system that allows us to iteratively estimate the depth profile of an object using both position and surface orientation. 1.1 Review: Where We Are and Shape From Shading Shape from Shading (SfS) is one of the class of problems we study in this course, and is part of a larger class of problems focusing on machine vision from image projection. In the previous lecture, we solved shape from shading for Hapke surfaces using the Image Irradiance Equation E(x, y) = R(p, q), where R(p, q) is the reflectance map that builds on the Bidirectional Reflectance Distribution Function (BRDF). Recall additionally that irradiance (brightness in the image) depends on: 1. Illumination 2. Surface material of object being imaged 3. Geometry of the object being imaged Recall from the previous lecture that for Hapke surfaces, we were able to solve the SfS problem in a particular direction, which we could then determine object surface orientation along. We can integrate along this direction (the profile direction) to find height z, but recall that we gain no information from the other direction. For the Hapke example, we have a rotated coordinate system governed by the following ODEs (note that xi will be the variable we parameterize our profiles along): 1. X-direction: dx dξ = ps 2. Y-direction: dy dξ = qs 3. By chain rule, Z-direction: dz dξ = ∂z dx ∂x dξ + ∂z dy ∂y dξ = pps + qqs = ps p + qs q Intuition: Infinitesimal steps in the image given by ξ gives which can be used for recovering surface orientation. dx dy dξ , dξ , and we are interested in finding the change in height dz dξ , Note: When dealing with brightness problems, e.g. SfS, we have implicitly shifted to orthographic projection (x = Zf0 X, y = f Z0 y). These methods can be applied to perspective projection as well, but the mathematics makes the intuition less clear. We can model orthographic projection by having a telecentric lens, which effectively places the object really far away from the image plane. 1 We can solve the 3 ODEs above using a forward Euler method with a given step size. For small steps, this approach will be accurate enough, and accuracy here is not too important anyway since we will have noise in our brightness measurements. Recall brightness of Hapke surfaces is given by: r n̂ · ŝ = n̂ · v̂ r cos θi = cos θe p 1 f rac1 + ps p + qs q 1 + p2 + q 2 rs · √ 1 + p s p + qs q E= Some cancellation and rearranging yields: ps p + qs q = rs E 2 − 1 Therefore, we have a direct relationship between our measured brightness E and our unknowns of interest: ∂z ∂z qs = pps + qqs = rs E 2 − 1 ps + ∂x ∂y Note here, however, that we do not know surface orientation based on this, since again for Hapke surfaces, we only know slope in one of our two directions in our rotated gradient space (p0 , q 0 ). A forward Euler approach will generate a set of independent profiles, and we do not have any information about the surface orientation at the interfaces of these independent profiles. This necessitates an initial curve containing initial conditions for each of these independent profiles. In 3D, this initial curve is parameterized by η, and is given by: (x(η), y(η), z(η)). Our goal is to find z(η, ξ), where: • η parameterizes the length of the initial curve • ξ parameterizes along the 1D profiles we carve out Therefore, we are able to determine the slope in a particular direction. 1.2 General Case Let us now generalize our Hapke case above to be any surface. We begin with our Image Irradiance Equation: E(x, y) = R(p, q) Let us now consider taking a small step δx, δy in the image, and our goal is to determine how z changes from this small step. We can do so using differentials and surface orientation: δz = ∂z ∂z δx + δy = pδx + qδy ∂x ∂y Therefore, if we know p and q, we can compute δz for a given δx and δy. But in addition to updating z using the equation above, we will also need to update p and q (intuitively, since we are still moving, the surface orientation can change): ∂p δx + ∂x ∂q δq = δx + ∂x δp = ∂p δy ∂y ∂q δy ∂y This notion of updating p and q provides motivation for keeping track/updating not only (x, y, z), but also p and q. We can think of solving our problem as constructing a characteristic strip of the ODEs above, composed of (x, y, z, p, q) ∈ R5 . In vector-matrix form, our updates to p and q become: ⎡ 2 ⎤ ⎡ ⎤ ⎡ ⎤ ∂ z ∂2z       r s   px py   2 ∂y∂x ∂x δp δx δx δx δx ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ = =⎣ =H ⎦ δy = δq δy δy δy ∂2z ∂2z s t q x qy 2 ∂x∂y ∂y Where H is the Hessian of second spatial derivatives (x and y) of z. Note that from Fubuni’s theorem and from our Taylor expansion from last lecture, py = qx . Intuition: 2 • First spatial derivatives ∂z ∂x • Second spatial derivatives and ∂z ∂y describe the surface orientation of an oject in the image. ∂2z ∂2z ∂x2 , ∂y 2 , and ∂2z ∂x∂y = ∂2z ∂y∂x describe curvature, or the change in surface orientation. • In 2D, the curvature of a surface is specified by the radius of curvature, or its inverse, curvature. • In 3D, however, the curvature of a surface is specified by a Hessian matrix H. This curvature matrix allows us to compute p and q. One issue with using this Hessian approach to update p and q: how do we update the second derivatives r, s, and t? Can we use 3rd order derivatives? It turns out that seeking to update lower-order derivatives with higher-order derivatives will just generate increasingly more unknowns, so we will seek alternative routes. Let’s try integrating our brightness measurements into what we already have so far: Image Irradiance Equation : E(x, y) = R(p, q) Chain Rule Gives : ∂E ∂R ∂p ∂R ∂q Ex = = + ∂x ∂p ∂x ∂q ∂x ∂R ∂R = r+ s ∂p ∂q ∂E ∂R ∂p ∂R ∂q = + Ey = ∂y ∂p ∂y ∂q ∂y ∂R ∂R = s+ t ∂p ∂q In matrix-vector form, we have:    Ex r = s Ey ⎡ ⎤ ⎡ ⎤ ∂R  ∂R ∂p s ⎢ ⎥ ⎢ ∂p ⎥ = H⎣ ⎦ t ⎣ ∂R ⎦ ∂R ∂q ∂q Notice that we have the same Hessian matrix that we had derived for our surface orientation update equation before! Intuition: These equations make sense intuitively - brightness will be constant for constant surface orientation in a model where brightness depends only on surface orientation. Therefore, changes in brightness correspond to changes in surface orientation. How do we solve for this Hessian H? We can: • Solve for Ex , Ey numerically from our brightness measurements. • Solve for Rp , Rq using our reflectance map. However, solving for H in this way presents us with 3 unknowns and only 2 equations (undetermined system). We cannot solve for H, but we can still compute δp, δq. We can pattern match δx and δy by using a (δx, δy)T vector in the same direction as the gradient of the reflectance map, i.e.     δx Rp = ξ δy Rq Where the vector on the lefthand side is our step in x and y, our vector on the righthand side is our gradient of the reflectance map in gradient space (p, q), and ξ is the step size. Intuitively, this is the direction where we can “make progress”. Substituting this equality into our update equation for p and q, we have:     δp R =H p ξ δq Rq   Ex = δξ Ey Therefore, we can formally write out a system of 5 first-order ordinary differential equations (ODEs) that generate our characteristic strip as desired: 3 1. dx dξ = Rp 2. dy dξ = Rq 3. dp dξ = Ex 4. dq dξ = Ey 5. dz dξ = pRp + qRq (“Output Rule”) Though we take partial derivatives on many of the righthand sides, we can think of these quantities as measurements or derived variations of our measurements, and therefore they do not correspond to partial derivatives that we actually need to solve for. Thus, this is why we claim this is a system of ODEs, and not PDEs. A few notes about this system of ODEs: • This system of 5 ODEs explores the surface along the characteristic strip generated by these equations. • Algorithmically, we (a) Look at/compute the brightness gradient, which helps us (b) Compute p and q, which (c) Informs us of Rp and Rq for computing the change in height z. • ODEs 1 & 2 and ODEs 3 & 4 are two systems of equations that each determine the update for the other system (i.e. there is dynamic behavior between these two systems). The 5th ODE is in turn updated from the results of updating the other 4 ODE quantities of interest. • The step in the image (x, y) depends on the gradient of the reflectance map in (p, q). • Analogously, the step in the reflectance map (p, q) depends on the gradient of the image in (x, y). • This system of equations necessitates that we cannot simply optimize with block gradient ascent or descent, but rather a process in which we dynamically update our variables of interest using our other updated variables of interest. • This approach holds generally for any reflectance map R(p, q). • We can express our image irradiance equation as a first order, nonlinear PDE:  ∂z ∂z  E(x, y) = R , ∂x ∂y Next, we will show how this general approach applied to a Hapke surface reduces to our previous SfS results for this problem. 1.2.1 Reducing General Form SfS for Hapke Surfaces Recall reflectance map for Hapke surfaces: r R(p, q) = 1 + ps p + qs q rs Taking derivatives using the system of 5 ODEs defined above, we have from the multivariate chain rule): Δ ∂R ∂p = 1 √ √1 rs 2 1+ps p+qs q ps Δ ∂R ∂q = 1 √ √1 rs 2 1+ps p+qs q qs 1. dx dξ : Rp = 2. dy dξ : Rq = 3. dz dξ : pRp + qRq = p p+qs q √ √s 2 rs 1+ps p+qs q Since the denominator is common in all three of these derivative equations, we can just conceptualize this factor as a speed/step size factor. Therefore, we can simply omit this factor when we update these variables after taking a step. With this, our updates become: 1. δx ← ps 2. δy ← qs 3. δz ← (ps p + qs q) = rs E 2 − 1 Which are consistent with our prior results using our Hapke surfaces. Next, we will apply this generalized approach to Scanning Electron Microscopes (SEMs): 4 1.2.2 Applying General Form SfS to SEMs For SEMs, the reflectance map is spherically symmetric around the origin: R(p, q) = f (p2 + q 2 ) for some f : R → R Using our system of 5 ODEs derived in the general case, we can compute our updates by first calculating derivatives dz dξ : Δ ∂R ∂p = 2f 0 (p2 + q 2 )p Δ ∂R ∂q = 2f 0 (p2 + q 2 )q 1. dx dξ : Rp = 2. dy dξ : Rq = 3. dz dξ : pRp + qRq = 2f 0 (p2 + q 2 )(p2 + q 2 ) dx dy dξ , dξ , and Again, here we can also simplify these updates by noting that the term 2f 0 (p2 + q 2 ) is common to all three derivatives, and therefore this factor can also be interpreted as a speed factor that only affects the step size. Our updates then become: 1. δx ← p 2. δy ← q 3. δz ← p2 + q 2 This tells us we will be taking steps along the brightness gradient. Our solution generates characteristic strips that contain information about the surface orientation. To continue our analysis of this SfS problem, in addition to defining characteristic strips, it will also be necessary to define the base characteristic. 1.3 Base Characteristics Recall that the characteristic strip that our system of ODEs solves for is given by a 5-dimensional surface composed of:  T x y z p q ∈ R5 Another important component when discussing our solution to this SfS problem is the base characteristic, which is the projection of the characteristic strip onto the x, y image plane:  T  T x(ξ) y(ξ) = projectionx,y {characteristic strip} = projectionx,y { x y z p q } A few notes/questions about base characteristics? • What are base characteristics? These can be thought of as the profiles in the (x, y) plane that originate (possibly in both directions) from a given point in the initial curve. • What are base characteristics used for? These are used to help ensure we explore as much of the image as possible. We can interpolate/average (depending on how dense/sparse a given region of the image is) different base characteristics to “fill out” the solution over the entire image space. • How are base characteristics used in practice? Like the solution profiles we have encountered before, base characteristics are independent of one another - this allows computation over them to be parallelizable. Intuition: These sets of base characteristics can be thought of like a wavefront propagating outward from the initial curve. We expect the solutions corresponding to these base characteristics to move at similar “speeds”. 1.4 Analyzing “Speed” When considering how the speeds of these different solution curves vary, let us consider ways in which we can set constant step sizes/speeds. Here are a few, to name: 1. Constant step size in z: This is effectively stepping “contour to contour” on a topographic map. Achieved by: Dividing by pRp + qRq =⇒ dz dξ = 1. 2. Constant step size inqimage: A few issues with this approach. First, curves may run at different rates. Second, methodology fails when Rp2 + Rq2 = 0. q p d Achieved by: Dividing by Rp2 + Rq2 =⇒ dξ ( δx2 + δy 2 ) = 1. 5 3. Constant Step Size in 3D/Object: Runs into issues when Rp = Rq = 0. q p Achieved by: Dividing by Rp2 + Rq2 + (pRp + qRq )2 =⇒ (δx)2 + (δy)2 + (δz)2 = 1. 4. Constant Step Size in Isophotes: Here, we are effectively taking constant steps in brightness. We will end up dividing by the dot product of the brightness gradient and the gradient of the reflectance map in (p, q) space. ∂R ∂E ∂R Achieved by: Dividing by ( ∂E ∂x ∂p + ∂y ∂q )δξ = ((Ex , Ey ) · (Rp , Rq ))δξ =⇒ δE = 1. 1.5 Generating an Initial Curve While only having to measure a single initial curve for profiling is far better than having to measure an entire surface, it is still undesirable and is not a caveat we are satisfied with for this solution. Below, we will discuss how we can generate an initial curve to avoid having to measure one for our SfS problem. As we discussed, it is undesirable to have to measure an initial curve/characteristic strip, a.k.a. measure:  T x(ξ) y(ξ) z(ξ) p(ξ) q(ξ) But we do not actually need to measure all of these. Using the image irradiance equation we can calculate ∂z ∂η : E(x, y) = R(p, q) ∂z ∂z ∂x ∂z ∂y = + ∂η ∂x ∂η ∂y ∂η ∂x ∂y =p +q ∂η ∂η Where p and q in this case are our unknowns. Therefore, in practice, we can get by with just initial curve, and do not need the full initial characteristic strip - i.e. if we have x(η), y(η), and z(η), then we can compute the orientation from the reflectance map using our brightness measurements. 1.5.1 Using Edge Points to Autogenerate an Initial Curve An even more sweeping question: do we even need an initial curve? Are there any special points on the object where we already know the orientation without a measurement? These points are along the edge, or occluding boundary, of our objects of interest. Here, we know the surface normal of each of these points. Could we use these edge points as “starting points” for our SfS solution? It turns out, unfortunately, that we cannot. This is due to the fact that as we infinisimally approach the edge of the boundary, ∂z ∂z we have that ∂x → ∞, ∂y → ∞. Even though the slope pq is known, we cannot use it numerically/iteratively with our step updates above. 1.5.2 Using Stationary Points to Autogenerate an Initial Curve Let us investigate another option for helping us to avoid having to generate an initial curve. Consider the brightest (or darkest) point on an image, i.e. the unique, global, isolated extremum. This point is a stationary point in (p, q) space, i.e. ∂R ∂R = = 0, granted that R(p, q) is differentiable. ∂p ∂q This would allow us to find the surface orientation (p, q) solely through optimizing the reflectance map function R(p, q). One immediate problem with this stationary point approach is that we will not be able to step x and y, since, from our system of ODEs, we have that at our stationary points, our derivatives of x and y are zero: dx = Rp = 0 dξ dy = Rq = 0 dξ Additionally, if we have stationary brightness points in image space, we encounter the “dual” of this problem. By definition stationary points in the image space imply that: ∂E ∂E = =0 ∂x ∂y 6 This in turn implies that p and q cannot be stepped: dp ∂E = =0 dξ ∂x dq ∂E = =0 dξ ∂y Intuition: Intuitively, what is happening with these two systems (each of 2 ODEs ((x, y) or p, q), as we saw in our 5 ODE system above) is that using stationary points in one domain amounts to updates in the other domain going to zero, which in turn prevents the other system’s quantities of interest from stepping. Since δz depends on δx, δy, δp, δq, and since δx, δy, δp, δq = 0 when we use stationary points as starting points, this means we cannot ascertain any changes in δz along these points. Note: The extremum corresponding to a stationary point can be maximum, as it is for Lambertian surfaces, or a minimum, as it is for Scanning Electron Microscopes (SEM). However, even though we cannot use stationary points themselves as starting points, could we use the area around them? If we stay close enough to the stationary point, we can approximate that these neighboring points have nearly the same surface orientation. • Approach 1: Construct a local tangent plane by extruding a circle in the plane centered at the stationary point with radius  - this means all points in this plane will have the same surface orientation as the stationary point. Note that mathematically, a local 2D plane on a 3D surface is equivalent to a 2-manifold [1]. This is good in the sense that we know the surface orientation of all these points already, but not so great in that we have a degenerate system - since all points have the same surface orientation, under this model they will all have the same brightness as well. This prevents us from obtaining a unique solution. • Approach 2: Rather than constructing a local planar surface, let us take a curved surface with non-constant surface orientation and therefore, under this model, non-constant brightness. 1.5.3 Example Suppose we have a reflectance map and surface function given by: Reflectance Map : R(p, q) = p2 + q 2 Surface Function : z(x, y) = x2 + y 2 Then, we can compute the derivatives of these as: ∂z = 2x ∂x ∂z q= = 4y ∂y E(x, y) = R(p, q) = p2 + q 2 = 4x2 + 16y 2 ∂E = 8x Ex = ∂x ∂E = 32y Ey = ∂y p= ∂E The brightness gradient ( ∂E ∂x , ∂y ) = (8x, 32y) = (0, 0) at the origin x = 0, y = 0. Can we use the brightness gradient to estimate local shape? It turns out the answer is no, again because of stationary points. But if we look at the second derivatives of brightness: ∂2E ∂E = (8x) = 8 2 ∂x ∂x ∂E ∂2E = = (32y) = 32 ∂y ∂y 2 ∂2E ∂ ∂ = = (32y) = (8y) = 0 ∂x∂y ∂x ∂y Exx = Eyy Exy 7 These second derivatives, as we will discuss more in the next lecture, will tell us some information about the object’s shape. As we have seen in previous lectures, these second derivatives can be computed by applying computational molecules to our brightness measurements. High-Level Algorithm (NOTE: This will be covered in greater detail next lecture) Let’s walk through the steps to ensure we can autogenerate an initial condition curve without the need to measure it: 1. Compute stationary points of brightness. 2. Use 2nd derivatives of brightness to estimate local shape. 3. Construct small cap (non-planar to avoid degeneracy) around the stationary point. 4. Begin solutions from these points on the edge of the cap surrounding the stationary point. 1.6 References 1. Manifold, https://en.wikipedia.org/wiki/Manifold 8