L5 - OPTI 280: Computer Programming Workshop Robin Palit...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: OPTI 280: Computer Programming Workshop Robin Palit College of Optical Sciences University of Arizona [email protected] OPTI 280 Spring 2010 Lecture 5 R. Palit, 1 Lesson Plan •  Matrices •  Element by element operations •  Matrix operations •  Applications of Matrices. •  •  •  •  Simultaneous Linear Equations. Imaging. Paraxial Raytrace. Rotation. NOTE: You will not be quizzed on material from this section. •  Advanced Linear Algebra Tools. •  Miscellaneous Commands. OPTI 280 Spring 2010 Lecture 5 R. Palit, 2 Debugging Tips •  Check typos, improper initializations, incorrect variables. •  Use variables that have one purpose and meaningful names. •  clear all variables prior to executing code. •  Run code in step mode to localize the error. •  Use breakpoint and check variables •  Make one change at a time. OPTI 280 Spring 2010 Lecture 5 R. Palit, 3 Matrices •  A matrix is an N-dimensional array of numbers (N ≥ 2). ȹ 0.95 0.48 0.45 0.44ȹ ȹ ȹ A = ȹ 0.23 0.89 0.01 0.61ȹ ȹ ȹ ȹ 0.60 0.0 0.82 0.79Ⱥ •  The elements of a matrix are accessed using (row,column) indices. € •  e.g. A(2,3) = 0.01, A(3,1) = 0.60, etc. •  Symbolically, A(mm,nn) = Am,n OPTI 280 Spring 2010 Lecture 5 R. Palit, 4 MatLab Element by Element Operations ȹ 0.2 -0.3 0.2 ȹ ȹ ȹ A = ȹ -0.1 0.6 0.5 ȹ ȹ ȹ ȹ 0.3 0.6 -0.3Ⱥ ȹ 0.9 0.2 -0.1ȹ ȹ ȹ B = ȹ -0.9 0.4 -0.3ȹ ȹ ȹ ȹ 0.2 -0.5 0.8 Ⱥ € •  Addition, +! € ȹ 1.1 -0.1 0.1ȹ ȹ ȹ C = A + B = ȹ -1.0 1.0 0.2ȹ ȹ ȹ ȹ 0.5 0.1 0.5Ⱥ •  Subtraction, -! € ȹ -0.7 -0.5 0.3 ȹ ȹ ȹ C = A - B = ȹ 0.8 0.2 0.8 ȹ ȹ ȹ ȹ 0.1 1.1 -1.1Ⱥ € OPTI 280 Spring 2010 Lecture 5 R. Palit, 5 MatLab Element by Element Operations ȹ 0.2 -0.3 0.2 ȹ ȹ ȹ A = ȹ -0.1 0.6 0.5 ȹ ȹ ȹ ȹ 0.3 0.6 -0.3Ⱥ ȹ 0.9 0.2 -0.1ȹ ȹ ȹ B = ȹ -0.9 0.4 -0.3ȹ ȹ ȹ ȹ 0.2 -0.5 0.8 Ⱥ •  Multiplication, .*! € C = A.* € ȹ 0.18 -0.06 -0.02ȹ ȹ ȹ = ȹ 0.09 0.24 -0.15ȹ ȹ ȹ ȹ 0.06 -0.3 -0.24Ⱥ •  Power, .^! ȹ 0.04 0.09 0.04ȹ ȹ ȹ C = A.^2 = ȹ 0.01 0.36 0.25ȹ ȹ ȹ ȹ 0.09 0.36 0.09Ⱥ •  Right Division, ./! € •  Left Division, .\! ȹ 4.5 -0.667 -0.5 ȹ ȹ ȹ C = A.\ B = ȹ 9.0 0.667 -0.6 ȹ ȹ ȹ ȹ 0.667 -0.833 -2.667Ⱥ R. Palit, 6 ȹ 0.222 -1.5 -2.0 ȹ ȹ ȹ C = A./B = ȹ 0.111 1.5 -1.667ȹ ȹ ȹ ȹ 1.5 -1.2 -0.375Ⱥ OPTI 280 Spring 2010 Lecture 5 € MatLab Element by Element Operations ȹ 0.2 -0.3 0.2 ȹ ȹ ȹ A = ȹ -0.1 0.6 0.5 ȹ ȹ ȹ ȹ 0.3 0.6 -0.3Ⱥ ȹ 0.9 + 0.5i 0.2 -0.1ȹ ȹ ȹ Z = ȹ -0.9 - 0.3i 0.4 -0.3ȹ ȹ ȹ 0.2 -0.5 + 0.1i 0.8 Ⱥ ȹ switch indices of rows and columns € •  Transpose, .’! •  Definition: € (A t ) mn = A nm ȹ 0.2 -0.1 0.3 ȹ ȹ ȹ C = A.' = ȹ -0.3 0.6 0.6 ȹ ȹ ȹ ȹ 0.2 0.5 -0.3Ⱥ € •  Conjugate Transpose, ‘! € •  Definition: switch indices of rows and columns AND take complex conjugate of the data (A ) mn = A † * nm ȹ 0.9 - 0.5i -0.9 + 0.3i ȹ 0.2 ȹ ȹ C = Z' = ȹ 0.2 0.4 -0.5 - 0.1iȹ ȹ ȹ -0.3 0.8 ȹ -0.1 Ⱥ OPTI 280 Spring 2010 Lecture 5 R. Palit, 7 € Transpose of Lena Lena - Original Lena - Transpose OPTI 280 Spring 2010 Lecture 5 R. Palit, 8 Matrix-Vector Multiplication •  An N-dimensional column vector f can be multiplied by a M x N matrix H to give a new vector g. N g m = ∑ H mn fn n =1 •  Example: € ȹ 1 2 3ȹ H = ȹ ȹ ȹ 4 5 6Ⱥ ȹ 7ȹ ȹ ȹ f = ȹ 8ȹ ȹ ȹ ȹ 9Ⱥ ȹ H11 * f1 + H12 * f2 + H13 * f3 ȹ ȹ 1 * 7 + 2 * 8 + 3 * 9 ȹ ȹ 50 ȹ g = ȹ ȹ = ȹ ȹ = ȹ ȹ €ȹ H 21 * f1 + H 22 * f2 + H 23 * f3 Ⱥ ȹ 4 * 7 + 5 * 8 + 6 * 9Ⱥ ȹ 122Ⱥ € OPTI 280 Spring 2010 Lecture 5 R. Palit, 9 Matrix Multiplication •  The MxK matrix A can be multiplied by the KxN matrix B to give the MxN matrix C. N Cmn = ∑ A mkB kn •  Example: k =1 ȹ 1 2 3ȹ A ȹ € = ȹ ȹ 4 5 6Ⱥ ȹ 7 8 9 10ȹ ȹ ȹ B = ȹ 11 12 13 14 ȹ ȹ ȹ ȹ 15 16 17 18Ⱥ ȹ 1 * 7 + 2 * 11 + 3 * 15 1 * 8 + 2 * 12 + 3 * 16 1 * 9 + 2 * 13 + 3 * 17 1 * 10 + 2 * 14 + 3 * 18 ȹ ȹ 74 80 86 92 ȹ C = ȹ ȹ = ȹ ȹ ȹ 4 * 7 + 5 * 11 + 6 * 15 4 * 8 + 5 * 12 + 6 * 16 4 * 9 + 5 * 13 + 6 * 17 4 * 10 + 5 * 14 + 6 * 18Ⱥ ȹ 173 188 203 218Ⱥ € € € OPTI 280 Spring 2010 Lecture 5 R. Palit, 10 Inverse of a Matrix •  Identity matrix •  An identity matrix can be created in MatLab using the eye() command. ȹ ȹ 1 ȹ ȹ 0 I = ȹ 0 ȹ ȹ ȹ ȹ 0 0 1 0 0 0 0 1 0 0 ȹ 0ȹ 0ȹ ȹ ȹ ȹ 1Ⱥ •  Not all matrices have inverse. Matrices that do not € have an inverse are called singular. •  The inverse of a matrix can be computed in MatLab using the inv() command. •  inv(A) ≠ A.^(-1) !!! •  Multiplication of a matrix by its inverse gives identity, A A-1 = I OPTI 280 Spring 2010 Lecture 5 R. Palit, 11 Matrix Division & Matrix Power •  Right division, /! •  A/B is roughly same as: A*inv(B) •  Exact computation: A/B = (B’\A’)’ •  Left division, \! •  A\B is roughly same as: inv(A)*B •  Exact computation: X=A\B if AX = B •  Power, ^! •  Matrix must be square. •  If N > 1 then A^N is the same as: •  A*A*…A, repeated N times. •  If N < 0 then the calculation involves eigenvectors and eigenvalues. OPTI 280 Spring 2010 Lecture 5 R. Palit, 12 Linear Equations •  Simultaneous linear equations are used in physics, optics and engineering. A r c ȹ 2 −3 1 ȹȹ x ȹ ȹ −2ȹ ȹ ȹȹ ȹ ȹ ȹ ȹ 1 −6 3 ȹȹ y ȹ = ȹ −2ȹ ȹ ȹȹ ȹ ȹ ȹ ȹ 3 € 3 −2Ⱥȹ z Ⱥ ȹ 2 Ⱥ 2 x − 3 y + z = −2 x − 6 y + 3z = −2 3 x + 3y − 2 z = 2 € Ar = c A −1Ar = A −1c r = A −1c OPTI 280 Spring 2010 Lecture 5 R. Palit, 13 Matrix Imaging Operator •  In a computer simulation of a system the imaging operator H is an MxN matrix that maps the Nx1 discrete representation of the object f to an Mx1 data vector g. Ⱥ h (1,1) h (1, 2) h (1, N ) Ⱥ Ⱥ Ⱥ H = Ⱥ Ⱥ Ⱥh ( M ,1) h ( M , 2) h ( M , N ) Ⱥ Ⱥ Ⱥ € The row index corresponds to a pixel in image space. •  •  •  row index => spatial location in image space (xd,yd) € The column index corresponds to a voxel in object space. column index => spatial location in object space (x,y,z) •  •  Values of H, h(m,n), are the response of pixel m as a result of imaging voxel n through the system. OPTI 280 Spring 2010 Lecture 5 R. Palit, 14 Matrix Imaging Operator •  One way to quantify the elements of H is to input an object f that is 1 in a single voxel and 0 elsewhere, model the physics of the system, and determine the response on the detector. This response fills one column of H. € Ⱥ0 Ⱥ Ⱥ0 Ⱥ Ⱥ h (1,1) h (1, 2) h (1, N ) Ⱥ Ⱥ Ⱥ Ⱥ Ⱥ Ⱥ H = Ⱥ Ⱥ Ⱥ Ⱥ1 Ⱥ Algorithm Ⱥ1 Ⱥ Ⱥh ( M ,1) h ( M , 2) h ( M , N ) Ⱥ Ⱥ Ⱥ Ⱥ0 Ⱥ Ⱥ Ⱥ f= g = h ( m, 2) = Ⱥ Ⱥ Ⱥ Ⱥ Ⱥ Ⱥ Ⱥ1 Ⱥ Ⱥ0 Ⱥ Nx1 Ⱥ Ⱥ Ⱥ Ⱥ € Ⱥ0 Ⱥ Mx1 OPTI 280 Spring 2010 Lecture 5 R. Palit, 15 € € Paraxial Optics •  Matrix multiplication is one way to perform paraxial ray tracing. Transfer: ω = nu y z t ȹ y ' ȹ ȹ 1 n ȹȹ y ȹ ȹȹ ȹ ȹ ȹ = ȹ ȹ ω 'Ⱥ ȹ 0 1Ⱥȹ ω Ⱥ φ € Refraction: z t ȹ 1 n ȹ T = ȹ ȹ 0 1Ⱥ ȹ € ȹ y ' ȹ ȹ 1 0ȹȹ y ȹ ȹ ȹ = ȹ ȹȹ ȹ ω 'Ⱥ ȹ −φ 1Ⱥȹ ω Ⱥ ȹ ȹ 1 0ȹ R = ȹ ȹ −φ 1Ⱥ ȹ R. Palit, 16 OPTI 280 Spring 2010 Lecture 5 € Paraxial Raytrace Through a System ȹ y ' ȹ ȹ y ȹ ȹ ȹ = Tk Rk … T3R2T2R1T ȹ ȹ 1 ȹ ω 'Ⱥ ȹ ω Ⱥ n=1.0 φ1 n=1.33 φ2 n=1.0 € z t1 t2 t3 ȹ y ' ȹ ȹ 1 t 3 ȹȹ 1 ȹ ȹ = ȹ ȹȹ ȹ ω 'Ⱥ ȹ 0 1 Ⱥȹ −φ 2 0ȹȹ 1 t 2 ȹȹ 1 0ȹȹ 1 t1 ȹȹ y ȹ ȹȹ 1.33 ȹȹ ȹ ȹȹ −φ1 1ȹȹ 0 1 ȹȹ ȹ 1Ⱥ 0 Ⱥȹ Ⱥȹ Ⱥ 1 Ⱥ ȹ € OPTI 280 Spring 2010 Lecture 5 R. Palit, 17 Linear Transformation •  Rotation of a vector •  In N-dimensions, rotation of a vector is equivalent to multiplication by a NxN matrix. Rotation about the x-axis (roll angle, α) •  In Cartesian space: ȹ x 'ȹ ȹ a11 a12 ȹ ȹ ȹ ȹ y 'ȹ = ȹ a21 a22 ȹ ȹ ȹ ȹ z' Ⱥ ȹ a31 a32 a13 ȹȹ x ȹ ȹȹ ȹ a23 ȹȹ y ȹ ȹȹ ȹ a33 Ⱥȹ z Ⱥ € v' = Rz (γ ) v ȹ 1 0 0 ȹ ȹ ȹ Rx (α ) = ȹ 0 cos(α ) − sin(α )ȹ ȹ ȹ ȹ 0 sin(α ) cos(α ) Ⱥ € Rotation Matrix € Rotation about the y-axis (pitch angle, β) ȹ cos( β) 0 sin( β) ȹ ȹ ȹ Ry ( β) = ȹ 0 1 0 ȹ ȹ ȹ ȹ − sin( β) 0 cos( β)Ⱥ Rotation about the z-axis (yaw angle, γ) € ȹ cos(γ ) − sin(γ ) 0ȹ ȹ ȹ Rz (γ ) = ȹ sin(γ ) cos(γ ) 0ȹ ȹ ȹ 0 1Ⱥ ȹ 0 R. Palit, 18 OPTI 280 Spring 2010 Lecture 5 Advanced Linear Algebra Commands ȹ 0.2 -0.3 0.2 ȹ ȹ ȹ A = ȹ -0.1 0.6 0.5 ȹ ȹ ȹ ȹ 0.3 0.6 -0.3Ⱥ ȹ 1 2 3ȹ ȹ ȹ Y = ȹ 4 5 6ȹ ȹ ȹ ȹ 7 8 9Ⱥ •  Determinant, det() € •  For a 2x2 matrix ȹ a b ȹ detȹ ȹ = ad − bc c d Ⱥ ȹ € •  For higher dimensions the determinant has a complicated analytic form but can still be computed. € •  A matrix that has a non-zero determinant is invertible. •  A matrix that has a determinant equal to zero is singular (i.e. does not have an inverse). •  Det(A) = -0.18 •  Det(Y) = 0! OPTI 280 Spring 2010 Lecture 5 R. Palit, 19 Advanced Linear Algebra Commands ȹ 0.2 -0.3 0.2 ȹ ȹ ȹ A = ȹ -0.1 0.6 0.5 ȹ ȹ ȹ ȹ 0.3 0.6 -0.3Ⱥ ȹ 1 2 3ȹ ȹ ȹ Y = ȹ 4 5 6ȹ ȹ ȹ ȹ 7 8 9Ⱥ •  Rank, rank() € •  The rank of a matrix is the number of linearly independent columns € (or rows) in the matrix. •  A matrix that has a rank that is as large as possible is said to be full rank.! •  rank(A) = 3! •  All columns are independent of each other. •  rank(Y) = 2! •  Superposition of columns 1 and 3 results in a scaled vector of column 2. •  Is either A or Y full rank? OPTI 280 Spring 2010 Lecture 5 R. Palit, 20 Very Advanced Linear Algebra Commands •  Eigenvalue Decomposition, eig()! •  If H is a square matrix, a set of vectors may exist such that Hun = λn un •  Vectors un are the eigenvectors of H and the scalar parameters λn are the eigenvalues of H. •  Singular Value Decomposition, svd()! •  If H is an MxN matrix (i.e. not square) an alternate decomposition called Singular Value Decomposition can be computed. € H Hun = µn un † HH v n = µn v n † •  Vectors un and vn are the singular values of H and the scalar parameters µn are the singular values of H. OPTI 280 Spring 2010 Lecture 5 R. Palit, 21 Very Advanced Linear Algebra Commands •  For a square matrix, eigenvalue decomposition is analogous to singular value decomposition.! •  The vectors and values from these decomposition techniques have numerous applications in imaging. •  One interesting application of SVD data is to construct an approximation of the inverse of a singular matrix. This approximation is called the pseudoinverse. R H =∑ + n =1 1 † un v n µn •  MatLab has a command to compute the pseudoinverse, pinv(). •  Computation of the pseudoinverse may be useful in € reconstructing an object from a measured image. OPTI 280 Spring 2010 Lecture 5 R. Palit, 22 round, fix, floor, ceil •  round(X) – rounds the elements of X to the nearest integers •  fix(X) – rounds the elements of X toward zero •  floor(X) – rounds the elements of X to the nearest integers less than or equal to X •  ceil(X) – rounds the elements of X to nearest integers greater than or equal to X •  clock – returns a 6 element date vector •  [year month day hour minute seconds] •  fix(clock) OPTI 280 Spring 2010 Lecture 5 R. Palit, 23 round, fix, floor, ceil round([-1.9 -0.2 1.25]) fix([-1.9 -0.2 1.25]) floor([-1.9 -0.2 1.25]) ceil([-1.9 -0.2 1.25]) [-2.0 0 1.0] [-1.0 0 1.0] [-2.0 -1.0 1.0] [-1.0 0 2.0] -2 -1 0 1 -1.9 -0.2 1.25 OPTI 280 Spring 2010 Lecture 5 R. Palit, 24 Summary & Homework •  •  •  •  •  Element by Element Operators. Matrix Operators. Applications of Matrices. MatLab Linear Algebra tools. Miscellaneous Commands. •  Read Section 20 to 23 in the Matlab Tutorial •  Read Section 4.1- 4.4 and 4.6 in the textbook (Herniter) •  Assignment 5 OPTI 280 Spring 2010 Lecture 5 R. Palit, 25 ...
View Full Document

This note was uploaded on 05/23/2010 for the course OPTI 280 taught by Professor Pau during the Spring '10 term at Arizona.

Ask a homework question - tutors are online