在MATLAB的Computer Vision System Toolbox™工具箱中，有FAST，Harris和Shi & Tomasi 角点检测子和SURF、MSER斑点检测子。这个工具箱还有SURF,FREAK,BRISK,LBP以及HOG描述子。

局部点的检测和提取：

一般用于：1）用于定位图像拼接或三维重建的锚点。2）在不需要图像分割的情况下，紧凑地表示图像内容以进行检测或分类。

1，图像配准与拼接MATLAB案例

第一步：加载需要拼接的图像；

% Load images.
buildingDir = fullfile(toolboxdir('vision'), 'visiondata', 'building');
buildingScene = imageSet(buildingDir);

% Display images to be stitched
montage(buildingScene.ImageLocation)

第二步：存储图像对；在I(n)和I(n-1)之间检测和匹配特征点；估计几何变换T(n)，从I(n)映射到I(n-1)；计算和转换映射I(n)成全景图像T(1)*...T(n-1)*T(n).

% Read the first image from the image set.
I = read(buildingScene, 1);

% Initialize features for I(1)
grayImage = rgb2gray(I);
points = detectSURFFeatures(grayImage);
[features, points] = extractFeatures(grayImage, points);

% Initialize all the transforms to the identity matrix. Note that the
% projective transform is used here because the building images are fairly
% close to the camera. Had the scene been captured from a further distance,
% an affine transform would suffice.
tforms(buildingScene.Count) = projective2d(eye(3));

% Iterate over remaining image pairs
for n = 2:buildingScene.Count

    % Store points and features for I(n-1).
    pointsPrevious = points;
    featuresPrevious = features;

    % Read I(n).
    I = read(buildingScene, n);

    % Detect and extract SURF features for I(n).
    grayImage = rgb2gray(I);
    points = detectSURFFeatures(grayImage);
    [features, points] = extractFeatures(grayImage, points);

    % Find correspondences between I(n) and I(n-1).
    indexPairs = matchFeatures(features, featuresPrevious, 'Unique', true);

    matchedPoints = points(indexPairs(:,1), :);
    matchedPointsPrev = pointsPrevious(indexPairs(:,2), :);

    % Estimate the transformation between I(n) and I(n-1).
    tforms(n) = estimateGeometricTransform(matchedPoints, matchedPointsPrev,...
        'projective', 'Confidence', 99.9, 'MaxNumTrials', 2000);

    % Compute T(1) * ... * T(n-1) * T(n)
    tforms(n).T = tforms(n-1).T * tforms(n).T;
end

开始使用projective2d outputLimits 方法来找到每个变换的输出限制。然后，使用输出限制来自动找到大致位于场景中心的图像。

imageSize = size(I);  % all the images are the same size
% Compute the output limits  for each transform
for i = 1:numel(tforms)
    [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(2)], [1 imageSize(1)]);
end

接下来，计算每个变换的平均极限X，并找到位于中心的图像。这里只使用限制X，因为场景是水平的。如果使用另一组图像，则可能需要使用X和Y限制来找到中心图像。

avgXLim = mean(xlim, 2);

[~, idx] = sort(avgXLim);

centerIdx = floor((numel(tforms)+1)/2);

centerImageIdx = idx(centerIdx);
Tinv = invert(tforms(centerImageIdx));

for i = 1:numel(tforms)
    tforms(i).T = Tinv.T * tforms(i).T;
end

第三步：初始化全景图像

for i = 1:numel(tforms)
    [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(2)], [1 imageSize(1)]);
end

% Find the minimum and maximum output limits
xMin = min([1; xlim(:)]);
xMax = max([imageSize(2); xlim(:)]);

yMin = min([1; ylim(:)]);
yMax = max([imageSize(1); ylim(:)]);

% Width and height of panorama.
width  = round(xMax - xMin);
height = round(yMax - yMin);

% Initialize the "empty" panorama.
panorama = zeros([height width 3], 'like', I);

第四步：创造全景图像

blender = vision.AlphaBlender('Operation', 'Binary mask', ...
    'MaskSource', 'Input port');

% Create a 2-D spatial reference object defining the size of the panorama.
xLimits = [xMin xMax];
yLimits = [yMin yMax];
panoramaView = imref2d([height width], xLimits, yLimits);

% Create the panorama.
for i = 1:buildingScene.Count

    I = read(buildingScene, i);

    % Transform I into the panorama.
    warpedImage = imwarp(I, tforms(i), 'OutputView', panoramaView);

    % Overlay the warpedImage onto the panorama.
    panorama = step(blender, panorama, warpedImage, warpedImage(:,:,1));
end

figure
imshow(panorama)

 

特征提取之灰度游程（行程）矩阵-GLRLM
Matlab特征提取之灰度游程（行程）矩阵-GLRLM

Matlab特征提取之灰度游程（行程）矩阵-GLRLM
前几天参加数学建模比赛，发现全网图形识别检测领域与灰度游程矩阵特征提取的matlab的参考材料很少，下面给出了求灰度游程矩阵以及部分特征提取的代码，供大家交流。
灰度游程（行程）矩阵
一幅图像的灰度游程矩阵反映了图像灰度关于方向，相邻间隔和变化幅度等

综合信息。灰度游程矩阵是对分析影像的局部模式和其排列规则基础之一。灰

度游程矩阵可以实现对一幅图像中同一方向同一灰度值连续出现个数的统计。在

一幅图像上，在某一方向上连续的像素点具有相同的灰度值，灰度游程矩阵就是

通过对这些像素点的分布进行统计得到纹理特征。本文中，我们使用P(i,j|θ)代表

灰度游程矩阵，(i，j)点代表灰度为i的像素在图像θ方向连续出现j次的计数。此

P(i,j|θ)代表灰度游程矩阵，θ一般有0度、45度、90度和135度。

$\mathrm{N}_{\mathrm{g}}$表示图像上的灰度级数目

$N_{\mathrm{r}}$表示图像上不同游程的数目

$\mathrm{N}_{\mathrm{p}}$表示图像上像素点的数目：
举个例子：矩阵A为：
$\left( \begin{array}{cccccccc}{0} & {0} & {1} & {1} & {2} & {2} & {3} & {3} \\ {1} & {1} & {2} & {2} & {3} & {3} & {0} & {0} \\ {2} & {2} & {3} & {3} & {0} & {0} & {1} & {1} \\ {3} & {3} & {0} & {0} & {1} & {1} & {2} & {2}\end{array}\right)$
那么P(i,j|00)（θ=00时），水平统计各像素出现的次数：
P(i,j|0)=$\left( \begin{array}{cccc}{0} & {4} & {0} & {0} \\ {0} & {4} & {0} & {0} \\ {0} & {4} & {0} & {0} \\ {0} & {4} & {0} & {0}\end{array}\right)$
0出现1次的次数为0,0出现两次的次数为4,0出现3次次数为0，以此类推…
同理θ=450：
P(i,j|450)=$\left( \begin{array}{cccc}{4} & {2} & {0} & {0} \\ {4} & {2} & {0} & {0} \\ {4} & {2} & {0} & {0} \\ {2} & {3} & {0} & {0}\end{array}\right)$
同理θ=900：
P(i,j|900)=$\left( \begin{array}{cccc}{8} & {0} & {0} & {0} \\ {8} & {0} & {0} & {0} \\ {8} & {0} & {0} & {0} \\ {8} & {0} & {0} & {0}\end{array}\right)$
代码：
// An highlighted block
%-----------------------------------
% author：Courage
% name：tdcup.model2.灰度游程矩阵
%-----------------------------------
function [glrlms,si]= grayrlmatrix(varargin)

[i, offset, nl, gl] = parseinputs(varargin{:});

if gl(2) == gl(1)
    si = ones(size(i));
else
    slope = (nl-1) / (gl(2) - gl(1));
    intercept = 1 - (slope*(gl(1)));
    si = round(imlincomb(slope,i,intercept,'double'));
end

si(si > nl) = nl;
si(si < 1) = 1;
numoffsets = size(offset,1);

if nl ~= 0
   
    for k = 1 : numoffsets
        glrlms{k} = computeglrlm(si,offset(k),nl);
    end
else
    glrlms = [];
end

% --------------------------------------------------------------------
function oneglrlm = computeglrlm(si,offset,nl)

switch offset
    case 1
        % 0 degree
        oneglrlm = rle_0(si,nl);
    case 2
        % 45 degree
        seq = zigzag(si);
        oneglrlm  = rle_45(seq,nl);
    case 3
        % 90 degree
        oneglrlm = rle_0(si',nl);
    case 4
       
        seq = zigzag(fliplr(si));
        oneglrlm = rle_45(seq,nl);
    otherwise
        error('only 4 directions supported')
end


% --------------------------------------------------------------------
function [i, offset, nl, gl] = parseinputs(varargin)

iptchecknargin(1,7,nargin,mfilename);

i = varargin{1};
iptcheckinput(i,{'logical','numeric'},{'2d','real','nonsparse'}, ...
    mfilename,'i',1);

offset = [1;2;3;4];

if islogical(i)
    nl = 2;
else
    nl = 8;
end
gl = getrangefromclass(i);


if nargin ~= 1

    paramstrings = {'offset','numlevels','graylimits'};

    for k = 2:2:nargin

        param = lower(varargin{k});
        inputstr = iptcheckstrs(param, paramstrings, mfilename, 'param', k);
        idx = k + 1;  %advance index to the value portion of the input.
        if idx > nargin
            eid = sprintf('images:%s:missingparametervalue', mfilename);
            msg = sprintf('parameter ''%s'' must be followed by a value.', inputstr);
            error(eid,'%s', msg);
        end

        switch (inputstr)

            case 'offset'

                offset = varargin{idx};
                iptcheckinput(offset,{'logical','numeric'},...
                    {'d','nonempty','integer','real'},...
                    mfilename, 'offset', idx);
                % must be row vector 
                if size(offset,2) ~= 1
                    eid = sprintf('images:%s:invalidoffsetsize',mfilename);
                    msg = 'offset must be an n x 1 array.';
                    error(eid,'%s',msg);
                end
                offset = double(offset);

            case 'numlevels'

                nl = varargin{idx};
                iptcheckinput(nl,{'logical','numeric'},...
                    {'real','integer','nonnegative','nonempty','nonsparse'},...
                    mfilename, 'nl', idx);
                if numel(nl) > 1
                    eid = sprintf('images:%s:invalidnumlevels',mfilename);
                    msg = 'nl cannot contain more than one element.';
                    error(eid,'%s',msg);
                elseif islogical(i) && nl ~= 2
                    eid = sprintf('images:%s:invalidnumlevelsforbinary',mfilename);
                    msg = 'nl must be two for a binary image.';
                    error(eid,'%s',msg);
                end
                nl = double(nl);

            case 'graylimits'

                gl = varargin{idx};
                iptcheckinput(gl,{'logical','numeric'},{'vector','real'},...
                    mfilename, 'gl', idx);
                if isempty(gl)
                    gl = [min(i(:)) max(i(:))];
                elseif numel(gl) ~= 2
                    eid = sprintf('images:%s:invalidgraylimitssize',mfilename);
                    msg = 'gl must be a two-element vector.';
                    error(eid,'%s',msg);
                end
                gl = double(gl);
        end
    end
end



十类游程矩阵相关的纹理统计特征：
1、短游程优势：Short Run Emphasis(SRE)
$\mathrm{SRE}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}}\left[\frac{\rho(\mathrm{i}, \mathrm{i} | \theta)}{\mathrm{j}^{2}}\right]}{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, \mathrm{j} | \theta)}$
2、长游程优势：Long Run Emphasis(LRE)
$\mathrm{LRE}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \mathrm{j}^{2} \rho(\mathrm{i}, \mathrm{j} | \theta)}{\sum_{i=1}^{\mathrm{N}} \sum_{j=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, \mathrm{j} | \theta)}$
3、灰度不均匀性：Gray Level Non．Uniformity(GLN)
$\operatorname{GLN}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}}\left[\sum_{\mathrm{j}=1}^{N_{\mathrm{r}}} \rho(\mathrm{i}, j | \theta)\right]^{2}}{\sum_{\mathrm{i}=1}^{\mathrm{N}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, \mathrm{j} | \theta)}$
4、长游程不均匀性：Run Length Non—Uniformity(RLN)
$\mathrm{RLN}=\frac{\sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}}\left[\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \rho(\mathrm{i}, \mathrm{j} | \theta)\right]^{2}}{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, j | \theta)}$
5、游程百分比：Run Percentage(RP)
$R P=\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}} \frac{\rho(i, j | \theta)}{N_{\rho}}$
6、低灰度级游程优势：Low Gray Level Run Emphasis(LGLRE)
$\mathrm{LGLRE}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}}\left[\frac{\rho(\mathrm{i}, j | \theta)}{\mathrm{i}^{2}}\right]}{\sum_{\mathrm{i}=1}^{\mathrm{N}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, j | \theta)}$
7、高灰度级游程优势：High Gray Level Run Emphasis(HGLRE)
$\mathrm{HGLRE}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \mathrm{i}^{2} \rho(\mathrm{i}, j | \theta)}{\sum_{\mathrm{i}=1}^{\mathrm{N}_{\mathrm{g}}} \sum_{\mathrm{j}=1}^{\mathrm{N}_{\mathrm{r}}} \rho(\mathrm{i}, j | \theta)}$
8、短游程低灰度级优势：Short Run Low Gray Level Emphasis(SRLGLE)
$S R L G L E=\frac{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}}\left[\frac{p(i, j | \theta)}{i^{2} j^{2}}\right]}{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}} p(i, j | \theta)}$
9、短游程高灰度级优势：Short Run High Gray Level Emphasis(SRHGLE)
$S R H G L E=\frac{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}}\left[\frac{p(i, j | \theta) i^{2}}{j^{2}}\right]}{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}} p(i, j | \theta)}$
10、长游程低灰度级优势：Long Run Low Gray Level Emphasis(LRLGLE)
$L R L G L E=\frac{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}}\left[\frac{p(i, j | \theta) j^{2}}{i^{2}}\right]}{\sum_{i=1}^{N_{g}} \sum_{j=1}^{N_{r}} p(i, j | \theta)}$
代码：
// An highlighted block

%--------------------------------------------------------------------------
% this program select a roi, qunatize to lower bit level and computing 
% gray level run length matrix and seven texture parameters viz., 
%    1. short run emphasis (sre) 
%    2. long run emphasis(lre)
%    3. gray level non-uniformity (gln)
%    4. run percentage (rp)
%    5. run length non-uniformity (rln)
%    6. low gray level run emphasis (lgre)
%    7. high gray level run emphasis (hgre)
%--------------------------------------------------------------------------
% author: courage
%--------------------------------------------------------------------------
clc, clear, close all
im=imread('file name with path');
figure
imshow(im)
im1=imcrop(im);
im2=im1(1:128,1:128);
im2=double(im2);
[m,n]=size(im2);
% --------- image quantization to 4 bits (16 gray levels)------------------
imax=max(max(im2));
imin=min(min(im2));
newim=im2-imin;
nmax=max(max(newim));
nmin=min(min(newim));
q=round(nmax/16);
[m,n]=size(newim);
quant=0;
for i=1:m
    for j=1:n
        i = newim(i,j);
        for b = 1:16
            if (i>quant)&(i<=quant+q)
                newim(i,j)=b/16;
                quant=quant+q;
            end            
        end
    end
end
newmax=max(max(newim));
newim1=newim/newmax;
newim2=round(newim1*16)+1;
dir=0; 
dist1=1;
if (dir == 1)
    newim2=newim2';
end
mx = max(max(newim2));
mn = min(min(newim2));
gl = (mx-mn)+1;
[p,q] = size(newim2);
n=p*q;
count=1;
c=1;
col=1;
grl(mx,p)=0;
maxcount(p*q)=0;
mc=0;
%---------------------computing gray level run length matrix---------------
for j=1:p
    for k=1:q-dist1
        mc=mc+1;
        g=newim2(j,k);
        f=newim2(j,k+dist1);
        if (g==f)&(g~=0)
            count=count+1;
            c=count;
            col=count;
            maxcount(mc)=count;
        else grl(g,c)=grl(g,c)+1;col=1;
            count=1;
            c=1;
        end
    end
    grl(f,col)=grl(f,col)+1;
    count=1;
    c=1;
end
i=(mx:mn);
m=grl(mn:mx,:);
m1=m';
maxrun=max(max(maxcount));
s=0;
g(gl)=0;
r(q)=0;
for u=1:gl
    for v=1:q
        g(u)=g(u)+m(u,v);
        s=s+m(u,v);
    end
end
for u1=1:q
    for v1=1:gl
        r(u1)=r(u1)+m1(u1,v1);
    end
end
[dim,dim1]=size(g);
sre=0; lre=0; gln=0; rln=0; rp=0; lgre=0; hgre=0;
for h1= 1:maxrun
    sre=sre+(r(h1)/(h1*h1));
    lre=lre+(r(h1)*(h1*h1));
    rln=rln+(r(h1)*r(h1));
    rp=rp+r(h1);
end
sre1=sre/s;
lre1=lre/s;
rln1=rln/s;
rp1=rp/n;
for h2=1:gl
    gln=(gln+g(h2)^2);
    lgre=lgre+(g(h2)/(h2*h2));
    hgre=hgre+(h2*h2)*g(h2);
end
gln1=gln/s;
lgre1=lgre/s;
hgre1=hgre/s;
clc
% ---------------------------display the parameters------------------------
disp(sprintf('%6.4f',sre1))
disp(sprintf('%6.4f',lre1))
disp(sprintf('%6.4f',gln1))
disp(sprintf('%6.4f',rp1))
disp(sprintf('%6.4f',rln1))
disp(sprintf('%6.4f',lgre1))
disp(sprintf('%6.4f',hgre1))

一.图像的边缘检测
图像的边缘是指其周围像素灰度急剧变化的那些像素的集合，它是图像最基本的特征。边缘存在于目标、背景、和区域之间，所以它是图像分割所依赖的最重要的依据。由于边缘是位置的标志，对灰度的变化不敏感，因此边缘也是图像匹配的重要特征。


二.角点特征检测
角点所在的领域通常也是图像中稳定的，信息丰富的区域，这些领域看能具有某些特性如旋转不变性，尺度不变性，仿射不变性和光照亮度不变性。基于图像边缘的检测方法和基于图像灰度的检测方法，前者往往需要对图像边缘进行编码，这在很大程度上依赖于图像的分割和边缘提取，具有较大的计算量，且一旦待检测目标局部发生变化，很有可能导致操作失败。且具有准确性，定位性，稳定性，实时性，鲁棒性。


三.SURF特征提取检测
SURF将DoH中的高斯二阶微分模板进行了近似简化，使得模板对图像的滤波只需要进行几个简单的加减法运算，并且，这种预算与滤波模板的尺寸无关，从而极大地提高了尺度不变特征的检测速度。

我喜欢的智慧石资源：基于Gabor特征提取和人工智能神经网络的人脸检测matlab代码




AMH8周年庆