ar = This problem develops a tool that will be used in later homework assignments for an application called photometric stereo. To approximate the derivatives of a function f(x) that is sampled on a grid 21,..., In where li+1 = Li +8, a typical finite difference approach is: af (:) f(Ri+1) - f(x) x=r; 8 When the sample spacing is 8 = 1, this approximation simplifies to af (2) : f(Li+1)-f(xi). We can express this relation for all Xsamples via the matrix-vector product [f'(21)] Tf(21)] f'(x2) f(22) Dn f'(X) = ax =fi f(an) where D, is the so-called first-difference matrix -1 1 -1 1 Dn -1 1 1 Here we choose to set Dn(n. 1) = 1, which corresponds to the (perhaps unexpected) approximation f'(In) f(21) - f(In). This choice is called a periodic boundary condition because essentially we are assuming that the domain wraps around. We make this assumption because the resulting Dn is a circulant matrix, so its eigenvectors can be computed in closed form! The goal of this problem is for you to derive and implement the analog of Dn for two-dimensional differentiation. Let f(x,y) be a function of two variables. We can approximate its partial derivatives using finite differences as follows: af (x,y) f(x+1, y) - f(x,y) = f(x+1,y)-f(x,y) (1) ax (x + 1)-2 af(x,y) f(x,y+1) - f(,y) f(x,y+1)-f(x,y) (2) ay (y +1) - y To simplify notation, define the mx n matrices FXY, DFDX, and DFDY having elements as follows: FXY(i,j) = f(ij) DFDX(i, j) af (1,1) a.r afij) DFDY(ij) The coordinate is along the column of FXY and the y coordinate is along the row of FXY, so we can think in terms of FXY(x,y). Define the corresponding mn x 1 vectors fxy, dfdx, and dfdy to be the vectorized versions obtianed using the vec(-) operation (HW1, Problem 9). With this notation, we can succinctly express equations (I) and (2) as afdx = A fxy, dfdy where A is a 2mn x mn matrix. (a) Find an expression for A in terms of the first difference matrices Dn and Dm, appropriately sized identity matrices, and appropriate Kronecker products of these matrices. Use periodic boundary conditions. (b) Once you have determined A, write the function first diffs 2d matrix that takes as input the dimen- sions m and n of FXY and returns the appropriate A matrix, stored in sparse format, since the resulting matrix will be prohibively large to store as a full double-precision matrix. PYTHON users should make use of the scipy.sparse library. (c) In this problem your function is designed to compute finite-difference approximations to derivatives along I and y. If you create an mx n array X that is a picture of a disk, then the finite derivatives will be mostly zero except near the edges of the disk. The code for generating such a disk is on the course website. Test your code using your function and turn in plots of the resulting images, which should match those below. Note that the PYTHON reshape requires setting order='F'. Fix,y) dix.ch df(xyl dy 5 5 5 10 10 10 15 15 15 20 20 ( ) 20 25 25 25 30 30 4 6 10 12 14 16 18 20 B a 10 12 14 15 18 20 2 6 B 10 12 14 16 18 20 clear % problem size m= 30; n = 20; [cols.rows] =meshgrid(1:n 1.m); % disk parameters cx = 8; cy=20; I= 5; % create matrix with disk X = (rows-cy). 2 +(cols - cx)^2 =r^2; figure imagesc(X) colormap gray title(f(x,y)) fig = gcf; fout="/fxy_matlab'; print(fig fout.-depsc') % test your first diffs matrix A = first_diffs_2d_matrix(m,n); x=XC.); DFDX = A(1:m*n,:)*x; DFDY = A(m*n+1:2*m*n.:)*x; figure imagesc(reshape(DFDX mn)) colormap gray title(dfxy) dx') fig = gcf; fout="/dfdx_matlab'; print(fig.fout.'-depsc') figure imagesc(reshape(DFDY.mn)) colormap gray title(df(x,y)/dy) fig = gcf; fout="/dfdy_matlab'; print(fig fout, -depsc')
