相机标定实践1-获取内参

相机标定，获取内参
至少3幅影响（求B，6参数，每幅影像对应求2个），每幅影像至少4个特征点（求H，8参数，每个特征点对应求2个）
参考资料：https://blog.csdn.net/zkl99999/article/details/48372203?depth_1-utm_source=distribute.pc_relevant.none-task&utm_source=distribute.pc_relevant.none-task

https://blog.csdn.net/a6333230/article/details/83478064

理论1
相机标定分：内参标定、外参标定、c矩阵
内参标定：采用棋盘图或者点阵图，将相机固定（焦距、光圈等内部元素）。用相机拍摄多组不同角度的棋盘图，这点和普通的matlab的棋盘图标定类似。matlab利用棋盘图标定，它也给出了标定的外参，但是给定的并非我们所求，而是将世界坐标系定在了棋盘图上，我们需要的世界坐标系是固定不变的或者至少是与相机的相对位置不会发生改变。
内参标定输入为：像素坐标u、v和世界坐标x、y、z，该z没有定义具体的数值，因为在此处的求解不会应用到，并且我们此步也无法求得。那么这些参数怎么来的呢？u、v为图像中棋盘格的角点坐标，与其对应的x、y表示该角点在整个棋盘图下的坐标，该坐标系定在棋盘图上，因此我们拍摄5次图像，这些图像相应的坐标应该相同。为了计算时数据的变换较为方便，这里z也给出，值给0或1都行，不影响计算。
内参标定输出为：内部参数+相机畸变参数。内部参数：fu、fv、u、v；相机畸变参数：k1，k2
内参求解步骤：
1.求单应矩阵




这里的h33设置为1……课本中给出的理由是它可以设置为任意数，二郎也和同事讨论过，没有定论，有一种理解有道理：计算中h33能都提到右边，那么在根据公式进行求解时，它变成几，都是对所有公式进行操作的，因此公式之间的关系不会发生改变，所求出来代表关系的h参数也不会受到影响，只是按h33进行了缩放。
比如：两边先同时除以h33（由于没有找到这里的h33=1的解释，因此二郎做了这种假设）
每个特征点，可以求得两个方程。这里有8个未知数，因此可以用4个特征点完成求解。
这里用到了棋盘格，在一张棋盘格图中，我们得到了很多（u，v）与（x，y）的对应，因此可以解出我们的单应矩阵。
这里每个棋盘格图对应了一个单应矩阵，这些单应矩阵主要的不同在于（R，t），每个棋盘图都有自己的世界坐标系，原点在棋盘图的一个角点上，这个matlab内部程序已经指定，一般情况下为左上角。此处单应矩阵的求解用到的最小二乘法。下面的推导和求解，均是建立在单应矩阵已知的情况下进行的。

实践1
1.拍摄标定板照片

function CreatePlate()
J = (checkerboard(300,3,4)>0.5);    %生成黑白棋盘图像
figure, imshow(J)   %显示黑白棋盘图像
imwrite(J,'plate.jpg'); %保存黑白棋盘图像
end 

然后手机拍照存储作为相机拍摄的标定板照片
设照片名称为“MyPlate1.jpg”

2.标定板角点检测
harris算法

function Myharris()
X = imread('MyPlate1.jpg');         % 读取图像
% figure,imshow(X);                       %在窗口1中打开原始图像
Info = imfinfo('MyPlate1.jpg');     %获取图像相关信息
if (Info.BitDepth > 8)
    f = rgb2gray(X);                    %若图像的位深大于8，则用rgb2gray将彩色图像转换为灰度图像 
end
%计算图像亮度f(x,y)在点(x,y)处的梯度-----------------------------------------------
ori_im = double(f) / 255;               %unit8转化为64位双精度double64，取值在0~1
fx = [-2 -1 0 1 2];                     % x方向梯度算子(用于Harris角点提取算法)
fy = [-2; -1; 0; 1; 2];                 % y方向梯度算子(用于Harris角点提取算法)
gx = filter2(fx, ori_im);               % x方向滤波
gy = filter2(fy, ori_im);               % y方向滤波
%构造自相关矩阵---------------------------------------------------------------
gx2 = gx .^ 2;
gy2 = gy .^ 2;
gxy = gx .* gy;
% clear gx;
% clear gy;
h= fspecial('gaussian', [7 7], 0.6);    % 产生5*5的高斯窗函数，sigma取值在0.3~0.9，此处取0.6
gx2 = filter2(h,gx2);
gy2 = filter2(h,gy2);
gxy = filter2(h,gxy);
%提取特征点---------------------------------------------------------------
[height,width]=size(ori_im);            %读取图像的高和宽
result = zeros(height, width);          % 纪录角点位置，角点处值为1
I = zeros(height, width);
Imax = 0;                               % 图像中最大的R值
k = 0.2;                               %k为常系数，经验取值范围为0.04~0.06
for i = 1 : height
    for j = 1 : width
        M = [gx2(i, j) gxy(i, j); gxy(i, j) gy2(i, j)];         % auto correlation matrix
        I(i,j) = det(M) - k * (trace(M)) ^ 2;                   % 计算R
        if I(i,j) > Imax
            Imax = I(i, j);
        end;
    end;
end;
T = 0.06 * Imax;%固定阈值，当R(i, j) > T时，则被判定为候选角点
%在计算完各点的值后，进行局部非极大值抑制-------------------------------------
cnt = 0;
for i = 2 : height-1
    for j = 2 : width-1
        % 进行非极大抑制，窗口大小3*3
        if (