2020.10.12-2020.10.18工作周报

摘要

本周主要做了一下四项工作：

• 跑鱼目论文的代码
• 吴恩达深度学习课程第三周
• 《视觉slam十四讲》第三讲、第四讲
• 读MICCAI上那篇手术机器人论文

一、读鱼目论文的代码，看能读懂多少

环境配好，跑起来了。过程中的三张截图如下：

运行完成之后，窗口全部关闭，生成了一个MKFTrajectory.txt文档，不知道是用来干嘛的。Readme下面提示说：Important: To evaluate the trajectory you will need to transform the result into the rigid body coordinate system. The transformation matrix MCS_to_Rigid_body can be found here. 下载了那个转换矩阵之后，还没明白怎么用的，不知道怎么把这个结果转换到rigid body corrdinate system。。

二、吴恩达的深度学习课程第三周

第三周的课程主要讲了有一层隐藏层的人工神经网络。比第二周Logistic回归更加复杂。实现了具有一个隐藏层的平面数据分类的神经网络。
代码整理如下：

1. Packages

#package imports
import numpy as np
import matplotlib.pyplot as plt
from testCases_v2 import *  # 这是Coursera平台自带的
import sklearn
import sklearn.datasets
import sklearn.linear_model
from planar_utils import plot_decision_boundary, sigmoid, load_planar_dataset,load_extra_datasets  # 这些也是Coursera自带的

%matplotlib inline

np.random.seed(1) # set a seed so that the results are consistent

2. Dataset

First, let’s get the dataset you will work on. The following code will load a “flower” 2-class dataset into variables X and Y.

X, Y = load_planar_dataset()

将这个数据集可视化：

# visualize the data
plt.scatter(X[0,:], X[1,:], c=Y, s=40, cmap=plt.com.Spectral);

3. Simple Logistic Regression

在建立一个完整的神经网络之前，我们先看看使用第二周的Logsitic回归在这个问题上的表现怎么样。

精确度为47%

4. 神经网络模型

我们的模型如下：

相关的数学公式：

注意：通常构建一个神经网络的方法步骤如下：

1. 定义神经网络的结构（输入层单元，隐藏层单元等）
2. 初始化模型的参数
3. 循环迭代：
   • 实现前向传播算法
   • 计算损失
   • 实现反向传播算法获得梯度
   • 更新参数（梯度下降算法）

4.1 定义神经网络的结构

# GRADED FUNCTION: layer_sizes

def layer_sizes(X, Y):
    """
    Arguments:
    X -- input dataset of shape (input size, number of examples)
    Y -- labels of shape (output size, number of examples)
    
    Returns:
    n_x -- the size of the input layer
    n_h -- the size of the hidden layer
    n_y -- the size of the output layer
    """
    ### START CODE HERE ### (≈ 3 lines of code)
    n_x = X.shape[0] # size of input layer
    n_h = 4  # 默认隐藏层4个单元
    n_y = Y.shape[0] # size of output layer
    ### END CODE HERE ###
    return (n_x, n_h, n_y)

4.2 初始化模型的参数

# GRADED FUNCTION: initialize_parameters

def initialize_parameters(n_x, n_h, n_y):
    """
    Argument:
    n_x -- size of the input layer
    n_h -- size of the hidden layer
    n_y -- size of the output layer
    
    Returns:
    params -- python dictionary containing your parameters:
                    W1 -- weight matrix of shape (n_h, n_x)
                    b1 -- bias vector of shape (n_h, 1)
                    W2 -- weight matrix of shape (n_y, n_h)
                    b2 -- bias vector of shape (n_y, 1)
    """
    
    np.random.seed(2) # we set up a seed so that your output matches ours although the initialization is random.
    
    ### START CODE HERE ### (≈ 4 lines of code)
    W1 = np.random.randn(n_h, n_x) * 0.01
    b1 = np.zeros((n_h, 1))
    W2 = np.random.randn(n_y, n_h) * 0.01
    b2 = np.zeros((n_y, 1))
    ### END CODE HERE ###
    
    assert (W1.shape == (n_h, n_x))
    assert (b1.shape == (n_h, 1))
    assert (W2.shape == (n_y, n_h))
    assert (b2.shape == (n_y, 1))
    
    parameters = {"W1": W1,
                  "b1": b1,
                  "W2": W2,
                  "b2": b2}
    
    return parameters

4.3 循环迭代

前向传播

# GRADED FUNCTION: forward_propagation

def forward_propagation(X, parameters):
    """
    Argument:
    X -- input data of size (n_x, m)
    parameters -- python dictionary containing your parameters (output of initialization function)
    
    Returns:
    A2 -- The sigmoid output of the second activation
    cache -- a dictionary containing "Z1", "A1", "Z2" and "A2"
    """
    # Retrieve each parameter from the dictionary "parameters"
    ### START CODE HERE ### (≈ 4 lines of code)
    W1 = parameters["W1"]
    b1 = parameters["b1"]
    W2 = parameters["W2"]
    b2 = parameters["b2"]
    ### END CODE HERE ###
    
    # Implement Forward Propagation to calculate A2 (probabilities)
    ### test
    #print(W1.shape)
    #print(X.shape)
    #print(b1.shape)
    #print(W2.shape)
    #print(b2.shape)
    
    ### START CODE HERE ### (≈ 4 lines of code)
    Z1 = np.matmul(W1, X) + b1  # W1: (n_h, n_x) X: (n_x, m)   Z1:(n_h, m)
    A1 = np.tanh(Z1) # A1: (n_h, m)
    Z2 = np.matmul(W2, A1) + b2 # W2:(n_y, n_h)  b2:(n_y, 1)   Z2: (n_y, m)
    A2 = sigmoid(Z2)   # A2:(n_y, m)
    ### END CODE HERE ###
    
    assert(A2.shape == (1, X.shape[1]))
    
    cache = {"Z1": Z1,
             "A1": A1,
             "Z2": Z2,
             "A2": A2}
    
    return A2, cache

计算loss

# GRADED FUNCTION: compute_cost

def compute_cost(A2, Y, parameters):
    """
    Computes the cross-entropy cost given in equation (13)
    
    Arguments:
    A2 -- The sigmoid output of the second activation, of shape (1, number of examples)
    Y -- "true" labels vector of shape (1, number of examples)
    parameters -- python dictionary containing your parameters W1, b1, W2 and b2
    [Note that the parameters argument is not used in this function, 
    but the auto-grader currently expects this parameter.
    Future version of this notebook will fix both the notebook 
    and the auto-grader so that `parameters` is not needed.
    For now, please include `parameters` in the function signature,
    and also when invoking this function.]
    
    Returns:
    cost -- cross-entropy cost given equation (13)
    
    """
    
    m = Y.shape[1] # number of example

    # Compute the cross-entropy cost
    ### START CODE HERE ### (≈ 2 lines of code)
    logprobs = np.multiply(np.log(A2),Y) + np.multiply(np.log(1-A2), 1-Y)
    cost = - (1/m) * np.sum(logprobs)
    ### END CODE HERE ###
    
    cost = float(np.squeeze(cost))  # makes sure cost is the dimension we expect. 
                                    # E.g., turns [[17]] into 17 
    assert(isinstance(cost, float))
    
    return cost

反向传播

# GRADED FUNCTION: backward_propagation

def backward_propagation(parameters, cache, X, Y):
    """
    Implement the backward propagation using the instructions above.
    
    Arguments:
    parameters -- python dictionary containing our parameters 
    cache -- a dictionary containing "Z1", "A1", "Z2" and "A2".
    X -- input data of shape (2, number of examples)
    Y -- "true" labels vector of shape (1, number of examples)
    
    Returns:
    grads -- python dictionary containing your gradients with respect to different parameters
    """
    m = X.shape[1]
    
    # First, retrieve W1 and W2 from the dictionary "parameters".
    ### START CODE HERE ### (≈ 2 lines of code)
    W1 = parameters["W1"]
    W2 = parameters["W2"]
    ### END CODE HERE ###
        
    # Retrieve also A1 and A2 from dictionary "cache".
    ### START CODE HERE ### (≈ 2 lines of code)
    A1 = cache["A1"]
    A2 = cache["A2"]
    ### END CODE HERE ###
    
    # Backward propagation: calculate dW1, db1, dW2, db2. 
    ### START CODE HERE ### (≈ 6 lines of code, corresponding to 6 equations on slide above)
    dZ2 = A2-Y
    dW2 = (1/m)*np.dot(dZ2,A1.T)
    db2 = (1/m)*np.sum(dZ2, axis=1, keepdims=True)
    dZ1 = np.dot(W2.T, dZ2) * (1 - np.power(A1, 2))
    dW1 = (1/m)*np.dot(dZ1, X.T)
    db1 = (1/m)*np.sum(dZ1, axis=1, keepdims=True)
    ### END CODE HERE ###
    
    grads = {"dW1": dW1,
             "db1": db1,
             "dW2": dW2,
             "db2": db2}
    
    return grads

更新参数

# GRADED FUNCTION: update_parameters

def update_parameters(parameters, grads, learning_rate = 1.2):
    """
    Updates parameters using the gradient descent update rule given above
    
    Arguments:
    parameters -- python dictionary containing your parameters 
    grads -- python dictionary containing your gradients 
    
    Returns:
    parameters -- python dictionary containing your updated parameters 
    """
    # Retrieve each parameter from the dictionary "parameters"
    ### START CODE HERE ### (≈ 4 lines of code)
    W1 = parameters["W1"]
    b1 = parameters["b1"]
    W2 = parameters["W2"]
    b2 = parameters["b2"]
    ### END CODE HERE ###
    
    # Retrieve each gradient from the dictionary "grads"
    ### START CODE HERE ### (≈ 4 lines of code)
    dW1 = grads["dW1"]
    db1 = grads["db1"]
    dW2 = grads["dW2"]
    db2 = grads["db2"]
    ## END CODE HERE ###
    
    # Update rule for each parameter
    ### START CODE HERE ### (≈ 4 lines of code)
    W1 = W1 - learning_rate * dW1
    b1 = b1 - learning_rate * db1
    W2 = W2 - learning_rate * dW2
    b2 = b2 - learning_rate * db2
    ### END CODE HERE ###
    
    parameters = {"W1": W1,
                  "b1": b1,
                  "W2": W2,
                  "b2": b2}
    
    return parameters

4.4 整合前面4.1-4.3三个部分

# GRADED FUNCTION: nn_model

def nn_model(X, Y, n_h, num_iterations = 10000, print_cost=False):
    """
    Arguments:
    X -- dataset of shape (2, number of examples)
    Y -- labels of shape (1, number of examples)
    n_h -- size of the hidden layer
    num_iterations -- Number of iterations in gradient descent loop
    print_cost -- if True, print the cost every 1000 iterations
    
    Returns:
    parameters -- parameters learnt by the model. They can then be used to predict.
    """
    
    np.random.seed(3)
    n_x = layer_sizes(X, Y)[0]
    n_y = layer_sizes(X, Y)[2]
    
    # Initialize parameters
    ### START CODE HERE ### (≈ 1 line of code)
    parameters = initialize_parameters(n_x, n_h, n_y)
    ### END CODE HERE ###
    
    # Loop (gradient descent)

    for i in range(0, num_iterations):
         
        ### START CODE HERE ### (≈ 4 lines of code)
        # Forward propagation. Inputs: "X, parameters". Outputs: "A2, cache".
        A2, cache = forward_propagation(X, parameters)
        
        # Cost function. Inputs: "A2, Y, parameters". Outputs: "cost".
        cost = compute_cost(A2, Y, parameters)
 
        # Backpropagation. Inputs: "parameters, cache, X, Y". Outputs: "grads".
        grads = backward_propagation(parameters, cache, X, Y)
 
        # Gradient descent parameter update. Inputs: "parameters, grads". Outputs: "parameters".
        parameters = update_parameters(parameters, grads)
        
        ### END CODE HERE ###
        
        # Print the cost every 1000 iterations
        if print_cost and i % 1000 == 0:
            print ("Cost after iteration %i: %f" %(i, cost))

    return parameters

4.5 预测

# GRADED FUNCTION: predict

def predict(parameters, X):
    """
    Using the learned parameters, predicts a class for each example in X
    
    Arguments:
    parameters -- python dictionary containing your parameters 
    X -- input data of size (n_x, m)
    
    Returns
    predictions -- vector of predictions of our model (red: 0 / blue: 1)
    """
    
    # Computes probabilities using forward propagation, and classifies to 0/1 using 0.5 as the threshold.
    ### START CODE HERE ### (≈ 2 lines of code)
    A2, cache = forward_propagation(X, parameters)
    predictions = (A2>0.5)
    ### END CODE HERE ###
    
    return predictions

# Build a model with a n_h-dimensional hidden layer
parameters = nn_model(X, Y, n_h = 4, num_iterations = 10000, print_cost=True)

# Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), X, Y)
plt.title("Decision Boundary for hidden layer size " + str(4))

# Build a model with a n_h-dimensional hidden layer
parameters = nn_model(X, Y, n_h = 4, num_iterations = 10000, print_cost=True)

# Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), X, Y)
plt.title("Decision Boundary for hidden layer size " + str(4))

输出：

精度为：98%

三、《视觉slam十四讲》——第三讲 三维空间刚体运动

本节目标：理解三维空间的刚体运动描述方式：旋转矩阵、变换矩阵、四元数和欧拉角

3.1 旋转矩阵

3.1.1 点和向量，坐标系

可以用外积表示向量的旋转。考虑两个不平行的向量$\mathbf{a}$,$\mathbf{b}$,我们要描述$\mathbf{a}$到$\mathbf{b}$之间是如何旋转的。我们可以用一个向量来描述三维空间中两个向量的旋转关系。在右手法则下，我们用右手的四个指头从$\mathbf{a}$转向$\mathbf{b}$，其大拇指朝向就是旋转向量的方向，事实上也是$\mathbf{a}\times \mathbf{b}$的方向。它的大小由$\mathbf{a}$和$\mathbf{b}$的夹角决定。通过这种方式，我们构造了从$\mathbf{a}$到$\mathbf{b}$的一个旋转向量。这个向量同样位于三维空间中，在此坐标系下，可以用三个实数来描述它。

3.1.2 坐标系间的欧式变换

与向量间的旋转类似，我们同样可以描述两个坐标系之间的旋转关系，再加上平移，统称为坐标系之间的变换关系。在机器人运动过程中，常见的做法是设定一个惯性坐标系（或者叫世界坐标系），可以认为它是固定不动的。同时，相机或机器人则是一个移动坐标系。我们会问：相机视野中某个向量$\mathbf{p}$，它的坐标为$p_c$，而从世界坐标系下看，它的坐标是$p_w$。这两个坐标之间是如何转换的呢？这时，就需要先得到该点针对机器人坐标系坐标值，再根据机器人位姿转换到世界坐标系中，这个转换关系由一个矩阵$T$来描述，如图3-2所示。
.

相机运动是一个刚体运动，它保证了同一个向量在各个坐标系下的长度和夹角都不会发生变化。这种变换称为欧式变换。一个欧式变换一般由一个旋转和一个平移两部分组成。首先来考虑旋转。我们设某个单位正交基$\left(\begin{array}{c}{e}_1,e_2,e_3\end{array}\right)$经过一次旋转，变成了$\left(\begin{array}{c}{e}_1^{\prime },e_2^{\prime },e_3^{\prime }\end{array}\right)$。那么，对于同一个向量$\mathbf{a}$（注意该向量并没有随着坐标系的旋转而发生改动），它在两个坐标系下的坐标为${\left[\begin{array}{c}{a}_1,a_2,a_3\end{array}\right]}^T$和${\left[\begin{array}{c}{a}_1^{\prime },a_2^{\prime },a_3^{\prime }\end{array}\right]}^T$。根据坐标的定义，有：
$$\left[\begin{array}{c}{e}_1,e_2,e_3\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]=\left[\begin{array}{c}{e}_1^{\prime },e_2^{\prime },e_3^{\prime }\end{array}\right]\left[\begin{array}{c}{a}_1^{\prime }\\ a_2^{\prime }\\ a_3^{\prime }\end{array}\right](3.4)$$
为了描述两个坐标之间的关系，我们对上面等式左右同时左乘$\left[\begin{array}{c}{e}_1^T\\ e_2^T\\ e_3^T\end{array}\right]$，那么左边的系数变成了单位矩阵，所以：

我们把中间的矩阵拿出来，定义成一个矩阵$R$。这个矩阵由两组基之间的内积组成，刻画了旋转前后同一个向量的坐标变换关系。只要旋转是一样的，那么这个矩阵也是一样的。可以说，矩阵$R$描述了旋转本身。因此它又称为旋转矩阵。

旋转矩阵有一些特别的性质。事实上，它是一个行列式为1的正交矩阵。反之，行列式为1的正交矩阵也是一个旋转矩阵。所以，我们可以把旋转矩阵的集合定义如下：
S O ( n ) = { R ∈ R n × n ∣ R R T = I , d e t ( R ) = 1 } (3.6) SO(n)=\{R\in\mathbb{R}^{n\times n}|RR^T=I,det(R)=1\} \quad \text{(3.6)} SO(n)={R∈Rn×n∣RRT=I,det(R)=1}(3.6)
$SO(n)$是特殊正交群（Special Orthogonal Group）的意思。这个集合由n维空间的旋转矩阵组成，特别的，SO(3)就是三维空间的旋转矩阵。通过旋转矩阵，我们可以直接谈论两个坐标系之间的旋转变换，而不用再从基开始谈起了。换句话说，旋转矩阵可以描述相机的旋转

由于旋转矩阵为正交阵，它的逆（即转置）描述了一个相反的旋转。按照上面的定义方式，有：
$${\mathbf{a}}^{\prime }=R^{-1}\mathbf{a}=R^T\mathbf{a}.\text{(3.7)}$$
显然$R^T$刻画了一个相反的旋转

在欧式变换中，除了旋转之外还有一个平移。考虑世界坐标系中的向量$\mathbf{a}$，经过一次旋转（用$R$描述）和一次平移$\mathbf{t}$后，得到了${\mathbf{a}}^{\prime }$，那么把旋转和平移合到一起，有：
$${\mathbf{a}}^{\prime }=R\mathbf{a}+\mathbf{t}.\text{(3.8)}$$
其中，$\mathbf{t}$称为平移向量。通过上式，我们用一个旋转矩阵$R$和一个平移向量$\mathbf{t}$完整地描述了一个欧式空间的坐标变换关系。

3.1.3 变换矩阵与齐次坐标

式（3.8）完整地表达了欧式空间的旋转与平移，不过还存在一个小问题：这里的变换关系不是一个线性关系。这样的形式在变换多次后会过于复杂。因此，我们要引入齐次坐标和变换矩阵重写式（3.8）：

这是一个数学技巧：我们把一个三维向量的末尾添加$1$，变成了四维向量，称为齐次坐标。对于这个四维向量，我们可以把旋转和平移写在一个矩阵里面，使得整个关系变成了线性关系。该式中，矩阵$T$称为变换矩阵（Transform Matrix）。

关于变换矩阵$T$，它具有比较特别的结构：左上角为旋转矩阵，右侧为平移向量，左下角为$0$向量，右下角为$1$。这种矩阵又称特殊欧式群（Special Euclidean Group）：
S E ( 3 ) = { T = [ R t 0 T 1 ] ∈ R 4 × 4 ∣ R ∈ S O ( 3 ) , t ∈ R 3 } . (3.12) SE(3) = \{T=\begin{bmatrix}R&t\\ \mathbf{0}^T & 1 \end{bmatrix} \in \mathbb{R}^{4\times 4} | R \in SO(3), \mathbf{t} \in \mathbb{R}^3 \}. \quad \text{(3.12)} SE(3)={T=[R0T​t1​]∈R4×4∣R∈SO(3),t∈R3}.(3.12)
与SO(3)一样，求解该矩阵的逆表示一个反向的变换：
$$T^{-1}=\left[\begin{array}{cc}{R}^T& R^T\mathbf{t}\\ 0^T& 1\end{array}\right].\text{(3.13)}$$

回顾一下这部分的内容：首先，我们说了向量和它的坐标表示，并介绍了向量间的运算；然后，坐标系之间的运动由欧式变换描述，它由平移和旋转组成。旋转可以由旋转矩阵$SO(3)$描述，而平移直接由一个${\mathbb{R}}^3$向量描述。最后，如果将平移和旋转放在一个矩阵中，就形成了变换矩阵$SE(3)$。

3.2 实践：Eigen

Eigen是一个C++开源线性代数库。它提供了快速的有关矩阵的线性代数运算，还包括解方程等功能。本节学习了如何使用Eigen来表示矩阵、向量，随后引申至旋转矩阵和变换矩阵的计算。

3.3 旋转向量和欧拉角

3.4 四元数

通过对上面小节的学习，基本了解了刚体运动在空间中的位姿变换。3.3和3.4节有点深奥~，暂时先不读了，等下周或者将来用到时再回过头来看吧。

四、读MICCAI的那篇手术机器人论文——机器人辅助微创手术中的实时立体重建

这是一篇比较短的论文，主要内容总结如下：
这篇论文主要提出了一种新颖的立体重建算法，该算法在一组候选特征匹配周围传播视差信息。这具有避免镜面高光，仪器遮挡以及与视图有关的照明偏差的问题。此外，该算法可以与允许深度分布在非常分散的视图中传播的任何特征匹配策略一起使用。为具有已知几何体的幻影模型提供了验证，该数据可在线获得，以便在现场建立结构化验证方案。通过在机器人辅助程序的各种体内图像上进行重建，进一步证明了所提出方法的实用价值，这些方法也可供社区使用。

篇末总结了该方法的优劣。为了在体内图像上定性评估所提出方法的性能，本文使用了使用daVinci®手术系统从不同手术程序中获取的几个数据集。 视差图的结果和相应的重建3D渲染结果如下图所示。很明显，所提出的技术有效地捕获了手术部位的3D几何形状。 它很好地解决了由于手术器械而导致的较大视差不连续以及镜面反射较大的问题。 但是，存在一些错误，尤其是在遮挡边界处，需要进一步解决。 重要的是要注意，本文没有明确地迎合镜面反射或模型遮挡，例如通过使用检测，并且没有采用任何手术器械跟踪。 通过将这种策略纳入本文的方法，应该可以显着改善结果。