计算机视觉（三）

一 单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表示。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，我们将频繁的使用单应性变换。本质上，单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：

或者

对于图像平面内的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标可以简单地使用一个矩阵来表示变换。
下面的函数可以实现对点进行归一化和转换齐次坐标的功能：

from numpy import *

def normalize(points):
    for row in points:
        row /= points[-1]
    return points


def make_homog(points):
    return vstack((points ,ones((1 ,points.shape[1]))))

二 直接线性变换算法

单应性矩阵可以有两幅图像（或者平面）中对应点对计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要4个对应点对。

DLT（Direct Linear Transformation，直接线性变换）是给定4个点或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点上，重新写出该方程，可以得到下面的方程：

def H_from_points(fp, tp):
    """使用线性DLT方法，计算单应性矩阵H，使fp映射到tp。点自动进行归一化"""

    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')

    # 对点进行归一化（对数值计算很重要）
    # --- 映射起始点 ---
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1)) + 1e-9
    C1 = diag([1/maxstd, 1/maxstd, 1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp = dot(C1,fp)
    
    # --- 映射对应点 ---
    m = mean(tp[:2], axis=1)
    maxstd = max(std(tp[:2], axis=1)) + 1e-9
    C2 = diag([1 / maxstd, 1 / maxstd, 1])
    C2[0][2] = -m[0] / maxstd
    C2[1][2] = -m[1] / maxstd
    tp = dot(C2, tp)
    
    # 创建用于线性方法的矩阵，对于每个对应对，在矩阵中会出现两行数值
    nbr_correspondences = fp.shape[1]
    A = zeros((2 * nbr_correspondences, 9))
    for i in range(nbr_correspondences):
        A[2*i] = [-fp[0][i], -fp[1][i],-1,0,0,0,
                  tp[0][i]*fp[0][i],tp[0][i]*fp[1][i],tp[0][i]]
        A[2*i+1] = [0,0,0,-fp[0][i],-fp[1][i],-1,
                    tp[1][i]*fp[0][i],tp[1][i]*fp[1][i],tp[1][i]]
        
    U,S,V = linalg.svd(A)
    H = V[8].reshape((3,3))
    
    #反归一化
    H = dot(linalg.inv(C2),dot(H,C1))
    
    #归一化，然后返回
    return H / H[2,2]

三 仿射变换

由于仿射变换具有6个自由度，因此我们需要三个对应点来估计矩阵H。通过将最后两个元素设置为0，即h7=h8=0，仿射变换可以用上面的DLT算法估计得出。

def Haffine_from_points(fp,tp):
    
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
        
    # 映射起始点
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1)) + 1e-9
    C1 = diag([1/maxstd, 1/maxstd, 1]) 
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp_cond = dot(C1,fp)
    
    # 映射对应点
    m = mean(tp[:2], axis=1)
    C2 = C1.copy() #must use same scaling for both point sets
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp_cond = dot(C2,tp)
    
    # 归一化后点的均值为0，所以平移量为0
    A = concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = linalg.svd(A.T)
    
    # 创建B.C矩阵
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]
    
    tmp2 = concatenate((dot(C,linalg.pinv(B)),zeros((2,1))), axis=1) 
    H = vstack((tmp2,[0,0,1]))
    
    # 反归一化
    H = dot(linalg.inv(C2),dot(H,C1))
    
    return H / H[2,2]

四 图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常在计算机图形学中，而且经常出现在计算机视觉算法总。扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。

from pylab import *
from numpy import *
from PIL import Image
from scipy import ndimage

im = array(Image.open('/Users/apple/Desktop/集大/集大.jpeg').convert('L'))
H = array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
im2 = ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]))
figure()
gray()
imshow(im2)
show()

 

五 图像中的图像

仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，是的他们能够和指定的区域或者标记物对齐。

def image_in_image(im1 ,im2 ,tp):

    m ,n = im1.shape[:2]
    fp = array([[0 ,m ,m ,0] ,[0 ,0 ,n ,n] ,[1 ,1 ,1 ,1]])

    H = homography.Haffine_from_points(tp ,fp)
    im1_t = ndimage.affine_transform(im1 ,H[:2 ,:2],
                                     (H[0 ,2] ,H[1 ,2]) ,im2.shape[:2])
    alpha = (im1_t > 0)

    return ( 1 -alpha ) *im2 + alpha *im1_t

创建alpha图像，该图像定义了每个像素从各个图像中获取的像素值成分多少，扭曲的图像实在扭曲区域边界之外以0来填充的图像，来创建一个alpah图像。

什么是alpha图像
将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。将扭曲的图像和第二幅图像融合，我们就创建了alpha图像
Alpha通道指的是特别的通道，意思是“非彩色”通道，主要是用来保存选区和编辑选区
在实验图片中，Alpha通道使用8位二进制数，就可以表示256级灰度，即256级的透明度。白色（值为255）的Alpha像素用以定义不透明的彩色像素，而黑色（值为0）的Alpha通道像素用以定义透明像素，介于黑白之间的灰度（值为30-255）的Alpha像素用以定义不同程度的半透明像素。
而以上代码的计算过程是：
两幅图像分别为A和B，由这两幅图像组合而成的图像称为C，则可用如下四元组表示图A和B，三元组表示图像C：
A：（Ra，Ga，Ba，AlphaA）
B：（Rb，Gb，Bb，AlphaB）
C：（Rc，Gc，Bc）
将图A和图B融合为图C
根据上述算法，则：
Rc=RaAlphaA+RbAlphaB
Gc=GaAlphaA+GbAlphaB
Bc=BaAlphaA+BbAlphaB

在此段代码中
将扭曲的图像和第二幅图像
融合，我们就创建了 alpha 图像 。该图像定义每个像素从各个图像中获取的像素值成分多少。这里我们基于以下事实，扭曲的图像是在扭曲区域边界之外以 0 来填充的图像，来创建一个二值的alpha 图像。
 

仿射变换实现

from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image

im1 = array(Image.open('/Users/apple/Desktop/集大/集大.jpeg').convert('L'))
im2 = array(Image.open('/Users/apple/Desktop/集大/集大2.jpeg').convert('L'))
tp = array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
im3 = image_in_image(im1,im2,tp)
figure()
gray()
imshow(im3)
axis('equal')
axis('off')
show()

分段仿射扭曲

给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，可以将图像和另一幅图像的对应标记点扭曲对应。
在Matplotlib 中有狄洛克三角剖分，但是使用jupyter notebook的时候报错了，改成了如下代码

from scipy.spatial import Delaunay

def triangulate_points(x,y):
 """ 二维点的Delaunay 三角剖分"""
 tri = Delaunay(np.c_[x,y]).simplices
 return tri


x,y = array(random.standard_normal((2,100)))
tri= triangulate_points(x,y)
figure()
for t in tri:
 t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
 plot(x[t_ext],y[t_ext],'r')
plot(x,y,'*')
axis('off')
show()

六.ransac全景拼接

from pylab import *
from numpy import *
from PIL import Image

# If you have PCV installed, these imports should work
from PCV.geometry import homography, warp
from PCV.localdescriptors import sift

"""
This is the panorama example from section 3.3.
"""

# set paths to data folder
featname = ['/Users/apple/Desktop/集大/'+str(i+1)+'.sift' for i in range(5)]
imname = ['/Users/apple/Desktop/集大/'+str(i+1)+'.jpeg' for i in range(5)]

# extract features and match
l = {}
d = {}
for i in range(5): 
    sift.process_image(imname[i],featname[i])
    l[i],d[i] = sift.read_features_from_file(featname[i])

matches = {}
for i in range(4):
    matches[i] = sift.match(d[i+1],d[i])

# visualize the matches (Figure 3-11 in the book)
'''
for i in range(4):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i+1]))
    figure()
    sift.plot_matches(im2,im1,l[i+1],l[i],matches[i],show_below=True)
'''

# function to convert the matches to hom. points
def convert_points(j):
    ndx = matches[j].nonzero()[0]
    fp = homography.make_homog(l[j+1][ndx,:2].T) 
    ndx2 = [int(matches[j][i]) for i in ndx]
    tp = homography.make_homog(l[j][ndx2,:2].T) 
    
    # switch x and y - TODO this should move elsewhere
    fp = vstack([fp[1],fp[0],fp[2]])
    tp = vstack([tp[1],tp[0],tp[2]])
    return fp,tp


# estimate the homographies
model = homography.RansacModel() 

fp,tp = convert_points(1)
H_12 = homography.H_from_ransac(fp,tp,model)[0] #im 1 to 2 

fp,tp = convert_points(0)
H_01 = homography.H_from_ransac(fp,tp,model)[0] #im 0 to 1 

tp,fp = convert_points(2) #NB: reverse order
H_32 = homography.H_from_ransac(fp,tp,model)[0] #im 3 to 2 

tp,fp = convert_points(3) #NB: reverse order
H_43 = homography.H_from_ransac(fp,tp,model)[0] #im 4 to 3    


# warp the images
delta = 500 # for padding and translation

im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12,im1,im2,delta,delta)

im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12,H_01),im1,im_12,delta,delta)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32,im1,im_02,delta,delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)


figure()
imshow(array(im_42, "uint8"))
axis('off')
show()

这里因为PCV因为版本问题需要修改print，安装百度的方法，从github下载并修改了它的代码后，setup.py好像因为路径问题找不到包，于是再次百度，将代码改为调用find_packages

提示安装成功 

 但是重启后运行依旧找不到包

同学说把pcv换成小写可以运行

 

 但是依旧提示失败

程序都是在tensorflow2的环境下执行的，mac有自带的python2版本，但是我安装和使用时都是指定的python3，理论上不应该存在装错位置找不到包才对

这里真的不知道怎么解决，如果找到解决办法的话会把这个实验继续完善