射影变换

IP2空间的图像校正与模拟贴图

Posted by Tower on January 4, 2018

======== 2018 新 年 分 割 线 ========

为了加快自己写周总结的速度,只能刻意练习了╮(╯▽╰)╭

射影变换

由照相机的小孔成像原理,当成像平面与所摄物体所在平面不平行时,产生的就是IP2中的射影变换。 二维射影平面内的一切射影变换本质都为3x3可逆阵,同时考虑到单应性,变换矩阵的自由度为8,只需4组对应点的x,y坐标值即可求。求解变换矩阵时可应用矩阵的SVD(奇异值分解)。图像校正的程序如下,需要注意的是:

  • Matlab和Python读取图像的原点都在左上角,在figure中用鼠标取得的坐标(x,y)中x表示到左边界的距离,y表示到上边界的距离。
  • 但读取的图像为矩阵,(i,j)中i表示所在的行数即可理解为到上边界的距离y,j表示所在的列数即到左边界的距离x。

Python3 使用现成函数校正图像

本程序使用的是python3.5中已经封装好的函数

cv2.getPerspectiveTransform(pts1,pts2) 得到从点集pts1到pts2的射影变换矩阵。

cv2.warpPerspective(img,M, (width, height)) 是从原图img通过矩阵M到新图的变换函数,其中width,height是原图的宽度和高度。

校正所需的结果已事先在程序中固定好了,需要的是用鼠标按照左上右上左下右下的顺序标定需要校正的图像中与校正结果中固定位置的4点相对应的4点。

  • 原图

 import cv2
 import numpy as np
 from matplotlib import pyplot as plt
 from pylab import *
 from PIL import Image

 def perspectiveTransform():
     img = Image.open("./project.jpg")
     im = np.array(img)
     imshow(im)
     print('Please click 4 points,order by lu,lr,dl,dr,(u=up,d=down,l=left,r=right)')
     # 点击顺序左上,右上,左下,右下,获取坐标保存在[x,y]列表中
     x = ginput(4)
     x = np.array(x)
     print('you clicked:', x)
     pts1 = np.float32(x)
     plt.close()
     img = cv2.imread("./project.jpg")
     height, width, channels = img.shape
     # 此处坐标(x,y)为x是到左边界的距离,而读出的图像矩阵x表示所在行数
     pts2 = np.float32([[20, 20], [940, 20], [20, 660], [940, 660]])
     pts3 = np.float32([[250, 180], [710, 180], [250, 500], [710, 500]])
     M = cv2.getPerspectiveTransform(pts1, pts2)
     dst = cv2.warpPerspective(img, M, (width, height))
     M3 = cv2.getPerspectiveTransform(pts1, pts3)
     dst3 = cv2.warpPerspective(img, M3, (width, height))
     print("completed,show result now!")
     plt.subplot(131), plt.imshow(img), plt.title("raw input")
     plt.subplot(132), plt.imshow(dst3), plt.title("projection regulate")
     plt.subplot(133), plt.imshow(dst), plt.title("enlarge output")
     plt.show()
     cv2.imwrite('./refine_result.jpg', dst3)
     cv2.imwrite('./refine_enlarge.jpg', dst)

 if __name__ == "__main__":
     perspectiveTransform()
     
  • 结果

Matlab射影贴图

Matlab中也有现成的可以直接做射影变换的函数:

imwarp(raw_img,H),其中raw_img为原始图像,H为射影矩阵

贴图的原图是整张图,直接按左上右上左下右下取原图的4角坐标即可。 运行时,在要贴的图上按对应顺序标定相应的4点。

将3张数学家的正面画像贴到魔方上

 warning off%不显示警告
 close all;clc;clear;
 img_path='./';%原图路径
 to_path='./';%贴图路径
 new_img='./results/';%新图路径
 img_list=dir([img_path,'*.jpg']);
 to_list=dir([img_path,'*.png']);
 if ~exist(new_img)
     mkdir(new_img);%本程序直接加在原图上了,没有用到new_img文件夹
 end
 num_img=length(img_list);
 %% 显示原图及bbox坐标
 for i=1:3
     tmp_i=i;
     img=imread([img_path,img_list(i).name]);
     to_img=imread('cubic.jpg');
     disp(['====',to_list(i).name,'====']);
     figure;imshow(to_img);
     disp('请以左上右上左下右下标定4个点');
     height=size(img,1);
     width=size(img,2);
     x=[1;width;1;width];
     y=[1;1;height;height];
     [xn,yn]=ginput(4);
    for tt=1:4
        xn(tt)=ceil(xn(tt));
        yn(tt)=ceil(yn(tt));
    end
     bound_u=min(yn);
     bound_d=max(yn);
     bound_l=min(xn);
     bound_r=max(xn);
     A=[x(1),y(1),1,0,0,0,-x(1)*xn(1),-xn(1)*y(1),-xn(1);
          0,0,0,x(1),y(1),1,-x(1)*yn(1),-y(1)*yn(1),-yn(1);
          x(2),y(2),1,0,0,0,-x(2)*xn(2),-xn(2)*y(2),-xn(2);
          0,0,0,x(2),y(2),1,-x(2)*yn(2),-y(2)*yn(2),-yn(2);
          x(3),y(3),1,0,0,0,-x(3)*xn(3),-xn(3)*y(3),-xn(3);
          0,0,0,x(3),y(3),1,-x(3)*yn(3),-y(3)*yn(3),-yn(3);
          x(4),y(4),1,0,0,0,-x(4)*xn(4),-xn(4)*y(4),-xn(4);
          0,0,0,x(4),y(4),1,-x(4)*yn(4),-y(4)*yn(4),-yn(4);];
     [~,~,V]=svd(A);
     h = V(:,9) ./ V(9,9);
     H= reshape(h,3,3);
     H=projective2d(H);
     new_im=imwarp(img,H);
     figure;imshow(new_im);
     to_img(bound_u:bound_d,bound_l:bound_r,:)=to_img(bound_u:bound_d,bound_l:bound_r,:)+new_im;
     figure;imshow(to_img);
     imwrite(to_img,'cubic.jpg');
     close all;
 end
  • 结果图 PS:像素有点渣╮(╯_╰)╭

射影变换作业小结

  • Python中cv2.imread()默认为bgr色彩模式所以显示有色差
  • 直接把原图(矩形)映射为平行四边形是不对的,平行四边形是仿射变换,相当于成像平面要平行于物面,但是此时物面是正朝外的,不应该得到这样倾斜的平行四边形,所以应该先映射为梯形,再截取。

  • 此处注意!

    无论求得的是原图到新图的变换矩阵还是新图到原图的变换矩阵, 最好遍历新图中的每一个点,求得其在原图中的对应点。新图大小若小于原图其实无所谓。 但是若新图大小大于原图,若是遍历原图,新图中必然有一些点没有对应到,就出现下图中左边的情况。 而遍历新图则可以保证新图中射影区域内的每一个点都有值,可能多个新图中的点对应原图中的一个点。  > 不过此时,找到原图中的对应点是要加一个是否超出原图边界的判断。当然,过程中算得的坐标都要取整。

  • 截取时可以将与原图等大的黑图变换为最后所要的平行四边形,黑色是0保留加上的颜色,白色255覆盖任意颜色。

最后的结果如下图:

具体程序分为3个文件:2个函数,1个主函数

  • SVD_H.m

    function [ H ]=SVD_H(x,y,xn,yn)
    %输入(y,x)为原图左上右上左下右下4点坐标
    %输入(yn.xn)为新图对应点坐标
    %此处在figure中取得的坐标要转置才为在矩阵中的行列坐标
    %从原图到新图的射影平面变换的3*3可逆阵
        A=[x(1),y(1),1,0,0,0,-x(1)*xn(1),-xn(1)*y(1),-xn(1);
             0,0,0,x(1),y(1),1,-x(1)*yn(1),-y(1)*yn(1),-yn(1);
             x(2),y(2),1,0,0,0,-x(2)*xn(2),-xn(2)*y(2),-xn(2);
             0,0,0,x(2),y(2),1,-x(2)*yn(2),-y(2)*yn(2),-yn(2);
             x(3),y(3),1,0,0,0,-x(3)*xn(3),-xn(3)*y(3),-xn(3);
             0,0,0,x(3),y(3),1,-x(3)*yn(3),-y(3)*yn(3),-yn(3);
             x(4),y(4),1,0,0,0,-x(4)*xn(4),-xn(4)*y(4),-xn(4);
             0,0,0,x(4),y(4),1,-x(4)*yn(4),-y(4)*yn(4),-yn(4);];
         [~,~,V]=svd(A);
         %reshape是将向量按照顺序一列一列的放进去,此处我们需要一行一行的放,所以要转置
         v=reshape(V(:,9)/V(9,9),3,3);
         H=v';
    end
    
  • Project_H.m

    function [ new_img ]=Project_H(img,H,new_img)
    %输入img为原图,H为变换矩阵
    %输入输出new_img为变换的新图
    %注意:H为原图到新图的,先求逆,遍历新图的点在原图中找到对应点,
    %当新图大于原图,很可能新图中多个点对原图中的一个点
        H2=inv(H);
        for m=1:size(new_img,1)
            for n=1:size(new_img,2)
                   tmp=H2*[m;n;1];
                   tmp=round(tmp/tmp(3));
                   if(tmp(1)<=size(img,1) &&tmp(1)>=1&&tmp(2)>=1&& tmp(2)<=size(img,2))
                       new_img(m,n,:)=img(tmp(1),tmp(2),:);
                   end
            end
        end
        new_img=uint8(new_img);
    end
    
  • project_final.m

    close all;clc;clear;
    img_path='.\FileRecv\作业\';%原图路径
    to_path='.\FileRecv\作业\';%原图路径
    new_img='.\FileRecv\作业\results';%新图路径
    img_list=dir([img_path,'*.jpg']);
    to_list=dir([img_path,'*.png']);
    if ~exist(new_img,'dir')
        mkdir(new_img);
    end
    num_img=length(img_list);
    for i=1:3
        img=imread([img_path,img_list(i).name]);
        black=img;
        black(:,:,:)=0;
        figure;imshow(img);
        to_img=imread('white_bg.jpg');
        disp(['====',to_list(i).name,'====']);
        figure;imshow(to_img);
        disp('请以左上右上左下右下标定4个点');
        height=size(img,1);
        width=size(img,2);
        y=[1;width;1;width];
        x=[1;1;height;height];
        %注意此处要加转置,读入的横纵坐标与矩阵中所在行列是相反的
        [yn,xn]=ginput(4);
       for tt=1:4
           xn(tt)=ceil(xn(tt));
           yn(tt)=ceil(yn(tt));
       end
       %把4个点都固定在左右两边上,至少为梯形
       yn=[1;size(to_img,2);1;size(to_img,2)];
       %原图x,y  新图xn,yn  切图xp,yn表示的都是在矩阵中的行列索引
       yp=yn;xp=xn;
       xp(4)=xp(2)-xp(1)+xp(3);
       %得到原图到新图的变换矩阵H
       H=SVD_H(x,y,xn,yn);
       raw_to=to_img;
       white_bg=to_img;
       close all;
       %% 法一:遍历原图的每个点
        for m=1:size(img,1)
            for n=1:size(img,2)
                   tmp=H*[m;n;1];
                   tmp=tmp/tmp(3);
                   tmp=min(max(round(tmp),[1;1;1]),[size(to_img,1);size(to_img,2);1]);
                   to_img(tmp(1),tmp(2),:)=img(m,n,:)
            end
        end
        to_img=uint8(to_img);
        figure;subplot(1,2,1);
        imshow(to_img);title('遍历原图中的每个点');
        %% 法二:遍历新图的每个点,找在原图中的位置
        raw_to=Project_H(img,H,raw_to);
        subplot(1,2,2);
        imshow(raw_to);title('遍历新图中的每个点');
        %% 利用黑色图像射影变换的白底覆盖梯形超出平行四边形的部分
        Hb=SVD_H(x,y,xp,yp);
        white_bg=Project_H(black,Hb,white_bg);
        white_bg=uint8(white_bg);
        figure;imshow(white_bg);title('从梯形汇中截取平行四边形的模板')
        %% 截取平行四边形
        raw_to=raw_to+white_bg;
        final_img=raw_to(1:xp(4),:,:);
        figure;
        subplot(1,2,1);
        imshow(img);title('原图');
        subplot(1,2,2);
        imshow(final_img);title('结果图');
        imwrite(final_img, [new_img,'\', to_list(1).name(1:end-7),'.jpg']);
    end