图像特征提取

% 读入图像并灰度化
I = imread('test.jpg');
I = rgb2gray(I);
% 生成高斯滤波器的核
w = fspecial('gaussian',3,0.5);
size_a = size(I);
% 进行高斯滤波
g = imfilter(I,w,'conv','symmetric','same');
% 降采样
t = g(1:2:size_a(1),1:2:size_a(2));
% 显示结果
figure(1);imshow(I);
figure(2);imshow(t);

function [ pyr ] = genPyr( img, type, level )
%   功能： 产生图像的高斯金字塔或拉普拉斯金字塔
%   输入：img-灰度图像；
%         type-变换的类型：'gauss' or 'laplace'；
%         level-分解层数
%   输出：PYR-1*LEVEL 元胞数组

pyr = cell(1,level);
pyr{1} = im2double(img);
%%%%%%图像的金字塔化%%%%%%

%计算图像的高斯金字塔
for p = 2:level
    pyr{p} = pyr_reduce(pyr{p-1});
end
if strcmp(type,'gauss'), return; end
 
% 调整图像尺寸
for p = level-1:-1:1
    osz = size(pyr{p+1})*2-1;
    pyr{p} = pyr{p}(1:osz(1),1:osz(2),:);
end

% 计算图像的拉普拉斯金字塔
for p = 1:level-1
    pyr{p} = pyr{p}-pyr_expand(pyr{p+1});
end
 
end

function [ imgout ] = pyr_reduce( img )
% 功能：计算图像的高斯金字塔
% 输入：img-灰度图像；
% 输出：imgout-分解后的金子塔图像数组

% 生成高斯核
kernelWidth = 5; 
cw = .375; 
ker1d = [.25-cw/2 .25 cw .25 .25-cw/2];
kernel = kron(ker1d,ker1d');
 
img = im2double(img);
sz = size(img);
imgout = [];
 
% 产生图像的高斯金字塔
for p = 1:size(img,3)
    img1 = img(:,:,p);
    imgFiltered = imfilter(img1,kernel,'replicate','same');
    imgout(:,:,p) = imgFiltered(1:2:sz(1),1:2:sz(2));
end
end

function [ imgout ] = pyr_expand( img )
% 功能：图像金字塔扩张
% 输入：img-待扩张的图像
% 输出：imgout-扩张后的图像

kw = 5; 
cw = .375; 
ker1d = [.25-cw/2 .25 cw .25 .25-cw/2];
kernel = kron(ker1d,ker1d')*4;
 
ker00 = kernel(1:2:kw,1:2:kw); 
ker01 = kernel(1:2:kw,2:2:kw); 
ker10 = kernel(2:2:kw,1:2:kw); 
ker11 = kernel(2:2:kw,2:2:kw); 
 
img = im2double(img);
sz = size(img(:,:,1));
osz = sz*2-1;
imgout = zeros(osz(1),osz(2),size(img,3));
 
for p = 1:size(img,3)
    img1 = img(:,:,p);
img1ph = padarray(img1,[0 1],'replicate','both');
img1pv = padarray(img1,[1 0],'replicate','both');    
    img00 = imfilter(img1,ker00,'replicate','same');
    img01 = conv2(img1pv,ker01,'valid'); 
    img10 = conv2(img1ph,ker10,'valid');
    img11 = conv2(img1,ker11,'valid');
    
    imgout(1:2:osz(1),1:2:osz(2),p) = img00;
    imgout(2:2:osz(1),1:2:osz(2),p) = img10;
    imgout(1:2:osz(1),2:2:osz(2),p) = img01;
    imgout(2:2:osz(1),2:2:osz(2),p) = img11;
end
end

%% 金字塔分解函数实现
close all;clear all;clc;
img=imread('qingdao.jpg');
n=3;
%[ pyr ] = genPyr( img, 'gauss', n );%高斯金字塔分解
[ pyr ] = genPyr( img, 'laplace', n );%拉普拉斯金字塔分解
for i=1:n
    figure(i);
    imshow(pyr{i});
end

function inv_m7 = invariable_moment(in_image)
% 功能：计算图像的Hu的七个不变矩
% 输入：in_image-RGB图像
% 输出：inv_m7-七个不变矩

% 将输入的RGB图像转换为灰度图像   
image=rgb2gray(in_image);     
%将图像矩阵的数据类型转换成双精度型
image=double(image);      
%%%=================计算 、 、 =========================
%计算灰度图像的零阶几何矩 
m00=sum(sum(image));     
m10=0;
m01=0;
[row,col]=size(image);
for i=1:row
    for j=1:col
        m10=m10+i*image(i,j);
        m01=m01+j*image(i,j);
    end
end
%%%=================计算 、 ================================
u10=m10/m00;
u01=m01/m00;
%%%=================计算图像的二阶几何矩、三阶几何矩============
m20 = 0;m02 = 0;m11 = 0;m30 = 0;m12 = 0;m21 = 0;m03 = 0;
for i=1:row
    for j=1:col
        m20=m20+i^2*image(i,j);
        m02=m02+j^2*image(i,j);
        m11=m11+i*j*image(i,j);
        m30=m30+i^3*image(i,j);
        m03=m03+j^3*image(i,j);
        m12=m12+i*j^2*image(i,j);
        m21=m21+i^2*j*image(i,j);
    end
end
%%%=================计算图像的二阶中心矩、三阶中心矩============
y00=m00;
y10=0;
y01=0;
y11=m11-u01*m10;
y20=m20-u10*m10;
y02=m02-u01*m01;
y30=m30-3*u10*m20+2*u10^2*m10;
y12=m12-2*u01*m11-u10*m02+2*u01^2*m10;
y21=m21-2*u10*m11-u01*m20+2*u10^2*m01;
y03=m03-3*u01*m02+2*u01^2*m01;
%%%=================计算图像的归格化中心矩====================
        n20=y20/m00^2;
        n02=y02/m00^2;
        n11=y11/m00^2;
        n30=y30/m00^2.5;
        n03=y03/m00^2.5;
        n12=y12/m00^2.5;
        n21=y21/m00^2.5;
%%%=================计算图像的七个不变矩======================
h1 = n20 + n02;
h2 = (n20-n02)^2 + 4*(n11)^2;
h3 = (n30-3*n12)^2 + (3*n21-n03)^2;  
h4 = (n30+n12)^2 + (n21+n03)^2;
h5 = (n30-3*n12)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n21-n03)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
h6 = (n20-n02)*((n30+n12)^2-(n21+n03)^2)+4*n11*(n30+n12)*(n21+n03);
h7 = (3*n21-n03)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n12-n30)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
 
inv_m7= [h1 h2 h3 h4 h5 h6 h7];   

function        [A_nm,zmlist,cidx,V_nm]        = zernike(img,n,m)
% 功能：计算输入图像的Zernike矩
% 输入：img-灰度图像
%       n-阶数
% 输出：V_nm-n阶的Zernike多项式，定义为在极坐标系中p，theta的函数
%       cidx-表示虚部值
%       A_nm-Zernike矩

if nargin>0
if nargin==1
        n=0;
end
d=size(img);
img=double(img);
% 取步长
xstep=2/(d(1)-1); 
ystep=2/(d(2)-1);
% 画方格
[x,y]=meshgrid(-1:xstep:1,-1:ystep:1); 
circle1= x.^2 + y.^2; 
% 提取符合circle1<=1的数
inside=find(circle1<=1); 
% 构造size（d）*size(d)的矩阵
mask=zeros(d);
% 构造size(inside)*size(inside)的全为1的矩阵赋值给mask（inside）
mask(inside)=ones(size(inside)); 
% 计算图像的复数表示
[cimg,cidx]=clipimg(img,mask); 
% 计算Z的实部和虚部
z=clipimg(x+i*y,mask); 
% 计算复数的模，sqrt(x,y),z=x+iy
p=0.9*abs(z);   ;
% 计算复数z的辐角值（tanz）
theta=angle(z); 
c=1;
for order=1:length(n)
        n1=n(order);
           if nargin<3
                m=zpossible(n1);
        end
        for r=1:length(m)
  V_nmt=zpoly(n1,m(r),p,theta); 
 % conj是求复数的共轭
zprod=cimg.*conj(V_nmt); 
 % (n1+1)/π*∑∑(zprod); 对于图像而言求和代替了求积分
A_nm(c)=(n1+1)*sum(sum(zprod))/pi; 
                zmlist(c,1:2)=[n1 m(r)]; 
                if nargout==4
                        V_nm(:,c)=V_nmt;
                end
                c=c+1;
        end
end
else
end

%%%%%子函数%%%%%
function [cimg,cindex,dim]=clipimg(img,mask)
%功能：计算复数的实部和虚部
dim=size(img);
cindex=find(mask~=0);
cimg=img(cindex);
return;

%%%%%子函数%%%%%
function  [m]=zpossible(n)
% 功能：判断n是偶数还是奇数，是偶数时，m取0,2,4,6等,否则取奇数赋值m
if iseven(n)
        m=0:2:n;
else
        m=1:2:n;
end
return;

%%%%%子函数%%%%%
function        [V_nm,mag,phase]=zpoly(n,m,p,theta)
%功能：计算Zernike矩多项式
R_nm=zeros(size(p)); % 产生size(p)*size(p)的零矩阵赋给R_nm
a=(n+abs(m))/2;
b=(n-abs(m))/2;
total=b;
for s=0:total
        num=((-1)^s)*fac(n-s)*(p.^(n-2*s)); % (-1).-1*(n-s)!r.^(n-2*s)
        den=fac(s)*fac(a-s)*fac(b-s); % s!*(a-s)!*(b-s)!
        R_nm=R_nm + num/den; % R_nm是一个实数值的径向多项式
end
mag=R_nm; % 赋值
phase=m*theta; 
V_nm=mag.*exp(i*phase); % V_nm为n阶的Zernike多项式，定义为在极坐标系中p，theta的函数
return;

%%%%%子函数%%%%%
function [factorial]=fac(n)
%功能：求n的阶乘
maxno=max(max(n));
zerosi=find(n<=0); %取n小于等于0的数
n(zerosi)=ones(size(zerosi));
factorial=n;
findex=n;
for i=maxno:-1:2
        cand=find(findex>2);
        candidates=findex(cand);
        findex(cand)=candidates-1;
        factorial(cand)=factorial(cand).*findex(cand);
end
return;

function [verdict]=iseven(candy)
verdict=zeros(size(candy));
isint=find(isint(candy)==1);
divided2=candy(isint)/2;
evens=(divided2==floor(divided2));
verdict(isint)=evens;
return;

function [verdict]=isint(candy)
verdict  = double(round(candy))==candy;
return;

%% 函数实现
close all;clear all;clc;
%% 
in_image=imread('qingdao.jpg');
inv_m7 = invariable_moment(in_image)    %%Hu矩
%% 
img=imread('qingdao.jpg');
img=rgb2gray(img);
[A_nm,zmlist,cidx,V_nm] = zernike(img);
A_nm    %%Zernike矩

close all;clear all;clc;
f=imread('xinglong.jpg');
f=rgb2gray(f);
%% 调用函数
e=log_edge(f,0.5);
figure(1);imshow(f);
figure(2);imshow(e);

function e=log_edge(a, sigma)
%功能：实现LoG算子提取边缘点
%输入：a-灰度图像
%      sigma-滤波器参数
%输出：e-边缘图像

%产生同样大小的边缘图像e，初始化为0.
[m,n]=size(a);
e=repmat(logical(uint8(0)),m,n);
rr=2:m-1;cc=2:n-1;
%选择点数为奇数的滤波器的尺寸fsize>6*sigma;
fsize=ceil(sigma*3)*2+1;
%产生LoG滤波器
op=fspecial('log',fsize,sigma);
%将LoG滤波器的均值变为0.
op=op-sum(op(:))/prod(size(op));
%利用LoG算子对图像滤波
b=filter2(op,a);
%设置过零检测的门限
%寻找滤波后的过零点，+-和-+表示水平方向从左到右和从右到左过零
%[+-]'和[-+]'表示垂直方向从上到下和从下到上过零
%这里我们选择边缘点为值为负的点
thresh=.75*mean2(abs(b(rr,cc)));
%[- +]的情况
[rx,cx]=find(b(rr,cc)<0&b(rr,cc+1)>0&abs(b(rr,cc)-b(rr,cc+1))>thresh); 
e((rx+1)+cx*m)=1;
%[- +]的情况
[rx,cx]=find(b(rr,cc-1)>0&b(rr,cc)<0&abs(b(rr,cc-1)-b(rr,cc))>thresh); 
e((rx+1)+cx*m)=1;
%[- +]的情况
[rx,cx]=find(b(rr,cc)<0&b(rr+1,cc)>0&abs(b(rr,cc)-b(rr+1,cc))>thresh); 
e((rx+1)+cx*m)=1;
%[- +]的情况
[rx,cx]=find(b(rr-1,cc)>0&b(rr,cc)<0&abs(b(rr-1,cc)-b(rr,cc))>thresh); 
e((rx+1)+cx*m)=1;
%某些情况下LoG滤波结果可能正好为0，下面考虑这种情况：
 %寻找滤波后的过零
%+0-和-0+表示水平方向从左到右和从右到左过零
%[+0-]'和[-0+]'表示垂直方向从上到下和从下到上过零
%边缘正好位于滤波值为零点上
[rz,cz]=find(b(rr,cc)==0);
if~isempty(rz)
%零点的线性坐标
zero=(rz+1)+cz*m;   
%[-0+]的情况
zz=find(b(zero-1)<0&b(zero+1)>0&abs(b(zero-1)-b(zero+1))>2*thresh);
e(zero(zz))=1;
%[+0-]'的情况
zz=find(b(zero-1)>0&b(zero+1)<0&abs(b(zero-1)-b(zero+1))>2*thresh);
e(zero(zz))=1;
%[-0+]的情况
zz=find(b(zero-m)<0&b(zero+m)>0&abs(b(zero-1)-b(zero+1))>2*thresh);
e(zero(zz))=1;
%[+0-]'的情况
zz=find(b(zero-m)>0&b(zero+m)<0&abs(b(zero-1)-b(zero+1))>2*thresh);
e(zero(zz))=1;
end

close all;clear all;clc;
f=imread('diaohua.jpg');
f=rgb2gray(f);
e = susan(f,25);
figure(1);imshow(f);
figure(2);imshow(e);

function image_out = susan(im,threshold)
% 功能：实现运用SUNSAN算子进行边缘检测
% 输入：image_in-输入的待检测的图像
%       threshold-阈值
% 输出：image_out-检测边缘出的二值图像

% 将输入的图像矩阵转换成double型
d = length(size(im));
if d==3
    image=double(rgb2gray(im));
elseif d==2
    image=double(im);
end
 
% 建立SUSAN模板
 
mask = ([ 0 0 1 1 1 0 0 ;0 1 1 1 1 1 0;1 1 1 1 1 1 1;1 1 1 1 1 1 1;1 1 1 1 1 1 1;0 1 1 1 1 1 0;0 0 1 1 1 0 0]);  
 
R=zeros(size(image));
% 定义USAN 区域
nmax = 3*37/4;

 [a b]=size(image);
new=zeros(a+7,b+7);
[c d]=size(new);
new(4:c-4,4:d-4)=image;
  
for i=4:c-4
    
    for j=4:d-4
        
        current_image = new(i-3:i+3,j-3:j+3);
        current_masked_image = mask.*current_image;
   
%   调用susan_threshold函数进行阈值比较处理
                
        current_thresholded = susan_threshold(current_masked_image,threshold);
        g=sum(current_thresholded(:));
        
        if nmax<g
            R(i,j) = g-nmax;
        else
            R(i,j) = 0;
        end
    end
end
 
image_out=R(4:c-4,4:d-4);

function thresholded = susan_threshold(image,threshold)
% 功能：设定SUSAN算法的阈值

[a b]=size(image);
intensity_center = image((a+1)/2,(b+1)/2);
 
temp1 = (image-intensity_center)/threshold;
temp2 = temp1.^6;
thresholded = exp(-1*temp2);

%% 基于小波变换模极大值的边缘检测
clear all;
load wbarb;
I = ind2gray(X,map);imshow(I);
I1 = imadjust(I,stretchlim(I),[0,1]);figure;imshow(I1);
[N,M] = size(I);
 
h = [0.125,0.375,0.375,0.125];
g = [0.5,-0.5];
delta = [1,0,0];
 
J = 3;
 
a(1:N,1:M,1,1:J+1) = 0;
dx(1:N,1:M,1,1:J+1) = 0;
dy(1:N,1:M,1,1:J+1) = 0;
d(1:N,1:M,1,1:J+1) = 0;
 
a(:,:,1,1) = conv2(h,h,I,'same');
dx(:,:,1,1) = conv2(delta,g,I,'same');
dy(:,:,1,1) = conv2(g,delta,I,'same');
 
x = dx(:,:,1,1);
y = dy(:,:,1,1);
d(:,:,1,1) = sqrt(x.^2+y.^2);
I1 = imadjust(d(:,:,1,1),stretchlim(d(:,:,1,1)),[0 1]);figure;imshow(I1);
 
lh = length(h);
lg = length(g);
 
for j = 1:J+1
  lhj = 2^j*(lh-1)+1;
  lgj = 2^j*(lg-1)+1;
  hj(1:lhj)=0;
  gj(1:lgj)=0;
  for n = 1:lh
    hj(2^j*(n-1)+1)=h(n);
  end
 
  for n = 1:lg
    gj(2^j*(n-1)+1)=g(n);
  end
  
  a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
  dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same');
  dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same');
 
  x = dx(:,:,1,j+1);
  y = dy(:,:,1,j+1);
  dj(:,:,1,j+1) = sqrt(x.^2+y.^2);
 
  I1 = imadjust(dj(:,:,1,j+1),stretchlim(dj(:,:,1,j+1)),[0 1]);figure;imshow(I1);
end  

% 输入图像，并将其转化成灰度图像
I=imread('qipan.jpg');
I=rgb2gray(I);

% 构造卷积核
%% 提取右侧上方边缘
F1=[-1 0 0 -1 -1
     1 1 1 0 0];
% 进行卷积运算
A1=conv2(double(I),double(F1));
% 转换成8位无符号整型并显示
A1=uint8(A1);
figure(1);imshow(A1);

% 构造卷积核
%% 提取左侧上方边缘
F2=[-1 -1 0 0 -1
     0 0 1 1 1];
% 进行卷积运算
A2=conv2(double(I),double(F2));
% 转换成8位无符号整型并显示
A2=uint8(A2);
figure(2);imshow(A2);

function e = straightline(f,alpha)
%% 功能：检测特定角度的边缘
%% 输入：f-待检测的图像
%%       alpha-检测角度
%%             1表示H垂直方向
%%             2表示V水平方向
%%             3表示45度
%%             4表示-45度
f = im2double(f);
%% 构造卷积核
H = [-1 -1 -1;2 2 2;-1 -1 -1];
V = [-1 2 -1;-1 2 -1;-1 2 -1];
P45 = [-1 -1 2;-1 2 -1;2 -1 -1];
M45 = [2 -1 -1;-1 2 -1;-1 -1 2];
switch(alpha)
    case 1
        e = imfilter(f,H);
    case 2
        e = imfilter(f,V);
    case 3
        e = imfilter(f,P45);
    case 4
        e = imfilter(f,M45);
    otherwise
        e = ones(size(f));
        fprintf('请检测输入：\nalpha-检测角度\n1表示H垂直方向\n2表示V水平方向\n3表示45度\n4表示-45度\n');
end

close all;clear all;clc;
f=imread('qipan.jpg');
f=rgb2gray(f);
%% 输入：f-待检测的图像
%%       alpha-检测角度
%%             1表示H垂直方向
%%             2表示V水平方向
%%             3表示45度
%%             4表示-45度
e1 = straightline(f,1);
e2 = straightline(f,2);
e3 = straightline(f,3);
e4 = straightline(f,4);
figure(1);imshow(f);
figure(2);imshow(e1);
figure(3);imshow(e2);
figure(4);imshow(e3);
figure(5);imshow(e4);

% 读入并显示原始图像
I=imread('hehua.jpg');
grayI=rgb2gray(I);
figure,imshow(grayI)

% 利用单尺度形态学梯度进行边缘检测
se=strel('square',3);
grad=imdilate(grayI,se)-imerode(grayI,se);
figure,imshow(grad)

% 利用多尺度形态学梯度进行边缘检测
se1=strel('square',1);
se2=strel('square',3);
se3=strel('square',5);
se4=strel('square',7);
grad1=imerode((imdilate(grayI,se2)-imerode(grayI,se2)),se1);
grad2=imerode((imdilate(grayI,se3)-imerode(grayI,se3)),se2);
grad3=imerode((imdilate(grayI,se4)-imerode(grayI,se4)),se3);
multiscaleGrad=(grad1+grad2+grad3)/3;
figure,imshow(multiscaleGrad)

close all;clear all;clc;
%% 运行时间较长
img=imread('sunflower.jpg');
imshow(img)
pt=log_Blob(rgb2gray(img));
draw(img,pt,'LoG Lindeberg');

function [points] = log_Blob(img,o_nb_blobs)
    % 功能：提取LoG斑点
    % 输入： 
    %      img –输入的图像
    %       o_nb_blobs -需要检测的斑点区域的数量
    % 输出：
    %       points -检测出的斑点
    % 参考文献：
    % Lindeberg, T. Feature Detection with Automatic Scale Selection
    % IEEE Transactions Pattern Analysis Machine Intelligence, 1998, 30, 
    % 77-116
    
    % 输入图像
    img = double(img(:,:,1));
        
    % 设定检测到斑点的数量
    if nargin==1
        nb_blobs = 120;
    else
        nb_blobs = o_nb_blobs;
    end
    
    % 设定LoG参数
    sigma_begin = 2;
    sigma_end   = 15;
    sigma_step  = 1;
    sigma_array = sigma_begin:sigma_step:sigma_end;
    sigma_nb    = numel(sigma_array);
        
    % 变量
    img_height  = size(img,1);
    img_width   = size(img,2);
        
    % 计算尺度规范化高斯拉普拉斯算子
    snlo = zeros(img_height,img_width,sigma_nb);
    for i=1:sigma_nb
        sigma       = sigma_array(i);

snlo(:,:,i) = sigma*sigma*imfilter(img,fspecial('log', floor(6*sigma+1), sigma),'replicate');

    end
        
    % 搜索局部极值
    snlo_dil             = imdilate(snlo,ones(3,3,3));
    blob_candidate_index = find(snlo==snlo_dil);
    blob_candidate_value = snlo(blob_candidate_index);
    [tmp,index]          = sort(blob_candidate_value,'descend');
    
blob_index           = blob_candidate_index( index(1:min(nb_blobs,numel(index))) );

[lig,col,sca]        = ind2sub([img_height,img_width,sigma_nb],blob_index);

    points               = [lig,col,3*reshape(sigma_array(sca),[size(lig,1),1])];
    
end

function draw(img,pt,str)
% 功能：在图像中绘制出特征点
% 输入： 
%      img –输入的图像
%       pt-检测出的特征点的坐标
%       str-在图上显示的名称

    figure('Name',str);
    imshow(img);
    hold on;
    axis off;
    switch size(pt,2)
        case 2
            s = 2;
            for i=1:size(pt,1)
                rectangle('Position',[pt(i,2)-s,pt(i,1)-s,2*s,2*s],'Curvature',[0 0],'EdgeColor','b','LineWidth',2);
            end
        case 3
            for i=1:size(pt,1)
                rectangle('Position',[pt(i,2)-pt(i,3),pt(i,1)-pt(i,3),2*pt(i,3),2*pt(i,3)],'Curvature',[1,1],'EdgeColor','w','LineWidth',2);
            end
    end
end

function points = gilles(im,o_radius)
    % 功能：提取Gilles斑点
    % 输入：
    %      im –输入的图像，灰度图
    %      o_radius -斑点检测半径（可选）
    % 输出：
    %      points -检测出的斑点
    %
    % 参考文献：
    % S. Gilles, Robust Description and Matching of Images. PhD thesis,
    % Oxford University, Ph.D. thesis, Oxford University, 1988.
   
    im = im(:,:,1);
 
    % 变量
    if nargin==1
        radius = 10;
    else
        radius = o_radius;
end
% 建立掩模(mask)区域
    mask = fspecial('disk',radius)>0;
 
    % 计算掩模区域的局部熵值
    loc_ent = entropyfilt(im,mask);%% 高版本Matlab
 
    % 寻找局部最值
    [l,c,tmp] = findLocalMaximum(loc_ent,radius);
 
    % 超过阈值的斑点确定为所要提取的斑点
    [l,c]     = find(tmp>0.95*max(tmp(:)));
    points    = [l,c,repmat(radius,[size(l,1),1])];
 
end

function [row,col,max_local] = findLocalMaximum(val,radius)
    % 功能：查找邻域极大值 
    % 输入：
    % val -NxM 矩阵；
    % radius –邻域半径；
    % 输出：
    % row –邻域极大值的行坐标；
    % col – 邻域极大值的列坐标；
    % max_local-邻域极大值。
  
    mask  = fspecial('disk',radius)>0;
    nb    = sum(mask(:));
    highest          = ordfilt2(val, nb, mask);
    second_highest   = ordfilt2(val, nb-1, mask);
    index            = highest==val & highest~=second_highest;
    max_local        = zeros(size(val));
    max_local(index) = val(index);
    [row,col]        = find(index==1);
 
end

close all;clear all;clc;
in_img=imread('long.jpg');
[posr,posc]=Harris1(in_img,0.04);

function [posr, posc]=Harris1(in_image,a)
% 功能：检测图像的Harris角点
% 输入：in_image－待检测的RGB图像数组
%       a－角点参数响应，取值范围为：0.04～0.06
% 输出：posr－所检测出角点的行坐标向量
%       posc－所检测出角点的列坐标向量

% 将RGB图像转化成灰度图像
in_image=rgb2gray(in_image);        
% unit8型转化为双精度double64型
ori_im=double(in_image); 
           
%%%%%%计算图像在x、y 两个方向的梯度%%%%%%
% x方向梯度算子模板
fx = [-1 0 1];                       
% x方向滤波
Ix = filter2(fx,ori_im);
% y方向梯度算子
fy = [-1;0;1];
% y方向滤波
Iy = filter2(fy,ori_im);  
              
%%%%%%计算两个方向的梯度乘积%%%%%%
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;

%%%%%%使用高斯函数对梯度乘积进行加权%%%%%%
% 产生7*7的高斯窗函数，sigma=2
h= fspecial('gaussian',[7 7],2);
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);

%%%%%%计算每个像元的Harris响应值%%%%%%
[height,width]=size(ori_im);
R = zeros(height,width);
% 像素(i,j)处的Harris响应值
for i = 1:height
    for j = 1:width
        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
        R(i,j) = det(M)-a*(trace(M))^2;          
    end
end

%%%%%%去掉小于阈值的Harris响应值%%%%%%
Rmax=max(max(R));
%  阈值
t=0.01* Rmax;
for i = 1:height
for j = 1:width
if R(i,j)<t
  R(i,j) = 0;
  end
end
end

%%%%%%进行3×3邻域非极大值抑制%%%%%%      
% 进行非极大抑制，窗口大小3*3
corner_peaks=imregionalmax(R);
countnum=sum(sum(corner_peaks));

%%%%%%显示所提取的Harris角点%%%%%% 
% posr是用于存放行坐标的向量
[posr, posc] = find(corner_peaks== 1);
% posc是用于存放列坐标的向量
figure(1);imshow(in_image);
hold on
for i = 1 : length(posr)
    plot(posc(i),posr(i),'r+');
end

function [row,col,max_local] = findLocalMaximum(val,radius)
    % 功能：查找邻域极大值 
    % 输入：
    % val -NxM 矩阵；
    % radius –邻域半径；
    % 输出：
    % row –邻域极大值的行坐标；
    % col – 邻域极大值的列坐标；
    % max_local-邻域极大值。
  
    mask  = fspecial('disk',radius)>0;
    nb    = sum(mask(:));
    highest          = ordfilt2(val, nb, mask);
    second_highest   = ordfilt2(val, nb-1, mask);
    index            = highest==val & highest~=second_highest;
    max_local        = zeros(size(val));
    max_local(index) = val(index);
    [row,col]        = find(index==1);
 
end

C = cornermetric(I,'Harris');

close all;clear all;clc;
%% 运行时间长
img=imread('door.jpg');
img=rgb2gray(img);
points=harrislaplace(img);
draw(img,points,'角点检测');

function points = harrislaplace(img)
% 功能：提取Harris-Laplace角点
% 输入：img-输入的RGB图像
% 输出：points-检测出的角点矩阵
% 参考文献：
% ==========
% K. Mikolajczyk and C. Schmid. Scale & affine invariant interest point % detectors.
% International Journal of Computer Vision, 2004
    
% 图像参数

    img         = double(img(:,:,1));
    img_height  = size(img,1);
    img_width   = size(img,2);
 
    % 尺度参数
    sigma_begin = 1.5;
    sigma_step  = 1.2;
    sigma_nb    = 13;
    sigma_array = (sigma_step.^(0:sigma_nb-1))*sigma_begin;
 
 
    % 第一部分：提取Harris角点
    harris_pts = zeros(0,3);
    for i=1:sigma_nb
 
        % 尺度
        % 积分尺度
s_I = sigma_array(i);   
        % 微分尺度
s_D = 0.7*s_I;
 
        
% 微分掩模
        x  = -round(3*s_D):round(3*s_D);
        dx = x .* exp(-x.*x/(2*s_D*s_D)) ./ (s_D*s_D*s_D*sqrt(2*pi));
        dy = dx';
 
        % 图像微分
        Ix = conv2(img, dx, 'same');
        Iy = conv2(img, dy, 'same');
 
        % 自相关矩阵
        g   = fspecial('gaussian',max(1,fix(6*s_I+1)), s_I);
        Ix2 = conv2(Ix.^2, g,  'same');
        Iy2 = conv2(Iy.^2, g,  'same');
        Ixy = conv2(Ix.*Iy, g, 'same');
 
        k = 0.06; cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2; 
       
 % 查询邻域极值
        [l,c,max_local] = findLocalMaximum(cim,3*s_I); 
        
% 设定局部邻域极值阈值
        t = 0.2*max(max_local(:));
 
        % 查找大于邻域极值阈值的点
        [l,c] = find(max_local>=t);
 
        n = size(l,1);
        harris_pts(end+1:end+n,:) = [l,c,repmat(i,[n,1])];
    end
 
 
    %%%%%%第二部分: LAPLACE变换%%%%%%
    % 计算尺度归一化Laplace算子
    laplace_snlo = zeros(img_height,img_width,sigma_nb);
    for i=1:sigma_nb
        % 尺度
s_L = sigma_array(i);
        laplace_snlo(:,:,i) = s_L*s_L*imfilter(img,fspecial('log', floor(6*s_L+1), s_L),'replicate');
    end
% 检测每个特征点在某一尺度LoG相应是否达到最大
    n   = size(harris_pts,1);
    cpt = 0;
    points = zeros(n,3);
    for i=1:n
        l = harris_pts(i,1);
        c = harris_pts(i,2);
        s = harris_pts(i,3);
        val = laplace_snlo(l,c,s);
        if s>1 && s<sigma_nb
            if val>laplace_snlo(l,c,s-1) && val>laplace_snlo(l,c,s+1)
                cpt = cpt+1;
                points(cpt,:) = harris_pts(i,:);
            end
        elseif s==1
            if val>laplace_snlo(l,c,2)
                cpt = cpt+1;
                points(cpt,:) = harris_pts(i,:);
            end
        elseif s==sigma_nb
            if val>laplace_snlo(l,c,s-1)
                cpt = cpt+1;
                points(cpt,:) = harris_pts(i,:);
            end
        end
    end
    points(cpt+1:end,:) = [];
    
       points(:,3) = 3*sigma_array(points(:,3));
       end

function [row,col,max_local] = findLocalMaximum(val,radius)
    % 功能：查找邻域极大值 
    % 输入：
    % val -NxM 矩阵；
    % radius –邻域半径；
    % 输出：
    % row –邻域极大值的行坐标；
    % col – 邻域极大值的列坐标；
    % max_local-邻域极大值。
  
    mask  = fspecial('disk',radius)>0;
    nb    = sum(mask(:));
    highest          = ordfilt2(val, nb, mask);
    second_highest   = ordfilt2(val, nb-1, mask);
    index            = highest==val & highest~=second_highest;
    max_local        = zeros(size(val));
    max_local(index) = val(index);
    [row,col]        = find(index==1);
 
end

function draw(img,pt,str)
% 功能：在图像中绘制出特征点
% 输入： 
%      img –输入的图像
%       pt-检测出的特征点的坐标
%       str-在图上显示的名称

figure('Name',str);
    imshow(img);
    hold on;
    axis off;
    switch size(pt,2)
        case 2
            s = 2;
            for i=1:size(pt,1)
                rectangle('Position',[pt(i,2)-s,pt(i,1)-s,2*s,2*s],'Curvature',[0 0],'EdgeColor','b','LineWidth',2);
            end
        case 3
            for i=1:size(pt,1)
                rectangle('Position',[pt(i,2)-pt(i,3),pt(i,1)-pt(i,3),2*pt(i,3),2*pt(i,3)],'Curvature',[1,1],'EdgeColor','w','LineWidth',2);
            end
    end
end