<jupyter_start><jupyter_text>CS217 : Assignment 1 ---Please edit the cell below to include your name and student ID  **name: <NAME>**SID: 45167563  Jupyter Notebook Tutorials We will make extensive use of Python's numerical arrays (NumPy) and interactive plotting (Matplotlib) in Jupyter notebooks for the course assignments. This first assignment is intended as a gentle warm up in case you haven't used these tools before. Start by reading through the following tutorials:https://jupyter-notebook.readthedocs.io/en/stable/notebook.htmlstarting-the-notebook-serverhttps://nbviewer.jupyter.org/github/jupyter/notebook/blob/master/docs/source/examples/Notebook/Notebook%20Basics.ipynbhttps://nbviewer.jupyter.org/github/jupyter/notebook/blob/master/docs/source/examples/Notebook/Running%20Code.ipynbThis page gives a good introduction to NumPy and many examples of using NumPy along with Matplotlib:http://www.scipy-lectures.org/intro/numpy/numpy.htmlYou should also get comfortable with searching through the documentation as neededhttps://docs.scipy.org/doc/numpy-1.13.0/reference/index.htmlhttps://matplotlib.org/api/_as_gen/matplotlib.pyplot.html  1. Thin lense approximationDerive the thin lens equation we discussed in class in a geometric way using only the constraints that (1) rays entering parallel to the lense converge on the focal point (2) the lens is symmetric as shown in the diagram below.   Prove that $\frac{1}{D} + \frac{1}{D'} = \frac{1}{f}$. You will likely want to draw your own diagram and introduce some additional variables in order to make your argument. Please add additional images and equations to the cell below as needed to support your argument.  Answer ![Len_1.png](attachment:Len_1.png) The above pair of similar triangles shows that:$\frac{y}{y^{'}} = \frac{D}{D^{'}}$ ![Len_2.png](attachment:Len_2.png) The above pair of similar triangles shows that:$\frac{S}{y^{'}} = \frac{f}{T}$As $T = D^{'} - f$, and $S = y$ from the same rectangle:$\frac{y}{y^{'}} = \frac{f}{D^{'} - f}$Therefore:$ \frac{D}{D^{'}} = \frac{f}{D^{'} - f}$$D(D^{'} - f) = D^{'}f$Divide $DD^{'}$ at both sides:$1 - \frac{f}{D^{'}} = \frac{f}{D}$$\frac{1}{f} - \frac{1}{D^{'}} = \frac{1}{D}$$\frac{1}{f} = \frac{1}{D^{'}} + \frac{1}{D}$  2. Cameras & Triangulation *** 2.1 Field of View Suppose your camera has an image sensor that is 640x480 pixels, a physical resolution of 10pixels/mm and a focal length of f=50mm. What is the horizontal field of view (in degrees)? What is the vertical field of view? Suppose you adjust the zoom on the camera, changing the focal length to 100mm. What is the new horizontal field of view?  Answer:The vertical length of image sensor is $$\frac{480}{10} = 48(mm)$$The tangent of half of the field of view is $$\frac{48}{2} \frac{1}{50} = \frac{12}{25}$$The vertical field of view is $2 arctan(\frac{12}{25})$The horizontal length of image sensor is $$\frac{640}{10} = 64(mm)$$The tangent of half of the field of view is $$\frac{64}{2} \frac{1}{50} = \frac{16}{25}$$The horizontal field of view is $2 arctan(\frac{16}{25})$Changing the focal length to $100mm$The new horizontal field of view is $2 arctan(\frac{8}{25})$<jupyter_code>import numpy as np

print("The vertical field of view is %.6f degrees." % (2 * np.arctan(12. / 25.) * 180 / np.pi))
print("The horizontal field of view is %.6f degrees." % (2 * np.arctan(16. / 25.) * 180 / np.pi))
print("The new horizontal field of view is %.6f degrees." % (2 * np.arctan(8. / 25.) * 180 / np.pi))<jupyter_output>The vertical field of view is 51.282012 degrees.
The horizontal field of view is 65.238486 degrees.
The new horizontal field of view is 35.489343 degrees.<jupyter_text>*** 2.2  Camera motionsYour camera starts out at the origin of the world coordinate system. You rotate it to the left about the y-axis by 45 degrees and then translate it right by 1 meter. Describe this motion of the camera concisely in terms of rotation matrices and translation vectors. Suppose there is a point with coordinates (1,1,1) meters in the world coordinate system. What will its coordinates be relative to the standard camera coordinate system after the camera has been moved?  Answer: ![coord.png](attachment:coord.png) Rotate the camera to the left about the y-axis by 45 degrees, the points rotate right about the y-axis by 45 degrees to the camera coordinate system. Therefore, the rotation matrix is :$$R = \left[\begin{array}{ccc} cos(\frac{\pi}{4}) & 0 & sin(\frac{\pi}{4})\\0 & 1 & 0\\-sin(\frac{\pi}{4}) & 0 & cos(\frac{\pi}{4})\end{array}\right] = \left[\begin{array}{ccc}\frac{\sqrt{2}}{2} & 0 & \frac{\sqrt{2}}{2}\\0 & 1 & 0\\-\frac{\sqrt{2}}{2} & 0 & \frac{\sqrt{2}}{2}\end{array}\right]$$Here assume the "translation to the right by 1 meter" is __examined under the world coordinates__.The translation vector is $$t = \left[\begin{array}{ccc} cos(\frac{\pi}{4}) & 0 & sin(\frac{\pi}{4})\\0 & 1 & 0\\-sin(\frac{\pi}{4}) & 0 & cos(\frac{\pi}{4})\end{array}\right]\left[\begin{array}{c} 0\\0\\1\end{array}\right] = \left[\begin{array}{c} \frac{\sqrt{2}}{2}\\0\\\frac{\sqrt{2}}{2}\end{array}\right]$$The camera coordinate is $$p = RP + Rt = R(P + t)$$When $$P = \left[\begin{array}{c} 1\\1\\1\end{array}\right]$$The camera coordinate $$p = \left[\begin{array}{ccc} cos(\frac{\pi}{4}) & 0 & sin(\frac{\pi}{4})\\0 & 1 & 0\\-sin(\frac{\pi}{4}) & 0 & cos(\frac{\pi}{4})\end{array}\right] \left[\begin{array}{c} 1\\1\\1\end{array}\right] + \left[\begin{array}{c} \frac{\sqrt{2}}{2}\\0\\\frac{\sqrt{2}}{2}\end{array}\right]$$$$p = \left[\begin{array}{c} \frac{3\sqrt{2}}{2}\\1\\\frac{\sqrt{2}}{2}\end{array}\right]$$  2.3 Estimate the extrinsic parameters for your eyesWork relative to a world coordinate system centered on the bridge of your nose with the z-axis pointed out of your face and x-axis to the right. Since eyes rotate in their sockets depending on how far away you focus (this phenomenon is called vergence) assume you are looking at a computer monitor that is 40 cm away from your face. You may find it interesting to find some red/blue glasses and modify the projection code to render anaglyphs of various 3D shapes overlaying the left and right eye views in red and blue. The following parameters should be passed into the project function mapping 3D point on computer monitor to perceived image of eyes.The distance between two eyes, $2b$: the translation matrix for left eye will be $[-b, 0, 0]$, the translation matrix for right eye will be $[b, 0, 0]$.The distance between face and the actual object, $z$ (given by 40cm): the rotated angle of right camera will be $-\theta = arctan(\frac{z}{b}) - \frac{\pi}{2}$, the rotated angle of left camera will be $\theta = \frac{\pi}{2} - arctan(\frac{z}{b})$.the right rotation matrix will be $$R_{right} = \left[\begin{array}{ccc}cos(\theta) & -sin(\theta) & 0 \\sin(\theta) & cos(\theta) & 0 \\0 & 0 & 1\end{array}\right]$$the left rotation matrix will be $$R_{left} = \left[\begin{array}{ccc}cos(\theta) & sin(\theta) & 0 \\-sin(\theta) & cos(\theta) & 0 \\0 & 0 & 1\end{array}\right]$$The focal length of both eyes $fL, fR$.The difference in height of both eyes to computer monitor focus point, so that the eye offset to world coordinate system can bu computed.  2.4 TriangulationWrite a function called **triangulate** that takes the coordinates of points in two images along with the camera parameters and returns the 3D coordinates of the points in world coordinates.  We discussed a couple different approaches to trianguation in class.  You are welcome to use whichever you prefer.  The provided functions **project** and **test_triangulation** outline a reasonable way to represent the camera parameters and should be useful in experimenting with your code.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
class Camera:
    """
    A simple data structure describing camera parameters 
    
    The parameters describing the camera
    cam.f : float   --- camera focal length (in units of pixels)
    cam.c : 2x1 vector  --- offset of principle point
    cam.R : 3x3 matrix --- camera rotation
    cam.t : 3x1 vector --- camera translation 

    
    """    
    def __init__(self,f,c,R,t):
        self.f = f
        self.c = c
        self.R = R
        self.t = t

        
    def project(self,pts3):
        """
        Project the given 3D points in world coordinates into the specified camera    

        Parameters
        ----------
        pts : 2D numpy.array (dtype=float)
            Coordinates of N points stored in a array of shape (3,N)

        Returns
        -------
        ipts : 2D numpy.array (dtype=float)
            Image coordinates of N points stored in a array of shape (N,2)

        """

        assert(pts3.shape[0]==3)

        # get point location relative to camera
        pcam = self.R.transpose() @ (pts3 - self.t)
         
        # project
        p = self.f * (pcam / pcam[2,:])
        
        # offset principal point
        pts2 = p[0:2,:] + self.c
        
        assert(pts2.shape[1]==pts3.shape[1])
        assert(pts2.shape[0]==2)
    
        return pts2
    
def generate_hemisphere(radius,center,npts):
    """
    Generate a set of 3D points which are randomly distributed on the
    surface of a hemisphere
    
    
    Parameters
    ----------
    radius : float
        Hemisphere radius

    center : numpy.array (dtype=float)
        3x1 vector specifying the center of the hemisphere
        
    npts : int
        number of points to generate

    Returns
    -------
    x : 2D numpy.array (dtype=float)
        3xnpts array containing coordinates of the points
        
    """
    
    assert(center.shape==(3,1))
    
    
    #generate randomly distributed points
    x = np.random.standard_normal((3,npts))
    
    #scale points to the surface of a sphere with given radius
    nx = np.sqrt(np.sum(x*x,axis=0))
    x = radius * x / nx
    
    # make points with positive z-coordinates negative
    # so that points are all on a half-sphere 
    x[2,:] = -np.abs(x[2,:])
    
    # translate to desired position
    x = x + center
    
    return x


# demo the generate_hemisphere function
if False:
    x = generate_hemisphere(1,np.array([[0,0,0],]).T,2000)
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(x[0,:],x[1,:],x[2,:],'.')
    ax.axis('square')
    ax.set_zlim(-1,1)
def visualize(a, b):
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1,projection='3d')
    ax.plot(a[:, 0], a[:, 1], a[:, 2],'ro',fillstyle='none')
    ax.plot(b[:, 0], b[:, 1], b[:, 2],'ro',fillstyle='none')
    plt.axis('equal')
    plt.xlabel('x')
    plt.ylabel('y')

def triangulate(pts2L,camL,pts2R,camR):
    """
    Triangulate the set of points seen at location pts2L / pts2R in the
    corresponding pair of cameras. Return the 3D coordinates relative
    to this global coordinate system


    Parameters
    ----------
    pts2L : 2D numpy.array (dtype=float)
        Coordinates of N points stored in a array of shape (2, N) seen from camL camera

    pts2R : 2D numpy.array (dtype=float)
        Coordinates of N points stored in a array of shape (2, N) seen from camR camera

    camL : Camera
        The first "left" camera view

    camR : Camera
        The second "right" camera view

    Returns
    -------
    pts3 : 2D numpy.array (dtype=float)
        array containing 3D coordinates of the points in global coordinates

    """
    N = pts2L.shape[1]
    invRL = np.linalg.inv(camL.R)
    R = invRL.dot(camR.R)
    t = invRL.dot(camR.t - camL.t)
    
    pL = np.concatenate((pts2L - camL.c, np.full((1, N), camL.f)), axis=0) / camL.f
    pR = np.concatenate((pts2R - camR.c, np.full((1, N), camR.f)), axis=0) / camR.f
    
    camera = np.concatenate((pL[:, np.newaxis, :], (-R.dot(pR))[:, np.newaxis, :]), axis=1)
    for i in range(N):
        z = np.linalg.lstsq(camera[:, :, i], t)[0]
        pL[:, i] = z[0] * pL[:, i]
        pR[:, i] = z[1] * pR[:, i]

    PL = camL.R.dot(pL) + camL.t
    PR = camR.R.dot(pR) + camR.t

    pts = (PL + PR) / 2
    return pts
%matplotlib notebook
#
# test your camera and triangulate function
#


# create a rotation matrix representing rotation around y-axis by amount theta
def roty(theta):
    st = np.sin(theta)
    ct = np.cos(theta)
    R = np.array([[ct,0,st],[0,1,0],[-st,0,ct]])
    return R


#compute rotation angle so that the camera is looking directly at the sphere

b = 5  #baseline between cameras
d = 10   #distance to object
theta = np.arctan(b/d)  
tL = np.array([[-b,0,0]]).T
tR = np.array([[b,0,0]]).T
camL = Camera(f=100,c=np.array([[50,50]]).T,t=tL,R=roty(theta))
camR = Camera(f=100,c=np.array([[50,50]]).T,t=tR,R=roty(-theta))

#generate 3D points
pts3 = generate_hemisphere(2, np.array([[0,0,d]]).T,500)

#project into each camera
pts2L = camL.project(pts3)
pts2R = camR.project(pts3)

#triangulate to recover 3d position
pts3t = triangulate(pts2L, camL, pts2R, camR)



#
# visualize results
#

# generate coordinates of a line segment running from the center
# of the camera to 3 units in front of the camera
lookL = np.hstack((tL,tL+camL.R @ np.array([[0,0,3]]).T))
lookR = np.hstack((tR,tR+camR.R @ np.array([[0,0,3]]).T))

# visualize the left and right image overlaid
fig = plt.figure()
ax = fig.add_subplot(2,2,1)
ax.plot(pts2L[0,:],pts2L[1,:],'b.')
ax.plot(pts2R[0,:],pts2R[1,:],'r.')
plt.axis('equal')

#visualize 3D layout of points, camera positions
# and the direction the camera is pointing
ax = fig.add_subplot(2,2,2,projection='3d')
ax.plot(pts3[0,:],pts3[1,:],pts3[2,:],'.')
ax.plot(tR[0],tR[1],tR[2],'ro')
ax.plot(tL[0],tL[1],tL[2],'bo')
ax.plot(lookL[0,:],lookL[1,:],lookL[2,:],'b')
ax.plot(lookR[0,:],lookR[1,:],lookR[2,:],'r')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')

# overhead view showing points, camera
# positions, and direction camera is pointed
ax = fig.add_subplot(2,2,3)
ax.plot(pts3[0,:],pts3[2,:],'.')
ax.plot(tL[0],tL[2],'bo')
ax.plot(lookL[0,:],lookL[2,:],'b')
ax.plot(tR[0],tR[2],'ro')
ax.plot(lookR[0,:],lookR[2,:],'r')
plt.axis('equal')
plt.grid()
plt.xlabel('x')
plt.ylabel('z')

# compare reconstruction
ax = fig.add_subplot(2,2,4,projection='3d')
ax.plot(pts3[0,:],pts3[1,:],pts3[2,:],'b.')
ax.plot(pts3t[0,:],pts3t[1,:],pts3t[2,:],'ro',fillstyle='none')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')<jupyter_output>//anaconda/envs/py3/lib/python3.6/site-packages/ipykernel_launcher.py:47: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.<jupyter_text>2.5 How sensitive is the 3D reconstruction to errors in the 2D point locations? You can approach this question empirically by generating a pair of test images and then adding random Gaussian noise to the 2D point locations. Plot out the error in the 3D recovery as a function of the noise in the 2D point locations for several different levels of noise. You can measure error as the sum of squared distances between the true locations and the recovered noisy locations.  Insert additional cells as needed below in order to carry out your experiments and describe your result.<jupyter_code>%matplotlib inline
x = np.arange(0, 1, 0.1).tolist()
y = []
fig = plt.figure(figsize=(10,10))
M = 5
for mu in range(M):
    for sigma in x:
        error2D = np.random.normal(mu, sigma, (pts2L.shape))
        pts2Le = pts2L + error2D
        pts2Re = pts2R + error2D
        pts3te = triangulate(pts2Le,camL,pts2Re,camR)
        error3D = np.sum((pts3 - pts3te)**2)
        y.append(error3D)
    ax = fig.add_subplot(M,1,mu + 1)
    plt.plot(x, y)
    plt.xlabel('sigma')
    plt.ylabel('error')
    plt.title("mu = %.2f" % (mu))
    y = []
fig.tight_layout()<jupyter_output>//anaconda/envs/py3/lib/python3.6/site-packages/ipykernel_launcher.py:47: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.<jupyter_text>2.6 How sensitive is the 3D reconstruction procedure to errors in the camera parameters? Compute reconstructions of the hemisphere where you vary (a) the image center (+/-10 pixels), (b) the focal length (+/-10 pixels), (c) the translation vector (+/-100mm), and (d) the rotation angle (+/-5 degrees). Vary each parameter of one of the cameras passed into triangulate while keeping the other camera and the images (xL,xR) fixed. Report in a plot or table, the reconstruction error as a function of these different parameters. What parameter is the reconstruction result quality most dependent on? Discuss how you might frame this question in a more general way that doesn’t depend on the units used to measure the camera parameters.  Insert additional cells as needed below in order to carry out your experiments and describe your result.<jupyter_code>def cmpL(wrong_cam, pts3, results, pts2L, pts2R, camR):
    wrong_pts3t = triangulate(pts2L, wrong_cam, pts2R, camR)
    tmp = np.sum((pts3 - wrong_pts3t)**2)
    print("Change Left: %.4f" % tmp)
    results.append(tmp)
    return results

def cmpR(wrong_cam, pts3, results, pts2L, pts2R, camL):
    wrong_pts3t = triangulate(pts2L, camL, pts2R, wrong_cam)
    tmp = np.sum((pts3 - wrong_pts3t)**2)
    print("Change Right: %.4f" % tmp)
    results.append(tmp)
    return results

def cmpBoth(wrong_camL, wrong_camR, pts3, results, pts2L, pts2R):
    wrong_pts3t = triangulate(pts2L, wrong_camL, pts2R, wrong_camR)
    tmp = np.sum((pts3 - wrong_pts3t)**2)
    print("Change Both: %.4f" % tmp)
    results.append(tmp)
    return results

resultL = []
resultR = []
resultB = []

#0
print(0)
wrong_camL = Camera(f=100, c=np.array([[60, 50]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[60, 50]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#1
print(1)
wrong_camL = Camera(f=100, c=np.array([[50, 60]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 60]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#2
print(2)
wrong_camL = Camera(f=100, c=np.array([[40, 50]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[40, 50]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#3
print(3)
wrong_camL = Camera(f=100, c=np.array([[50, 40]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 40]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#4
print(4)
wrong_camL = Camera(f=110, c=np.array([[50, 50]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=110, c=np.array([[50, 50]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#5
print(5)
wrong_camL = Camera(f=90, c=np.array([[50, 50]]).T, t=tL, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=90, c=np.array([[50, 50]]).T, t=tR, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#6
print(6)
wrong_camL = Camera(f=100, c=np.array([[50, 50]]).T, t=tL + np.array([[-100,0,0]]).T, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 50]]).T, t=tR + np.array([[100,0,0]]).T, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#7
print(7)
wrong_camL = Camera(f=100, c=np.array([[50, 50]]).T, t=tL + np.array([[100,0,0]]).T, R=roty(theta))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 50]]).T, t=tR + np.array([[-100,0,0]]).T, R=roty(-theta))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#8
print(8)
wrong_camL = Camera(f=100, c=np.array([[50, 50]]).T, t=tL, R=roty(theta + 5 * np.pi / 180))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 50]]).T, t=tR, R=roty(-theta - 5 * np.pi / 180))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
#9
print(9)
wrong_camL = Camera(f=100, c=np.array([[50, 50]]).T, t=tL, R=roty(theta - 5 * np.pi / 180))
resultL = cmpL(wrong_camL, pts3, resultL, pts2L, pts2R, camR)
wrong_camR = Camera(f=100, c=np.array([[50, 50]]).T, t=tR, R=roty(-theta + 5 * np.pi / 180))
resultR = cmpR(wrong_camR, pts3, resultR, pts2L, pts2R, camL)
resultB = cmpBoth(wrong_camL, wrong_camR, pts3, resultB, pts2L, pts2R)
fig = plt.figure(figsize=(10,10))
x = np.arange(len(resultL))
plt.title("Change Left Camera")
plt.plot(x, resultL)
plt.show()
fig = plt.figure(figsize=(10,10))
x = np.arange(len(resultR))
plt.title("Change Right Camera")
plt.plot(x, resultR)
plt.show()
fig = plt.figure(figsize=(10,10))
x = np.arange(len(resultB))
plt.title("Change Both Camera")
plt.plot(x, resultB)
plt.show()<jupyter_output>0
Change Left: 863.7305
Change Right: 661.0927
Change Both: 676.1580
1
Change Left: 134.0450
Change Right: 135.4428
Change Both: 531.0574
2<jupyter_text>The camera parameter that the reconstruction quality most depend on is translation vector, for which it has most error. The most stable camera parameter is the focal length.To frame the question independent from the camera parameter units, simply provide the intrinsic matrix of $K$ that contains the pixel magnification factors, so that pysical distance can be converted to pixels and vice versa.  3. Camera CalibrationWrite a function **calibrate** that takes the coordinates of points in an image and in 3D and resturns estimated camera parameters.  Your function should first use the linear technique to estimate the camera matrix and then decompose the matrix into intrinsic and extrinsic components using RQ matrix decomposition. Your function should output the same data structure used previously to represent a camera.<jupyter_code>def calibrate(pts2, pts3):
    """
    Give a set of 3D points and their correpsonding 2D image
    coordinates, estimate the parameters of the camera.

    Parameters
    ----------
    pts2 : 2D numpy.array (dtype=float)
        Coordinates of N points stored in a array of shape (N,2) seen from camL camera
        
    pts3 : 2D numpy.array (dtype=float)
        array containing 3D coordinates of the points relative to this camera's origin

    Returns
    -------
     cam : Camera
        The estimated camera parameters

    """
    N = pts2.shape[1]

    y = np.concatenate((np.full((1, N), 0), np.full((1, N), -1)), axis=0)
    y = np.concatenate((y, pts2[1, :][np.newaxis, :]), axis = 0).T
    x = np.concatenate((np.full((1, N), 1), np.full((1, N), 0)), axis=0)
    x = np.concatenate((x, -pts2[0, :][np.newaxis, :]), axis = 0).T
    p3 = np.concatenate((pts3, np.full((1, N), 1)), axis=0).T[:, np.newaxis, :]
    p3 = np.repeat(p3, 3, axis=1)
    p3 = np.repeat(p3, 2, axis=0)
    p = np.ones(p3.shape)
    p[::2, :, :] = np.multiply(p[::2, :, :], y[:, :, np.newaxis])
    p[1::2, :, :] = np.multiply(p[1::2, :, :], x[:, :, np.newaxis])
    A = np.multiply(p3, p).reshape(2 * N, 12)
    U, S, Vh = np.linalg.svd(A)
    M = Vh[np.argmin(S), :].reshape(3, 4)
    C1 = M[:, :-1]
    R, K = np.linalg.qr(np.linalg.inv(C1))
    K = np.linalg.inv(K)
    t = np.dot(np.linalg.inv(C1), M[:, -1][:, np.newaxis])
    K = K / K[-1, -1]
    D = np.diag(np.sign(np.diag(K)))
    K = np.dot(K, D)
    R = np.dot(D, R)
    t = np.dot(D, t)
    
    t = np.dot(np.linalg.det(R), t)
    R = np.dot(np.linalg.det(R), R).T
    return Camera(f=K[0, 0], c=np.array([[K[0, -1], K[1, -1]]]).T, t=t, R=R)<jupyter_output><empty_output><jupyter_text>3.1 Synthetic calibration testTest your method on a synthetic dataset using the cameras and points from the previous problem on triangulation. Attempt to recover camL and camR and compare your recovered values to the true parameters used to synthesize the test image. If you add Gaussian noise to the 2D point locations (but no the 3D locations) which camera parameter(s) are most affected?  Which are most stable? Please justify your answer with some simple experiment.<jupyter_code>def cmp(pts2, pts3, real_cam):
    cam = calibrate(pts2, pts3)
    print("Difference in focal length: %.4f" % (cam.f - real_cam.f))
    print("Distance in camera center: %.4f" % np.sqrt(np.sum((cam.c - real_cam.c)**2)))
    print("Difference in translation vector: \n", cam.t - real_cam.t)
    print("Difference in rotation matrix: \n", cam.R - real_cam.R)

print("Left Camera")
cmp(pts2L, pts3, camL)
print("Right Camera")
cmp(pts2R, pts3, camR)

# x = np.arange(0, 1, 0.1).tolist()
# for mu in range(3, 5):
#     for sigma in x:
#         error2D = np.random.normal(mu, sigma, (pts2L.shape))
#         pts2Le = pts2L + error2D
#         pts2Re = pts2R + error2D
#         print("Noise mu = %.4f, sigma = %.4f" % (mu, sigma))
#         print("Left Camera")
#         cmp(pts2Le, pts3, camL)
#         print("Right Camera")
#         cmp(pts2Re, pts3, camR)<jupyter_output>Left Camera
[[-1.00000000e+02 -2.21900569e-12  5.00000000e+01]
 [-0.00000000e+00 -1.00000000e+02  5.00000000e+01]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
Difference in focal length: 0.0000
Distance in camera center: 0.0000
Difference in translation vector: 
 [[-2.40696352e-13]
 [-2.57571742e-14]
 [-1.50546242e-13]]
Difference in rotation matrix: 
 [[ 1.44328993e-15  3.52495810e-15 -1.94289029e-15]
 [-1.08369572e-14  2.22044605e-16  1.37918684e-14]
 [ 1.94289029e-15 -1.71822567e-14  1.22124533e-15]]
Right Camera
[[-1.00000000e+02  2.57950417e-12  5.00000000e+01]
 [-0.00000000e+00 -1.00000000e+02  5.00000000e+01]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
Difference in focal length: -0.0000
Distance in camera center: 0.0000
Difference in translation vector: 
 [[-1.77635684e-14]
 [-6.43929354e-15]
 [ 1.79634085e-13]]
Difference in rotation matrix: 
 [[ 1.76525461e-14  1.93733918e-14  3.49720253e-14]
 [-1.17917136e-14  0.00000000e+00  1.97367937e-14]
 [-3.49720253e-14 [...]<jupyter_text>The focal length is the most affected parameter in calibration.The rotation matrix is the most stable parameter in calibration.  3.2 Planar calibration testTest your calibration method on one of the calibration images provided on the course website. The corner points of the checkerboard have physical coordinates given by a grid with 8 vertices along one axis and 6 along the other axis with a spacing of 2.8cm.What intrinsic camera parameters do you recover? Are the camera pixels square? What ambiguities arise from calibrating using a single image of a planar calibration object? To answer this question you may need to consult one of the recommended textbooks for the class or other online resources.<jupyter_code># This code is a modification of the following tutorial:
# https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_calib3d/py_calibration/py_calibration.html
%matplotlib inline
import cv2
import os

nx = 8
ny = 6
objp = np.zeros((nx * ny, 3), np.float32)
objp[:, :2] = np.mgrid[0:nx, 0:ny].T.reshape(-1, 2)
objp[:, :2] *= 2.8
objp[:, -1] += 1

objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
dirname = 'data/calib'
for a, b, c in os.walk(dirname):
    for fname in c:
        img = cv2.imread(os.path.join(dirname, fname))
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        ret, corners = cv2.findChessboardCorners(gray, (nx, ny), None)
        if ret == True:
            objpoints.append(objp)
            corners2 = cv2.cornerSubPix(gray,corners, (11, 11), (-1, -1), criteria)
            imgpoints.append(corners2)
            # Draw and display the corners
            img = cv2.drawChessboardCorners(img, (nx, ny), corners2,ret)
            plt.imshow(img)
            plt.show()
        break
    break
objpoints = np.array(objpoints)
imgpoints = np.array(imgpoints)
objpoints = objpoints[0].T
imgpoints = imgpoints[0].reshape(imgpoints.shape[1], 2).T
camera = calibrate(imgpoints, objpoints)
print(camera.c)
print(camera.R)
print(camera.f)
print(camera.t)<jupyter_output><empty_output><jupyter_text>![plane.png](attachment:plane.png)Textbook 22.1.2If all calibration points lie on the same plane, the camera center has infinite many solutions lying on a line except for the intersection point of the object plane, with image-equivalent 2D projection results, which means that the distance from camera to the object is lost due to normalization.  3.3 Geometric recovery challengeIn the course files, you will find a pair of images of an object with unknown geometry taken with the same camera as the calibration images. Your goal is to estimate the dimensions of the imaged box.1. You will want to get a good estimate of the intrinsic camera parameters.  Since our simple **calibrate** function doesn't work well for planar objects, I recommend using an existing calibration tool.  One option for python is to use the opencv calibration functionality  (see links and examples below).  Another alternative is to use the MATLAB camera calibration toolbox and then import your estimates back into this notebook.  2. For the pair of images, you need to estimate the extrinsic camera parameters (R, t) that relate the camera views. You can do this using the locations of the outside corners of the checkerboard and known intrinsic parameters and peform a factorization of the Essential matrix, or alternately using non-linear optimization (e.g. using **scipy.optimize.least_squares**) to directly search for R,t which projects the points to the observed location. 3. Finally, use your triangulation routine to estimate the dimensions of the object imaged. You can select correspondences in the two images to triangulate by hand. There will be a prize for the person whose estimate the box geometry is closest to the true solution. Please document your approach and include any code you used with your homework submission.<jupyter_code># This code is a modification of the following tutorial:
# https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_calib3d/py_calibration/py_calibration.html
%matplotlib notebook
import cv2
import os
import selectpoints

nx = 8
ny = 6

objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
dirname = 'data/calib'
fname = '13.jpg'
img = cv2.imread(os.path.join(dirname, fname))
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(img)

#selectpoints will show the axis and wait for the specified 
# number of user clicks.
k = 4
spoints = selectpoints.select_k_points(ax,k)<jupyter_output><empty_output><jupyter_text>--- Notes and resources: CalibrationOpenCV provides a set of tools for camera calibration and even a demo script which should work with a bit of modification.  You will need to install the package (e.g., by running ***conda install opencv***) in order to use it in your notebook.  Consult these references to get some idea of the functionality.  The code snippet below demonstrates automatic detection of the chessboard corners in the image.https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_calib3d/py_calibration/py_calibration.htmlhttps://github.com/opencv/opencv/blob/master/samples/python/calibrate.py<jupyter_code>import cv2
nx = 8
ny = 6

# Make a list of calibration images
fname = 'data/calib/c1.jpg'
img = cv2.imread(fname)
plt.imshow(img)

# Convert to grayscale
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

# Find the chessboard corners
ret, corners = cv2.findChessboardCorners(gray, (nx, ny), None)

# If found, draw corners
if ret == True:
    # Draw and display the corners
    cv2.drawChessboardCorners(img, (nx, ny), corners, ret)
    plt.imshow(img)
    plt.show()<jupyter_output><empty_output><jupyter_text>Interactive Point SelectionFor interactively selecting points in an image by hand (e.g. clicking on the corners of the box) I have provided a bit of code **selectpoints.py** to assist.  Note: this will only work with the **%matplotlib notebook** directive which gives interactive plots<jupyter_code>import selectpoints

image = plt.imread('data/challenge/13.jpg')

#display the image in the notebook using a grayscale colormap
# here we keep track of the image axis variable ax, which we
# need to pass along to the select_k_points function
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(image)

#selectpoints will show the axis and wait for the specified 
# number of user clicks.
k = 4
spoints = selectpoints.select_k_points(ax,k)
# once you have finished clicking, you can then execute this cell to get
# the point coordinates

# get the x,y coordinates of the selected points
xp = spoints.xs
yp = spoints.ys

print((xp,yp))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13833446015461948, 'fraction_numerical': 0.04585672712307828, 'mean_word_length': 3.9758666081614744, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 9, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 17, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '531680', 'license': 'MIT', 'n_tokens_mistral': 12952, 'n_tokens_neox': 11961, 'n_words': 3773, 'path': 'HW1/assignment1.ipynb', 'repo_name': 'SIGUSRz/LightGeo', 'size': '34021'}