CS131 Homework

Filters-Instagram

This notebook includes both coding and written questions. Please hand in this notebook file with all the outputs and your answers to the written questions.

# Setup

Convolutions

@(Commutative Property)
Recall that the convolution of an image $f: ℝ^2→ℝ$ and a kernel $h: ℝ^2→ℝ$ is defined as follows:

$\begin{aligned}(f∗h)[m,n] &=∑_{i=−∞}^∞∑_{j=−∞}^∞f[i,j]⋅h[m−i,n−j] \tag{1.0.1}\\ &= (h∗f)[m,n] \\ &= ∑_{i=−∞}^∞∑_{j=−∞}^∞h[i,j]⋅f[m−i,n−j]\end{aligned}$

将式中积分变量 i 和 j 置换为 m−x 和 n−y，即可证明

@(Linear and Shift Invariance)
Let $f$ be a function $ℝ^2→ℝ$. Consider a system $f\xrightarrow sg$, where $g=(f∗h)$ with some kernel $h:ℝ^2→ℝ$. Show that $S$ defined by any kernel $h$ is a Linear Shift Invariant (LSI) system. In other words, for any $h$, show that $S$ satisfies both of the following:

$S[a⋅f1+b⋅f2]=a⋅S[f1]+b⋅S[f2]$

$\text{If } f[m,n]\xrightarrow sg[m,n] \text{ then } f[m−m_0,n−n_0]\xrightarrow sg[m−m_0,n−n_0]$

向量化

向量化可以提高计算效率，实现更快的卷积计算。

【算法+图像处理】2D卷积与快速卷积算法C语言实现

向量化的思路就是：

1. 将卷积核进行向量化
2. 将迭代过程中的计算块（计算中心及其领域）进行向量化，由于运算的对象都是核，此处可以将向量化的计算块拼成矩阵
3. 两者进行内积计算，并reshape成原有图像大小

Cross-correlation

Cross-correlation of two 2D signals $f$ and $g$ is defined as follows:

$(f⋆g)[m,n]=∑_{i=−∞}^∞∑_{j=−∞}^∞f[i,j]⋅g[i−m,j−n]$

Template Matching with Cross-correlation

使用互相关性实现模板匹配
Suppose that you are a clerk at a grocery store. One of your responsibilites is to check the shelves periodically and stock them up whenever there are sold-out items. You got tired of this laborious task and decided to build a computer vision system that keeps track of the items on the shelf.
Luckily, you have learned in CS131 that cross-correlation can be used for template matching: a template $g$ is multiplied with regions of a larger image $f$ to measure how similar each region is to the template.

The template of a product (template.jpg) and the image of shelf (shelf.jpg) is provided. We will use cross-correlation to find the product in the shelf.
Implement cross_correlation function in filters.py and run the code below.

Hint: you may use the conv_fast function you implemented in the previous question.

from filters import cross_correlation

Interpretation

How does the output of cross-correlation filter look like? Was it able to detect the product correctly? Explain what might be the problem with using raw template as a filter.

Zero-mean cross-correlation

A solution to this problem is to subtract off the mean value of the template so that it has zero mean.
Implement zero_mean_cross_correlation function in filters.py and run the code below.

Normalized Cross-correlation

One day the light near the shelf goes out and the product tracker starts to malfunction. The zero_mean_cross_correlation is not robust to change in lighting condition. The code below demonstrates this.
A solution is to normalize the pixels of the image and template at every step before comparing them. This is called normalized cross-correlation.

The mathematical definition for normalized cross-correlation of $f$ and template $g$ is:

$(f⋆g)[m,n]=∑_{i,j} \frac{f[i,j] - \overline {f_{m,n}}}{\sigma_{f_{m,n}}} \cdot \frac{g[i-m, j-n] - \bar g}{\sigma_g}$

where

$f_{m,n}$is the patch image at position (m,n)
$\overline{f_{m,n}}$ is the mean of the patch image $f_{m,n}$
$\sigma_{f_{m,n}}$ is the standard deviation of the patch image $f_{m,n}$
$\bar g$ is the mean of the template $g$
$σ_g$ is the standard deviation of the template $g$

不可以直接对整幅待匹配的图做归一化操作。

正确的图像区域，如果是做整幅图的归一化，明显受到了其他区域诸如亮度颜色等的影响，干扰了匹配。
应该对每个待匹配块做归一化操作，然后比较匹配程度。

Separable Filters

Consider a $M_1×N_1$ image I and a $M_2×N_2$ filter F. A filter F is separable if it can be written as a product of two 1D filters: $F=F_1F_2$.

For example,

$F = \begin{bmatrix}[1 & -1\\ 1& - 1 \end{bmatrix}$

can be written as a matrix product of

$F_1 = \begin{bmatrix}1\\1\end{bmatrix}, F_2 = [1,-1]$

Therefore $F$ is a separable filter.

For any separable filter $F=F_1F_2 \to I∗F=(I∗F_1)∗F_2$

Complexity comparison

1. $I∗F$乘法操作的次数为 $M_1⋅N_1$
2. $(I∗F_1)∗F_2$乘法操作的次数为 $M_1+N_1$

Edges-Smart Car Lane Detection

Canny Edge Detector

In this part, you are going to implment Canny edge detector. The Canny edge detection algorithm can be broken down in to five steps:

1. Smoothing：平滑处理（去躁）
2. Finding gradients：寻找梯度
3. Non-maximum suppression：非极大值抑制
4. Double thresholding：双阈值算法
5. Edge tracking by hysteresis：滞后边缘追踪

Smoothing

We first smooth the input image by convolving it with a Gaussian kernel. The equation for a Gaussian kernel of size (2k+1)×(2k+1) is given by:

$h_{ij}=\frac{1}{2πσ^2}exp(−\frac{(i−k)^2+(j−k)^2}{2σ^2}),0≤i,j<2k+1$

Implement gaussian_kernel in edge.py and run the code below.

from edge import conv, gaussian_kernel

需要先实现edge.py中的conv卷积函数，可以使用作业hw1中的conv_nested，conv_fast,conv_faster等函数

# Test with different kernel_size and sigma

What is the effect of the kernel_size and sigma?

高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点权值是随该点与中心点的距离单调增减的，因为边缘是一种图像局部特征，如果平滑运算对离算子中心很远的像素点仍然有很大作用，则平滑运算会使图像失真。

高斯滤波器宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的．σ越大，高斯滤波器的频带就越宽，平滑程度就越好．通过调节平滑程度参数σ，可在图像特征过分模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷。

Finding gradients

The gradient of a 2D scalar function $I:ℝ^2→ℝ$ in Cartesian coordinate is defined by:

$\begin{aligned}∇I(x,y) &= [\frac{∂I}{∂x},\frac{∂I}{∂y}]\\ &= [\lim_{\Delta x \to 0} \frac{I(x+\Delta x, y) - I(x, y)}{\Delta x}, \lim_{\Delta y \to 0} \frac{I(x, y + +\Delta y) - I(x, y)}{\Delta y}] \\ &\approx [\frac{I(x + 1, y) - I(x-1, y)}{2}, \frac{I(x, y+1) - I(x,y-1)}{2}] \end{aligned}$

Note that the partial derivatives can be computed by convolving the image $I$ with some appropriate kernels $D_x$ and $D_y$:

$\begin{aligned}\frac{\partial I}{\partial x} &\approx I * D_x = G_x \\\frac{\partial I}{\partial y} &\approx I * D_y = G_y \end{aligned}$

Find the kernels $D_x$ and $D_y$ and implement partial_x and partial_y using conv defined in edge.py.

-Hint: Remeber that convolution flips the kernel.

from edge import partial_x, partial_y

# Compute partial derivatives of smoothed image

What is the reason for performing smoothing prior to computing the gradients?

平滑是为了去噪，噪声会使得计算梯度时得到幅度较大的异常值。

Now, we can compute the magnitude and direction of gradient with the two partial derivatives:

$G = \sqrt{G_x^2 + G_y^2}\\\theta = arctan(\frac{G_y}{G_x})$

Implement gradient in edge.py which takes in an image and outputs $G$ and $Θ$.

from edge import gradient

Non-maximum suppression

You should be able to note that the edges extracted from the gradient of the smoothed image is quite thick and blurry. The purpose of this step is to convert the “blurred” edges into “sharp” edges. Basically, this is done by preserving all local maxima in the gradient image and discarding everything else. The algorithm is for each pixel $(x,y)$ in the gradient image:

1. Round the gradient direction $Θ[y,x]$ to the nearest 45 degrees, corresponding to the use of an 8-connected neighbourhood.

当然这种做法其实是简化版本的，虽然Canny的论文上是这么写的，但是后来也有提出，梯度肯定不止这四个，针对别的梯度方向的比较，使用插值等方法来进行非极大值抑制. Canny算子中的非极大值抑制（Non-Maximum Suppression）分析

2. Compare the edge strength of the current pixel with the edge strength of the pixel in the positive and negative gradient direction. For example, if the gradient direction is south (theta=90), compare with the pixels to the north and south.
3. If the edge strength of the current pixel is the largest; preserve the value of the edge strength. If not, suppress (i.e. remove) the value.

 

非极大值抑制属于一种边缘细化的方法，梯度大的位置有可能为边缘，在这些位置沿着梯度方向，找到像素点的局部最大值，并将其非最大值抑制。

形象的说就是一条粗粗的区域都根据梯度强度判断成了边缘，但目标只需要一条细细的边，这个时候就找区域内的最大的那条作为边（注意是沿着梯度方向），其他的或抛弃或作为备选区域

from edge import non_maximum_suppression

Double Thresholding

The edge-pixels remaining after the non-maximum suppression step are (still) marked with their strength pixel-by-pixel. Many of these will probably be true edges in the image, but some may be caused by noise or color variations, for instance, due to rough surfaces. The simplest way to discern between these would be to use a threshold, so that only edges stronger that a certain value would be preserved. The Canny edge detection algorithm uses double thresholding. Edge pixels stronger than the high threshold are marked as strong; edge pixels weaker than the low threshold are suppressed and edge pixels between the two thresholds are marked as weak.

Implement double_thresholding in edge.py

双阀值方法，设置一个maxval，以及minval，梯度大于maxval则为强边缘，梯度值介于maxval与minval则为弱边缘点，小于minval为抑制点。

双阈值的结果在第五步可以用来作为边缘追踪的依据。

from edge import double_thresholding

Edge tracking

Strong edges are interpreted as “certain edges”, and can immediately be included in the final edge image. Weak edges are included if and only if they are connected to strong edges. The logic is of course that noise and other small variations are unlikely to result in a strong edge (with proper adjustment of the threshold levels). Thus strong edges will (almost) only be due to true edges in the original image. The weak edges can either be due to true edges or noise/color variations. The latter type will probably be distributed in dependently of edges on the entire image, and thus only a small amount will be located adjacent to strong edges. Weak edges due to true edges are much more likely to be connected directly to strong edges.

Implement link_edges in edge.py

由于边缘是连续的，因此可以认为弱边缘如果为真实边缘则和强边缘是联通的，可由此判断其是否为真实边缘。

Canny edge detector

Implement canny in edge.py using the functions you have implemented so far. Test edge detector with different parameters.

from edge import canny

Suppose that the Canny edge detector successfully detects an edge in an image. The edge is then rotated by $θ$, where the relationship between a point on the original edge $(x,y)$ and a point on the rotated edge $(x′,y′)$ is defined as

$x′=xcosθ \\ y′=xsinθ$

Will the rotated edge be detected using the same Canny edge detector? Provide either a mathematical proof or a counter example.

-Hint: The detection of an edge by the Canny edge detector depends only on the magnitude of its derivative. The derivative at point $(x, y)$ is determined by its components along the $x$ and $y$ directions. Think about how these magnitudes have changed because of the rotation.

$Mag_{origin}=\sqrt{dx^2+dy^2}=|dx|$

$Mag_{rotate} = \sqrt{dx'^2+dy'^2} = \sqrt{(dxcos\theta)^2 + (dysin\theta)^2} = |dx|$

幅度没变，所以边缘检测不变。

 

After running the Canny edge detector on an image, you notice that long edges are broken into short segments separated by gaps. In addition, some spurious edges appear. For each of the two thresholds (low and high) used in hysteresis thresholding, explain how you would adjust the threshold (up or down) to address both problems. Assume that a setting exists for the two thresholds that produces the desired result. Briefly explain your answer.

调参方法：

如果长直线被断成小段，说明weak_edges的阈值太大，weak_edges的数量较少，此时应当调低weak_edges的阈值。

如果spurious (伪直线)出现，说明strong_edges的阈值太小，strong_edges的数量较多，此时应当调高strong_edges的阈值。

Extra Credit: Optimizing Edge Detector

One way of evaluating an edge detector is to compare detected edges with manually specified ground truth edges. Here, we use precision, recall and F1 score as evaluation metrics. We provide you 40 images of objects with ground truth edge annotations. Run the code below to compute precision, recall and F1 score over the entire set of images. Then, tweak the parameters of the Canny edge detector to get as high F1 score as possible. You should be able to achieve F1 score higher than 0.31 by carefully setting the parameters.

from os import listdir

Lane Detection

In this section we will implement a simple lane detection application using Canny edge detector and Hough transform.
Here are some example images of how your final lane detector will look like.

The algorithm can broken down into the following steps:

1. Detect edges using the Canny edge detector.
2. Extract the edges in the region of interest (a triangle covering the bottom corners and the center of the image).
3. Run Hough transform to detect lanes.

Edge detection

Lanes on the roads are usually thin and long lines with bright colors. Our edge detection algorithm by itself should be able to find the lanes pretty well. Run the code cell below to load the example image and detect edges from the image.

from edge import canny

Extracting region of interest (ROI)

We can see that the Canny edge detector could find the edges of the lanes. However, we can also see that there are edges of other objects that we are not interested in. Given the position and orientation of the camera, we know that the lanes will be located in the lower half of the image. The code below defines a binary mask for the ROI and extract the edges within the region.

H, W = img.shape

Fitting lines using Hough transform

The output from the edge detector is still a collection of connected points. However, it would be more natural to represent a lane as a line parameterized as $y=ax+b$, with a slope $a$ and $y$-intercept $b$. We will use Hough transform to find parameterized lines that represent the detected edges.

In general, a straight line $y=ax+b$ can be represented as a point $(a,b)$ in the parameter space. However, this cannot represent vertical lines as the slope parameter will be unbounded. Alternatively, we parameterize a line using $θ∈[−π,π]$ and $ρ∈ℝ$ as follows:

$ρ=x⋅cosθ+y⋅sinθ$

Using this parameterization, we can map everypoint in $xy$-space to a sine-like line in $θρ$-space (or Hough space). We then accumulate the parameterized points in the Hough space and choose points (in Hough space) with highest accumulated values. A point in Hough space then can be transformed back into a line in $xy$-space.

See lectures on Hough transform.

Implement hough_transform in edge.py.

from edge import hough_transform

之所以只检测出一条直线，是因为上面的边缘检测部分不是很完美，尤其是左边车道线的边缘检测出现了问题，才会导致左边车道线做霍夫变换也得不到理想的结果

Panorama-ImageStiching

This assignment covers Harris corner detector(角点检测器), RANSAC(随机抽样一致) and HOG descriptor(HOG描述子).

Panorama stitching is an early success of computer vision. Matthew Brown and David G. Lowe published a famous panoramic image stitching paper in 2007. Since then, automatic panorama stitching technology has been widely adopted in many applications such as Google Street View, panorama photos on smartphones,and stitching software such as Photosynth and AutoStitch.

In this assignment, we will detect and match keypoints from multiple images to build a single panoramic image. This will involve several tasks:

1. Use Harris corner detector to find keypoints.
2. Build a descriptor to describe each point in an image.
   Compare two sets of descriptors coming from two different images and find matching keypoints.
3. Given a list of matching keypoints, use least-squares method to find the affine transformation matrix that maps points in one image to another.
4. Use RANSAC to give a more robust estimate of affine transformation matrix.
   Given the transformation matrix, use it to transform the second image and overlay it on the first image, forming a panorama.
5. Implement a different descriptor (HOG descriptor) and get another stitching result.

Harris Corner Detector

In this section, you are going to implement Harris corner detector for keypoint localization. Review the lecture slides on Harris corner detector to understand how it works. The Harris detection algorithm can be divide into the following steps:

1. Compute $x$ and $y$ derivatives $(I_x,I_y)$ of an image
2. Compute products of derivatives ($I_x^2, I_y^2, I_{xy}$) at each pixel
3. Compute matrix M at each pixel, where $M = \sum_{x,y}w(x,y)\begin{bmatrix}I_x^2 & I_xI_y \\ I_xI_y & I_y^2\end{bmatrix}$
4. Compute corner response $R=Det(M)−k(Trace(M)^2)$ at each pixel
5. Output corner response map $R(x,y)$

Step 1 is already done for you in the function harris_corners in panorama.py. Complete the function implementation and run the code below.

-Hint: You may use the function scipy.ndimage.filters.convolve, which is already imported in panoramy.py
关于Harris角点检测，这里有篇博客总结的很好。Harris角点

from panorama import harris_corners

Once you implement the Harris detector correctly, you will be able to see small bright blobs around the corners of the sudoku grids and letters in the output corner response image. The function corner_peaks from skimage.feature performs non-maximum suppression to take local maxima of the response map and localize keypoints.

对响应输出进行非极大抑制。

# Perform non-maximum suppression in response map

Describing and Matching Keypoints

We are now able to localize keypoints in two images by running the Harris corner detector independently on them. Next question is, how do we determine which pair of keypoints come from corresponding locations in those two images? In order to match the detected keypoints, we must come up with a way to describe the keypoints based on their local appearance. Generally, each region around detected keypoint locations is converted into a fixed-size vectors called descriptors.

Creating Descriptors

In this section, you are going to implement a simple_descriptor; each keypoint is described by normalized intensity in a small patch around it.

为了进行两幅图像中的关键点的点与点之间的匹配，需要对关键点进行描述。
此时，提出描述符对关键点的局部特征进行描述。
一般来说，描述符就是将检测到的关键点的附近区域转换成固定大小的向量。

from panorama import harris_corners

Matching Descriptors

Then, implement match_descriptors function to find good matches in two sets of descriptors. First, calculate Euclidean distance between all pairs of descriptors from image 1 and image 2. Then use this to determine if there is a good match: if the distance to the closest vector is significantly (by a factor which is given) smaller than the distance to the second-closest, we call it a match. The output of the function is an array where each row holds the indices of one pair of matching descriptors.

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

from panorama import simple_descriptor, match_descriptors, describe_keypoints

Transformation Estimation

We now have a list of matched keypoints across the two images. We will use this to find a transformation matrix that maps points in the second image to the corresponding coordinates in the first image. In other words, if the point $p_1=[y_1,x_1]$ in image 1 matches with $p_2=[y_2,x_2]$ in image 2, we need to find an affine transformation matrix H such that

\tilde {p_2} H = \tilde{p_1}$$ where $\tilde {p_1}, \tilde{p_2}$ are homogenous coordinates of $p_1$ and $p_2$.

Note that it may be impossible to find the transformation $H$ that maps every point in image 2 exactly to the corresponding point in image 1. However, we can estimate the transformation matrix with least squares. Given N matched keypoint pairs, let $X_1$ and $X_2$ be $N×3$ matrices whose rows are homogenous coordinates of corresponding keypoints in image 1 and image 2 respectively. Then, we can estimate H by solving the least squares problem

$X_2H = X_1$

Implement fit_affine_matrix in panorama.py

-Hint: read the documentation about np.linalg.lstsq

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

from panorama import fit_affine_matrix

After checking that your fit_affine_matrix function is running correctly, run the following code to apply it to images.
Images will be warped and image 2 will be mapped to image 1. Then, the two images are merged to get a panorama. Your panorama may not look good at this point, but we will later use other techniques to get a better result.

from utils import get_output_space, warp_image

merged = img1_warped + img2_warped

RANSAC

Rather than directly feeding all our keypoint matches into fit_affine_matrix function, we can instead use RANSAC (“RANdom SAmple Consensus”) to select only “inliers” to use to compute the transformation matrix.
The steps of RANSAC are:

1. Select random set of matches
2. Compute affine transformation matrix
3. Find inliers using the given threshold
4. Repeat and keep the largest set of inliers
5. Re-compute least-squares estimate on all of the inliers

Implement ransac in panorama.py, run through the following code to get a panorama. You can see the difference from the result we get without RANSAC.

from panorama import ransac

output_shape, offset = get_output_space(img1, [img2], [H])

merged = img1_warped + img2_warped

Histogram of Oriented Gradients (HOG)

In the above code, you are using the simple_descriptor, and in this section, you are going to implement a simplified version of HOG descriptor.

HOG stands for Histogram of Oriented Gradients. In HOG descriptor, the distribution ( histograms ) of directions of gradients ( oriented gradients ) are used as features. Gradients ( x and y derivatives ) of an image are useful because the magnitude of gradients is large around edges and corners ( regions of abrupt intensity changes ) and we know that edges and corners pack in a lot more information about object shape than flat regions.

The steps of HOG are:

1. compute the gradient image in $x$ and $y$
   Use the sobel filter provided by skimage.filters
2. compute gradient histograms
   Divide image into cells, and calculate histogram of gradient in each cell.
3. normalize across block
   Normalize the histogram so that they
4. flattening block into a feature vector

Implement hog_descriptor in panorama.py, and run through the following code to get a panorama image.

关于HOG特征的算法流程，可以阅读此文

【特征检测】HOG特征算法

from panorama import hog_descriptor

# Extract features from the corners

from panorama import ransac

output_shape, offset = get_output_space(img1, [img2], [H])

merged = img1_warped + img2_warped

实现了HOG特征描述符之后，可以替代掉之前的归一化密度作为特征的匹配。效果不错。

Extra Credit: Better Image Merging

You will notice the blurry region and unpleasant lines in the middle of the final panoramic image. In the cell below, come up with a better merging scheme to make the panorama look more natural. Be creative!

怎么消除拼接处的两条拼接痕迹呢？

1. OpenCV探索之路（二十四）图像拼接和图像融合技术
2. 《Automatic Panoramic Image Stitching using Invariant Features》

加权融合

博客上的处理思路是加权融合，在重叠部分由前一幅图像慢慢过渡到第二幅图像，即将图像的重叠区域的像素值按一定的权值相加合成新的图像。

我的理解就是两块重叠的区域，从左向右每个像素由图1和图2融合，最开始图1的权重最大，图2的权重最小，向右融合过程中图1权重逐渐减小，图2权重逐渐增大。

融合的公式可以简单理解成：

$dst[i,j]=Img_1[i,j]∗α_{[i,j]}+Img_2[i,j]∗(1−α_{[i,j]})$

而权重$α$的值则与像素点在重叠区域的位置相关，可以定义成：

$\alpha_{[row, col]} = (col - Left_{region})/(Width_{region})$

论文《Automatic Panoramic Image Stitching using Invariant Features》论文针对拼接步骤，提出了三点：

1. 自动拼接校直 Automatic Panorama Straightening

校直是用BA做的相机参数的估计，用来解决拼接过程中的图像的旋转

2. 增益补偿 Gain Compensation

增益补偿是用过使用最小二乘法计算每个拼接块的增益，用来均衡整幅图像中各个拼接块的亮度不均问题

3. 多频带混合 Multi-Band Blending

# Modify the code below

OpenPano

Github上宇哥开源的拼接库OpenPano

Extra Credit: Stitching Multiple Images

Work in the cell below to complete the code to stitch an ordered chain of images.

Given a sequence of m images $(I_1,I_2,...,I_m)$, take every neighboring pair of images and compute the transformation matrix which converts points from the coordinate frame of $I_{i+1}$ to the frame of $I_i$. Then, select a reference image $I_{ref}$, which is in the middle of the chain. We want our final panorama image to be in the coordinate frame of $I_{ref}$. So, for each $I_i$ that is not the reference image, we need a transformation matrix that will convert points in frame $i$ to frame $ref$.

-Hint:

• If you are confused, you may want to review the Linear Algebra slides on how to combine the effects of multiple transformation matrices.
• The inverse of transformation matrix has the reverse effect. Please use numpy.linalg.inv function whenever you want to compute matrix inverse.

img1 = imread('yosemite1.jpg', as_grey=True)

#  进行拼接，拼接方式一

多幅图的拼接就可以看出来简单方式的问题了：

1. 拼接痕迹的存在，可以使用加权平均，
2. 各个拼接块亮度不均匀，需要做增益补偿，
3. 存在未知的图像的旋转，可能是由于拍摄时相机发生了位姿的变动，需要做优化处理（校直），
4. 重叠部分的拼接并不能直接实现所有像素的对应叠加，可以看出来存在模糊的情况，可以使用多频带混合的方式去做处理。

Resizing-SeamCarving

复现Seam Carving（接缝拼接）算法，实现基于内容感知的图像智能缩放，并且实现图像中目标的移除，了解图像能量图、动态编程。

@Reference
CS131: Computer Vision Foundations and Applications

The material presented here is inspired from:

1. paper on seam carving: http://graphics.cs.cmu.edu/courses/15-463/2007_fall/hw/proj2/imret.pdf
2. tutorial: http://cs.brown.edu/courses/cs129/results/proj3/taox
3. tutorial: http://www.cs.cmu.edu/afs/andrew/scs/cs/15-463/f07/proj2/www/wwedler/

The whole seam carving process was covered in lecture 8, please refer to the slides for more details to the different concepts introduced here.

这份作业的主题是实现SeamCarving(接缝裁剪)的算法，即基于图像内容感知的图像缩放算法。

1. 实现图像裁剪
2. 实现图像扩大
3. 其他实验

也可以参考这篇专栏讲了如何用Python实现Seam carving算法，原理讲解的挺清楚的。
如何用Python实现神奇切图算法Seam Carving ？

Sean Carving 的基本思路

1. 计算图像能量图
2. 寻找最小能量线
3. 移除得到的最小能量线，让图片的宽度缩小一个像素

Seam carving is an algorithm for content-aware image resizing.
To understand all the concepts in this homework, make sure to read again the slides from lecture 8
http://vision.stanford.edu/teaching/cs131_fall1718/files/08_seam_carving.pdf

from skimage import io, util

Energy function

We will now implemented the energy_function to compute the energy of the image.
The energy at each pixel is the sum of:

1. absolute value of the gradient in the $x$ direction
2. absolute value of the gradient in the $y$ direction

The function should take around 0.01 to 0.1 seconds to compute.

图像能量图此处定义为图像每个像素的两方向梯度幅值之和，当然能量也可以是别的定义方式：梯度幅度、熵图、显著图等等。计算梯度要注意转换为灰度图再进行。

from seam_carving import energy_function

# Compute energy function

Compute cost

Now implement the function compute_cost.

Starting from the energy map, we’ll go from the first row of the image to the bottom and compute the minimal cost at each pixel.
We’ll use dynamic programming to compute the cost line by line starting from the first row.
The function should take around 0.05 seconds to complete.

此处要求使用动态规划计算整幅图的cost图，用于第三步寻找最小能量线。
此处动态规划的意思是从图像的最上面一行开始，计算能量的最小累加路线，cost图中每一行的每个像素从它八相邻的上一行的三个像素中取最小cost值的元素进行累加。
直到累加到最后一行，同时在计算过程中保留每个像素的累加路径，即像素的cost值是从上一行哪一个元素相加得到的（用-1,0,1表示上左，上中，上右）。
题目中要求只用一个循环实现，则需要用向量化的形式进行计算。

from seam_carving import compute_cost

# Vertical Cost Map

# Horizontal Cost Map

Finding optimal seams

Using the cost maps we found above, we can determine the seam with the lowest energy in the image.
We can then remove this optimal seam, and repeat the process until we obtain a desired width.

由于上面计算能量cost图时，已经保存了每个像素点的叠加路径，因此可以在最后一行找到最小的cost，并从该点开始，向上溯源，找到一条最佳的接缝，即最小的能量线。

Backtrack seam

Implement function backtrack_seam

from seam_carving import backtrack_seam

vcost, vpaths = compute_cost(img, energy)

到此，就实现了一条接缝的查找，接下来移除它就可以实现一个像素宽度的缩小。

In the image above, the optimal vertical seam (minimal cost) goes through the portion of sky without any cloud, which yields the lowest energy.

Reduce

We can now use the function backtrack and remove_seam iteratively to reduce the size of the image through seam carving.

Each reduce can take around 10 seconds to compute, depending on your implementation.
If it’s too long, try to vectorize your code in compute_cost to only use one loop.

from seam_carving import reduce

# Reduce image width

裁剪部分我是按照每一行删除一个元素的方式去实现的。
能够达到题目要求的10s左右的时间要求。

We observe that resizing from width 640 to width 400 conserves almost all the important part of the image (the person and the castle), where a standard resizing would have compressed everything.

All the vertical seams removed avoid the person and the castle.

程序同样可以缩减高度，程序中对于另一个维度的裁剪做了判断，
只需要seam carving流程走之前旋转图像，走完流程之后再将裁剪宽度的图像旋转回来就成了裁剪高度了。

# Reduce image height

For reducing the height, we observe that the result does not look as nice.

The issue here is that the castle is on all the height of the image, so most horizontal seams will go through it.
Interestingly, we observe that most of the grass is not removed. This is because the grass has small variation between neighboring pixels (in a kind of noisy pattern) that make it high energy.
The seams removed go through the sky on the left, go under the castle to remove some grass and then back up in the low energy blue sky.

虽然实现了高度的裁剪，但是图像看起来怪怪的，这是可以从图像上做分析的。
城堡明显具有很强的能量，详见第一步能量图计算，因此裁剪高度时，每一条最小能量线都会越过它。
草坪看起来其实跟噪声一样，每个像素点与相邻像素都有或大或小的差异，因此梯度也是较大的，因此能量也较高，不会首先裁剪草坪而是更平滑的蓝天，
而城堡下面的草坪凹陷是因为城堡上方的蓝天已被裁剪光了，最小能量线无法从上方穿过时才从城堡下方的草坪开始裁剪的。

Image Enlarging

Enlarge naive

We now want to tackle the reverse problem of enlarging an image.
One naive way to approach the problem would be to duplicate the optimal seam iteratively until we reach the desired size.

图像扩大的一种方式是直接在最小能量线的地方复制它。
在enlarge_naive函数中实现它。

from seam_carving import enlarge_naive

W_new = 800

The issue with enlarge_naive is that the same seam will be selected again and again, so this low energy seam will be the only to be duplicated.

Another way to get k different seams is to apply the process we used in function reduce, and keeping track of the seams we delete progressively.
The function find_seams(image, k) will find the top k seams for removal iteratively.

The inner workings of the function are a bit tricky so we’ve implemented it for you, but you should go into the code and understand how it works.
This should also help you for the implementation of enlarge.

作业已经帮忙实现了：查找k条不同的接缝，这样就可以根据k条接缝来轮流进行删除了。

from seam_carving import find_seams

不同的颜色是因为按照排序，给每条路线的标记的值不同。

Enlarge

We can see that all the seams found are different, and they avoid the castle and the person.

One issue we can mention is that we cannot enlarge more than we can reduce. Because of our process, the maximum enlargement is the width of the image W because we first need to find W different seams in the image.

One effect we can see on this image is that the blue sky at the right of the castle can only be enlarged x2. The concentration of seams in this area is very strong.
We can also note that the seams at the right of the castle have a blue color, which means they have low value and were removed in priority in the seam selection process.

值得注意的是，对于图像能放大的像素宽度肯定不会超过图像能缩减的像素宽度，毕竟最小能量线就只有那么多条。

实现enlarge函数。

from seam_carving import enlarge

W_new = 800

# Map of the seams for horizontal seams.

H_new = 600

As you can see in the example above, the sky above the castle has doubled in size, the grass below has doubled in size but we still can’t reach a height of 600.
The algorithm then needs to enlarge the castle itself, while trying to avoid enlarging the windows for instance.

Other experiments on the image

Feel free to experiment more on this image, try different sizes to enlarge or reduce, or check what seams are chosen…

Reducing by a 2x factor often leads to weird patterns.
Enlarging by more than 2x is impossible since we only duplicate seams. One solution is to enlarge in mutliple steps (enlarge x1.4, enlarge again x1.4…)

放大两倍或者缩小二分之一图像会变得很奇怪。
这里提出一个方式就是进行多次Resize步骤，每次x1.4倍，x1.4倍这样的方式来放大，缩小同理。

# Reduce image width

Faster reduce (extra credit 0.5%)

Implement a faster version of reduce called reduce_fast in the file seam_carving.py.

We will have a leaderboard on gradescope with the performance of students.

The autograder tests will check that the outputs match, and run the reduce_fast function on a set of images with varying shapes (say between 200 and 800).

怎么能实现更快地裁剪呢？

题目的提示是：do we really need to compute the whole cost map again at each iteration?
分析了一下reduce的过程，的确有一块是可以改进的。

while out.shape[1] > size:

就是计算energy能量图这一步，因为能量图的计算其实就是梯度的计算。

如果一个区域内的图像像素不变动的话，那这个区域内的梯度应该也不会放生改变，基于这一点，可以对计算能量图的步骤进行改进。

只需要对上一次optical seam最小能量线覆盖宽度内的区域进行梯度的重新更新，两侧区域的能量值不需要变动。
比如接缝最左边为left，则left-1处像素不变，梯度改变，left-2处及之前像素不变梯度不变。
比如接缝最右边为right，则right+1处像素不变，梯度改变，left+2处及之后像素不变梯度不变。

然后考虑特殊情况。

from seam_carving import reduce_fast

Reducing and enlarging on another image

Also check these outputs with another image.

# Load image

out = reduce(img2, 300)

out = enlarge(img2, 800)

Forward Energy

Forward energy is a solution to some artifacts that appear when images have curves for instance.
Implement the function compute_forward_cost. This function will replace the compute_cost we have been using until now.

基本的seam carving 会影响物体的原有结构，比如曲线啊。
因此提出了Forward Energy这个方法。

详见 CS131 Lecture8 的PPT P76-P86。

介绍上说的是 Instead of removing the seam of least energy, remove the seam that inserts the least energy to the image !
之前的做法是去除掉具有最小能量的接缝，而Forward Energy的做法是去除掉给整张图像带来最小能量的接缝。
具体的可视化解释在PPT中，这里就放公式了。
之前的能量图定义为

$M(i,j)=E(i,j)+min\begin{cases}M(i-1, j-1) \\ M(i-1, j) \\ M(i-1, j+1)\end{cases}$

Forward能量图的定义为

$M(i,j)=E(i,j)+min\begin{cases}M(i-1, j-1) + C_L(i,j) & C_L(i,j) = |I(i,j+1) - I(i,j-1)| + |I(i-1,j) - I(i,j-1)|\\ M(i-1, j) + C_V(i,j) & C_V(i,j) = |I(i,j + 1) - I(i, j - 1)| \\ M(i-1, j+1) + C_R(i,j) & C_R(i,j) = |I(i,j+1) - I(i,j-1)| + |I(i-1,j) - I(i,j+1)|\end{cases}$

# Load image

from seam_carving import compute_forward_cost

from seam_carving import compute_forward_cost

We observe that the forward energy insists more on the curved edges of the image.
可以看到，Forward Energy方法在保留边缘曲线的结构上比原方法好。

from seam_carving import reduce

The issue with our standard reduce function is that it removes vertical seams without any concern for the energy introduced in the image.

In the case of the dinosaur above, the continuity of the shape is broken. The head is totally wrong for instance, and the back of the dinosaure lacks continuity.

Forward energy will solve this issue by explicitely putting high energy on a seam that breaks this continuity and introduces energy.

# This step can take a very long time depending on your implementation.

Object removal (extra credit 0.5%)

Object removal uses a binary mask of the object to be removed.

Using the reduce and enlarge functions you wrote before, complete the function remove_object to output an image of the same shape but without the object to remove.

移除目标和保留目标

移除目标的思路应该很清晰，就是将图像中目标处的像素视为能量最小或最大，即可将目标去除了。

两个问题：

1. 如何进行标记，从而每次删除能量线的时候，都能删到并且尽可能多的删到带有标记的能量线？同时坐标下标不会混乱。
2. 如何在删除目标之后保持和原图同样的尺寸

解决方式：

1. 判断最小能量线的时候，将mask区域乘上负数，自然找到的接缝能量会最小，并且穿过mask区域的像素点会更多。
2. 做裁剪的时候，将mask根据接缝和图像一起进行裁剪，自然不会引起下标混乱。
3. 删除目标之后再调用之前的enlarge函数，进行图像的放大，放大到原图尺寸。

# Load image

# Use your function to remove the object

Segmentation-Clustering

实现聚类算法：K-Means和HAC(凝聚聚类算法)；
构建图像的特征序列，进行聚类分割；
实现分割结果与GroundTruth的量化评估，根据指标评估聚类算法。

This notebook includes both coding and written questions. Please hand in this notebook file with all the outputs and your answers to the written questions.

This assignment covers K-Means and HAC methods for clustering.

这份作业包含的主题：用于聚类的K-means 和 HAC（层次凝聚聚类）方法，以及聚类算法的量化评估。

Introduction

In this assignment, you will use clustering algorithms to segment images. You will then use these segmentations to identify foreground and background objects.

Your assignment will involve the following subtasks:

1. Clustering algorithms: Implement K-Means clustering and Hierarchical Agglomerative Clustering.
2. Pixel-level features: Implement a feature vector that combines color and position information and implement feature normalization.
3. Quantitative Evaluation: Evaluate segmentation algorithms with a variety of parameter settings by comparing your computed segmentations against a dataset of ground-truth segmentations.

Clustering Algorithms

# Generate random data points for clustering

K-Means Clustering

As discussed in class, K-Means is one of the most popular clustering algorithms. We have provided skeleton code for K-Means clustering in the file segmentation.py. Your first task is to finish implementing kmeans in segmentation.py. This version uses nested for loops to assign points to the closest centroid and compute new mean of each cluster.

目前题图允许使用循环实现点到最近中心的分配即可。
K-means的步骤：

1. 随机初始化聚类中心
2. 将点分配给最近的中心
3. 计算每个聚类的新的中心
4. 如果聚类的分配不再变动则停止
5. 从第二步继续循环

from segmentation import kmeans

简单的实现了K-Means算法，第三步判断终止条件很重要，添加终止条件后只迭代了几次就结束了，大大降低了运行时间。

We can use numpy functions and broadcasting to make kmeans faster. Implement kmeans_fast in segmentation.py. This should run at least 10 times faster than the previous implementation.

思路：

题目的提示是使用numpy函数和broadcasting机制
np.repeat and np.argmin 有帮助
还有一种思路就是向量化了，不做迭代而是放在矩阵里去计算。

from segmentation import kmeans_fast

快了十六倍 ，666，向量化是个好东西！！！

Hierarchical Agglomerative Clustering

Another simple clustering algorithm is Hieararchical Agglomerative Clustering, which is somtimes abbreviated as HAC. In this algorithm, each point is initially assigned to its own cluster. Then cluster pairs are merged until we are left with the desired number of predetermined clusters (see Algorithm 1).

Implement hiererachical_clustering in segmentation.py.

按照上图的算法定义是：

1. 初始化每个点为一个簇
2. 迭代直到只剩下K个簇：
   1. 计算所有簇与簇之间（簇对）的距离
   2. 合并最近的两个簇

这个算法的关键，主要是如何定义簇与簇之间的距离。定义有下面三种：

1. 单链（Single-link）：不同两个聚类中离得最近的两个点之间的距离，即MIN；
2. 全链（Complete-link）：不同两个聚类中离得最远的两个点之间的距离，即MAX；
3. 平均链（Average-link）：不同两个聚类中所有点对距离的平均值，即AVERAGE；

题目代码中提示的做法是，利用簇中心的欧几里得距离来定义簇与簇之间的距离。
hierarchical_clustering函数将按照这个思路实现。

Numpy中的距离度量函数: 距离度量以及python实现(一)

from segmentation import hierarchical_clustering

运行时间比预期的块，总结一下就是这类算法的提速方式：

1. 第一个是善用向量化，用矩阵来计算问题
2. 第二个就是善用numpy提供的函数，解决一些需要搜索遍历的问题

Pixel-Level Features

Before we can use a clustering algorithm to segment an image, we must compute some feature vectore for each pixel. The feature vector for each pixel should encode the qualities that we care about in a good segmentation. More concretely, for a pair of pixels $p_i$ and $p_j$ with corresponding feature vectors $f_i$ and $f_j$, the distance between $f_i$ and $f_j$ should be small if we believe that $p_i$ and $p_j$ should be placed in the same segment and large otherwise.

上一步的聚类算法实现其实是流程的第二步，第一步其实是特征的提取，有了特征才能进行像素（或者说样本）的聚类。只不过第一步中，拿了点的位置作为特征来进行聚类了。

# Load and display image

Color Features

One of the simplest possible feature vectors for a pixel is simply the vector of colors for that pixel. Implement color_features in segmentation.py. Output should look like the following:

在color_features实现简单的颜色特征的聚类。

from segmentation import color_features

In the cell below, we visualize each segment as the mean color of pixels in the segment.

将聚类结果用颜色进行显示（颜色是指源图各个分割区域内的颜色均值）

from utils import visualize_mean_color_image

Color and Position Features

Another simple feature vector for a pixel is to concatenate its color and position within the image. In other words, for a pixel of color $(r,g,b)$ located at position $(x,y)$ in the image, its feature vector would be $(r,g,b,x,y)$. However, the color and position features may have drastically different ranges; for example each color channel of an image may be in the range $[0,1)$, while the position of each pixel may have a much wider range. Uneven scaling between different features in the feature vector may cause clustering algorithms to behave poorly.

One way to correct for uneven scaling between different features is to apply some sort of normalization to the feature vector. One of the simplest types of normalization is to force each feature to have zero mean and unit variance.

Implement color_position_features in segmentation.py.

Output segmentation should look like the following:

结合颜色特征和位置特征的聚类

其实就是生成特征为 $(r,g,b,x,y)$ 的序列，但需要考虑的问题是尺度的不同，解决的方案也很简单，就是进行归一化，对每一个维度进行归一化。公式如下：

$feature_{new} = \frac{feature - Mean(feature)}{Std(feature)}$

from segmentation import color_position_features

visualize_mean_color_image(img, segments)

Extra Credit: Implement Your Own Feature

For this programming assignment we have asked you to implement a very simple feature transform for each pixel. While it is not required, you should feel free to experiment with other feature transforms. Could your final segmentations be improved by adding gradients, edges, SIFT descriptors, or other information to your feature vectors? Could a different type of normalization give better results?

Implement your feature extractor my_features in segmentation.py

Depending on the creativity of your approach and the quality of your writeup, implementing extra feature vectors can be worth extra credit (up to 5% of total points).

Describe your approach: (YOUR APPROACH)

实现自己定义的特征

可以实现一些简单的特征，比如添加梯度、边缘、SIFT描述子或者其他信息。
我添加了一个梯度强度信息，结合颜色特征，位置特征。

from segmentation import my_features

Quantitative Evaluation

Looking at images is a good way to get an idea for how well an algorithm is working, but the best way to evaluate an algorithm is to have some quantitative measure of its performance.

For this project we have supplied a small dataset of cat images and ground truth segmentations of these images into foreground (cats) and background (everything else). We will quantitatively evaluate different segmentation methods (features and clustering methods) on this dataset.

We can cast the segmentation task into a binary classification problem, where we need to classify each pixel in an image into either foreground (positive) or background (negative). Given the ground-truth labels, the accuracy of a segmentation is (TP+TN)/(P+N).

Implement compute_accuracy in segmentation.py.
实现准确率计算就是正确个数与总个数之比。

from segmentation import compute_accuracy

You can use the script below to evaluate a segmentation method’s ability to separate foreground from background on the entire provided dataset. Use this script as a starting point to evaluate a variety of segmentation parameters.

下面这段程序是已经完成的，主要逻辑就是

1. 加载数据集，包括原图和分割结果真值
2. 调用之前实现的聚类分割算法进行分割
3. 将真值和分割结果进行比对，计算准确率

可以在下面进行调参，更改参数，查看效果。

from utils import load_dataset, compute_segmentation

# Visualize segmentation results

Include a detailed evaluation of the effect of varying segmentation parameters (feature transform, clustering method, number of clusters, resize) on the mean accuracy of foreground-background segmentations on the provided dataset. You should test a minimum of 10 combinations of parameters. To present your results, add rows to the table below (you may delete the first row).

Observe your results carefully and try to answer the following question:

1. Based on your quantitative experiments, how do each of the segmentation parameters affect the quality of the final foreground-background segmentation?
2. Are some images simply more difficult to segment correctly than others? If so, what are the qualities of these images that cause the segmentation algorithms to perform poorly?
3. Also feel free to point out or discuss any other interesting observations that you made.

基于对结果的观察：

1. 基于颜色-位置特征的聚类效果比基于颜色特征的聚类效果差
2. 聚类数目3个就是最好的了，增大减小都影响效果
3. 尺度是越接近原图，也就是越接近1，准确度越高

的确存在一些图更难识别。这些图的一些特性包括

1. 背景较为复杂
2. 前景颜色、亮度都与背景较为接近

同时基于颜色-位置特征的聚类中，尺度并不是越大越好，可能是我做的实验次数不够多。

Recognition-Classification

实现分类算法。
使用SVD进行图像压缩。
使用KNN进行图像识别（分类）。
使用PCA（主成分分析）和LDA（线性判别分析）对数据进行降维。

This notebook includes both coding and written questions. Please hand in this notebook file with all the outputs as a pdf on gradescope for “HW6 pdf”. Upload the three files of code (compression.py, k_nearest_neighbor.py and features.py) on gradescope for “HW6 code”.

This assignment covers:

1. image compression using SVD
2. kNN methods for image recognition.
3. PCA and LDA to improve kNN

Image Compression

Image compression is used to reduce the cost of storage and transmission of images (or videos).
One lossy compression method is to apply Singular Value Decomposition (SVD) to an image, and only keep the top n singular values.

机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用

image = io.imread('pitbull.jpg', as_grey=True)

Let’s implement image compression using SVD.
We first compute the SVD of the image, and as seen in class we keep the $n$ largest singular values and singular vectors to reconstruct the image.

Implement function compress_image in compression.py.

在compress_image函数中实现奇异值分解，并保留前N个较大的奇异值和对应的特征向量。
题目中的压缩率在我理解就是压缩后有效数据与压缩之前图像数据的比值。

from compression import compress_image

Original image shape: (1704, 1280)

# Number of singular values to keep

Data size (original): 2181120

Data size (original): 2181120

Data size (original): 2181120

实现了用SVD实现了图像的压缩，同时可以发现：n=100时，压缩率0.13时，图像质量就已经很接近原图了。
这样看来，SVD的确用来做降维是很合适的。毕竟能够大大提高运行速度。

Face Dataset

We will use a dataset of faces of celebrities. Download the dataset using the following command:

sh get_dataset.sh

The face dataset for CS131 assignment.
The directory containing the dataset has the following structure:

faces/

Each class has 50 training images and 10 testing images.

命令行执行sh get_dataset.sh就可以下载了。

# 加载数据集

Class names: ['angelina jolie', 'anne hathaway', 'barack obama', 'brad pitt', 'cristiano ronaldo', 'emma watson', 'george clooney', 'hillary clinton', 'jennifer aniston', 'johnny depp', 'justin timberlake', 'leonardo dicaprio', 'natalie portman', 'nicole kidman', 'scarlett johansson', 'tom cruise']

# 可视化数据集中的一些样本

# 将图像数据展开为行，即每个样本为一行

Training data shape: (800, 4096)

k-Nearest Neighbor

We’re now going to try to classify the test images using the k-nearest neighbors algorithm on the raw features of the images (i.e. the pixel values themselves). We will see later how we can use kNN on better features.

Here are the steps that we will follow:

1. We compute the L2 distances between every element of X_test and every element of X_train in compute_distances
2. We split the dataset into 5 folds for cross-validation in split_folds.
3. For each fold, and for different values of k, we predict the labels and measure accuracy.
4. Using the best k found through cross-validation, we measure accuracy on the test set.

kNN（k最近邻算法）的步骤如下：

1. 计算测试数据X_test中每个样本与训练数据X_train中每个样本之间的距离，距离计算方式为L2距离。
2. 将数据集分割为5Folder（5折）
3. 对每个子集和每个不同的参数k，预测结果并计算准确率
4. 通过交叉验证找到最佳的参数k，并在测试集上评估准确率

from k_nearest_neighbor import compute_distances

dists shape: (160, 800)

预测结果分两步：

1. 找到预测点的最近的k个点
2. 判断k个点中频率最高的train_label

from k_nearest_neighbor import predict_labels

Got 38 / 160 correct => accuracy: 0.237500

Cross-Validation

We don’t know the best value for our parameter k.
There is no theory on how to choose an optimal k, and the way to choose it is through cross-validation.

We cannot compute any metric on the test set to choose the best k, because we want our final test accuracy to reflect a real use case. This real use case would be a setting where we have new examples come and we classify them on the go. There is no way to check the accuracy beforehand on that set of test examples to determine k.

Cross-validation will make use split the data into different fold (5 here).
For each fold, if we have a total of 5 folds we will have:

• 80% of the data as training data
• 20% of the data as validation data

We will compute the accuracy on the validation accuracy for each fold, and use the mean of these 5 accuracies to determine the best parameter k.

from k_nearest_neighbor import split_folds

# Step 3: Measure the mean accuracy for each value of `k`

# plot the raw observations

从这幅图看，前30之内有两个峰值就达到了最高的准确率和较好的标准差。
可以进一步细化step进行测试。

# Step 3: Measure the mean accuracy for each value of `k`

# plot the raw observations

# Based on the cross-validation results above, choose the best value for k,   

For k = 2, got 42 / 160 correct => accuracy: 0.262500

既然k=2和k=17准确率相同，应该是选择k较小的，这样可以提高运行速率吧。

PCA

Principal Component Analysis (PCA) is a simple yet popular and useful linear transformation technique that is used in numerous applications, such as stock market predictions, the analysis of gene expression data, and many more. In this tutorial, we will see that PCA is not just a “black box”, and we are going to unravel its internals in 3 basic steps.

Introduction
The sheer size of data in the modern age is not only a challenge for computer hardware but also a main bottleneck for the performance of many machine learning algorithms. The main goal of a PCA analysis is to identify patterns in data; PCA aims to detect the correlation between variables. If a strong correlation between variables exists, the attempt to reduce the dimensionality only makes sense. In a nutshell, this is what PCA is all about: Finding the directions of maximum variance in high-dimensional data and project it onto a smaller dimensional subspace while retaining most of the information.

A Summary of the PCA Approach

• Standardize the data.
• Obtain the Eigenvectors and Eigenvalues from the covariance matrix or correlation matrix, or perform Singular Vector Decomposition.
• Sort eigenvalues in descending order and choose the k eigenvectors that correspond to the k largest eigenvalues where k is the number of dimensions of the new feature subspace (k≤d).
• Construct the projection matrix $W$ from the selected k eigenvectors.
• Transform the original dataset $X$ via $W$ to obtain a $k$-dimensional feature subspace Y.

浅谈对主成分分析（PCA）算法的理解

PCA的主旨就是：在高维数据中找到最大方差的方向，并将高维数据投射到低维的子空间，这样能够保留尽可能多的信息。

1. 标准化数据
2. 从协方差矩阵或互相关矩阵或者用SVD奇异值分解，来获得特征向量和特征值
3. 将特征值降序排序，根据新特征空间的维度k，选择前k个最大的特征值和对应的特征向量（k <= d）
4. 根据k个特征向量构建投影矩阵W
5. 将数据集X通过W转换到k维的特征子空间

from features import PCA

Eigendecomposition

The eigenvectors and eigenvalues of a covariance (or correlation) matrix represent the “core” of a PCA: The eigenvectors (principal components) determine the directions of the new feature space, and the eigenvalues determine their magnitude. In other words, the eigenvalues explain the variance of the data along the new feature axes.

Implement _eigen_decomp in pca.py.

特征向量和特征值代表了PCA的核心：

1. 特征向量（主成分）决定了新的特征空间的方向
2. 特征值决定了特征向量的幅度，换言之，特征值表示了数据在新特征方向的方差。方差越大，信息量越大。

 

协方差矩阵计算方法

# Perform eigenvalue decomposition on the covariance matrix of training data.

(4096,)

直接解特征向量太慢了。

Singular Value Decomposition

Doing an eigendecomposition of the covariance matrix is very expensive, especially when the number of features (D = 4096 here) gets very high.

To obtain the same eigenvalues and eigenvectors in a more efficient way, we can use Singular Value Decomposition (SVD). If we perform SVD on matrix $X$, we obtain $U$, $S$ and $V$ such that:

$X = USV^T$

• the columns of U are the eigenvectors of $XX^T$
• the columns of VT are the eigenvectors of $X^TX$
• the values of S are the square roots of the eigenvalues of $X^TX (or XX^T)$

Therefore, we can find out the top k eigenvectors of the covariance matrix $X^TX$ using SVD.

Implement _svd in pca.py.

这里讲到了直接利用协方差矩阵计算特征分解是非常耗时的事情，特别是维度很高的情况下。
其实更高效的方式是使用SVD（奇异值分解）。
具体的左奇异向量和右奇异向量与特征向量的对应关系在上面描述了，所以利用SVD就能得到特征向量和特征值

# Perform SVD on directly on the training data.

(800,)

# Test whether the square of singular values and eigenvalues are the same.

Dimensionality Reduction

The top k principal components explain most of the variance of the underlying data.

By projecting our initial data (the images) onto the subspace spanned by the top k principal components,
we can reduce the dimension of our inputs while keeping most of the information.

In the example below, we can see that using the first two components in PCA is not enough to allow us to see pattern in the data. All the classes seem placed at random in the 2D plane.

维度缩减
通过将前k个主成分投影在子空间，就可以达到缩减维度的目的，同时能够保留最大的信息。
需要在feature.py中实现 transform 函数，同时在fit函数中实现X的均值计算

# Dimensionality reduction by projecting the data onto

一脸懵逼，这里我们可以看到，仅仅用两个最大的主成分是无法看出数据的特征来的。

Visualizing Eigenfaces

The columns of the PCA projection matrix pca.W_pca represent the eigenvectors of $X^TX$.

We can visualize the biggest singular values as well as the corresponding vectors to get a sense of what the PCA algorithm is extracting.

If we visualize the top 10 eigenfaces, we can see tht the algorithm focuses on the different shades of the faces. For instance, in face n°2 the light seems to come from the left.

for i in range(10):

# Reconstruct data with principal components

(800, 4096)

100个主成分能够分清男女了，200个主成分大概能够分清是哪个人了。

Written Question 1

Question: Consider a dataset of N face images, each with shape (h,w). Then, we need O(Nhw) of memory to store the data. Suppose we perform dimensionality reduction on the dataset with p principal components, and use the components as bases to represent images. Calculate how much memory we need to store the images and the matrix used to get back to the original space.

Said in another way, how much memory does storing the compressed images and the uncompresser cost.

1. 降维后的图像维度为 X_proj.shape = (N, n_components)，所以实际压缩的图像占用空间即为O(Np)
2. 为了实现压缩和复原的其他资源（姑且称之为解压器），需要保存的变量包括self.W_pca和self.mean，实际上内存大头就是样本的所有主成分表，所以解压器占用的内存空间为O(pD) = O(phw)

Reconstruction error and captured variance

We can plot the reconstruction error with respect to the dimension of the projected space.
The reconstruction gets better with more components.
We can see in the plot that the inflexion point is around dimension 200 or 300. This means that using this number of components is a good compromise between good reconstruction and low dimension.

# Plot reconstruction errors for different k

从上面这幅图可以获取到的有用信息就是：

1. 1000以内，误差随着主成分数目增加而减小
2. 大概在200-400间，出现了下降速率的拐角，此处可以来个权衡，因为主成分数目并不是越多越好，它影响了运行的效率。选择拐角处的值，既能保证误差率较低，还能兼顾计算时间。

We can do the same process to see how much variance is captured by the projection.
Again, we see that the inflexion point is around 200 or 300 dimensions.

同样，方差的拐点也在200-400之间

# Plot captured variance

kNN with PCA

Performing kNN on raw features (the pixels of the image) does not yield very good results.
Computing the distance between images in the image space is not a very good metric for actual proximity of images.
For instance, an image of person A with a dark background will be close to an image of B with a dark background, although these people are not the same.

Using a technique like PCA can help discover the real interesting features and perform kNN on them could give better accuracy.

However, we observe here that PCA doesn’t really help to disentangle the features and obtain useful distance metrics between the different classes. We basically obtain the same performance as with raw features.

上面这段话讲的挺好，直接使用原始的特征（即图像的像素）并不能表现的很好，因为在计算图像之间的距离时，两个不同物体具有相同的黑色背景时，背景明显在计算上减小了两者之间的距离，虽然我们要判断的是前景的类别。
使用PCA能够得到有趣的特征，能够使kNN能够达到更高的准确率。
但是，在这里，PCA并不能在不同类别里区分开特征，获取到有用的距离。这里看起来和使用像素是一样的效果。

num_test = X_test.shape[0]

Got 45 / 160 correct => accuracy: 0.281250

PCA有用，但就涨了这么点准确率。

因为，在这里，PCA并不能在不同类别里区分开特征，获取到有用的距离。这里看起来和使用像素是一样的效果。

Fisherface: Linear Discriminant Analysis

LDA is a linear transformation method like PCA, but with a different goal.
The main difference is that LDA takes information from the labels of the examples to maximize the separation of the different classes in the transformed space.

Therefore, LDA is not totally unsupervised since it requires labels. PCA is fully unsupervised.

In summary:

• PCA perserves maximum variance in the projected space.
• LDA preserves discrimination between classes in the project space. We want to maximize scatter between classes and minimize scatter intra class.

1. 机器学习中的数学(4)-线性判别分析（LDA）, 主成分分析(PCA)
2. 如何理解LDA算法？能够简洁明了地说明一下LDA算法的中心思想吗？

LDA和PCA的区别

LDA的主旨就是从标签中获取信息，在转换的目标空间最大化不同种类之间的隔离。

from features import LDA

Dimensionality Reduction via PCA

To apply LDA, we need D<N. Since in our dataset, N=800 and D=4096, we first need to reduce the number of dimensions of the images using PCA.
More information at: http://www.scholarpedia.org/article/Fisherfaces

通过PCA实现维度的缩减
这里提到数据特征维度要比样本数少，所以需要进行一步特征的维度缩减，就直接拿PCA做降维好了。

N = X_train.shape[0]

Scatter matrices

We first need to compute the within-class scatter matrix:

$S_W = \sum_{i=1}^cS_i$

where $S_i = \sum{x_k \in Y_i} (x_k - \mu_i)(x_k - \mu_i)^T$ is the scatter of class $i$.

We then need to compute the between-class scatter matrix:

$S_B=∑_{i=1}^cN_i(μ_i−μ)(μ_i−μ)^T$

where $N_i$ is the number of examples in class $i$.

散点图矩阵

上面的公式偏向结论的，具体的推导在Lecture13 的PPT 第42页开始。

# Compute within-class scatter matrix

(784, 784)

# Compute between-class scatter matrix

(784, 784)

Solving generalized Eigenvalue problem

Implement methods fit and transform of the LDA class.

根据Lecture 第55页，最终将优化问题泛化为了一个求特征值的问题。编程时直接使用线性代数库求解特征值即可。

lda.fit(X_train_pca, y_train)

# Dimensionality reduction by projecting the data onto

第一幅图的效果好是因为已经fit过了，即进行过了训练。
第二幅图的效果不好，泛化能力明显不行。

实现维度缩减的最后一步，投射到子空间的方式和PCA是一样的，
那么PCA和LDA唯一的不同就是寻找主成分的理念不同了，

PCA强调的是在实现维度缩减的同时最大化保留的信息量
LDA强调的则是最小化类内的区别，最大化类间的区别

kNN with LDA

Thanks to having the information from the labels, LDA gives a discriminant space where the classes are far apart from each other.
This should help kNN a lot, as the job should just be to find the obvious 10 clusters.

However, as we’ve seen in the previous plot (section 4.3), the training data gets clustered pretty well, but the test data isn’t as nicely clustered as the training data (overfitting?).

Perform cross validation following the code below (you can change the values of k_choices and n_choices to search). Using the best result from cross validation, obtain the test accuracy.

LDA具有最大化类间差距最小化类内差距的功效，所以在分类问题能够有很大帮助

Object-Detection

1. 实现人脸检测
2. 使用Hog表征和滑窗法进行人脸检测
3. 使用图像金字塔对尺度问题进行改进
4. 使用DPM(可变形组件模型)进行人脸检测

In this homework, we will implement a simplified version of object detection process. Note that the tests on the notebook are not comprehensive, autograder will contain more tests.

from __future__ import print_function

Hog Representation

In this section, we will compute the average hog representation of human faces.

There are 31 aligned face images provided in the \face folder. They are all aligned and have the same size. We will get an average face from these images and compute a hog feature representation for the averaged face.

Use the hog function provided by skimage library, and implement a hog representation of objects.
Implement hog_feature function in detection.py

Hog表征
这里要计算人脸的平均hog表征。可以直接调用skimage包中feature模块内的hog函数进行hog特征的提取。

【特征检测】HOG特征算法

image_paths = fnmatch.filter(os.listdir('./face'), '*.jpg')

Sliding Window

Implement sliding_window function to have windows slide across an image with a specific window size. The window slides through the image and check if an object is detected with a high score at every location. These scores will generate a response map and you will be able to find the location of the window with the highest hog score.

滑窗法
滑窗法就是定义一个窗体对象（匹配对象），遍历图像中的图像块，计算卷积响应（此处对两者Hog特征序列进行内积计算），输出值最高的地方就是图像中与匹配对象最相似的地方。

image_path = 'image_0001.jpg'

Sliding window successfully found the human face in the above example. However, in the cell below, we are only changing the scale of the image, and you can see that sliding window does not work once the scale of the image is changed.

上面的图像能够检测到人脸的正确区域，但是下面这张图却出现了错误。说明滑窗这种方法对于图像的尺度十分敏感。

image_path = 'image_0001.jpg'

Image Pyramids

In order to make sliding window work for different scales of images, you need to implement image pyramids where you resize the image to different scales and run the sliding window method on each resized image. This way you scale the objects and can detect both small and large objects.

针对上面的问题，需要采用图像金字塔的方案来进行匹配。

图像金字塔分层算法

代码中直接是使用resize函数，生成了不同尺度的图像，具体的图像金字塔实现可以看博客，上采样和下采样是利用高斯滤波器和插值滤波器实现的。

Image Pyramid

Implement pyramid function in detection.py, this will create pyramid of images at different scales. Run the following code, and you will see the shape of the original image gets smaller until it reaches a minimum size.

image_path = 'image_0001.jpg'

Pyramid Score

After getting the image pyramid, we will run sliding window on all the images to find a place that gets the highest score. Implement pyramid_score function in detection.py. It will return the highest score and its related information in the image pyramids.

image_path = 'image_0001.jpg'

From the above example, we can see that image pyramid has fixed the problem of scaling. Then in the example below, we will try another image and implement deformable part model.

通过调整image = rescale(image, 1.2)中的缩放比例，可以看到，图像金字塔的方案还是能解决掉上一步中尺度引起的问题的。

image_path = 'image_0338.jpg'

Deformable Parts Detection

In order to solve the problem above, you will implement deformable parts model in this section, and apply it on human faces.

The first step is to get a detector for each part of the face, including left eye, right eye, nose and mouth.

For example for the left eye, we have provided the groundtruth location of left eyes for each image in the \face directory. This is stored in the lefteyes array with shape (n,2), each row is the (r,c) location of the center of left eye. You will then find the average hog representation of the left eyes in the images.

DPM（Deformable Part Model）原理详解

下面四块分别是左眼、右眼、鼻子和眼镜的检测器。
代码是作业本身就提供了的，不过可以仔细研究一下具体实现。

Run through the following code to get a detector for left eyes.

image_paths = fnmatch.filter(os.listdir('./face'), '*.jpg')

Run through the following code to get a detector for right eye.

righteye_h = 10

Run through the following code to get a detector for nose.

nose_h = 30

Run through the following code to get a detector for mouth

mouth_h = 20

Human Parts Location

Implement compute_displacement to get an average shift vector mu and standard deviation sigma for each part of the face. The vector mu is the distance from the main center, i.e the center of the face, to the center of the part.

找到各个部位的位置
这里计算的是相对位置，即各个组件的中心与脸部中心的距离。

# test for compute_displacement

Your implementation is correct!

lefteye_mu, lefteye_std = compute_displacement(lefteyes, face_shape)

After getting the shift vectors, we can run our detector on a test image. We will first run the following code to detect each part of left eye, right eye, nose and mouth in the image. You will see a response map for each of them.

分别使用脸部和各个组件的检测器，对整张图像进行匹配计算，得到每个组件的响应图。

image_path = 'image_0338.jpg'

max_score, lefteye_r, lefteye_c, lefteye_scale, lefteye_response_map = \

max_score, righteye_r, righteye_c, righteye_scale, righteye_response_map = \

max_score, nose_r, nose_c, nose_scale, nose_response_map = \

max_score, mouth_r, mouth_c, mouth_scale, mouth_response_map =\

After getting the response maps for each part of the face, we will shift these maps so that they all have the same center as the face. We have calculated the shift vector mu in compute_displacement, so we are shifting based on vector mu. Implement shift_heatmap function in detection.py.

得到各个组件的响应图之后，移动这些响应图，使得它们的中心和脸部的中心相同，方便求和。

face_heatmap_shifted = shift_heatmap(face_response_map, [0,0])

Gaussian Filter

In this part, apply gaussian filter convolution to each heatmap. Blur by kernel of standard deviation sigma, and then add the heatmaps of the parts with the heatmap of the face. On the combined heatmap, find the maximum value and its location. You can use function provided by skimage to implement gaussian_heatmap.

这个部分呢，就是将之前的各个组件的响应图给叠加起来，然后判断图像区域内哪个点的响应最大，则该处判断为人脸区域。

heatmap_face= face_heatmap_shifted

Result Analysis

Does your DPM work on detecting human faces? Can you think of a case where DPM may work better than the detector we had in part 3 (sliding window + image pyramid)? You can also have examples that are not faces.

DPM的原理是很能理解的，就是将各个组件的响应值累加起来，判断哪个区域拥有最高的响应值。
上图的检测出现了问题，怀疑是脸部组件的检测不咋滴。

Extra Credit

You have tried detecting one face from the image, and the next step is to extend it to detecting multiple occurences of the object. For example in the following image, how do you detect more than one face from your response map? Implement the function detect_multiple, and write code to visualize your detected faces in the cell below.

Tracking-OpticalFlow

使用光流实现目标追踪。
实现LK光流法。
使用高斯迭代法求解LK光流。
使用金字塔方法优化LK光流。

This assignment covers Lukas-Kanade tracking method. Please hand in motion.py and this notebook file to Gradescope.
这份作业是实现Lukas-Kanade光流追踪算法。

Displaying Video

We have done some cool stuff with static images in past assignemnts. Now, let’s turn our attention to videos! For this assignment, the videos are provided as time series of images. We also provide utility functions to load the image frames and visualize them as a short video clip.

Note: You may need to install video codec like FFmpeg. For Linux/Mac, you will be able to install ffmpeg using apt-get or brew. For Windows, you can find the installation instructions here.

from utils import animated_frames, load_frames

Lucas-Kanade Method for Optical Flow

Deriving optical flow equation

Optical flow methods are used to estimate motion at each pixel location in two consecutive image frames. The optical flow equation at a spatio-temporal point $p=(x,y,t)$ is given by:

$I_x(p)v_x + I_y(p)v_y + I_t(p) = 0$

where $Ix,I_y$ and $I_t$ are partial derivatives of pixel intensity $I$, and $v_x={\Delta x}/{\Delta t}$ and $v_{y}={\Delta y}/{\Delta t}$ are flow vectors.

Let us derive the equation in order to understand what it actually means. First, we make a reasonable assumption (a.k.a. brightness constancy) that the pixel intensity of a moving point stays the same between two consecutive frames with small time difference. Consider pixel intensity $I(x,y,t)$ of a point $(x,y)$ in the first frame $t$. Suppose that the point has moved to $(x+Δx,y+Δy)$ after $Δt$. According to the brightness constancy constraint, we can relate intensities of the point in the two frames using the following equation:

$I(x,y,t)=I(x+Δx,y+Δy,t+Δt)$

• a. Derive the optical flow equation from the brightness constancy equation. Clearly state any assumption you make during derivation.

光流公式的推导

利用灰度不变假设，我们有： $$I(x,y,t)=I(x+Δx,y+Δy,t+Δt)$$
对右边进行泰勒展开，保留一阶项，得：$$I(x+Δx,y+Δy,t+Δt)≈I(x,y,t)+I_xΔx+I_yΔy+I_tΔt$$
因为我们假设了灰度不变，所以下个时刻的灰度等于之前的灰度，从而：$$I_xΔx+I_yΔy+I_tΔt=0$$
两边除以$Δt$ ，得：$$I_x\frac{\Delta x}{\Delta t} + I_y\frac{\Delta y}{\Delta t} + I_t = 0$$
其中，$Δx/Δt$表示像素在x轴上的运动速度，$Δy/Δt$表示像素在y轴上的运动速度，用$v_x$和$v_y$表示，最终得$$I_x§v_x + I_y§v_y + I_t§ = 0$$

• b. Can the optical flow equation be solved given two consecutive frames without further assumption? Which values can be computed directly given two consecutive frames? Which values cannot be computed without additional information?

不能求解。

上式中，$I_x、I_y、I_t$可以经由连续的两幅图像求解，$I_x$表示图像在点$P$处$x$方向的梯度，$I_y$表示图像在点$P$处$y$方向的梯度，$I_t$表示图像灰度对时间的变化量。
但是仍有两个未知变量$v_x$和$v_y$，最终得到的是一个两个变量的一次方程。
无法求解，需要引入其他约束。

Overview of Lucas-Kanade method

The Lucas–Kanade method assumes that the motion of the image contents between two frames is approximately constant within a neighborhood of the point p under consideration (spatial coherence).

Consider a neighborhood of $p$, $N(p)=p_1,…,p_n$ (e.g. 3x3 window around $p$). According to the optical flow equation and spatial coherence assumption, the following should be satisfied:

For every $p_i∈N(p)$, $$I_x(p_i)v_x+I_y(p_i)v_y=−I_t(p_i)$$

These equations can be written in matrix form $Av=b$, where $$A = \begin{bmatrix}I_x(p_1) & I_y(p_1) \ I_x(p_2) & I_y(p_2) \ \vdots & \vdots \ I_x(p_n) & I_y(p_n) \end{bmatrix} \hspace{1em} v = \begin{bmatrix}v_x \ v_y\end{bmatrix} \hspace{1em} b = \begin{bmatrix}-I_t(p_1) \ -I_t(p_2) \ -I_t(p_n)\end{bmatrix}$$

Note that this linear system may not have solution for v as it is usually over-determined. Instead, Lucas-Kanade method estimates the flow vector by solving the least-squares problem, $A^TAv=A^Tb$.

• a. What is the condition for this equation to be solvable?
• b. Reason about why Harris corners might be good features to track using Lucas-Kanade method.

Lucas-Kanade算法的思想
上一个光流公式的推导可以得知：求解光流公式需要引入其他约束，而LK光流法的约束就是，假设某一个窗口内的像素具有相同的运动。

如果是采用3x3的窗口，则会另有9个公式，此时采用最小二乘来估计最优的$v_x$和$v_y$即可。

问题a：方程求解的条件是什么？参见lecture17 PPT

1. $A^TA$可逆
2. 因为噪声的存在, $A^TA$不能太小

$A^TA$的特征值$λ_1$和$λ_2$都不能太小

3. $A^TA$应当是well-conditioned

$\lambda_1/\lambda_2$不能太大($\lambda_1$是较大的特征值)

问题b：为何Harris角点会是LK光流法追踪的好特征？参见lecture17 PPT第25页
因为$M=A^TA$是二阶矩矩阵，$A^TA$的特征向量和特征值与边缘的方向和强度相关。

• 较大的特征值对应的特征向量表示强度变化最快的方向
• 另外一个特征向量则与它正交

Implementation of Lucas-Kanade method

In this section, we are going to implement basic Lucas-Kanade method for feature tracking. In order to do so, we first need to find keypoints to track. Harris corner detector is commonly used to initialize the keypoints to track with Lucas-Kanade method. For this assignment, we are going to use skimage implementation of Harris corner detector.

实现基础版本LK光流
第一步需要找到关键点来进行追踪。
Lucas-Kanade方法通常是使用Harris角点检测器来进行关键点的初始化。
skimage中提供了Harris角点检测的实现。

skimage中的Harris角点检测可以参考：Programming Computer Vision with Python （学习笔记九）

from skimage import filters

Implement function lucas_kanade in motion.py and run the code cell below. You will be able to see small arrows pointing towards the directions where keypoints are moving.

实现两帧之间的LK光流追踪

from motion import lucas_kanade

We can estimate the position of the keypoints in the next frame by adding the flow vectors to the keypoints.
通过当前帧关键点的位置和光流的方向，就可以预测下一帧中的关键点的位置

# Plot tracked kepoints

上图可以看到，光流法预测的第二帧中的特征点地位置还是基本正确的。

Feature Tracking in multiple frames

Now we can use Lucas-Kanade method to track keypoints across multiple frames. The idea is simple: compute flow vectors at keypoints in $i$-th frame, and add the flow vectors to the points to keep track of the points in $i+1$-th frame. We have provided the function track_features for you. First, run the code cell below. You will notice that some of the points just drift away and are not tracked very well.

Instead of keeping these ‘bad’ tracks, we would want to somehow declare some points are ‘lost’ and just discard them. One simple way to is to compare the patches around tracked points in two subsequent frames. If the patch around a point is NOT similar to the patch around the corresponding point in the next frame, then we declare the point to be lost. Here, we are going to use mean squared error between two normalized patches as the criterion for lost tracks.

Implement compute_error in motion.py, and re-run the code cell below. You will see many of the points disappearing in later frames.

作业中已经提供了追踪连续图像帧中的特征点的功能，在track_features函数中。
其中思路就是根据初始的关键点位置和相邻两帧LK光流法得到的光流方向，计算预测后续每一帧中的关键点位置。
但是由于出现一些关键点的漂移（计算精确度以及关键点的消失），整体的追踪效果会很差。
为了去除这些“坏点” 的影响，提出了一种思路就是：判断关键点的丢失。

原理就是：判断相邻两帧关键点区域图像块的相似程度，如果不相似则判断为该点丢。它的依据就是Lucas-Kanade的假设，某一窗口内的像素具有相同的移动。

计算方式：计算两个归一化的图像块的均方差MSE，如果大于阈值，则抛弃该点。

from utils import animated_scatter

从上面的运行效果看，保留下来的关键点都能够准确的被预测到。
但是也可以看到就是，当运动幅度较大的时候，好多关键点都被丢弃了。

Pyramidal Lucas-Kanade Feature Tracker

In this section, we are going to implement a simpler version of the method described in “Pyramidal Implementation of the Lucas Kanade Feature Tracker”.

金字塔光流法

TLD算法学习之L-K光流法理论篇二

Iterative Lucas-Kanade method

One limitation of the naive Lucas-Kanade method is that it cannot track large motions between frames. You might have noticed that the resulting flow vectors (blue arrows) in the previous section are too small that the tracked keypoints are slightly off from where they should be. In order to address this problem, we can iteratively refine the estimated optical flow vectors. Below is the step-by-step description of the algorithm:

Let $p=[p_x\space p_y]^T$ be a point on frame $I$. The goal is to find flow vector $v=[v_x\space v_y]^T$ such that $p+v$ is the corresponding point of $p$ on the next frame $J$.

• Initialize flow vector: $$v = \begin{bmatrix}0\0\end{bmatrix}$$
• Compute spatial gradient matrix:

$G = \sum_{x=p_x-w}^{p_x + w}\sum_{y=p_y -w}^{p_y+w}\begin{bmatrix}I_x^2(x,y) & I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y) & I_y^2(x,y)\end{bmatrix}$

• for $k=1$ to $K$

• Compute temporal difference: $δI_k(x,y)=I(x,y)−J(x+g_x+v_x,y+g_y+v_y)$
• Compute image mismatch vector

$b_k = \sum_{x=p_x-w}^{p_x + w}\sum_{y=p_y -w}^{p_y+w}\begin{bmatrix}δI_k(x,y)I_x(x,y) \\ δI_k(x,y)I_y(x,y) \end{bmatrix}$

Compute optical flow: $v^k=G^{−1}b_k$
Update flow vector for next iteration: $v:=v+v^k$

• Return $v$

Implement iterative_lucas_kanade method in motion.py and run the code cell below. You should be able to see slightly longer arrows in the visualization.

from motion import iterative_lucas_kanade

可以看到，迭代式光流法可以检测到大运动的光流。

# Plot tracked kepoints

# Detect keypoints to track in the first frame

运行到最后，相比于之前LK光流，迭代式光流法追踪并保留了更多的关键点。

Coarse-to-Fine Optical Flow

The iterative method still could not track larger motions. If we downscaled the images, larger displacements would become easier to track. On the otherhand, smaller motions would become more difficult to track as we lose details in the images. To address this problem, we can represent images in multi-scale, and compute flow vectors from coarse to fine scale.

Run the following code cell to visualize image pyramid.

的迭代金字塔方式依旧不能追踪更大的运动。

原因是因为，如果进行下采样得到的图像，大的运行是更容易追踪到了，但是另一方面，较小的运动却被忽略掉了，因为这部分信息被下采样给丢失了。

为了解决这个问题，这里使用多尺度来表征图像，并且从粗粒度到细粒度来计算光流。

#下面这段函数是使用skimage提供的函数来生成不同尺度的图像。

Following is the description of pyramidal Lucas-Kanade algorithm:

Let $p$ be a point on image $I$ and $s$ be the scale of pyramid representation.

• Build pyramid representations of $I$ and $J$: ${I^L}\{L=0,…,L_m\}$and${J^L}\{L=0,…,L_m\}$
• Initialize pyramidal guess $$g^{L_m}=[g{L_m}_xg^{L_m}_y]T=[0\space 0]^T$$
• for $L=L_m$ to $0$ with step of $-1$
  • Compute location of $p$ on $I^L$: $p^L=p/s^L$
  • Let $d^L$ be the optical flow vector at level $L$:$$d^{L:=IterativeLucasKanade(I}L,J^L,pL,g^L)$$
  • Guess for next level L−1: $g^{L−1}=s(g^L+d^L)$
• Return $d=g^0+d^0$

Implement pyramid_lucas_kanade.

上面的是金字塔法Lukas-Kanade光流的算法描述

from motion import pyramid_lucas_kanade

上图可知，金字塔方式LK光流既能检测到大运行，也能检测到小运动。但是幅度也太大了点。

# Plot tracked kepoints

from utils import animated_scatter

总体比较下来，金字塔法光流耗时最多。
效果也不稳定，貌似第二帧关键点都追丢了，怀疑是计算得到的光流幅度太大引起的。

Object Tracking

Let us build a simple object tracker using the Lucas-Kanade method we have implemented in previous sections. In order to test the object tracker, we provide you a short face-tracking sequence. Each frame in the sequence is annotated with the ground-truth location (as bounding box) of face.

An object tracker is given an object bounding box in the first frame, and it has to track the object by predicting bounding boxes in subsequent frames.

目标追踪

在第一帧中提供了目标的外包围框，要在接下来的帧序列中预测目标运行。

from utils import animated_bbox, load_bboxes

In order to track the object, we first find keypoints to track inside the bounding box. Then, we track those points in each of the following frames and output a tight bounding box around the tracked points. In order to prevent all the keypoints being lost, we detect new keypoints within the bounding box every 20 frames.

找到包围盒中的关键点。
在接下来的帧中追踪这些关键点，并用包围盒框住这些关键点。
为了防止关键点的消失，每二十帧进行一次关键点的检测。

# Find features to track within the bounding box

from motion import compute_error

ani = animated_bbox(frames, bboxes)

金字塔光流法
到最后竟然包围盒消失了，两个原因吧

• 一个是特征点的检测算法对于光照敏感，
• 一个是光流追踪没有完善。
  其实，上一步金字塔光流法其实还没有解决问题哩。

迭代光流法
IOU分数能达到0.32，但是明显的错误在于，在过程中更新关键点时，将复杂背景下的一些脸部之外的关键点包括了进来，导致了最后的追踪不仅包含了最开始的目标，还包括了复杂背景的一些关键点。

Evaluating Object Tracker: intersection over union (IoU)

Intersection over union is a common metric for evaluating performance of an object tracker. Implement IoU in motion.py to evaluate our object tracker. With default parameters, you will get IoU score of ~0.32.

根据真值，评估自己的追踪效果。
计算方法看图：

Intersection over Union (IoU) for object detection

from motion import IoU

0.32643400622199703

Extra Credit: Optimizing Object Tracker

Optimize object tracker in the code cell below. You may modify the code and define new parameters. We will grant extra credit for IoU score > 0.45.

这里要求改进，提高准确率。

思路就是,从上一视频的问题入手，错误的一个明显原因就是引入了复杂背景的一些关键点，所以一个解决方案就是避免引入这些关键点。

方案一：

• 每隔20帧重新获取关键点
• 对新的关键点进行筛选，保留与初始关键点最相似的点。其实就是根据找特征点的匹配。

这就需要对关键点进行描述了。当然，其实也可以对Harris检测角点这一方法换成别的关键点检测算法来做。