CS131 HW3

最新推荐文章于 2024-05-20 20:26:39 发布

ZWHY.

最新推荐文章于 2024-05-20 20:26:39 发布

阅读量1.4k

点赞数 2

分类专栏： CS131 Homework 文章标签： python 线性代数

本文链接：https://blog.csdn.net/weixin_43843322/article/details/119984188

版权

CS131 Homework 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

Homework 3

该作业涵盖了用于全景拼接的 Harris 角点检测器、RANSAC 和 HOG 描述符。

在本次作业中，我们将从多张图像中检测和匹配关键点以构建单个全景图像。这将涉及几个任务：

1.使用 Harris 角点检测器查找关键点。
2.构建一个描述符来描述图像中的每个点。
比较来自两个不同图像的两组描述符并找到匹配的关键点。
3.给定匹配的关键点列表，使用最小二乘法找到将一个图像中的点映射到另一个图像的仿射变换矩阵。
4.使用 RANSAC 对仿射变换矩阵进行更稳健的估计。
给定变换矩阵，用它来变换第二幅图像并将其叠加在第一幅图像上，形成全景图。
5.实现不同的描述符（HOG 描述符）并获得另一个拼接结果。

Part 1 Harris 角检测器

在本节中，您将实现 Harris 角点检测器以进行关键点定位。查看有关 Harris 角检测器的讲座幻灯片以了解其工作原理。Harris检测算法可以分为以下几个步骤：
1.计算一张图片x和y方向的梯度
2.计算每个像素的导数乘积(I2x,I2y,Ixy)
3.在每个像素处计算矩阵M，其中在这里插入图片描述
4.计算每个像素点的角响应
5.输出角响应映射R(x,y）

步骤1已经在panorama.py中的harris_corners函数中完成了。我们使用了Sobel算子，它计算x和y方向上每个像素的平滑梯度。有关sobel内核和操作符的更多信息，请参阅skimage文档中的sobel_v和sobel_h。

对于步骤3，我们已经在开始代码中为您创建了一个统一的窗口函数w。你可以假设窗口的大小是奇数。

完成harris_corners的函数实现并运行下面的代码。

提示:有两种方法可以解决这个问题
**向量化:**如果你想真正高效，你可以使用scipy. nimage .filters.convolve函数。它已经导入panorama.py，并一次计算每个像素的响应映射R。您可以在在这里插入图片描述

中计算窗口梯度，然后计算响应映射，无需任何for循环!

**迭代:**更直观的解决方案是迭代图像的每个像素，根据Ix, Iy, Ixy中像素梯度的周围邻域计算M，然后计算响应映射像素R(x,y)。你可能会发现hw1中conv_nested和conv_fast的实现是有用的参考!

注意，您需要显式地指定零填充来匹配Harris响应映射定义，但是我们也接受scipy. nimage .filters.convolve的默认行为。如果您使用零填充，那么向量化和for循环实现都将得到相同的答案!

“Alternate Accepted Harris Corner Solution”图像显示了scipy. nimage .filters的结果。 convolve的默认反射填充，而“Harris Corner Solution”的图像呈现的是零填充的解决方案。类似地，“Alternate Accepted Detected Corners Solution”图像将显示scipy. nimage .filters的结果。 convolve的默认反射填充，而“检测到的角落解决方案”图像呈现的是零填充解决方案。两者都是可接受的解决方案!

harris_corners函数（角点检测）：

`def harris_corners(img, window_size=3, k=0.04):
	    H, W = img.shape
	    window = np.ones((window_size, window_size))
	
	    response = np.zeros((H, W))
	
	    # 1. Compute x and y derivatives (I_x, I_y) of an image
	    dx = filters.sobel_v(img)
	    dy = filters.sobel_h(img)
	
	    ### YOUR CODE HERE
	    #计算Ix2、Iy2、IxIy
	    Ix2 = np.multiply(dx,dx)
	    Iy2 = np.multiply(dy,dy)
	    IxIy = np.multiply(dx,dy)
	    #每个像素点的权重均为一
	    A = convolve(Ix2,window)
	    B = convolve(Iy2,window)
	    C = convolve(IxIy,window)
	    #计算行列式和迹
	    detM = np.multiply(A, B) - np.multiply(C, C)
	    traceM = A + B
	
	    response = detM - k * np.square(traceM)
	    ### END YOUR CODE
	
	    return response`

Part 2 重点描述与匹配

Part 2.1 创建描述符
实现simple_descriptor函数，其中每个关键点由其周围小块的归一化强度描述。

def simple_descriptor(patch):
    """
    Describe the patch by normalizing the image values into a standard
    normal distribution (having mean of 0 and standard deviation of 1)
    and then flattening into a 1D array.

    The normalization will make the descriptor more robust to change
    in lighting condition.

    Hint:
        In this case of normalization, if a denominator is zero, divide by 1 instead.

    Args:
        patch: grayscale image patch of shape (H, W)

    Returns:
        feature: 1D array of shape (H * W)
    """
    ##归一化，在HW1中用过，进行归一化处理后将使得描述符在光照变化的情况下更加鲁棒
    feature = []
    ### YOUR CODE HERE
    mean = np.mean(patch)
    #np.std:返回给定数组的标准差
    sigma = np.std(patch)
    if sigma == 0.0:
        sigma = 1
    #normalized.flatten()：将矩阵降维成一维行向量
    normalized = (patch - mean) / sigma
    feature = normalized.flatten()
    ### END YOUR CODE
    return feature

Part 2.2 匹配描述符
a.使用标准化的密度来作为描述子
b.使用欧几里得距离来对描述子进行匹配,当最短距离与第二短距离的比值小于阈值，则判定为匹配

def simple_descriptor(patch):
    """
    Describe the patch by normalizing the image values into a standard
    normal distribution (having mean of 0 and standard deviation of 1)
    and then flattening into a 1D array.

    The normalization will make the descriptor more robust to change
    in lighting condition.

    Hint:
        In this case of normalization, if a denominator is zero, divide by 1 instead.

    Args:
        patch: grayscale image patch of shape (H, W)

    Returns:
        feature: 1D array of shape (H * W)
    """
    #归一化，在HW1中用过，进行归一化处理后将使得描述符在光照变化的情况下更加鲁棒
    feature = []
    ### YOUR CODE HERE
    mean = np.mean(patch)
    #np.std:返回给定数组的标准差
    sigma = np.std(patch)
    if sigma == 0.0:
        sigma = 1
    
    normalized = (patch - mean) / sigma
    #normalized.flatten()：将矩阵降维成一维行向量
    feature = normalized.flatten()
    ### END YOUR CODE
    return feature
  
```python
def describe_keypoints(image, keypoints, desc_func, patch_size=16):
    """
    Args:
        image: grayscale image of shape (H, W)
        keypoints: 2D array containing a keypoint (y, x) in each row
        desc_func: function that takes in an image patch and outputs
            a 1D feature vector describing the patch
        #接收并输出图像补丁的函数，描述该补丁的一维特征向量 
        patch_size: size of a square patch at each keypoint

    Returns:
        desc: array of features describing the keypoints
    """
    #将图像转换成数组
    image.astype(np.float32)
    desc = []
    #枚举：从索引0开始枚举列表中所有元素
    for i, kp in enumerate(keypoints):
        y, x = kp
        patch = image[y-(patch_size//2):y+((patch_size+1)//2),
                      x-(patch_size//2):x+((patch_size+1)//2)]
        #append：往列表的最后位置添加括号里的内容
        desc.append(desc_func(patch))
    return np.array(desc)

def match_descriptors(desc1, desc2, threshold=0.5):
    """
    Match the feature descriptors by finding distances between them. A match is formed
    when the distance to the closest vector is much smaller than the distance to the
    second-closest, that is, the ratio of the distances should be strictly smaller
    than the threshold (not equal to). Return the matches as pairs of vector indices.

    Hint:
        The Numpy functions np.sort, np.argmin, np.asarray might be useful

        The Scipy function cdist calculates Euclidean distance between all
        pairs of inputs
    Args:
        desc1: an array of shape (M, P) holding descriptors of size P about M keypoints
        desc2: an array of shape (N, P) holding descriptors of size P about N keypoints

    Returns:
        matches: an array of shape (Q, 2) where each row holds the indices of one pair
        of matching descriptors
        通过寻找特征描述符之间的距离来匹配特征描述符。 匹配就形成了  
        当到最近向量的距离远小于到第二近距离的，也就是说，距离的比例应该严格地更小  
        于阈值(不等于)。 以向量索引对的形式返回匹配项。  

        提示:  
        Numpy函数np.sort, np.argmin, np.asarray可能有用 
        
        np.argmin():给出水平方向的最小值下标。
        np.asarray：转化为数组
        
        Scipy函数cdist计算所有物体之间的欧氏距离对输入
        
        cdist返回值：一个二维数组，第一行为desc1的第一个值和desc2的所有值的欧氏距离，第二行为
        desc1的第二个值和desc2的所有值的欧氏距离。。。
        参数:  
        desc1:一个形状(M, P)数组，包含M个关键点大小为P的描述符  
        desc2:一个形状数组(N, P)，包含大小为P的描述符，包含N个关键点  
        
        返回:  
        match:一个形状数组(Q, 2)，其中每一行保存着一对的索引  
        匹配的描述符  
    """
    matches = []
    M = desc1.shape[0]
    #cdist函数就是计算欧几里得距离的函数
    dists = cdist(desc1, desc2)
    ### YOUR CODE HERE
    for i in range(M):
        sort=np.sort(dists[i,])
        if sort[0]/sort[1]<threshold:
#matches里面存的不是数值，如果计算得到某个点匹配，就把该点的坐标记录下来。i横坐标。 np.argmin(dists[i,:])竖坐标
            matches.append([i, np.argmin(dists[i,:])])
    #将list类型转换为数组类型，因为后面对matches的使用要求其是数组类型，否则会报错
    matches=np.asarray(matches)
    ### END YOUR CODE

    return matches

Part 3 变换估计

计算图2到图1的仿射变换矩阵。由于匹配点对众多，这里采用最小二乘法估计变换矩阵。

def fit_affine_matrix(p1, p2):
    """ 
    Fit affine matrix such that p2 * H = p1. First, pad the descriptor vectors
    with a 1 using pad() to convert to homogeneous coordinates, then return
    the least squares fit affine matrix in homogeneous coordinates.

    Hint:
        You can use np.linalg.lstsq function to solve the problem. 

        Explicitly specify np.linalg.lstsq's new default parameter rcond=None 
        to suppress deprecation warnings, and match the autograder.

    Args:
        p1: an array of shape (M, P) holding descriptors of size P about M keypoints
        p2: an array of shape (M, P) holding descriptors of size P about M keypoints

    Return:
        H: a matrix of shape (P+1, P+1) that transforms p2 to p1 in homogeneous
        coordinates
        
       拟合仿射矩阵，使p2 * H = p1。 首先，使用pad()将描述符向量填充1，转换为齐次坐标，然后返        回齐次坐标下的最小二乘拟合仿射矩阵。
       pad():
        将[[ 0.3 -0.2]
         [-0.4 -0.9]
         [ 0.1  0.1]]
         填充为：
        [[ 0.3 -0.2  1. ]
         [-0.4 -0.9  1. ]
         [ 0.1  0.1  1. ]]
 
       提示:  
       您可以使用np. linalgl .lstsq函数来解决这个问题。  

       显式地指定np.linalg。 lstsq的新默认参数rcond=None  
       禁止弃用警告，并匹配自动分级器。  
        
       参数:  
       p1:一个形状(M, P)数组，包含M个关键点大小为P的描述符  
       p2:一个形状(M, P)数组，包含M个关键点大小为P的描述符  

       返回:  
       H:形状为(P+1, P+1)的矩阵，它使p2在齐次中变换为p1坐标  
    """

    assert (p1.shape[0] == p2.shape[0]),\
        'Different number of points in p1 and p2'
    p1 = pad(p1)
    p2 = pad(p2)

    ### YOUR CODE HERE
    #np.linalg.lstsq的返回值有四个，定义四个参数来接收它们
    #linalg.lstsq(a, b, rcond='warn')计算的是a×H=b中的矩阵H，所以要将p1与 p2换个位置
    H, residuals, rank, s = np.linalg.lstsq(p2, p1, rcond=None)
    ### END YOUR CODE

    # Sometimes numerical issues cause least-squares to produce the last
    # column which is not exactly [0, 0, 1]
    H[:,2] = np.array([0, 0, 1])
    return H

Part 4 RANSAC

def ransac(keypoints1, keypoints2, matches, n_iters=200, threshold=20):
    """
    Use RANSAC to find a robust affine transformation:

        1. Select random set of matches
        2. Compute affine transformation matrix
        3. Compute inliers via Euclidean distance
        4. Keep the largest set of inliers
        5. Re-compute least-squares estimate on all of the inliers

    Update max_inliers as a boolean array where True represents the keypoint
    at this index is an inlier, while False represents that it is not an inlier.

    Hint:
        You can use np.linalg.lstsq function to solve the problem. 

        Explicitly specify np.linalg.lstsq's new default parameter rcond=None 
        to suppress deprecation warnings, and match the autograder.

        You can compute elementwise boolean operations between two numpy arrays,
        and use boolean arrays to select array elements by index:
        https://numpy.org/doc/stable/reference/arrays.indexing.html#boolean-array-indexing 

    Args:
        keypoints1: M1 x 2 matrix, each row is a point
        keypoints2: M2 x 2 matrix, each row is a point
        matches: N x 2 matrix, each row represents a match
            [index of keypoint1, index of keypoint 2]
        n_iters: the number of iterations RANSAC will run
        threshold: the number of threshold to find inliers

    Returns:
        H: a robust estimation of affine transformation from keypoints2 to
        keypoints 1
        使用RANSAC找到一个健壮的仿射变换:  

        1.  随机选择一组匹配  
        2.  计算仿射变换矩阵  
        3.  通过欧几里得距离计算嵌线  
        4.  保留最大的嵌套  
        5.  重新计算所有嵌线的最小二乘估计  

        将max_inliers更新为一个布尔数组，其中True表示关键点  
        在这个索引处是一个嵌套器，而False表示它不是嵌套器。  

        提示:  
        您可以使用np. linalgl .lstsq函数来解决这个问题。  

        显式地指定np.linalg。 lstsq的新默认参数rcond=None  
        禁止弃用警告，并匹配自动分级器。  

        你可以计算两个numpy数组之间的elementwise布尔运算，  
        使用布尔数组通过索引选择数组元素:  
        https://numpy.org/doc/stable/reference/arrays.indexing.html#boolean-array-indexing  

        参数:  
        关键点1:M1 × 2矩阵，每一行是一个点  
        关键点2:M2 x2矩阵，每一行是一个点  
        匹配:N x2矩阵，每一行代表一个匹配  
        [index of keypoint t1, index of keypoint 2]  
        n_iter: RANSAC将运行的迭代次数  
        阈值:查找嵌套的阈值数目  

        返回:  
        H:从关键点2到的仿射变换的稳健估计  
        要点1  
    """
    # Copy matches array, to avoid overwriting it
    orig_matches = matches.copy()
    matches = matches.copy()

    N = matches.shape[0]
    n_samples = int(N * 0.2)

    matched1 = pad(keypoints1[matches[:,0]])
    matched2 = pad(keypoints2[matches[:,1]])
    #以matches数组的列数定义一个一维向量，默认为false（即不是关键点）
    max_inliers = np.zeros(N, dtype=bool)
    n_inliers = 0

    # RANSAC iteration start
    
    # Note: while there're many ways to do random sampling, please use
    # `np.random.shuffle()` followed by slicing out the first `n_samples`
    # matches here in order to align with the auto-grader.
    # Sample with this code:
    #shuffle直接在原来的数组上进行操作，改变原来数组的顺序，无返回值(是对列表x中的所有元素随机打乱顺序，若x不是列表，则报错)。
    '''
        np.random.shuffle(matches)
        samples = matches[:n_samples]
        sample1 = pad(keypoints1[samples[:,0]])
        sample2 = pad(keypoints2[samples[:,1]])
    '''
    ### YOUR CODE HERE
    for i in range(n_iters):
        inliersArr = np.zeros(N)
        #打乱matches顺序
        np.random.shuffle(matches)
        #拿出前五分之一作为一组
        samples = matches[:n_samples]
        #取出两张图的关键点坐标
        p1 = pad(keypoints1[samples[:,0]])
        p2 = pad(keypoints2[samples[:,1]])
        #计算仿射变换矩阵
        H, residuals, rank, s = np.linalg.lstsq(p2, p1, rcond=None)
        #将图二的点与仿射变换矩阵相乘
        output = np.dot(matched2, H)
        #将匹配结果小于阈值的进行计数
        #np.linalg.norm：以行向量的形式计算范数
        inliersArr = np.linalg.norm(output-matched1, axis=1)**2 < threshold

        inliersCount = np.sum(inliersArr)
        #保持max_inliers记录最大数目的匹配点
        if inliersCount > n_inliers:
            max_inliers = inliersArr.copy() #还是注意深拷贝和浅拷贝啊
            n_inliers = inliersCount
    
    # 迭代完成，拿最大数目的匹配点对进行估计变换矩阵
    H, residuals, rank, s = np.linalg.lstsq(matched2[max_inliers], matched1[max_inliers], rcond=None)
  
    ### END YOUR CODE
    return H, orig_matches[max_inliers]

Part 5 Histogram of Oriented Gradients (HOG)

HOG特征提取步骤（简化版本）：

1.计算图像中x和y方向上的梯度
2.计算梯度直方图将图像分成多个单元，计算每个单元内的梯度分布，梯度分布n个方向中，统计单元内每个方向的梯度强度之和
3.归一化图像块的直方图
4.将块直方图转换成向量
待续

ZWHY.

关注

2
点赞
踩
12

收藏

觉得还不错? 一键收藏
2
评论
CS131 HW3

Homework 3该作业涵盖了用于全景拼接的 Harris 角点检测器、RANSAC 和 HOG 描述符。在本次作业中，我们将从多张图像中检测和匹配关键点以构建单个全景图像。这将涉及几个任务：1.使用 Harris 角点检测器查找关键点。2.构建一个描述符来描述图像中的每个点。比较来自两个不同图像的两组描述符并找到匹配的关键点。3.给定匹配的关键点列表，使用最小二乘法找到将一个图像中的点映射到另一个图像的仿射变换矩阵。4.使用 RANSAC 对仿射变换矩阵进行更稳健的估计。给定变换矩阵，用
复制链接

扫一扫