Lecture 3 Notes: Time to Contact, Focus of Expansion, Direct Motion Vision Methods, Noise Gain

6.801/6.866: Machine Vision, Lecture 3 1 Lecture 3: Time to Contact, Focus of Expansion, Direct Motion Vision Methods, Noise Gain 1.1 Noise Gain Example/motivation: Indoor GPS. Rather than using localization of position with satellites, use indoor cellular signals. Fun fact: EM waves travel at 1 ns/foot. Dilution of Precision: • How far off GPS is w.r.t. your location. • Important to note that your dilution of precision can vary in different directions - e.g. horizontal precision is oftentimes greater than vertical position. 1.2 1.2.1 Forward and Inverse Problems of Machine Vision Scalar Case One way to conceptualize the goals of machine vision, as well as to highlight why noise gain is important, is by considering the following problems with some one-dimensional input x and an output y = f (x): Δ • The forward problem: x → y = f (x) Δ • The inverse problem∗ : y = f (x) → x *(this term comes up a lot in machine vision, computer graphics, and robotics) In machine vision, we oftentimes focus on solving the inverse problem, rather than the forward problem. Intuitively, we usually observe some y = f (x), and from this infer the latent parameters x using our model f . When it is possible to express the inverse of a function in closed form or via a matrix/coefficient, we can simply solve the inverse problem using: x = f −1 (y). More importantly, to build a robust machine vision system to solve this inverse problem, it is critical that small perturbations in y = f (x) do not need to large changes in x. Small perturbations need to be taken into account in machine vision problems because the sensors we use exhibit measurement noise. The concept of noise gain can come in to help deal with this uncertainty. Consider a perturbation δy that leads to a perturbation δx when we solve the inverse problem. In the limit, as δ ∈ R → 0, then we arrive at the definition of noise gain: noise gain = δx 1 = 0 = δy f (x) 1 df (x) dx 1 (1) Like other concepts/techniques we’ve studied so far, let’s understand when this system fails. Below are two cases; we encourage to consider why they fail from both a mathematical and intuitive perspective (hint: for the mathematical component, look at the formula above, and for the intuitive component, think about how the effect on x from a small change in y in a curve that is nearly flat): • f 0 (x) = 0 (flat curve) • f 0 (x) ≈ 0 (nearly flat curve) 1.2.2 Vector Case Now that we’ve analyzed this problem in the scalar case, let us now consider it in the vector/multi-dimensional case. Since images are inherently multidimensional, this more general vector case is where we will find ourselves. First, we can restate these problems: • Forward Problem: x = Mb, M ∈ Rm×n for m, n ∈ N • Inverse Problem: b = M−1 x, M ∈ Rm×n for m, n ∈ N But how good is this answer/approach? If x changes, how much does b change? We quantify noise gain in this case as follows: “noise gain” → ||δb|| ,δ ∈ R ||δx|| (2) *NOTE: This multidimensional problem is more nuanced because we may be in the presence of an anisotropic (spatially non-uniform) noise gain - e.g. there could be little noise gain in the x1 direction, but a lot of noise gain in the x2 direction. As in the previous case, let’s analyze when this approach fails. To do this, let’s consider M −1 to help build intuition for why this fails: ⎡ ⎤ · · · 1 ⎣· · ·⎦ M −1 = (3) det|M | · · · Let’s ignore the specific entries of M −1 for now, and focus on the fact that we need to compute the determinant of M . When is this determinant zero? This will be the case whenever there exists linear dependence in the columns of M. As we saw before, two cases that can yield to poor performance will be: • det|M | = 0: This corresponds to a non-invertible matrix, and also causes the noise term to blow up. • det|M | ≈ 0: Though this matrix may be invertible, it may cause numerical instability in the machine vision system, and can also cause the noise term to blow up. Let’s also revisit, just as a refresher, the inverse of a 2 × 2 matrix:  a A= c b d  (4) (5) A−1 =  1 d −c detA   1 −b d = a −c ad − bc  −b a Now let’s verify that this is indeed the inverse:     1 1 d −b a b ad − bc −1 A A= = c d ad − bc −c a ad − bc cd − cd 1.3   1 −ab + ab = ad − bc 0 (6)  0 = I2 1 Review from Lecture 2 Before we dive into the next set of concepts, let’s also revisit some of the concepts discussed in the previous lectures. 2 (7) 1.3.1 Two-Pixel Motion Estimation, Vector Form First, let’s recall the two pixel motion estimation set of equations: uEx1 + vEy1 + Et1 = 0 uEx2 + vEy2 + Et2 = 0 (8) (9) Rewritten in matrix form:    1 u Ey2 = v −E Ex1 Ey2 − Ey1 Ex2 x2 −Ey1 Ex1  (10) Take note of the denominator on the right-hand side. Does this term look familiar to the determinant term above? We were indeed solving an instance of the inverse problem. If this determinant-like quantity (Ex1 Ey2 − Ey1 Ex2 ) is small, the noise is greatly amplified. We saw that this happens when brightness gradients are similar to one another. 1.3.2 Constant Brightness Assumption, and Motion Equation Derivation Recall the constant brightness assumption’s mathematical formulation: Constant Brightness Assumption =⇒ dE =0 dt (11) *(Please note the quantity above is a total derivative.) Intuition Behind This: As the object/camera moves, the physical properties of the camera do not change and therefore the total derivative of the brightness w.r.t. time is 0. From the chain rule, we can rewrite this total derivative assumption: dE = 0 =⇒ uEx + vEy + Et = 0 dt (12) *(Recall this is for when x and y are parameterized w.r.t. time, i.e. x = x(t), y = y(t).) The above constraint is known as the Brightness Change Constraint Equation (BCCE). 1.3.3 Optical Mouse Problem Recall our motion estimation problem with the optical mouse, in which our objective is no longer to find the point where the BCCE is strictly zero (since images are frequently corrupted by noise through sensing), but to minimize the LHS of the BCCE, i.e: ZZ Δ min{J(u, v) = (uEx + vEy + Et )2 dxdy} (13) u,v We solve the above using unconstrained optimization and by taking a “least-squares” approach (hence why we square the LHS of the BCCE). Solve by setting the derivatives of the two optimizing variables to zero: dJ(u, v) dJ(u, v) = 0, =0 du dv 1.3.4 (14) Perspective Projection Recall the perspective projection equations in the scalar form: x X y Y = (x-component), = (y-component) f Z f Z (15) *(Note that capital coordinates are in the world space, and lowercase coordinates are in the image space.) What if these quantities are changing in the world w.r.t. time? Take time derivatives: 1 dx 1 dX 1 dZ = − 2X dt f dt Z dt Z Writing these for x and y: 3 (16) • x: 1 fu = 1 ZU − W X Z Z • y: 1 fv = 1 ZV − W Y Z Z Note again the following definitions: • u → image velocity in the x direction • v → image velocity in the y direction • U → world velocity in the X direction • V → world velocity in the Y direction When are these points in (u, v) space interesting? When u = v = 0 - this is the Focus of Expansion (FOE). The FOE given by (x0 , y0 ) in two dimensions: (x0 , y0 ) = ( 1.4 f U f V , ) ZW ZW (17) Time to Contact (TTC) In the previous lecture, we discussed the derivation of Time to Contact (TTC) in terms of meters: Z Δ Z meters == dZ = meters = seconds W seconds dt (18) Let us express the inverse of this Time to Contact (TTC) quantity as C: Δ C= W 1 = Z TTC (19) Let us now suppose we parameterize our spatial velocities u and v according to u = Cx, v = Cy, where C is the inverse TTC value we introduced above. Then, substituting these into the BCCE equation, we have: Recall BCCE: uEx + vEy + Et = 0 Substitute: C(xEx + yEy ) + Et = 0 Et Solve for C: C=− xEx + yEy (20) (21) (22) The denominator in the derivation of C is the “radial gradient”: g = xEx + yEy = (x, y) · (Ex , Ey ) (23) Building Intuition: If we conceptualize 2D images as topographic maps, where brightness is the third dimension of the surface (and the spatial dimensions x and y comprise the other two dimensions), then the brightness gradient is the direction of steepest ascent up the brightness surface. Another note: (x, y) is a radial vector, say in a polar coordinate system. Hence why the above dot product term is coined the name “radial gradient”. This gradient is typically normalized by its L2 /Euclidean norm to illustrate the multiplication of the brightness gradients with a radial unit vector): g= p x y x2 + y 2 ( p ,p ) · (Ex , Ey ) 2 2 2 x + y2 x +y (24) This g quantity can be thought of as: “How much brightness variation is in an outward direction from the center of the image?” For a more robust estimate, let us again employ the philosophy that estimating from more data points is better. We’ll again take a least-squares approach, and minimize across the entire image using the parameterized velocities we had before. In this case, since we are solving for inverse Time to Contact, we will minimize the error term over this quantity: ZZ Δ min{J(C) = (C(xEx + yEy ) + Et )2 dxdy} (25) C 4 Without the presence of measurement noise, the optimal value of C gives us an error of zero, i.e. perfect adherence to the BCCE. However, as we’ve seen with other cases, this is not the case in practice due to noise corruption. We again will use unconstrained optimization to solve this problem. Taking the derivative of the objective J(C) and setting it to zero, we obtain: ZZ dJ(C) = 0 =⇒ 2 (C(xEx + yEy ) + Et )(xEx + yEy )dxdy = 0 dC (26) This in turn gives us: ZZ C (xEx + yEy )2 dxdy + ZZ (xEx + yEy )Et dxdy = 0 RR (xEx + yEy )Et dxdy 1 = C = − RR TTC (xEx + yEy )2 dxdy (27) (28) A few notes here: • Et can be computed by taking differences between a pixel at different points in time. • To implement a framework like this, we can do so with accumulators. • This is an instance of “direct computation”. • When objects/important components of a scene are far away, we can combine image pixels and run these algorithms at multiple scales - this allows us to compute/analyze motion velocities over different timescales. This can also be helpful for building more computationally-tractable motion estimation implementations, since the number of pixels over which computations must occur can be reduced quadratically with the scale. • Note one problem we run into with this approach - each pixel we apply this estimation approach to introduces one equation, but two unknowns. • Note that neighboring pixels typically exhibit similar behavior to one another. Let’s briefly return to the concept of the radial gradient. When is this zero? • E = 0 everywhere (coal mine) • (x, y) · (Ex , Ey ) = 0 (radial gradient is zero) Now, let’s take a more general case when we have world motion, i.e. U = 6 0, V = 6 0. For our two components x and y: • x-component: u U XW = − f Z f Z fU W u= −X Z Z u = A − XC Δ fU Δ W Where: A = ,C = Z Z (29) (30) (31) (32) • y-component: v V Y W = − f Z f Z fV W v= −Y Z Z v =B−YC Δ fV Δ W Where: B = ,C = Z Z (33) (34) (35) (36) Note that for the world quantities A and B, we also have the following identities (note that the Focus of Expansion (FOE) is given by the point (x0 , y0 ): 5 • A= fU Z = Cx0 • B= fV Z = Cy0 Building Intuition: “As I approach the wall, it will loom outward and increase in size.” Now returning to our BCCE, this time with A, B, and C: AEx + BEy + C(xEx + yEy ) + Et = 0 (37) We can again use least-squares to minimize the following objective enforcing the BCCE. This time, our optimization aim is to minimize the objective function J(A, B, C) using the quantities A, B, and C: ZZ Δ min {J(A, B, C)} = (AEx + BEy + C(xEx + yEy ) + Et )2 dxdy} (38) A,B,C Use unconstrained optimzation with calculus and partial derivatives to solve. Since we have three variables to optimize over, we have three first-order conditions (FOCs): • dJ(A,B,C) dA =0 • dJ(A,B,C) dB =0 • dJ(A,B,C) dC =0 Using the Chain rule for each of these FOCs, we can derive and rewrite each of these conditions to obtain 3 equations and 3 Δ unknowns. Note that G = xEx + yEy . • A variable: ZZ 2 ZZ A (AEx + BEy + C(xEx + yEy ))Ex = 0 ZZ ZZ ZZ Ex2 + B Ex Ey + C GEx = − E x Et (39) (40) • B variable: ZZ (AEx + BEy + C(xEx + yEy ))Ex = 0 ZZ ZZ ZZ ZZ A Ey Ex + B Ey2 + C GEy = − Ey Et 2 (41) (42) • C variable: ZZ (AEx + BEy + C(xEx + yEy ))Ex = 0 ZZ ZZ ZZ 2 GEx + B GEy + C G =− GEt 2 ZZ A Putting all of these equations together: ZZ ZZ ZZ ZZ Ex Ey + C GEx = − Ex Et A Ex2 + B ZZ ZZ ZZ ZZ A E y Ex + B Ey2 + C GEy = − Ey Et ZZ ZZ ZZ ZZ A GEx + B GEy + C G2 = − GEt This can be compactly written as a matrix-vector product equation: ⎡ RR ⎤⎡ ⎤ RR RR ⎡RR ⎤ Ex2 Ex Ey RR Ex G A Ex Et RR RR RR ⎢ ⎥ ⎦ ⎣ ⎦ ⎣ ⎣ RR Ey Ex RR Ey2 RREy2G⎦ B = − RR Ey Et C GEt GEx GEy G As in the time-to-contact problem above, this can again be implemented using accumulators. Let’s end on a fun fact: Did you know that optical mice have frame rates of 1800 fps? 6 (43) (44)