cp_descriptor学习记录

输入参数(img,keypoints)

输出参数（descriptors）

参数的初始化

NBP = 4;
NBO = 8;
key_num = size(keypoints, 2);% key_num = 4
descriptors = zeros([NBP*NBP*NBO key_num]);

计算Image1的梯度幅值

        一维二维（非标准的）sobel算子：一维是x,y方向，二维是对角线方向） 

        更像是3x3的罗伯特交叉梯度算子

%======= compute the image gradient vertors =====%
[M, N] = size(img); % M is the height of image, N is the width of image
du_filter = [-1 0 1];
dv_filter = du_filter';
duv_filter = diag([-1 0 1]);
dvu_filter = flip(duv_filter,2);
gradient_u = imfilter(img, du_filter);
gradient_v = imfilter(img, dv_filter);
gradient_uv = imfilter(img,duv_filter);
gradient_vu = imfilter(img,dvu_filter);
gradient_x = 1.414*gradient_u + gradient_uv - gradient_vu;
gradient_y = 1.414*gradient_v + gradient_uv + gradient_vu;
% dx_filter = [-1 0 1];
% dy_filter = dx_filter';
% gradient_x = imfilter(img, dx_filter);
% gradient_y = imfilter(img, dy_filter);
magnitudes = sqrt( gradient_x.^2 + gradient_y.^2);

 计算梯度角度与高斯核

        x表示cout的行

        y表示cout的列

        s表示像素个数

        theta表示角度

        sintheta表示sin(theta) 分量      

        costheta表示cos(theta)  分量  

        高斯差分金字塔

        [xx,yy]表示选取高斯核的大小

        wincoef表示高斯核

%====== compute angle of gradient =======%
angles = zeros(size(img));
for i = 1:M
    for j = 1:N
        if gradient_x(i,j) <= 0
            angles(i,j) = atan(gradient_y(i,j)/gradient_x(i,j)) + pi;
        elseif gradient_x(i,j) > 0 && gradient_y(i,j) >= 0
            angles(i,j) = atan(gradient_y(i,j)/gradient_x(i,j));
        elseif gradient_x(i,j) > 0 && gradient_y(i,j) < 0
            angles(i,j) = atan(gradient_y(i,j)/gradient_x(i,j)) + 2*pi;
        end
    end
end
x = keypoints(1,:); %u
y = keypoints(2,:); %v
s = keypoints(3,:); %scale
theta = keypoints(4,:);
sintheta = sin(theta);
costheta = cos(theta);

[xx,yy] = meshgrid(-NBP/2 : NBP/2);
wincoef = exp( -(xx.^2 + yy.^2)/NBP^2 * 2);

下列代码用于计算描述符

%========= compute descriptors ==========%
parfor p = 1: key_num % 623s
    magnitude = magnitudes; % 幅值
    sp = s(p);
    xp= x(p); % u == x
    yp= y(p); % v == y
    sinth0 = sintheta(p) ;
    costh0 = costheta(p) ;
    W = sp; % scale = 36
    ss = W/2; % scale / 2 = 18

    descriptor = zeros(NBP, NBP, NBO); % 4 * 4 * 8

    % 稀疏矩阵（极大值-起始点-终止点）
    % 归一化
    pp = magnitudes(max(-W, 1-yp)+yp : min(+W, M-yp)+yp, max(-W, 1-xp)+xp: min(W, N-xp)+xp);
    pp = pp./max(pp(:));
    xx = sort(pp(:));
    % 采样4个值
    xind1 = round((1 - 1/5) * size(xx,1));
    xind2 = round((1 - 2/5) * size(xx,1));
    xind3 = round((1 - 3/5) * size(xx,1));
    xind4 = round((1 - 4/5) * size(xx,1));
    pp = (pp>=xx(xind1)) + (pp<xx(xind1) & (pp>=xx(xind2))) * 0.75 + (pp<xx(xind2) & (pp>=xx(xind3))) * 0.5...
        +(pp<xx(xind3) & (pp>=xx(xind4))) * 0.25;
    magnitude(max(-W, 1-yp) +yp: min(+W, M -yp)+yp, max(-W, 1-xp)+xp: min(W, N- xp)+xp) = pp;

    [dx,dy] = meshgrid(max(-W, 1-xp): min(W, N - xp),max(-W, 1-yp) : min(+W, M - yp));
    nx = (costh0 * dx + sinth0 * dy ) / ss;
    ny = (-sinth0 * dx + costh0 * dy ) / ss;
    for kk = 1:numel(dx) % 统计元素个数
        % 具体的梯度直方图需要2个数据（梯度+角度）
        mag = magnitude(yp + dy(kk), xp + dx(kk));
        angle = angles(yp + dy(kk), xp + dx(kk)) ;
        angle = mod(angle - theta(p), pi);
        % Cubic interpolation  三次插值
        nt = NBO * angle / pi ;
        binx = floor( nx(kk) - 0.5 ) ;
        biny = floor( ny(kk) - 0.5 ) ;
        bint = floor( nt );
        rbinx = nx(kk) - (binx+0.5) ;
        rbiny = ny(kk) - (biny+0.5) ;
        rbint = nt - bint ;
        for dbinx = 0:1
            for dbiny = 0:1
                for dbint = 0:1
                    if  binx+dbinx >= -(NBP/2) && ...
                            binx+dbinx <   (NBP/2) && ...
                            biny+dbiny >= -(NBP/2) && ...
                            biny+dbiny <   (NBP/2) &&  ~isnan(bint)

                        weight = wincoef(binx+dbinx + NBP/2 + 1, biny+dbiny + NBP/2+ 1)...
                            * mag * abs(1 - dbinx - rbinx) * abs(1 - dbiny - rbiny) ...
                            * abs(1 - dbint - rbint) ;

                        descriptor(binx+dbinx + NBP/2 + 1, biny+dbiny + NBP/2+ 1, mod((bint+dbint),NBO)+1) = ...
                            descriptor(binx+dbinx + NBP/2+ 1, biny+dbiny + NBP/2+ 1, mod((bint+dbint),NBO)+1 ) +  weight ;
                    end
                end
            end
        end
    end
    descriptor = reshape(descriptor, 1, NBP * NBP * NBO);
    descriptor = descriptor ./ norm(descriptor);

    descriptors(:,p) = descriptor'; % arrange based on column
end

初始化变量：（举例其中p=1：623）

        magnitude ：将幅值赋值给新变量，不改变原有变量

        sp = 像素数

        xp= x的第一个元素

        yp= y的第一个元素

        sinth0 = sintheta(p)的第一个元素

        costh0 = costheta(p)的第一个元素

        W = sp; % scale = 36 表示一个ceil的半径包含36个像素

        ss = W/2; % scale / 2 = 18

        descriptor = zeros(NBP, NBP, NBO); % 4 * 4 * 8 初始化描述符

        该特征描述符的维度为 128 是因为采用了 4×4×8 的小方块，即将整个关键点周围的区域分成了 16 个子区域（4×4），每个子区域又根据梯度方向划分成了 8 个方向的直方图，因此每个子区域的直方图维度为 8，那么整个关键点周围的特征描述符维度就是 16（子区域数）×8（每个子区域的直方图维度）=128。    

parfor p = 1: key_num % 623 
   magnitude = magnitudes; % 幅值
    sp = s(p);
    xp= x(p); % u == x
    yp= y(p); % v == y
    sinth0 = sintheta(p) ;
    costh0 = costheta(p) ;
    W = sp; % scale = 36
    ss = W/2; % scale / 2 = 18

    descriptor = zeros(NBP, NBP, NBO); % 4 * 4 * 8

        pp根据关键点位置和其所处的尺度 s 获取从原始梯度图中取出的用于描述该点的梯度信息。由于特征描述符对纹理变换不敏感，因此采用了类似于可变形卷积网络中空间采样方法的方式。具体来说，根据当前关键点的位置(xp​,yp​) 和其所处的尺度 sp​，取以该点为中心，宽度为 2sp​ 的矩形区域内的梯度信息。将该矩形区域内的梯度大小归一化，使其最大值为 1，方便后续直方图统计。将归一化后的梯度大小按升序排列，并将排列结果存储在向量 xx 中，方便后续进行采样。

        根据 xx 所述的排序信息，对每个元素进行采样并把新的梯度大小重新存储在 pp 中（这里通过对 pp 进行赋值实现，具体是表示每一个点的比重），这里采用了一种简单的方式，即在最大的梯度处将其设为 1，小于最大梯度、大于次最大梯度处的梯度设置为 0.75，以此类推。将新得到的 pp 存储到原始梯度图中，方便后续处理。

• xp 和 yp 是网格中心的坐标，即网格坐标系的原点
• W 是网格坐标系半边长
• N 和 M 是网格坐标系的宽度和高度
• costh0 和 sinth0 分别表示旋转角度 th0 的余弦值和正弦值（th0 是一个弧度值）
• ss 是网格坐标系的缩放比例

        在这个网格坐标系内，每个点由 (dx, dy) 表示，它们分别表示该点相对于网格中心的横向和纵向偏移量。在这里，我们使用了 MATLAB 中的 meshgrid 函数生成了一个以网格中心为原点的二维数组 [dx,dy]。然后，我们通过如下公式将该网格坐标系旋转和缩放：

• nx = (costh0 * dx + sinth0 * dy) / ss
• ny = (-sinth0 * dx + costh0 * dy) / ss

        其中 nx 和 ny 表示旋转和缩放后的坐标系中的点的横坐标和纵坐标。这些公式实际上是对于原坐标系中的每个点 [dx, dy]，先将其进行旋转，旋转后的点坐标为 [costh0 * dx + sinth0 * dy, -sinth0 * dx + costh0 * dy]。然后再将其进行缩放，缩放比例为 ss，得到最终点的坐标 [nx, ny]。这些操作使得生成的网格坐标系在旋转和缩放后与原坐标系不同，但仍然呈矩形的形状。

        在此代码中，爆改了新坐标的计算公式

        

    % 稀疏矩阵（极大值-起始点-终止点）
    % 归一化
    pp = magnitudes(max(-W, 1-yp)+yp : min(+W, M-yp)+yp, max(-W, 1-xp)+xp: min(W, N-xp)+xp);
    pp = pp./max(pp(:));
    xx = sort(pp(:));
    % 采样4个值
    xind1 = round((1 - 1/5) * size(xx,1));
    xind2 = round((1 - 2/5) * size(xx,1));
    xind3 = round((1 - 3/5) * size(xx,1));
    xind4 = round((1 - 4/5) * size(xx,1));
    pp = (pp>=xx(xind1)) + (pp<xx(xind1) & (pp>=xx(xind2))) * 0.75 + (pp<xx(xind2) & (pp>=xx(xind3))) * 0.5...
        +(pp<xx(xind3) & (pp>=xx(xind4))) * 0.25;
    magnitude(max(-W, 1-yp) +yp: min(+W, M -yp)+yp, max(-W, 1-xp)+xp: min(W, N- xp)+xp) = pp;

    [dx,dy] = meshgrid(max(-W, 1-xp): min(W, N - xp),max(-W, 1-yp) : min(+W, M - yp));
    nx = (costh0 * dx + sinth0 * dy ) / ss;
    ny = (-sinth0 * dx + costh0 * dy ) / ss;

 构建SIFT模型：

   dx 和 dy 是之前计算得到的相对于当前关键点位置 (xp​,yp​) 的水平和竖直方向上的偏移量，在这里使用 numel() 函数获取元素个数，对其进行遍历。

        计算以当前关键点为中心的 16×16 的像素块中每个像素的梯度幅值和梯度方向角。magnitude 和 angles 分别是之前计算得到的梯度幅值和梯度方向矩阵；yp 和 xp 分别是当前关键点在图像中的纵坐标和横坐标；theta(p) 是该关键点的主方向，用来进行方向归一化。使用 mod() 函数对梯度方向角进行归一化，保证其在 0−180∘0−180∘ 范围内。

        对梯度方向角进行三次插值并计算插值后的方向和坐标，其中 NBO 是描述符每个方向划分的子区间数目，nx(kk) 和 ny(kk) 分别是在当前关键点方向下的水平和竖直偏移量。binx、biny、rbinx、rbiny、bint 和 rbint 则分别对应空间的坐标和方向的偏差值。

   binx、biny、rbinx、rbiny、bint 和 rbint 分别代表了方向梯度直方图（HOG）特征计算中，以当前关键点为中心的 16×16 的像素块中，某个像素位置和梯度方向与该关键点的偏移量。具体解释如下：

• binx 表示当前像素相对于关键点在 x 轴方向上的偏移量，单位为像素。
• biny 表示当前像素相对于关键点在 y 轴方向上的偏移量，单位为像素。
• rbinx 表示当前像素与关键点在 x 轴方向上的相对距离（相对于当前像素所在的单元格），范围为[0,1)。
• rbiny 表示当前像素与关键点在 y 轴方向上的相对距离（相对于当前像素所在的单元格），范围为[0,1)。
• bint 表示当前像素的梯度方向角与该关键点的主方向之间的偏移量，单位为弧度。
• rbint 表示当前像素的梯度方向角与该关键点的主方向之间的相对偏移量（相对于当前像素所在的方向单元格），范围为[0,1)。

这些偏移量的计算是为了将局部特征描述符归一化，并实现描述符的旋转不变性。

    

    for kk = 1:numel(dx)
        mag = magnitude(yp + dy(kk), xp + dx(kk)); 
        angle = angles(yp + dy(kk), xp + dx(kk)) ;  
        angle = mod(angle - theta(p), pi); 
% Cubic interpolation           
        nt = NBO * angle / pi ;
        binx = floor( nx(kk) - 0.5 ) ;
        biny = floor( ny(kk) - 0.5 ) ;
        bint = floor( nt );
        rbinx = nx(kk) - (binx+0.5) ;
        rbiny = ny(kk) - (biny+0.5) ;
        rbint = nt - bint ;
        for dbinx = 0:1
           for dbiny = 0:1
               for dbint = 0:1
                    if  binx+dbinx >= -(NBP/2) && ...
                        binx+dbinx <   (NBP/2) && ...
                        biny+dbiny >= -(NBP/2) && ...
                        biny+dbiny <   (NBP/2) &&  ~isnan(bint)

                          weight = wincoef(binx+dbinx + NBP/2 + 1, biny+dbiny + NBP/2+ 1)...
                              * mag * abs(1 - dbinx - rbinx) * abs(1 - dbiny - rbiny) ...
                              * abs(1 - dbint - rbint) ;

                          descriptor(binx+dbinx + NBP/2 + 1, biny+dbiny + NBP/2+ 1, mod((bint+dbint),NBO)+1) = ...
                              descriptor(binx+dbinx + NBP/2+ 1, biny+dbiny + NBP/2+ 1, mod((bint+dbint),NBO)+1 ) +  weight ;
                    end
               end
           end
        end
    end