[Solved] CS148 Homework 5: Sampling and Geometry

Introduction to Computer Graphics and Imaging

Youll notice that this problem set is a few more pages than the usual CS 148 homework; part of this problem set is a reading assignment. The problems are intended to make sure you understand both what weve covered in class and what the narration below discusses: You are responsible for both on the final exam!

In the past few weeks, we have covered a diverse set of topics relevant to computer graphics, from sampling and Fourier theory to computer aided geometric design. In this homework, we will be exploring a number of these topics to help you get better acquainted with the issues we have discussed in class.

Well start by discussing some of the details of sampling theory laid out in class to help determine the resolution you need to display different continuous images. Recall that we stated the following theorem in class:

Theorem(-ish). Any function can be written as a combination of simple wiggles.

Of course, the word any here is very suspect, but as computer scientists well ignore the intricacies of Lebesgue integration and Schwartz distributions in favor of an entirely un-rigorous approach.^[1]

Before we define the Fourier transform, we should state a famous mathematical relationship that you likely have seen before:

eⁱ = cos + i sin (1)

This equation often is called Eulers Formula. Lets fill in some details:

• e = 718 . . . is the natural logarithm base.

• i = 1 is the imaginary unit of the complex plane C.
• R is any real number.

If you never have seen (1), you should make the following straightforward observation: This equation makes no sense. And, youre not wrong! In particular, raising any number to an imaginary exponent i isnt a reasonable operation in the usual sense of exponentiation as repeated multiplication.

In fact, in many ways we can think of Eulers formula as a notational convenience. We can prove that the right hand side of (1) behaves identically to exponentiation in almost every possible way. This observation allows us to compactify derivations of many facts we know from high school trigonometry. For instance, consider the following derivation:

cos(+)+ i sin(+) = eⁱ⁽⁺⁾ by (1)

= eⁱeⁱ by properties of exponentiation

= (cos + i sin )(cos + i sin ) by (1)

= (cos cos sin sin )+ i(sin cos + cos sin ) since i^[2] = 1

Make sure you understand each step here before moving on!

But, what did we just show? If you take the real and imaginary components of the left and right hand sides of the derivation above separately, we find:

cos(+) = cos cos sin sin sin(+) = sin cos + cos sin

You probably worked much harder to prove these statements in high school trigonometry.

Anyway, our goal today is to understand the Fourier transform of a function. We now have all the tools we need to do so.

Definition (Fourier Transform). The Fourier transform of a function f : R R is the function f() given by

f dx (2)

The variable here represents a frequency. So, f() represents the component of f that wiggles at frequency .

Problem 1 (10 points). Recall that in Lecture 11 of CS 148 we stated that the Fourier transform does the

simplest possible thing: multiplying a function against a wiggle at a fixed frequency and integrating. Use Eulers formula (1) to justify why the Fourier transform (2) does exactly this. Why does f() need real and imaginary components?

As tempting as it is for your mathematician instructor, we will not derive most of the (beautiful) relationships between properties of a function and properties of its Fourier transform. There is, however, one relationship that is fundamental to the sampling theory we discussed in lecture.

Problem 2 (15 points). Take any real number a > 0, and define h(x) = f(ax). Derive the following relationship:

h (3)

Explain how this relationship substantiates the claim made in class that increasing the rate at which you sample a (band-limited) function makes it less likely that you will see artifacts.²

Figure 1: For problem 3(a).

Hint: To derive (3), substitute h into the definition of the Fourier transform and make the substitution y = ax.

There are myriad other facts we could derive from the Fourier transform f that are relevant to sampling theory as it applies to computer graphics, but well abandon our discussion at this point in favor of moving to the topic of geometric modeling.

Problem 3 (10 points). (a) The control points of a quadratic Bezier curve g(t) and a cubic Bezier curve f(t) are shown in Figure 1. Construct the points g(3/2) and f(1/3). Youll need to turn in a copy of Figure 1 with your constructions drawn on top; you may wish to enlarge these structures and/or print them on separate pages for more space. (b) Find the (quadratic and cubic, resp.) blossoms of the functions g(t) = t² 2t + 1 and h(t) = (t 2)³.

One of the key ideas of computer-aided geometric design (CAGD) as it pertains to curves is that of finding a basis for polynomial functions of a given degree. As we discussed in lecture, the idea of a basis here is slightly more abstract than what you might have seen in Math 51, so its worth putting some thought into our construction.

We consider a basis of the set of polynomials of degree n to be a set of n + 1 polynomial functions p₁(t), . . . , p_n+1(t) such that any degree-n polynomial can be written as a linear combination of these p_is. For example, an obvious basis for the degree-three polynomials would be p₁(t) = 1, p₂(t) = t, p₃(t) = t², p₄(t) = t³. Why? Because the general form of a degree-three

polynomial is at³ + bt² + ct + d, or equivalently ap₄(t)+ bp₃(t)+ cp₂(t)+ dp₁(t).

Certain polynomial bases are convenient for certain problems, so a strategy in CAGD is to work in a basis that simplifies our math as much as possible. We explored this strategy in class when we derived the Hermite basis, which simplified the problem of finding a curve with given begin and end points/tangents. In our final written problem, we will derive an alternative basis for cubic polynomials simplifying math for Bezier curves.

Problem 4 (15 points). Suppose we label the control points~a,^~b,~c, and d of a cubic curve in Figure 2 and^~ use the structure drawn to find the point f(t) (assume the control points specify the curve between f(0) and f(1)). Label the constructed points as linear combinations of the control points; we provide you with one label to get you started. Use your label for f(t) to define the Bernstein basis for cubic polynomials, that is, a basis that makes it easy to write down the cubic curve associated with a set of four control points.

Figure 2: For problem 4.

Figure 3: Examples of MLS output.

For the coding part of this assignment, you will implement the Moving Least Squares (MLS) image deformation algorithm introduced in SIGGRAPH 2006.^[3] Examples of the output of MLS are shown in Figure 3. The purpose of this algorithm is fairly straightforward: The user marks a number of handles on an image and then starts moving them around. As the handles are moved, the image deforms smoothly while satisfying the constraints that the handles are mapped to their targets. MLS can be used to achieve lots of interesting image effects with relatively simple computations.

To motivate MLS, we first consider a slightly easier problem. Lets say we mark a number of points ~p₁, . . . ,~p_n R² on the plane, and we specify a list of targets ~q₁, . . . ,~q_n R². We wish to find an affine transformation of the plane given by ~pfor some M R²²,^~t R² such that M for each i {1, . . . , n}. When is n > 2, the is key: it might be the case that no affine transformation can map all the ~p_is to their corresponding~q_is exactly!

Thus, we can only hope to approximate such a relationship. To this end, well introduce one more set of inputs: a set of weight values w₁, . . . , w_n > 0. A large weight w_i says we really care that ~p_i maps to , while a small w_i 1 means we dont really care if M and^~t dont do a good job mapping p_i to q_i.

Note that if M ^~ then M 0. In other words, we desire that

~q_ik² 0. So, lets define an energy:

n

E(M, (4)

i=1

We can make a number of important observations about E:

• We always have E
• If there exist^~t and M such that ~q_i = Mexactly for all i {1, . . . , n}, then E will be minimized with value 0 at this point.
• When such a zero minimum is not possible, the minimum E > 0 occurs when M and ^~t give the best affine map approximating ~p_i 7 ~q_i in a least squares sense. The weights w_i prioritize with ) pairs matter the most.

Thus, if the user inputs lots of pairs ), we can frame an optimization problem to find M and

~t as follows:

min E(M,~t) (5)

MR22,~tR2

You probably spent much of calculus class solving problems exactly of this nature! Simply take derivatives of E with respect to the elements of M and^~t, set them to zero, and solve.

Well resist the urge to force you to derive these formulas yourself, although youre encouraged to check them on scratch paper; our discussion above plus a few tricks from Math 51 are all you need. Instead, well jump right to the solution of our minimization problem (5). To clean up our notation a bit, lets define some helper variables:

wi~pi

(6) w_i

wi~qi

(7) w_i

Its easy to differentiate (5) with respect to the elements of^~t (at least try this part at home!), yielding the following relationship when we set this derivative to zero:

~topt =~q Mopt~p (8)

That is, we can find the optimal translation^~t_opt easily enough once we know the optimal matrix

Mopt.

Finding M_opt takes a bit more bookkeeping. Well need two more sets of helper variables:

p_i (9) q_i =~q_i ~q (10) We can reduce finding M_opt to a least-squares problem, so when the smoke clears well get the following normal equation solution:

!1

Moptwjqjpwipipi> (11)

This finishes our story. From the inputs , . . . , (~p_n,~q_n, w_n) we can compute ~p and ~q. We can use these to compute p_i and q_i for all i. Then, we use (11) to compute M_opt. Finally, we find^~t_opt using (8). The formulas might look complicated, but theyre quite simple to evaluate; the 2 2 inverse in (11) can be found in closed-form with the usual formula (built into Eigen using the .inverse() method of Matrix2d).

We may wish to put additional constraints on M_opt. Its easy to see that regardless of our constraints on M_opt, the formula (8) for^~t_opt stays the same. This aside, if we want only to consider similarity matrices M that contain no shearthat is, they only can rotate and scalethen we constrain M^>M = ²I for some R (think about why!) and find:

similarity 1 n

^M_opt = w_i p_{i i} (12)

i=1 i

where

n

= w_ikp_ik² (13)

i=1

If we force M only to encode rotations without any scaling, then (12) still works, but we replace with

(14)

Here, for vector~a = (a₁, a₂) we denote~a = (a₂, a₁).

MLS simply applies the formulas we developed above with a clever choice of the w_i values. Were going to render an image by drawing a k k grid of quads with evenly-spaced texture coordinates. In our code the grid of texture coordinates is stored in the variable textureCoords; you should not change the values in this array. Well deform the grid to a new set of positions (stored in imageCoords) using the MLS algorithm and render the deformed image using textured triangles; see display() to see our basic setup.

The user shift-clicks the control points ~p_i to make red handles appear (shift-clicking those red handles again makes them disappear); the indices of these points in imageCoords are stored in the set active. The user then can drag the handles to new targets~q_i, the positions of which are stored in imageCoords. This much functionality is implemented for you in the skeleton code.

Your job is to implement deformImage(), which writes new values for elements of the array imageCoords. The basic algorithm goes like this:

• Iterate over each texture coordinate ~p in the array textureCoords.
• If ~p is one of the ~p_is, then its position in imageCoords is ~q_i and can be left alone. Otherwise proceed with the steps below.
• Compute weights w_i for the provided constant > Recall that ~p_i denotes the texture coordinate of the ith point in active. Notice what this choice of w_i is doing: Control points that are close to ~p get more weight than those that are far away.
• Compute M_opt and^~t_opt using the weights w_i you just computed and the formulas above.
• Replace the position in imageCoords with~q = M_opt~p .

In other words, each (sampled) point on the base texture computes its own least-squares map and uses that map to find a target.

Problem 5 (50 points). Implement MLS deformation starting from the skeleton code. You must implement the affine, similarity, and rotation-only versions of the algorithm, using command line options stored in deformationType to determine what to do in the application.

Problem 6 (5 points extra credit). When does our choice of w_i in MLS break down? Can you suggest a replacement? Describe your solution in the written homework; no need to add it to your code, as the implementation change might (hint!) require a using a structure other than a grid mesha considerable change!

[1] Please come to Justins office hours or post on Piazza if you would like work out the details of making this theory more rigorous. Some fascinating and deep mathematics are hiding in the straightforward-looking derivations we cover in CS 148!

[2] For students with a strong math background, note that in CS 148 you need not be concerned with whether your integrals converge or are well-defined: feel free to apply and reasonable-looking operation without regard for such analytical details!

[3] Youre encouraged to read the original paper for a more detailed discussion. Note that the paper uses row rather than column vectors, so their math is a bit different.