cp_registration.m的学习记录

七仔的鸽子

已于 2023-04-14 10:37:48 修改

阅读量220

点赞数

分类专栏：图像匹配算法文章标签：学习算法

于 2023-04-11 11:14:17 首次发布

本文链接：https://blog.csdn.net/weixin_47000169/article/details/130034242

版权

图像匹配算法专栏收录该内容

10 篇文章 1 订阅

订阅专栏

第一部分：

这一部分不用做过多解释，主要是对cp_cornerDetection.m的学习记录这一篇文章中的3个输出进行处理。cor1, orientation1, I1edge表示的是cout, angle, edge_map，cout表示的就是每一个元组中符合条件的极大值，起始点，终止点。angle表示的就是符合论文中所给条件所构成的νjfm的极角。edge_map表示的是元组中处理后的轮廓点。处理表示的是判断轮廓中的点数是否大于设定值，才会被认为为是一个轮廓。

之后cor12表示的是将红外图像与可见光图像经过cp_cornerDetection处理后的满足要求的极大值点保存（为什么要进行转化呢，因为图片的坐标原点再左上角，而现在使用的要将其逆时针旋转90°）

%% Corner detection and orientation computed
tic,
%                                                   I  C ang sigm H   L  endp gap
[cor1, orientation1, I1edge] = cp_cornerDetection(I1,[],[], [], [], 0,  1,  1,    Lc,iteration==0); % Length of curvature
[cor2, orientation2, I2edge] = cp_cornerDetection(I2,[],[], [], [], [], 1,  1,    Lc,iteration==0);
CFO_time = toc
cor12 = [cor1(:,2) cor1(:,1); 0 0; cor2(:,2) cor2(:,1)];
I1 = double(I1);
I2 = double(I2);
I2ori = double(I2ori);
% fprintf('\tNumber of corner1: %d || Number of corner2: %d \n',length(cor1),length(cor2));

第二部分：

画出x,y方向的侯选角

在 MATLAB 中，quiver 函数用于在 2D 平面上画出向量图。通过指定起点、长度和方向，quiver 函数会在图形窗口中绘制箭头代表向量，并可进行颜色和风格的自定义。

quiver 函数的基本语法如下所示：
quiver(x, y, u, v)
quiver(X,Y,U,V) 在由 X 和 Y 指定的笛卡尔坐标上绘制具有定向分量 U 和 V 的箭头。例如，第一个箭头源于点 X(1) 和 Y(1)，按 U(1) 水平延伸，按 V(1) 垂直延伸。默认情况下，quiver 函数缩放箭头长度，使其不重叠。

%% show orientation of each corner
if iteration == 0 && showflag
    xu1 = 4 * cos(orientation1);
    yv1 = 4 * sin(orientation1);
    xu2 = 4 * cos(orientation2);
    yv2 = 4 * sin(orientation2);
    image1 = I1/255;
    image2 = I2/255;
%     image1 = I1edge;
%     image2 = I2edge;
    figure,
    subplot(121),imshow(image1);hold on, plot(cor1(:,2),cor1(:,1),'y+'),
    title('红外图像角点')
    vect = quiver(cor1(:,2),cor1(:,1),xu1,yv1);
    vect.Color = 'r';
    vect.LineWidth = 0.8;
    subplot(122),imshow(image2);hold on, plot(cor2(:,2),cor2(:,1),'y+'),
    title('可见光图像角点')
    vect = quiver(cor2(:,2),cor2(:,1),xu2,yv2);
    vect.Color = 'r';
    vect.LineWidth = 0.8;
end

第三部分：

使用下面这段代码前的理论支持：

%% extract descriptor of Image 1
% 这里后期需要根据自己的图片进行画框
tic;
scale12 = 36; % square width is 6(pixels) * 6(squares) = 36 pixels totally
cols1=cor1(:,1); % swap row and col 交换x，y
rows1=cor1(:,2);  % row is x axis | col is y axis
s1 = scale12 * ones(size(cols1));
if iteration > 0% no need to caculate orientation
    o1 = zeros(size(cols1));
else
    o1 = orientation1;
end
key1 = [rows1';cols1';s1';o1'];
des1 = cp_descriptor(I1,key1);
%% extract descriptors of Image 2 under multi scales
cols2 =cor2(:,1); % swap row and col
rows2  =cor2(:,2);  % row is x axis | col is y axis
if iteration > 0  % no need to caculate orientation
    o2 = zeros(size(cols2));
else
    o2 = orientation2;
end
s2 = scale12 * ones(size(cols2));
key2 = [rows2';cols2';s2';o2'];
% zoom transformation parameters
zoomstepup = 1;
zoomstepdown = 0.5;
des2 = zeros([128 length(cor2) zoomascend+zoomdescend+1]);
des2(:,:,1) = cp_descriptor(I2,key2);
level = 1;
scale = size(I2ori,1) / size(I2,1);
if iteration == 0
    for i=1:zoomascend
        level = level + 1;
        key2zoom = key2;
        I2zoom = imresize(I2ori,(1+zoomstepup*i) / scale);
        key2zoom(1:2,:) = floor( (1+zoomstepup*i) * key2zoom(1:2,:) );
        des2(:,:,level) = cp_descriptor(I2zoom,key2zoom);
    end
    for i=1:zoomdescend
        level = level + 1;
        key2zoom = key2;
        I2zoom = imresize(I2ori,(1-zoomstepdown*i)/scale);
        key2zoom(1:2,:) = floor( (1-zoomstepdown*i) * key2zoom(1:2,:) );
        des2(:,:,level) = cp_descriptor(I2zoom,key2zoom);
    end
end
descriptortoc = toc

a部分：

对红外图像进行进行SIFT算法处理，提取特征描述符，其中呢，极值与角度都在cp_cornerDetection中计算得到。直接传入cp_descriptor中。

在每个像素处，梯度都有一个大小和一个方向

%% extract descriptor of Image 1
% 这里后期需要根据自己的图片进行画框
tic;
scale12 = 36; % square width is 6(pixels) * 6(squares) = 36 pixels totally
cols1=cor1(:,1); % swap row and col 交换x，y
rows1=cor1(:,2);  % row is x axis | col is y axis
s1 = scale12 * ones(size(cols1));
if iteration > 0% no need to caculate orientation
    o1 = zeros(size(cols1));
else
    o1 = orientation1;
end
key1 = [rows1';cols1';s1';o1'];
des1 = cp_descriptor(I1,key1);

b部分：

对可见光图像进行进行SIFT算法处理，提取特征描述符，其中呢，极值与角度都在cp_cornerDetection中计算得到。分为两部分传入cp_descriptor中，存入des2中，des2分为两部分，一部分存放原有的特征描述符，第二部分存放将其放大一倍后x,y之后输入到cp_descriptor中得到的描述符。下面两个for循环只在iteration == 0才执行。具体是提取缩放后的特征描述符。并存入到des2中。

%% extract descriptors of Image 2 under multi scales
cols2 =cor2(:,1); % swap row and col
rows2  =cor2(:,2);  % row is x axis | col is y axis
if iteration > 0  % no need to caculate orientation
    o2 = zeros(size(cols2));
else
    o2 = orientation2;
end
s2 = scale12 * ones(size(cols2));
key2 = [rows2';cols2';s2';o2'];
% zoom transformation parameters
zoomstepup = 1;
zoomstepdown = 0.5;
des2 = zeros([128 length(cor2) zoomascend+zoomdescend+1]);
des2(:,:,1) = cp_descriptor(I2,key2);
level = 1;
scale = size(I2ori,1) / size(I2,1);
if iteration == 0
    for i=1:zoomascend
        level = level + 1;
        key2zoom = key2;
        I2zoom = imresize(I2ori,(1+zoomstepup*i) / scale);
        key2zoom(1:2,:) = floor( (1+zoomstepup*i) * key2zoom(1:2,:) );
        des2(:,:,level) = cp_descriptor(I2zoom,key2zoom);
    end
    for i=1:zoomdescend
        level = level + 1;
        key2zoom = key2;
        I2zoom = imresize(I2ori,(1-zoomstepdown*i)/scale);
        key2zoom(1:2,:) = floor( (1-zoomstepdown*i) * key2zoom(1:2,:) );
        des2(:,:,level) = cp_descriptor(I2zoom,key2zoom);
    end
end
descriptortoc = toc

第四部分：

双边检测，先了解cp_match函数cp_match函数记录_七仔的鸽子的博客

cp_match函数就是对特征描述符进行双边匹配，选出符合匹配条件的key与value

%% *----Matchings the keys from two images---*
%% coarsely match with BBF method
[matchIndex1, matchIndex2, zoom] = cp_match(des1',des2, 0.97);
zoomscale = (zoom == 1)*1 + (1<zoom & zoom <= (1+zoomascend))*(1+(zoom-1)*zoomstepup) +...
    (zoom > (1+zoomascend))*(1-(zoom-1-zoomascend)*zoomstepdown)
regis_points111 = [cor1(matchIndex1,2) cor1(matchIndex1,1) o1(matchIndex1)]; % (u,v,orientation)
regis_points222 = [cor2(matchIndex2,2) cor2(matchIndex2,1) o2(matchIndex2)];
fprintf('\tNumber of features1: %d || Number of features2: %d \n',length(regis_points111),length(regis_points222));
if size(regis_points111,1) < 2
    error('ERROR: No suffcient points( < 2 ) ! Faild to registrate!');
end
if showflag
    cp_showMatch(I1,imresize(I2ori,zoomscale/scale),regis_points111(:,1:2),...
        zoomscale*regis_points222(:,1:2),[],'Putative matches');
end
% cp_showMatch(I1,I2,regis_points111,regis_points222,[],'BBF matching result');

其中这一段表示缩放因子：

(zoom == 1)*1：这个条件表示当zoom等于1时，缩放因子就是1，即不进行缩放。通过这个条件可以避免进行不必要的缩放操作。
(1<zoom & zoom <= (1+zoomascend))*(1+(zoom-1)*zoomstepup)：这个条件表示当zoom大于1小于等于(1+zoomascend)时，缩放因子是(1+(zoom-1)*zoomstepup)。这段代码可以实现缩放增加操作，使得当缩放比例在这个范围内时，缩放因子随着缩放比例增加而增加，从而达到更好的视觉效果。
(zoom > (1+zoomascend))*(1-(zoom-1-zoomascend)*zoomstepdown)：这个条件表示当zoom大于(1+zoomascend)时，缩放因子是(1-(zoom-1-zoomascend)*zoomstepdown)。这段代码可以实现缩放减小操作，使得当缩放比例超过这个范围时，缩放因子随着缩放比例减小而减小，从而达到更好的视觉效果。

zoomscale = (zoom == 1)*1 + (1<zoom & zoom <= (1+zoomascend))*(1+(zoom-1)*zoomstepup) +...
    (zoom > (1+zoomascend))*(1-(zoom-1-zoomascend)*zoomstepdown)

下面这段代码是将cp_cornerDetection.m的学习记录处理好的极大值点的x，y与对应的角度存入数组，regis_points111的大小为 [ N x 3 ]

regis_points111 = [cor1(matchIndex1,2) cor1(matchIndex1,1) o1(matchIndex1)]; % (u,v,orientation)
regis_points222 = [cor2(matchIndex2,2) cor2(matchIndex2,1) o2(matchIndex2)];

之后进入cp_showMatch将匹配点绘制出来。

第五部分： Scale invariant tilt angle consistency matching（比例不变倾角一致性匹配）

%% Scale invariant tilt angle consistency matching
tic;
if iteration == 0
    delta0 = mod(round(180 / pi*(regis_points222(:,3) - regis_points111(:,3))),360);
    % delta0 = delta0 - mod(delta0,2); % Take Δorientation = 2°as a various level
    dd = 5;
    d_delta = 0:dd:360;
    n_delta = histc(delta0,d_delta);
    [~,nindex] = sort(n_delta ,'descend' );
    n0 = nindex(1);
    nmat = [n0^2, n0, 1;(n0-1)^2, n0-1, 1;(n0+1)^2, n0+1, 1] \ ...
        [n_delta(n0) n_delta(n0-1+360/dd*(n0==1)) n_delta(n0+1-(n0==360/dd))]';
    Modetheta_discrete = -nmat(2)/ 2 /nmat(1); % -b/2a
    Modetheta = Modetheta_discrete * dd; % 抛物线插值
    % x_dd = (Modetheta_discrete*dd - 1.2*dd):0.01*dd:(Modetheta_discrete*dd +1.2*dd);
    % xx = x_dd / dd;
    % yy = xx.^2 * nmat(1)+xx * nmat(2) + nmat(3);
    % figure, bar(d_delta,n_delta,0.7,'Facecolor',[1 0.8 0.6]);
    % hold on, plot(x_dd-dd,yy,'LineWidth',2),
    %          plot(dd*[n0 n0-1+360/dd*(n0==1) n0+1-(n0==360/dd)]-dd,[n_delta(n0) n_delta(n0-1+360/dd*(n0==1)) n_delta(n0+1-(n0==360/dd))],'r*');
    % hold off;
else
    Modetheta = 0;
end
% image to the same direction then those correct matches line will be parallel

trans222 = [regis_points222(:,1)-size(I2,2)/2  regis_points222(:,2)-size(I2,1)/2] * ...
    [cos(Modetheta*pi/180) -sin(Modetheta*pi/180); sin(Modetheta*pi/180) cos(Modetheta*pi/180)];
trans222 = [trans222(:,1)+size(I2,2)/2 trans222(:,2)+size(I2,1)/2];
phi_uv = cp_atan(zoomscale*trans222(:,2)-regis_points111(:,2),...
    zoomscale*trans222(:,1)+size(I1,2)-regis_points111(:,1));
% % show result of bilateral matching
if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(:,1:2),...
        zoomscale*trans222(:,1:2),[],['After ' num2str(iteration) ' resized and rotated']);
end
dd = 5; d_phi = -90:dd:90;
n_phi = histc(phi_uv,d_phi);
[~,nindex] = sort(n_phi ,'descend' ); n0 = nindex(1);
interval_index = find(n_phi < (n_phi(n0) * 0.2));
interval_index = interval_index - n0;
left_phi = find(interval_index < 0); right_phi = find(interval_index > 0);
maxtheta1 = -dd * interval_index(left_phi(end)); maxtheta2 = dd * interval_index(right_phi(1));
nmat = [n0^2, n0, 1;(n0-1+180/dd*(n0==1))^2, n0-1+180/dd*(n0==1), 1;(n0+1-(n0==180/dd))^2, n0+1-(n0==180/dd), 1]\...
    [n_phi(n0) n_phi(n0-1+180/dd*(n0==1)) n_phi(n0+1-(n0==180/dd))]';
ModePhi_discrete = -nmat(2)/ 2 /nmat(1) - 90/dd; % -b/2a
ModePhi = ModePhi_discrete * dd;
delta1 = ModePhi - maxtheta1; delta2 = ModePhi + maxtheta2;
[valid0,~] = find( (phi_uv>= delta1) & (phi_uv <=delta2) );

Dist = sqrt((zoomscale*trans222(valid0,2)-regis_points111(valid0,2)).^2+(zoomscale*trans222(valid0,1)+size(I1,2)-regis_points111(valid0,1)).^2);
meandist = mean(Dist);
[valid1,~] = find( (Dist>= 0.5*meandist) & (Dist <=1.5*zoomscale*meandist) );

regis_points11 = regis_points111(valid0(valid1),:);
regis_points22 = regis_points222(valid0(valid1),:);

% % show result of Tile angle consistency matching
if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(valid0(valid1),1:2),...
        zoomscale*trans222(valid0(valid1),1:2),[],['After ' num2str(iteration) ' Tile angle consistency']);
end
SITAC_TIME = toc
correctindex = cp_mismatchRemoval(regis_points11,regis_points22,I1,I2,maxErr);
RANSAC_TIME = toc
regis_points1 = regis_points11(correctindex,1:2);
regis_points2 = regis_points22(correctindex,1:2);
Runtime = CFO_time + descriptortoc + RANSAC_TIME;

if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(valid0(valid1(correctindex)),1:2),...
        zoomscale*trans222(valid0(valid1(correctindex)),1:2),[],['After ' num2str(iteration) ' Tile angle consistency']);
    cp_showMatch(I1,I2,regis_points1,regis_points2,[],'Fine matching result');% show matches between two images
end
fprintf('【1】Completed image registration\n ');

a部分：

用匹配后的对应点的角度求差得到 $\Delta \varphi$ （弧度转化为角度，之后再对360度进行求模 $mod(x,y) = x - floor(x./y).*y$ ，为了将 $-\varphi$ 转换到(0°,359°）区间上）。

之后论文中使用6°作为间隔，而此处使用5°作为间隔来绘制直方图。

histc(delta0,d_delta);

举例：

a = [1,1,2,3,4,5,6]；
b = 0:1.5:6；
c = histc(a,b)
figure
bar(b,c,'histc')

对n_delta进行降序排列，选出具有最多点的区间的索引。

开始抛物线插值，计算系数矩阵

image to the same direction then those correct matches line will be parallel（图像的方向相同，那么那些正确的匹配线将是平行的）

对于trans222的解释，

regis_points222(:,1) 表示 regis_points222 矩阵中所有点的横坐标。size(I2,2) 表示图像 I2 的宽度，也就是图像在横轴上的尺寸。

将 regis_points222(:,1) 减去 size(I2,2)/2 的目的是将每个点的横坐标向左平移图像宽度的一半，这样可以将图像中心移动到坐标系的原点。

同理，代码中的 regis_points222(:,2)-size(I2,1)/2 将每个点的纵坐标向上平移图像高度的一半，也是为了将图像中心移动到坐标系的原点。

之后又将trans222返回到图像处理坐标轴中。

phi_uv就是公式16中的 $\varphi$ 。

展示旋转之后的图像。

%% Scale invariant tilt angle consistency matching
tic;
if iteration == 0
    delta0 = mod(round(180 / pi*(regis_points222(:,3) - regis_points111(:,3))),360);
    % delta0 = delta0 - mod(delta0,2); % Take Δorientation = 2°as a various level
    dd = 5;
    d_delta = 0:dd:360;
    n_delta = histc(delta0,d_delta);
    [~,nindex] = sort(n_delta ,'descend' );
    n0 = nindex(1);
    nmat = [n0^2, n0, 1;(n0-1)^2, n0-1, 1;(n0+1)^2, n0+1, 1] \ ...
        [n_delta(n0) n_delta(n0-1+360/dd*(n0==1)) n_delta(n0+1-(n0==360/dd))]';
    Modetheta_discrete = -nmat(2)/ 2 /nmat(1); % -b/2a
    Modetheta = Modetheta_discrete * dd; % 抛物线插值
    % x_dd = (Modetheta_discrete*dd - 1.2*dd):0.01*dd:(Modetheta_discrete*dd +1.2*dd);
    % xx = x_dd / dd;
    % yy = xx.^2 * nmat(1)+xx * nmat(2) + nmat(3);
    % figure, bar(d_delta,n_delta,0.7,'Facecolor',[1 0.8 0.6]);
    % hold on, plot(x_dd-dd,yy,'LineWidth',2),
    %          plot(dd*[n0 n0-1+360/dd*(n0==1) n0+1-(n0==360/dd)]-dd,[n_delta(n0) n_delta(n0-1+360/dd*(n0==1)) n_delta(n0+1-(n0==360/dd))],'r*');
    % hold off;
else
    Modetheta = 0;
end
% image to the same direction then those correct matches line will be parallel

trans222 = [regis_points222(:,1)-size(I2,2)/2  regis_points222(:,2)-size(I2,1)/2] * ...
    [cos(Modetheta*pi/180) -sin(Modetheta*pi/180); sin(Modetheta*pi/180) cos(Modetheta*pi/180)];
trans222 = [trans222(:,1)+size(I2,2)/2 trans222(:,2)+size(I2,1)/2];
phi_uv = cp_atan(zoomscale*trans222(:,2)-regis_points111(:,2),...
    zoomscale*trans222(:,1)+size(I1,2)-regis_points111(:,1));
% % show result of bilateral matching
if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(:,1:2),...
        zoomscale*trans222(:,1:2),[],['After ' num2str(iteration) ' resized and rotated']);
end

b部分：对角度与长度进行筛选

phi_uv就是计算两幅图之间匹配点连线的倾斜角度。

n_phi是对phi_uv进行分类，放入直方图中。

nmat如a部分的nmat计算思路一样，找到出现次数最多的点与左右两边的点索引，以及点出现的次数，进而来进行抛物线插值。

dd = 5; 
d_phi = -90:dd:90;
n_phi = histc(phi_uv,d_phi);
[~,nindex] = sort(n_phi ,'descend' ); 
n0 = nindex(1);
interval_index = find(n_phi < (n_phi(n0) * 0.2));
interval_index = interval_index - n0;
left_phi = find(interval_index < 0); 
right_phi = find(interval_index > 0);
maxtheta1 = -dd * interval_index(left_phi(end)); 
maxtheta2 = dd * interval_index(right_phi(1));
nmat = [n0^2, n0, 1;(n0-1+180/dd*(n0==1))^2, n0-1+180/dd*(n0==1), 1;(n0+1-(n0==180/dd))^2, n0+1-(n0==180/dd), 1]\...
    [n_phi(n0) n_phi(n0-1+180/dd*(n0==1)) n_phi(n0+1-(n0==180/dd))]';
ModePhi_discrete = -nmat(2)/ 2 /nmat(1) - 90/dd; % -b/2a
ModePhi = ModePhi_discrete * dd;
delta1 = ModePhi - maxtheta1; 
delta2 = ModePhi + maxtheta2;
[valid0,~] = find( (phi_uv>= delta1) & (phi_uv <=delta2) );

Dist = sqrt((zoomscale*trans222(valid0,2)-regis_points111(valid0,2)).^2+...
    (zoomscale*trans222(valid0,1)+size(I1,2)-regis_points111(valid0,1)).^2);
meandist = mean(Dist);
[valid1,~] = find( (Dist>= 0.5*meandist) & (Dist <=1.5*zoomscale*meandist) );

regis_points11 = regis_points111(valid0(valid1),:);
regis_points22 = regis_points222(valid0(valid1),:);

% % show result of Tile angle consistency matching
if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(valid0(valid1),1:2),...
        zoomscale*trans222(valid0(valid1),1:2),[],['After ' num2str(iteration) ' Tile angle consistency']);
end

c部分：随机样本一致性（RANSAC）算法

在cp_mismatchRemoval函数记录（使用修改后的RANSCA算法进一步筛选点之后，在b部分的基础上画图，之后再没有旋转的图上将最终选取的点绘制出来。

SITAC_TIME = toc
correctindex = cp_mismatchRemoval(regis_points11,regis_points22,I1,I2,maxErr);
RANSAC_TIME = toc
regis_points1 = regis_points11(correctindex,1:2);
regis_points2 = regis_points22(correctindex,1:2);
Runtime = CFO_time + descriptortoc + RANSAC_TIME;

if showflag
    cp_showMatch(I1,imresize(imrotate(I2,Modetheta,'crop'),zoomscale),regis_points111(valid0(valid1(correctindex)),1:2),...
        zoomscale*trans222(valid0(valid1(correctindex)),1:2),[],['After ' num2str(iteration) ' Tile angle consistency']);
    cp_showMatch(I1,I2,regis_points1,regis_points2,[],'Fine matching result');% show matches between two images
end
fprintf('【1】Completed image registration\n ');

这个函数大体解释完了，等后期继续补充细节。