计算机视觉：作业4图像匹配，缝合和单应性

16-720A计算机视觉：作业4图像匹配，缝合和单应性

讲师：Kris M. Kitani助教：Leonid，Mohit，Arjun，Rawal，Aashi，Tanya
截止日期：2018年3月27日星期二晚上11:59
版本1(3月8日，12pm）

总积分：100
额外积分：30

兴趣点检测器和描述符是大多数计算机视觉应用程序的核心。 这项任务的目的是通过研究一种确定单应性的应用程序来获得见识。 我们首先快速回顾一下，建立对平面单应性背后的理论的基本理解，然后继续展示如何在全景照片拼接和增强现实应用中常用的模板图像和源图像之间进行转换。

与这种情况下的先前作业相比，我们提供了使用Python或MATLAB进行作业的机会。 请选择一种语言，然后用该语言一致地完成整个作业。

图1：在分配任务结束时，您将能够将多张照片拼接成无缝的单张照片，其源图像来自M. Kazhdan

使用说明
1.诚信与协作：如果您是一个小组，请在您的文章中注明合作者的姓名。代码不应共享或复制。除非允许，否则请勿使用外部代码。严禁抄袭，否则可能会导致本课程失败。

2.早点开始！特别是如果您想最后使用AR应用程序！

3.广场：如有任何疑问，请先查看广场。我们为每个问题创建了文件夹，并在每个问题的顶部创建了一个粘性线程。在提交新问题并将其标签添加到即时贴之前，请仔细阅读所有在适当评论中标记的帖子。

4.书面记录：请注意，在此作业中，我们不接受您书面记录的手写扫描。请以电子方式输入对理论问题和实验讨论的答案。

5.代码：请坚持讲义中提到的功能原型。这使得TA的验证代码更加容易。

6.提交：请将所有结果包括在您在Gradescope上提交的pdf撰写文件中。对于代码提交，创建一个组成您的Matlab文件的zip文件<andrew-id> .zip。请确保删除您生成的所有临时文件。
最终上传的文件应按此布局排列。如果您是在MATLAB而不是Python中完成作业，则可以随意忽略python文件夹。
同样，如果您使用Python而不是MATLAB，那么。Pl

 <AndrewId>.zip
{ <AndrewId>/
 matlab/

makeTestPattern.m

computeBrief.m

briefLite.m

testMatch.m

briefRotTest.m

computeH.m

computeH norm.m

computeH ransac.m

HarryPotterize.m

compositeH.m

imageStitching.m

imageStitching noClip.m
 python/

q2.py

q3.py

q4.py

ec.py

run q2.py

run q3.py

run q4.py

run ec.py

7.负责此项任务的TA：Mohit Sharma(mohits1@andrew.cmu.edu）和Leonid Keselman(lkeselma@andrew.cmu.edu）。

1 同形(30分）

假设我们有两个摄像机C1和C2，它们看着3D空间中的公共平面Ⅱ。 Ⅱ上的任何3D点P都会在第一个摄像机C₁上的p =(u1; v1; 1）T处生成投影的2D点，而在第二个摄像机C₂上的q =(u2; v2; 1）T处会生成投影2D点。由于P限于平面Ⅱ，我们期望p和q之间存在关系。特别是，存在一个公共的3 * 3矩阵H，因此对于任何p和q，以下条件均成立：

我们称这种关系为平面单应性。回想一下p和q都在齐次坐标中，并且等式=意味着p与Hq成比例(回想齐次坐标）。事实证明，这种关系对于通过纯旋转而没有平面约束的摄像机也适用。

Q1.1同构(5分）
证明存在一个H，它满足给定的两个3 * 4摄像机投影矩阵M₁和M₂(分别对应于摄像机C₁，C₂和平面Ⅱ），满足方程1。不要为H生成实际的代数表达式。我们所要的只是证明H的存在。

注意：当平面Ⅱ包含两个相机的中心时，可能会发生退化的情况，在这种情况下，可以满足公式1的H的无限选择。您可以在答案中忽略这种情况。

Q1.2旋转下的单应性(5分）
证明给定两个由纯旋转分开的相机，存在满足x₁ = Hx₂的单应性H。也就是说，对于摄像机1，x₁ = K1 [I 0] X，对于摄像机2，x₂ = K₂ [R 0]X。请注意，K₁和K₂是两个摄像机的3 * 3本征矩阵，并且不同。 I是3 * 3单位矩阵，0是3 * 1零向量，X是3D空间中的一个点。 R是相机的3 * 3旋转矩阵。

Q1.3对应(10分）
令x1是图像中的一组点，而x2是另一台相机拍摄的图像中的相应点。假设存在一个单应性H：

其中x₁ⁱ = [x₁ⁱ y₁ⁱ 1] ^T在同构坐标中，x₁ⁱ ∈ x₁，H是3 * 3矩阵。对于每个点对，此关系可以重写为

A_ih = 0
其中h是从H重塑的列向量，而A_i是具有从点x₁ⁱ和x₂ⁱ派生的元素的矩阵。这可以帮助从给定的点对应关系计算H。

1. h有多少个自由度？ (3分）
2.求解h需要多少个点对？ (2分）
3.派生A_i(5分）

Q1.4了解旋转下的单应性(5分）

假设摄像机围绕其中心C旋转，则保持固有参数K不变。令H是将视图从一个相机方向映射到第二方向上的视图的单应性。令H为两者之间的旋转角度。证明H²是对应于2θ旋转的单应性。请把答案限制在两行之内。冗长的证据表明您做的事情太复杂(或错误）。

Q1.5平面单应性的局限性(2分）
为什么平面单应性还不足以将任意场景图像映射到另一个视点？用一两句话简洁地陈述您的答案。

Q1.6透视投影下的线的行为(3分）
我们在课堂上说过，透视投影保留直线(3D的直线投影为2D的直线）。代数验证情况是否如此，即验证x = PX中的投影P是否保留线条。

注意：兴趣点
在找到图像对之间的单应性之前，我们需要找到两个图像之间的对应点对。但是我们如何获得这些积分呢？一种方法是手动选择它们(使用cpselect），这很繁琐且效率低下。 CV方法是在图像对中找到兴趣点并自动进行匹配。为了能够做一些有趣的事情，我们在这里不会重新实现特征检测器，而是使用内置的MATLAB方法。兴趣点检测器(例如Harris，FAST，SIFT，SURF等）的目的是在图像中找到特定的显着点，在该显着点附近我们提取特征描述符(例如Brief，ORB等）。这些描述符尝试以简洁而又具有描述性的方式总结特征点周围的图像内容(对于许多计算机视觉任务，通常在表示和计算复杂度之间进行权衡；您可以使用一个非常高维度的特征描述符，可以确保您获得良好的匹配，但计算成本可能会高得惊人）。然后，匹配是一项任务，即尝试在与当前描述符最匹配的新图像上计算出的描述符列表中找到一个描述符。这可能就像两个描述符之间的欧几里德距离一样简单，也可能更复杂，这取决于描述符的组成方式。

出于此练习的目的，您可以将FAST或Harris检测器与Brief描述符一起使用。您可以使用Matlab内置函数detectFASTFeatures或detectHarrisFeatures来计算兴趣点。同样，在Python中，您可以快速使用skimage.feature形式的函数，例如skimage.feature.corner harris或skimage.feature.corner。

2 简述(25分）
现在我们有了兴趣点，告诉我们在哪里可以找到图像中最有用的点，我们可以计算出可用于与不同图像中同一点的其他视图匹配的描述符。 Brief描述符以兴趣点的特征尺度为中心，以兴趣点为中心的9 * 9块p编码信息。如果您需要刷新内存，请参见点特征检测器的讲义。

2.1创建一组简短测试(5分）
描述符本身是一个n位长的向量，其中每个位是以下简单测试的结果：


将n设置为256位。无需将测试结果编码为实际位。它是细将它们编码为256个元素的矢量。

用于计算t(p; x; y）的256个测试对(x; y）有很多选择(n位中的每一个）。作者在[1]中描述和测试了其中的一些。阅读该文件的第3.2节，并实现这些解决方案之一。您应该生成一组静态的测试对，并将该数据保存到文件中。您将使用这些对进行随后的Brief描述符的所有计算。

Q 2.1：编写函数以创建x和y对，我们将使用它们对进行比较以计算t：
[compareX，compareY] = makeTestPattern(patchWidth，nbits）

patchWidth是图像补丁的宽度(通常为9），nbits是Brief描述符中的测试次数n。 compareX和compareY是patchWidth * patchWidth图像补丁的线性索引，每个都是nbit * 1个向量。对给定参数patchWidth = 9和n = 256运行此例程，然后将结果保存在testPattern.mat中。在提交中包括此文件。

2.2计算简要描述符(5分）
现在是时候为检测到的关键点计算BRIEF描述符了。

问题2.2：编写函数：

[locs，desc] = computeBrief(im，locs，compareX，compareY）

其中im是值为0到1的灰度图像，locs是哈里斯(Harris）关键点检测器返回的关键点位置，compareX和compareY是在2.1节中计算的测试模式，并保存到testPattern.mat中。

该函数返回locs(一个m * 2向量），其中前两列是关键点的图像坐标，而desc则是一个堆叠的Brief描述符的m * n位矩阵。 m是图像中有效描述符的数量，并且会有所不同。您可能需要注意输入检测器的位置，因为它们可能位于图像的边缘，在该边缘我们无法提取宽度patchWidth的完整块。因此，输出位置的数量可以少于输入位置。

2.3全部放在一起(5分）
问题2.3：编写函数：
[locs，desc] = briefLite(im），

它接受值为0到1之间的灰度图像im并返回loc(m * 2矢量），其中前两列是关键点的图像坐标，desc是堆叠的Brief描述符的m * n位矩阵。 m是图像中有效描述符的数量，并且会有所不同。 n是Brief描述符的位数。

此功能应执行所有必要步骤以从图像中提取描述符，包括
Θ加载参数和测试模式
Θ获取关键点位置
Θ计算一组有效的Brief描述符

2.4检查点：描述符匹配(5分）
尽管视图，光照等发生了变化，但描述符的优势在于能够与同一世界点生成的其他描述符匹配。用于计算两个描述符之间相似度的距离度量至关重要。对于简而言之，该距离度量是汉明距离。汉明距离只是两个描述符中不同的位数。 (请注意，位的位置很重要。）

为了执行上述描述符匹配，我们在matlab / briefMatch.m中提供了以下功能：

[matches] = briefMatch(desc1, desc2, ratio)

它接受来自第一张图片的m₁ * n位摘要描述符的堆栈和来自第二张图片的m₂ * n位摘要描述符的堆栈，并返回一个 p * 2匹配矩阵，其中第一列是desc1的索引，第二列是是desc2的索引。请注意，m₁，m₂和p的大小可能不同，并且p <= min(m₁; m₂）。

Q 2.4.1：编写测试脚本testMatch来加载两个鸡汤图像，计算特征匹配。使用提供的plotMatches和briefMatch函数可视化结果。

plotMatches(im1, im2, matches, locs1, locs2)

其中im1和im2是从0到1的灰度图像，matches是briefMatch返回的匹配项列表，而locs1和locs2是briefLite中关键点的位置。

保存结果图并将其提交到PDF中。简要讨论任何效果较差或更好的情况。

图2是示例结果。调试建议：对代码的一个很好的测试是检查您是否可以将图像与其自身进行匹配。

2.5简介和轮换(5分）
您可能已经注意到旋转下的性能较差。让我们对此进行调查！

问2.5：拍摄模型ckenbroth.jpg测试图像并将其与自身匹配，同时以10度为增量旋转第二个图像(提示：旋转）。计算每次旋转时正确匹配的次数，并构建一个显示旋转角度与正确匹配次数的条形图。将其包含在PDF中，并解释为什么您认为描述符具有这种行为。创建执行此任务的脚本briefRotTest.m。

图2：模型chickenbroth.jpg和chickenbroth 01.jpg的简要匹配示例。

2.6提升绩效-(额外功劳）

下面介绍的额外信用机会是可选的，它们提供了探索计算机视觉并改善上面开发的技术性能的途径。

1.(5分）正如我们所看到的，“摘要”不是旋转不变的。使用到目前为止开发的工具设计一个简单的解决方案来解决此问题(回想边缘检测和/或Harris角的协方差矩阵）。在PDF中说明您的设计决策以及如何选择使用的任何参数。演示算法在与大旋转相关的图像对上的有效性。

2.(5分）此Brief的实现具有一定的规模不变性，但有其局限性。当您将一张图片与一半大小的图片匹配到同一张图片时会发生什么？请参阅[2]的第3节，该技术将使您的检测器对比例变化更加鲁棒。实施并通过几个测试图像进​​行演示。您可以简单地重新缩放我们提供给您的一些测试图像。

3单应性计算(35分）
Q3.1计算单应性(10分）

编写一个函数computeH，该函数根据一组匹配的点对来估计平面单应性。

function [H2to1] = computeH(x1, x2)

x₁和x₂是N * 2矩阵，其中包含两个图像之间的点对坐标(x; y）。 H2to1应该是3 * 3矩阵，以便以最小二乘的方式从图像2到图像1获得最佳单应性。这应该从1中派生的矩阵得出。对于求解器，请随意使用您喜欢的任何东西。尽管本讲义末尾的SVD注释是更直接的方法之一。

单应归一化
归一化可提高解决方案的数值稳定性，因此应始终对坐标数据进行归一化。标准化有两个步骤：
1.将点的平均值转换为原点。
2.缩放点，使到原点的最大距离为√2。
这是线性变换

Q3.2归一化的单应性(5分）
实现函数computeH规范：

function [H2to1] = computeH norm(x1, x2).

该函数应该对x₁和x₂中的坐标进行归一化，并如上所述调用calculateH(x₁，x₂）。

RANSAC
RANSAC算法通常可以将任何模型拟合到嘈杂的数据。您将为图像之间的(平面）单应图实现它。

Q3.3实施RANSAC计算单应性(10分）

编写一个函数：
function [bestH2to1, inliers] = computeH ransac(locs1, locs2)

其中bestH2to1应该是在RANSAC期间发现的大多数内点的单应性H。 H将是一个单应性，如果x₂是locs2中的一个点，而x₁是locs1中的一个对应的点，则x₁ = Hx₂。 locs1和locs2是包含匹配点的N * 2个矩阵。 inliers是长度为N的向量，在属于共识集的那些匹配项中为1，在其他情况下为0。使用computeH范数计算单应性。

自动化的全息术估计/扭曲以增强现实

Q3.4放在一起(10分）

编写一个脚本HarryPotterize.m
1.读取cv cover.jpg，cv desk.png和hp cover.jpg。
2.使用computeH ransac自动计算单应性。
3.使用提供的warpH函数将hp cover.jpg变形为cv desk.png图像的尺寸。
4.此时，您应注意，尽管图像被扭曲到正确的位置，但它并没有填充与书本相同的空间。您为什么认为这种情况正在发生？您将如何修改hp cover.jpg以解决此问题？
5.实现函数：

function [ composite img ] = compositeH( H2to1, template, img )

现在将这个扭曲的图像与桌面图像组成，如图4所示。
6.将您生成的图像包括在您的文章中。另请在您的文章中打印最终的H矩阵(已规范化，因此右下角的值为1）

图3：课本　　　　　　图4：HarryPotterized教科书

图5：问题4.1的示例输出：原始图像img1和img2(左和中）和img2变形为适合img1(右）。 请注意，变形的图像会裁剪出图像。 我们将在Q 4.2中解决此问题

4 拼接全景图(10分）

我们还可以使用同形异义图从同一场景的多个视图创建全景图像。例如，当视图之间没有摄像机平移时(例如，仅围绕摄像机中心旋转），这是可能的。首先，您将使用简短匹配使用图像之间的匹配点对应关系生成全景图。我们将假定图像之间的匹配点对应关系中没有错误(尽管可能存在一些错误，甚至很小的错误也会产生巨大的影响）。在下一部分中，您将扩展该技术以处理嘈杂的关键点匹配。

您将需要使用提供的 function warp im = warpH(im, H, out size),

该函数使用单应变换H扭曲图像im。warp_im中的像素在矩形(1; 1）到(out_size）的坐标处进行采样(2）; out_size(1））。将源图像中像素的坐标取为(1; 1）到(size(im，2）; size(im，1））并根据H进行变换。

Q 4.1(5分）

在此问题中，您将实现并使用以下功能：

[panoImg] = imageStitching(img1, img2, H2to1)

在Dusquesne斜面上的两个图像上。该功能接受两幅图像，并从单应性估计功能接收输出。该功能：
1.使用提供的warpH(）函数将img2变形为img1的参考框架
2.混合img1和变形的img2，然后输出全景图像。

对于此问题，请使用提供的图像pnc1作为img1和pnc0作为img2。点对应点pts由您的简要描述符匹配生成。

将ransacH(）应用于这些对应关系以计算H2to1，这是从pnc0到pnc1的单应性。然后使用warpH(）将单应性应用于pnc0。
注意：由于扭曲的图像将被向右平移，因此您将需要更大的目标图像。

可视化变形的图像。请在撰写的PDF文件中包含图片和H2to1矩阵(右下角索引为1）。 (以及说明您使用了哪个图像对）

Q 4.2(3分）

请注意Q 4.1的输出如何修剪到边缘？我们现在将解决此问题。实现功能

[panoImg] = imageStitching noClip(img1，img2，H2to1）

它采用与Q 4.1中相同的输入类型并产生相同的输出。

为了防止边缘出现剪切，我们需要将图像1和图像2都扭曲到一个公共的第三参考帧中，在该参考帧中可以显示两个图像而没有任何剪切。具体来说，我们想找到一个仅进行缩放和平移的矩阵M：

warp im1 = warpH(im1, M, out size);
warp im2 = warpH(im2, M*H2to1, out size);

这样会在公共参考系中生成变形图像，在该参考系中im1和im2中的所有点均可见。为此，我们将仅采用尺寸超出的宽度或高度作为输入，并根据给定图像计算另一个尺寸或高度，以使变形图像不会在全景图像中被挤压或拉长。现在，假设我们仅将图像的宽度(即out size(2））作为输入，因此应该计算正确的高度(即out size(1））。

提示：将根据H2to1和两个图像的极点(角）进行计算。确保M仅包含比例尺(查找全尺寸全景图像的纵横比）和平移。

可视化变形的图像。请将该图像包括在您的撰写PDF中。 (以及说明您使用了哪个图像对）

问题4.3(2分）：您现在拥有自动生成全景图所需的所有工具。编写一个函数，接受两个图像作为输入，为两个图像计算关键点和描述符，通过匹配关键点描述符找到推定的特征对应关系，使用RANSAC估计单应性，然后将其中一个图像与该应同性翘曲以使其对齐，然后对齐覆盖它们。

im3 = generatePanorama(im1，im2）

在图像对data / pnc1.jpg，data / pnc0.jpg或data / incline L.jpg，data / incline R.jpg上运行代码。但是，在调试期间，请尝试按比例缩小图像的版本，以保持运行时间短。将生成的全景图保存在全尺寸图像上，并在插图中包含图形和计算的单应性矩阵。

Q 4.4(3分）额外信用：收集一对您自己的图像(使用手机），并使用上一节中的代码将它们拼接在一起。在记录中包括图像对及其结果。

Q 4.5(2分）额外积分：收集至少6张图像并将它们缝合成一张noClip图像。您可以收集自己的图像，也可以使用Matt Uyttendaele的PNC公园图像。我们将PNC公园图像(约1/4大小），ORB关键点和描述符用于我们的参考解决方案。

图6：最终全景视图。 单应性估计使用RANSAC。

5 额外积分(15分）(15分）

为了获得更多荣誉，我们将实施一种更智能的图像混合形式，将接缝最小化。这可用于处理针迹图像之间的曝光差异，校正运动模糊，甚至可以用作photoshop修复工具。有关照片拼接的示例，请参见http://www.cs.jhu.edu/~misha/Code/SMG/ or http://www.cs.jhu.edu/~misha/Code/DMG/Version4.5/

图7：包含图像的泊松拼接示例。

5.1问题5.1(15分）

在本节中，我们将实现Poisson Image Stitching [3]，可用于创建无缝的照片蒙太奇。

编写一个名为poisson blend(background，foreground，mask）的函数，该函数获取3个大小相等的图像(背景和前景为RGB，掩模为二进制），并使用前景中的渐变和背景中的边界条件求解泊松方程。请在撰写本文时同时包含fg1，bg1，mask1和fg2，bg2，mask2图像的结果。

泊松图像拼接背后的基本思想是建立一个线性系统Ax = b，其中A的每一行都是图像像素蒙版的梯度算符(例如，每个像素4个，四个像素-1 +1），而b中的行包含从前景开始的已知渐变。 A中的其他行是边界上的像素，它们具有已知值，存储在b中。然后求解x会​​为您提供整个蒙版色块的像素值。分别对每个图像通道执行此操作。请参阅本文，本页或此页面(如果需要更多信息）。

参考文献
[1] Michael Calonder et al. BRIEF: Binary robust independent elementary features". In: ECCV. 2010. url: https://www.robots.ox.ac.uk/~vgg/rg/papers/CalonderLSF10.pdf.
[2] David G. Lowe. Distinctive Image Features from Scale-Invariant Keypoints".In: IJCV (2004). url: https://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf.
[3] Patrick Perez, Michel Gangnet, and Andrew Blake. Poisson image editing".In: TOG (2003). url: https://www.cs.virginia.edu/~connelly/class/2014/comp_photo/proj2/poisson.pdf.

6 SVD注意事项
矩阵A的奇异值分解(SVD）表示为：

这里，U是称为“左奇异矢量”的列矢量的矩阵；类似地，V被称为“右奇异矢量”；矩阵Σ是对角矩阵；每个对角元素σi被称为“奇异值”，这些按大小顺序排序。
在我们的例子中，它是一个9 * 9的矩阵。

Θ如果σ₉= 0，则系统是精确确定的，存在单应性，并且所有点都精确拟合。
Θ如果σ₉> = 0，则表示系统超定。存在单应性，但并非所有点都完全适合(它们在最小二乘误差意义上适合）。此值表示合身性。
Θ通常，您将至少有四个对应关系。如果不是，则系统未确定。我们不会在这里处理这些。

U的列是AA^T的特征向量。 V的列是A^TA的特征向量。我们可以使用该事实来求解方程式Ah = 0中的h。使用此知识，让我们重新构造求解Ax = 0的问题。我们希望最小二乘意义上的求解误差最小。理想情况下，乘积Ah应该为0。
因此，平方和误差可以写成：

相对于h最小化此错误，我们得到：

这意味着h的值等于对应于零特征值的特征向量(或在出现噪声的情况下最接近零）。因此，我们选择A^TA的最小特征值，即Σ中的σ₉，Ah = 0的最小二乘解是对应的特征向量(在矩阵V的第9列）。