【图像提取】基于matlab PCA-CSIFT feature图像特征提取【含Matlab源码 1174期】

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【图像提取】基于matlab PCA-CSIFT feature图像特征提取【含Matlab源码 1174期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、简介

三、部分源代码

%%initial image and create input H image.
clear;
clc
img=imread('test.png');
row=256;
colum=256;
img=imresize(img,[row,colum]);
transfer_matrix=[0.06,0.63,0.27;0.3,0.04,-0.35;0.34,-0.6,0.17];
R=img(:,:,1);
G=img(:,:,2);
B=img(:,:,3);
e=transfer_matrix(1,1)*R+transfer_matrix(1,2)*G+transfer_matrix(1,3)*B;
e1=transfer_matrix(2,1)*R+transfer_matrix(2,2)*G+transfer_matrix(2,3)*B;
e2=transfer_matrix(3,1)*R+transfer_matrix(3,2)*G+transfer_matrix(3,3)*B;
H=(double(e1)./double(e2));
origin=img;
img=H;
[m,n]=size(img);
%% Scale-Space Extrema Detection
% original sigma and the number of actave can be modified.
sigma=1.6;
octave=3;%6*sigma*k^(octave*level)<=min(m,n)/(2^(octave-2))
level=5;
k=2^(1/level);
D=cell(1,octave);
for i=1:octave
D(i)=mat2cell(zeros(row*2^(2-i)+2,colum*2^(2-i)+2,level),row*2^(2-i)+2,colum*2^(2-i)+2,level);
end
% first image in first octave is created by interpolating the original one.
temp_img=img(1:1/2:m,1:1/2:n);
temp_img=padarray(temp_img,[1,1],'pre');
temp_img=padarray(temp_img,[1,1],'replicate');
%create the DoG pyramid.
for i=1:octave
    temp_D=D{i};
    sigma=sigma*sqrt(2)^(i-1);
    for j=1:level
        p=(level)*(i-1);
        figure(1);
        subplot(octave,level,p+j);
        f=fspecial('gaussian',[1,floor(6*sigma*k^(p+j-1))],sigma*k^(p+j-1));
        L1=temp_img;
        if(i==1&&j==1)
        L2=conv2(temp_img,f,'same');
        L2=conv2(L2,f','same');
        temp_D(:,:,j)=L2-L1;
        imshow(uint8(255 * mat2gray(temp_D(:,:,j))));
        L1=L2;
        else
        L2=conv2(temp_img,f,'same');
        L2=conv2(L2,f','same');
        temp_D(:,:,j)=L2-L1;
        L1=L2;
        if(j==level)
            temp_img=L1(2:end-1,2:end-1);
        end
        imshow(uint8(255 * mat2gray(temp_D(:,:,j))));
        end
    end
    D{i}=temp_D;
    temp_img=temp_img(1:2:end,1:2:end);
    temp_img=padarray(temp_img,[1,1],'replicate');
end
%% Keypoint Localistaion
% search each pixel in the DoG map to find the extreme point
interval=level-1;
extrema=[];
for i=1:octave
    [m,n,~]=size(D{i});
    m=m-2;
    n=n-2;
    volume=m*n/(4^(i-1));
    for k=2:interval      
        for j=1:volume
            % starter=D(ceil(i/m)+1,mod(i-1,n)+1+1,j);
            sub=D{i}(ceil(j/m):ceil(j/m)+2,mod(j-1,n)+1:mod(j-1,n)+3,k-1:k+1);
            large=max(max(max(sub)));
            little=min(min(min(sub)));
            if(large==D{i}(ceil(j/m)+1,mod(j-1,n)+1+1,k))
                temp=[i,k,j,1];
                extrema=[extrema,temp];
            end
            if(little==D{i}(ceil(j/m)+1,mod(j-1,n)+1+1,k))
                temp=[i,k,j,-1];
                extrema=[extrema,temp];
            end
        end
    end
end
%% accurate keypoint localization 
%eliminate the point with low contrast or poorly localised on an edge
% x:|,y:-- x is for vertial and y is for horizontal
% value comes from the paper.
threhold=0.03;
r=10;
extr_volume=length(extrema)/4;
[m,n]=size(img);
for i=1:extr_volume
    x=floor((extrema(4*(i-1)+3)-1)/(m/(2^(extrema(4*(i-1)+1)-2))))+1;
    y=mod((extrema(4*(i-1)+3)-1),n/(2^(extrema(4*(i-1)+1)-2)))+1;
    rx=x+1;
    ry=y+1;
    rz=extrema(4*(i-1)+2);
    z=D{extrema(4*(i-1)+1)}(rx,ry,rz);
    if(abs(z)<threhold)
        extrema(4*(i-1)+4)=0;
    else 
        Dxx=D{extrema(4*(i-1)+1)}(rx-1,ry,rz)+D{extrema(4*(i-1)+1)}(rx+1,ry,rz)-2*D{extrema(4*(i-1)+1)}(rx,ry,rz);
        Dyy=D{extrema(4*(i-1)+1)}(rx,ry-1,rz)+D{extrema(4*(i-1)+1)}(rx,ry+1,rz)-2*D{extrema(4*(i-1)+1)}(rx,ry,rz);
        Dxy=D{extrema(4*(i-1)+1)}(rx-1,ry-1,rz)+D{extrema(4*(i-1)+1)}(rx+1,ry+1,rz)-D{extrema(4*(i-1)+1)}(rx-1,ry+1,rz)-D{extrema(4*(i-1)+1)}(rx+1,ry-1,rz);
        deter=Dxx*Dyy-Dxy*Dxy;
        R=(Dxx+Dyy)/deter;
        R_threshold=(r+1)^2/r;
        if(deter<0||R>R_threshold)
            extrema(4*(i-1)+4)=0;
        end
        
    end
end
flag=[];
for i=1:1:extr_volume
    if(extrema(4*(i-1)+4)==0)
        flag=[flag,(4*(i-1)+1:4*i)];
    end
end
extrema(flag)=[];
extr_volume=length(extrema)/4;
% eliminate the points on the image edge.
flag=[];
for i=1:extr_volume
    x=floor((extrema(4*(i-1)+3)-1)/(m/(2^(extrema(4*(i-1)+1)-2))))+1;
    y=mod((extrema(4*(i-1)+3)-1),n/(2^(extrema(4*(i-1)+1)-2)))+1;
    if(x==1||x==m/(2^(extrema(4*(i-1)+1)-2))||y==1||y==n/(2^(extrema(4*(i-1)+1)-2)))
        flag=[flag,(4*(i-1)+1:4*i)];
    end
end
extrema(flag)=[];
extr_volume=length(extrema)/4;
figure(2)
imshow(origin);
for i=1:extr_volume
    x=floor((extrema(4*(i-1)+3)-1)/(m/(2^(extrema(4*(i-1)+1)-2))))+1;
    y=mod((extrema(4*(i-1)+3)-1),n/(2^(extrema(4*(i-1)+1)-2)))+1;
    rx=x/2^(octave-1-extrema(4*(i-1)+1));
    ry=y/2^(octave-1-extrema(4*(i-1)+1));
    text(ry,rx,'*','Color','green');
end
hold on
%% pre-compute eigenspace of each extre point
flag=[];
A=[];
for i=1:extr_volume
    x=floor((extrema(4*(i-1)+3)-1)/(m/(2^(extrema(4*(i-1)+1)-2))))+1;
    y=mod((extrema(4*(i-1)+3)-1),n/(2^(extrema(4*(i-1)+1)-2)))+1;
    rx=x+1;
    ry=y+1;
    rz=extrema(4*(i-1)+2);
    width=41;   
    if(x>(width/2)&&y>(width/2)&&x<(m/2^(extrema(4*(i-1)+1)-2)-width/2-2)&&y<(n/2^(extrema(4*(i-1)+1)-2)-width/2-2))
    patch=D{extrema(4*(i-1)+1)}(rx-(width-1)/2:rx+(width-1)/2,ry-(width-1)/2:ry+(width-1)/2,rz);    
    [Dx,Dy]=gradient(patch);
    G=zeros(1,2*width*width);
        j=1;
        for p=1:width
            for q=1:width
                G(j)=Dx(p,q);
                G(j+1)=Dy(p,q);
                j=j+2;
            end
        end
        A=[A;G];
    else
        flag=[flag,(4*(i-1)+1:4*i)];
    end
end
extrema(flag)=[];
extr_volume=length(extrema)/4;
[p,q]=size(A);
M=A-repmat(mean(A,1),p,1);
conv=M'*M;
[V,~]=eig(conv);
N=20;
project=V(:,(end-N+1):end)';
%% Orientation Assignment(Main orientations assignment)
kpori=[];
flag=[];
for i=1:extr_volume
    %search in the certain scale
    sigma=1.6*sqrt(2)^(extrema(4*(i-1)+1));
    width=2*round(3*1.5*sigma);
    count=1;
    x=floor((extrema(4*(i-1)+3)-1)/(m/(2^(extrema(4*(i-1)+1)-2))))+1;
    y=mod((extrema(4*(i-1)+3)-1),n/(2^(extrema(4*(i-1)+1)-2)))+1;
    %make sure the point in the searchable area
    if(x>(width/2)&&y>(width/2)&&x<(m/2^(extrema(4*(i-1)+1)-2)-width/2-2)&&y<(n/2^(extrema(4*(i-1)+1)-2)-width/2-2))
        rx=x+1;
        ry=y+1;
        rz=extrema(4*(i-1)+2);
        reg_volume=width*width;%3? thereom
        % make weight matrix
        weight=zeros(width,width);
        for p=1:width
            for q=1:width
                weight(p,q)=exp(-((p-width/2)^2+(q-width/2)^2)/(1.5*sigma));
            end
        end
        %calculate region pixels' magnitude and region orientation
        reg_mag=zeros(1,count);
        reg_theta=zeros(1,count);
    for l=(rx-width/2):(rx+width/2-1)
        for k=(ry-width/2):(ry+width/2-1)
            reg_mag(count)=sqrt((D{extrema(4*(i-1)+1)}(l+1,k,rz)-D{extrema(4*(i-1)+1)}(l-1,k,rz))^2+(D{extrema(4*(i-1)+1)}(l,k+1,rz)-D{extrema(4*(i-1)+1)}(l,k-1,rz))^2);
            reg_theta(count)=atan2((D{extrema(4*(i-1)+1)}(l,k+1,rz)-D{extrema(4*(i-1)+1)}(l,k-1,rz)),(D{extrema(4*(i-1)+1)}(l+1,k,rz)-D{extrema(4*(i-1)+1)}(l-1,k,rz)))*(180/pi);
            count=count+1;
        end
    end

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/119280898