halvesize matlab,《数字图像处理原理与实践（MATLAB版）》一书之代码Part10（完结篇）...

本文系《数字图像处理原理与实践(MATLAB版)》一书之代码系列的Part10(完结篇)，辑录该书第416至第430页之代码，供有需要读者下载研究使用。至此全书代码发布已经完毕，希望这些源码能够对有需要的读者有所帮助。代码执行结果请参见原书配图，建议下载代码前阅读下文：

关于《数字图像处理原理与实践(MATLAB版)》一书代码发布的说明

http://blog.csdn.net/baimafujinji/article/details/40987807

附全书代码链接

http://blog.csdn.net/baimafujinji/article/details/44856089

http://blog.csdn.net/baimafujinji/article/details/44815599

http://blog.csdn.net/baimafujinji/article/details/42110831

http://blog.csdn.net/baimafujinji/article/details/42080713

http://blog.csdn.net/baimafujinji/article/details/41971057

http://blog.csdn.net/baimafujinji/article/details/41970367

http://blog.csdn.net/baimafujinji/article/details/41926571

http://blog.csdn.net/baimafujinji/article/details/41598455

http://blog.csdn.net/baimafujinji/article/details/41146381

http://blog.csdn.net/baimafujinji/article/details/41117641

首先给出的是原书P429所列之程序源码(注意这是一个测试程序)，详情请参见原书之相关内容。

===========================================================

I1=imreadbw('box.png');

I2=imreadbw('box_in_scene.png');

I1=I1-min(I1(:));

I1=I1/max(I1(:));

I2=I2-min(I2(:));

I2=I2/max(I2(:));

%SIFT特征检测

[frames1,descr1,gss1,dogss1] = do_sift(I1, 'Verbosity', 1, ...

'NumOctaves', 4, 'Threshold', 0.1/3/2);

[frames2,descr2,gss2,dogss2] = do_sift(I2, 'Verbosity', 1, ...

'NumOctaves', 4, 'Threshold', 0.1/3/2);

%计算匹配

descr1 = descr1';

descr2 = descr2';

do_match(I1, descr1, frames1',I2, descr2, frames2' ) ;

===========================================================

鉴于上述代码调用的函数所占篇幅过长，以下部分之代码(被调用函数)并未辑录于书中，在此详细列出，以供有需要的读者参阅。特别说明，下面的程序源码由Yantao Zheng和 NOEMIE PHULPIN编写，笔者略有修改(主要是剔除了算法执行计时的部分)，应用时请尊重原作者权利。

===========================================================

function [frames,descriptors,scalespace,difofg]=do_sift(I,varargin)

warning off all;

[M,N,C] = size(I) ;

% Lowe's choices

S=3 ;

omin= 0 ;

%O=floor(log2(min(M,N)))-omin-4 ; % Up to 16x16 images

O = 4;

sigma0=1.6*2^(1/S) ;

sigman=0.5 ;

thresh = 0.2 / S / 2 ; % 0.04 / S / 2 ;

r = 18 ;

NBP = 4 ;

NBO = 8 ;

magnif = 3.0 ;

% Parese input

compute_descriptor = 0 ;

discard_boundary_points = 1 ;

verb = 0 ;

% Arguments sanity check

if C > 1

error('I should be a grayscale image') ;

end

frames = [] ;

descriptors = [] ;

scalespace = do_gaussian(I,sigman,O,S,omin,-1,S+1,sigma0) ;

difofg = do_diffofg(scalespace) ;

for o=1:scalespace.O

% Local maxima of the DOG octave

oframes1 = do_localmax(  difofg.octave{o}, 0.8*thresh, difofg.smin  ) ;

oframes = [oframes1 , do_localmax( - difofg.octave{o}, 0.8*thresh, difofg.smin)] ;

if size(oframes, 2) == 0

continue;

end

% Remove points too close to the boundary

rad = magnif * scalespace.sigma0 * 2.^(oframes(3,:)/scalespace.S) * NBP / 2 ;

sel= ...

oframes(1,:)-rad >= 1                     & ...

oframes(1,:)+rad <= size(scalespace.octave{o},2) & ...

oframes(2,:)-rad >= 1                     & ...

oframes(2,:)+rad <= size(scalespace.octave{o},1)      ;

oframes=oframes(:,sel) ;

% Refine the location, threshold strength and remove points on edges

oframes = do_extrefine(...

oframes, ...

difofg.octave{o}, ...

difofg.smin, ...

thresh, ...

r) ;

% Compute the orientations

oframes = do_orientation(...

oframes, ...

scalespace.octave{o}, ...

scalespace.S, ...

scalespace.smin, ...

scalespace.sigma0 ) ;

% Store frames

x     = 2^(o-1+scalespace.omin) * oframes(1,:) ;

y     = 2^(o-1+scalespace.omin) * oframes(2,:) ;

sigma = 2^(o-1+scalespace.omin) * scalespace.sigma0 * 2.^(oframes(3,:)/scalespace.S) ;

frames = [frames, [x(:)' ; y(:)' ; sigma(:)' ; oframes(4,:)] ] ;

sh = do_descriptor(scalespace.octave{o}, ...

oframes, ...

scalespace.sigma0, ...

scalespace.S, ...

scalespace.smin, ...

'Magnif', magnif, ...

'NumSpatialBins', NBP, ...

'NumOrientBins', NBO) ;

descriptors = [descriptors, sh] ;

end

===========================================================

function num = do_match(im1, des1, loc1, im2, des2, loc2)

% This function matchings the SIFT keys from two iamges

% distRatio: Only keep matches in which the ratio of vector angles from the

%   nearest to second nearest neighbor is less than distRatio.

% Postprocessing: check each matching point, eliminate false matches by voting from

% neighbouring area

% Author: Yantao Zheng. Nov 2006.  For Project of CS5240

distRatio = 0.75;

matched_points_img1 =[];

% For each descriptor in the first image, select its match to second image.

des2t = des2';                          % Precompute matrix transpose

for i = 1 : size(des1,1)

dotprods = des1(i,:) * des2t;        % Computes vector of dot products

[vals,indx] = sort(acos(dotprods));  % Take inverse cosine and sort results

% Check if nearest neighbor has angle less than distRatio times 2nd.

if (vals(1) < distRatio * vals(2))

match(i) = indx(1);

matched_points_img1 = [matched_points_img1; i];

else

match(i) = 0;

end

end

%---------------------------------------------------------

% check each matching point, eliminate false matches by voting from

% neighbouring area

%------

dis_thres = 0.3;

orien_thres = 0.3;

num = sum(match > 0);

final_match = zeros(1,length(match));  % if == 1, the corresponding match(i) is accepted. if ==0, the corresponding match(i) is rejected

dis_img_1= zeros(num ,num); % each row contains the spatial distance vector that descript the correctness of the point

dis_img_2= zeros(num ,num);

orien_diff_img_1= zeros(num ,num); % each row contains the orientation difference vector that descript the correctness of the point

orien_diff_img_2= zeros(num ,num);

for k = 1: num

dis_img_1(k, k)  = 0;

dis_img_2(k, k)  = 0; % the distance between the point and distance is 0

orien_diff_img_1(k,k) = 0;

orien_diff_img_2(k,k) = 0;

for j = k+ 1: num

dis_img_1(k, j) = sqrt( (loc1(matched_points_img1(k), 1) - loc1(matched_points_img1(j), 1))^2 ...

+ (loc1(matched_points_img1(k), 2) - loc1(matched_points_img1(j), 2))^2 );

dis_img_1(j, k) =  dis_img_1(k, j); % dis_img_1 is a symmetric matrix

% compute the corresponding distances of the matching points in

% image 2

dis_img_2(k, j) = sqrt( ((loc2(match(matched_points_img1(k)), 1) - loc2(match(matched_points_img1(j)), 1))^2 + (loc2(match(matched_points_img1(k)), 2) - loc2(match(matched_points_img1(j)), 2))^2 ));

dis_img_2(j, k) =  dis_img_2(k, j); % dis_img_1 is a symmetric matrix

orien_diff_img_1(k, j) =  loc1(matched_points_img1(k), 4) - loc1(matched_points_img1(j), 4);

orien_diff_img_1(j, k) =  orien_diff_img_1(k, j); % dis_img_1 is a symmetric matrix

orien_diff_img_2(k, j) =  loc2(match(matched_points_img1(k)), 4) - loc2(match(matched_points_img1(j)), 4);

orien_diff_img_2(j, k) =  orien_diff_img_2(k, j); % dis_img_1 is a symmetric matrix

end

end

% normalize the distance and orein_diff vector

for ii =1: num

dis_img_1(ii, :) = dis_img_1(ii, :) ./ norm(dis_img_1(ii, :) );

dis_img_2(ii, :) = dis_img_2(ii, :) ./ norm(dis_img_2(ii, :) );

orien_diff_img_1(ii,:) = orien_diff_img_1(ii,:) ./( eps + norm(orien_diff_img_1(ii,:)));

orien_diff_img_2(ii,:) = orien_diff_img_2(ii,:) ./( eps +norm(orien_diff_img_2(ii,:)));

end

for m = 1: num

dis_coherence = dot(dis_img_1(m,:), dis_img_2(m,:));

orein_coh = dot(orien_diff_img_1(m,:), orien_diff_img_2(m,:));

if dis_coherence > dis_thres && (orein_coh > orien_thres  || ( sum(orien_diff_img_1(m,:)>0) == 0 && sum(orien_diff_img_2(m,:)>0) == 0 ) )

% if the orientations are the same, then orein_coh will be 0, so

% cope with this with another condition

final_match(matched_points_img1(m)) = 1;

end

end

num = sum(match > 0);

fprintf('>>>>RESULT: Found %d matches.\n', num);

% Create a new image showing the two images side by side.

im3 = appendimages(im1,im2);

% Show a figure with lines joining the accepted matches.

figure('Position', [100 100 size(im3,2) size(im3,1)]);

colormap('gray');

imagesc(im3);

hold on;

cols1 = size(im1,2);

for i = 1: size(des1,1)

if (match(i) > 0)

line([loc1(i,1) loc2(match(i),1)+cols1], ...

[loc1(i,2) loc2(match(i),2)], 'Color', 'c');

end

end

hold off;

===========================================================

function descriptors = do_descriptor(octave, ... % gaussian scale space of an octave

oframes, ...  % frames containing keypoint coordinates and scale, and orientation

sigma0, ...   % base sigma value

S, ...        % level of scales in the octave

smin, ...

varargin)

for k=1:2:length(varargin)

switch lower(varargin{k})

case 'magnif'

magnif = varargin{k+1} ;

case 'numspatialbins'

NBP = varargin{k+1} ;

case  'numorientbins'

NBO = varargin{k+1} ;

otherwise

error(['Unknown parameter ' varargin{k} '.']) ;

end

end

num_spacialBins = NBP;

num_orientBins = NBO;

key_num = size(oframes, 2);

% compute the image gradients

[M, N, s_num] = size(octave); % M is the height of image, N is the width of image; num_level is the number of scale level of the octave

descriptors = [];

magnitudes = zeros(M, N, s_num);

angles = zeros(M, N, s_num);

% compute image gradients

for si = 1: s_num

img = octave(:,:,si);

dx_filter = [-0.5 0 0.5];

dy_filter = dx_filter';

gradient_x = imfilter(img, dx_filter);

gradient_y = imfilter(img, dy_filter);

magnitudes(:,:,si) =sqrt( gradient_x.^2 + gradient_y.^2);

%     if sum( gradient_x == 0) > 0

%         fprintf('00');

%     end

angles(:,:,si) = mod(atan(gradient_y ./ (eps + gradient_x)) + 2*pi, 2*pi);

end

x = oframes(1,:);

y = oframes(2,:);

s = oframes(3,:);

% round off

x_round = floor(oframes(1,:) + 0.5);

y_round = floor(oframes(2,:) + 0.5);

scales =  floor(oframes(3,:) + 0.5) - smin;

for p = 1: key_num

s = scales(p);

xp= x_round(p);

yp= y_round(p);

theta0 = oframes(4,p);

sinth0 = sin(theta0) ;

costh0 = cos(theta0) ;

sigma = sigma0 * 2^(double (s / S)) ;

SBP = magnif * sigma;

%W =  floor( sqrt(2.0) * SBP * (NBP + 1) / 2.0 + 0.5);

W =   floor( 0.8 * SBP * (NBP + 1) / 2.0 + 0.5);

descriptor = zeros(NBP, NBP, NBO);

% within the big square, select the pixels with the circle and put into

% the histogram. no need to do rotation which is very expensive

for dxi = max(-W, 1-xp): min(W, N -2 - xp)

for dyi = max(-W, 1-yp) : min(+W, M-2-yp)

mag = magnitudes(yp + dyi, xp + dxi, s); % the gradient magnitude at current point(yp + dyi, xp + dxi)

angle = angles(yp + dyi, xp + dxi, s) ;  % the gradient angle at current point(yp + dyi, xp + dxi)

angle = mod(-angle + theta0, 2*pi);      % adjust the angle with the major orientation of the keypoint and mod it with 2*pi

dx = double(xp + dxi - x(p));            % x(p) is the exact keypoint location (floating number). dx is the relative location of the current pixel with respect to the keypoint

dy = double(yp + dyi - y(p));            % dy is the relative location of the current pixel with respect to the keypoint

nx = ( costh0 * dx + sinth0 * dy) / SBP ; % nx is the normalized location after rotation (dx, dy) with the major orientation angle. this tells which x-axis spatial bin the pixel falls in

ny = (-sinth0 * dx + costh0 * dy) / SBP ;

nt = NBO * angle / (2* pi) ;

wsigma = NBP/2 ;

wincoef =  exp(-(nx*nx + ny*ny)/(2.0 * wsigma * wsigma)) ;

binx = floor( nx - 0.5 ) ;

biny = floor( ny - 0.5 ) ;

bint = floor( nt );

rbinx = nx - (binx+0.5) ;

rbiny = ny - (biny+0.5) ;

rbint = nt - bint ;

for(dbinx = 0:1)

for(dbiny = 0:1)

for(dbint = 0:1)

% if condition limits the samples within the square

% width W. binx+dbinx is the rotated x-coordinate.

% therefore the sampling square is effectively a

% rotated one

if( binx+dbinx >= -(NBP/2) && ...

binx+dbinx

biny+dbiny >= -(NBP/2) && ...

biny+dbiny

weight = wincoef * 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

end

descriptor = reshape(descriptor, 1, NBP * NBP * NBO);

descriptor = descriptor ./ norm(descriptor);

%Truncate at 0.2

indx = find(descriptor > 0.2);

descriptor(indx) = 0.2;

descriptor = descriptor ./ norm(descriptor);

descriptors = [descriptors, descriptor'];

end

===========================================================

function J=extrefine(oframes,octave,smin,thres,r)

%% file:        extrefine.m

% author:      Noemie Phulpin

% description: Refine the location, threshold strength and remove points on edges

[M,N,S]=size(octave);

[L,K]=size(oframes);

comp=1;

for p = 1:K

b=zeros(1,3) ;

A=oframes(:,p);

x=A(1)+1;

y=A(2)+1;

s=A(3)+1-smin;

%Local maxima extracted from the DOG have coordinates 1<=x<=N-2, 1<=y<=M-2

% and 1<=s-mins<=S-2. This is also the range of the points that we can refine.

if(x < 2 || x > N-1 || y < 2 || y > M-1 || s < 2 || s > S-1)

continue ;

end

val=octave(y,x,s);

Dx=0;Dy=0;Ds=0;Dxx=0;Dyy=0;Dss=0;Dxy=0;Dxs=0;Dys=0 ;

dx = 0 ;

dy = 0 ;

for iter = 1:5

A = zeros(3,3) ;

x = x + dx ;

y = y + dy ;

if (x < 2 || x > N-1 || y < 2 || y > M-1 )  break ; end

% Compute the gradient.

Dx = 0.5 * (octave(y,x+1,s) - octave(y,x-1,s));

Dy = 0.5 * (octave(y+1,x,s) - octave(y-1,x,s)) ;

Ds = 0.5 * (octave(y,x,s+1) - octave(y,x,s-1)) ;

% Compute the Hessian.

Dxx = (octave(y,x+1,s) + octave(y,x-1,s) - 2.0 * octave(y,x,s)) ;

Dyy = (octave(y+1,x,s) + octave(y-1,x,s) - 2.0 * octave(y,x,s)) ;

Dss = (octave(y,x,s+1) + octave(y,x,s-1) - 2.0 * octave(y,x,s)) ;

Dys = 0.25 * ( octave(y+1,x,s+1) + octave(y-1,x,s-1) - octave(y-1,x,s+1) - octave(y+1,x,s-1) ) ;

Dxy = 0.25 * ( octave(y+1,x+1,s) + octave(y-1,x-1,s) - octave(y-1,x+1,s) - octave(y+1,x-1,s) ) ;

Dxs = 0.25 * ( octave(y,x+1,s+1) + octave(y,x-1,s-1) - octave(y,x-1,s+1) - octave(y,x+1,s-1) ) ;

% Solve linear system.

A(1,1) = Dxx ;

A(2,2) = Dyy ;

A(3,3) = Dss ;

A(1,2) = Dxy ;

A(1,3) = Dxs ;

A(2,3) = Dys ;

A(2,1) = Dxy ;

A(3,1) = Dxs ;

A(3,2) = Dys ;

b(1) = - Dx ;

b(2) = - Dy ;

b(3) = - Ds ;

c=b*inv(A);

% If the translation of the keypoint is big, move the keypoint and re-iterate the computation. Otherwise we are all set.

if (c(1) >  0.6 && x < N-2 )

if (c(1) < -0.6 && x > 1)

dx=0;

else

dx=1;

end

else

if (c(1) < -0.6 && x > 1)

dx=-1;

else

dx=0;

end

end

if (c(2) >  0.6 && y < N-2 )

if (c(2) < -0.6 && y > 1)

dy=0;

else

dy=1;

end

else

if (c(2) < -0.6 && y > 1)

dy=-1;

else

dy=0;

end

end

if( dx == 0 && dy == 0 ) break ; end

end

%we keep the value only of it verify the conditions

val = val + 0.5 * (Dx * c(1) + Dy * c(2) + Ds * c(3)) ;

score = (Dxx+Dyy)*(Dxx+Dyy) / (Dxx*Dyy - Dxy*Dxy) ;

xn = x + c(1) ;

yn = y + c(2) ;

sn = s + c(3) ;

if (abs(val) > thres) && ...

(score < (r+1)*(r+1)/r) && ...

(score >= 0) && ...

(abs(c(1)) < 1.5) && ...

(abs(c(2)) < 1.5) && ...

(abs(c(3)) < 1.5) && ...

(xn >= 0) && ...

(xn <= M-1) && ...

(yn >= 0) && ...

(yn <= N-1) && ...

(sn >= 0) && ...

(sn <= S-1)

J(1,comp)=xn-1;

J(2,comp)=yn-1;

J(3,comp)=sn-1+smin;

comp=comp+1;

end

end

return

===========================================================

function L = do_gaussian(I,sigman,O,S,omin,smin,smax,sigma0)

%% file:        do_gaussian.m

% author:      Noemie Phulpin

% description: gaussian scale space of image I

%%

if(nargin<7)

sigma0=1.6*k;

end

if omin<0

for o=1:-omin

I=doubleSize(I);

end

elseif omin>0

for o=1:-omin

I=halveSize(I);

end

end

[M,N] = size(I);                      %size of image

k = 2^(1/S);                          %scale space multiplicative step k

sigma0=1.6*k;                         %definition by Lowe

dsigma0 = sigma0*sqrt(1-1/k^2);       %scale step factor

sigmaN=0.5;                           %nominal smoothing of the image

so=-smin+1;                           %index offset

%scale space structure

L.O = O;

L.S = S;

L.sigma0 = sigma0;

L.omin = omin;

L.smin = smin;

L.smax = smax;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%First Octave

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%initilize the octave with S sub-levels

L.octave{1} = zeros(M,N,smax-smin+1);

%initilize the first sub-level

sig=sqrt( (sigma0*k^smin)^2 - (sigmaN/2^omin)^2 );

%b=smooth2(I,sig) ;

%[N1,M1]=size(b)

%b(1:4,1:4)

%c=imsmooth(I,sig) ;

%[N2,M2]=size(c)

%c(1:4,1:4)

L.octave{1}(:,:,1) = smooth(I,sig);

%other sub-levels

for s=smin+1:smax

dsigma = k^s * dsigma0;

L.octave{1}(:,:,s+so) = smooth( squeeze(L.octave{1}(:,:,s-1+so)) ,dsigma);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Folowing Octaves

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%convert all octaves

for o=2:O

sbest = min(smin+S,smax);

TMP = halvesize( squeeze(L.octave{o-1}(:,:,sbest+so)) );

sigma_next = sigma0*k^smin;

sigma_prev = sigma0*k^(sbest-S);

if (sigma_next>sigma_prev)

sig=sqrt(sigma_next^2-sigma_prev^2);

TMP= smooth( TMP,sig);

end

[M,N] = size(TMP);

L.octave{o} = zeros(M,N,smax-smin+1);

L.octave{o}(:,:,1) = TMP;

%other sub-levels

for s=smin+1:smax

dsigma = k^s * dsigma0;

L.octave{o}(:,:,s+so) = smooth( squeeze(L.octave{o}(:,:,s-1+so)) ,dsigma);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Auxiliary functions

function J = halvesize(I)

J=I(1:2:end,1:2:end);

function J = doubleSize(I)

[M,N]=size(I) ;

J = zeros(2*M,2*N) ;

J(1:2:end,1:2:end) = I ;

J(2:2:end-1,2:2:end-1) = ...

0.25*I(1:end-1,1:end-1) + ...

0.25*I(2:end,1:end-1) + ...

0.25*I(1:end-1,2:end) + ...

0.25*I(2:end,2:end) ;

J(2:2:end-1,1:2:end) = ...

0.5*I(1:end-1,:) + ...

0.5*I(2:end,:) ;

J(1:2:end,2:2:end-1) = ...

0.5*I(:,1:end-1) + ...

0.5*I(:,2:end) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

===========================================================

function oframes = do_orientation(...

oframes, ...

octave, ...

S, ...

smin, ...

sigma0 )

% this function computes the major orientation of the keypoint (oframes).

% Note that there can be multiple major orientations. In that case, the

% SIFT keys will be duplicated for each major orientation

% Author: Yantao Zheng. Nov 2006.  For Project of CS5240

frames = [];

win_factor = 1.5 ;

NBINS = 36;

histo = zeros(1, NBINS);

[M, N, s_num] = size(octave); % M is the height of image, N is the width of image; num_level is the number of scale level of the octave

key_num = size(oframes, 2);

magnitudes = zeros(M, N, s_num);

angles = zeros(M, N, s_num);

% compute image gradients

for si = 1: s_num

img = octave(:,:,si);

dx_filter = [-0.5 0 0.5];

dy_filter = dx_filter';

gradient_x = imfilter(img, dx_filter);

gradient_y = imfilter(img, dy_filter);

magnitudes(:,:,si) =sqrt( gradient_x.^2 + gradient_y.^2);

angles(:,:,si) = mod(atan(gradient_y ./ (eps + gradient_x)) + 2*pi, 2*pi);

end

% round off the cooridnates and

x = oframes(1,:);

y = oframes(2,:) ;

s = oframes(3,:);

x_round = floor(oframes(1,:) + 0.5);

y_round = floor(oframes(2,:) + 0.5);

scales = floor(oframes(3,:) + 0.5) - smin;

for p=1:key_num

s = scales(p);

xp= x_round(p);

yp= y_round(p);

sigmaw = win_factor * sigma0 * 2^(double (s / S)) ;

W = floor(3.0* sigmaw);

for xs = xp - max(W, xp-1): min((N - 2), xp + W)

for ys = yp - max(W, yp-1) : min((M-2), yp + W)

dx = (xs - x(p));

dy = (ys - y(p));

if dx^2 + dy^2 <= W^2 % the points are within the circle

wincoef = exp(-(dx^2 + dy^2)/(2*sigmaw^2));

bin = round( NBINS *  angles(ys, xs, s+ 1)/(2*pi) + 0.5);

histo(bin) = histo(bin) + wincoef * magnitudes(ys, xs, s+ 1);

end

end

end

theta_max = max(histo);

theta_indx = find(histo> 0.8 * theta_max);

for i = 1: size(theta_indx, 2)

theta = 2*pi * theta_indx(i) / NBINS;

frames = [frames, [x(p) y(p) s theta]'];

end

end

oframes = frames;

% for each keypoint

===========================================================

function J = do_localmax(octave, thresh,smin)

%% file:       do_localmax.m

% author:      Noemie Phulpin

% description: returns the indexes of the local maximizers of the octave.

%%

[N,M,S]=size(octave);

nb=1;

k=0.0002;

%for each point of this scale space

% we look for extrama bigger than thresh

J = [];

for s=2:S-1

for j=20:M-20

for i=20:N-20

a=octave(i,j,s);

%             V = [ thresh octave(i-1,j-1,s-1) octave(i-1,j,s-1) octave(i-1,j+1,s-1) octave(i,j-1,s-1) octave(i,j+1,s-1) octave(i+1,j-1,s-1) octave(i+1,j,s-1) octave(i+1,j+1,s-1) octave(i-1,j-1,s) octave(i-1,j,s) octave(i-1,j+1,s) octave(i,j-1,s) octave(i,j+1,s) octave(i+1,j-1,s) octave(i+1,j,s) octave(i+1,j+1,s) octave(i-1,j-1,s+1) octave(i-1,j,s+1) octave(i-1,j+1,s+1) octave(i,j-1,s+1) octave(i,j+1,s+1) octave(i+1,j-1,s+1) octave(i+1,j,s+1) octave(i+1,j+1,s+1)];

%             maximum=max(V);

if a>thresh+k ...

&& a>octave(i-1,j-1,s-1)+k && a>octave(i-1,j,s-1)+k && a>octave(i-1,j+1,s-1)+k ...

&& a>octave(i,j-1,s-1)+k && a>octave(i,j+1,s-1)+k && a>octave(i+1,j-1,s-1)+k ...

&& a>octave(i+1,j,s-1)+k && a>octave(i+1,j+1,s-1)+k && a>octave(i-1,j-1,s)+k ...

&& a>octave(i-1,j,s)+k && a>octave(i-1,j+1,s)+k && a>octave(i,j-1,s)+k ...

&& a>octave(i,j+1,s)+k && a>octave(i+1,j-1,s)+k && a>octave(i+1,j,s)+k ...

&& a>octave(i+1,j+1,s)+k && a>octave(i-1,j-1,s+1)+k && a>octave(i-1,j,s+1)+k ...

&& a>octave(i-1,j+1,s+1)+k && a>octave(i,j-1,s+1)+k && a>octave(i,j+1,s+1)+k ...

&& a>octave(i+1,j-1,s+1)+k && a>octave(i+1,j,s+1)+k && a>octave(i+1,j+1,s+1)+k

%if (a-maximum>0.00004)

%if (a-maximum>0.001)

J(1,nb)=j-1;

J(2,nb)=i-1;

J(3,nb)=s+smin-1;

nb=nb+1;

end

end

end

end

===========================================================

function D = do_diffofg(L)

%% file:        do_diffofg.m

% author:      Noemie Phulpin

% description: substraction of consecutive levels of the scale space SS.

%%

D.smin = L.smin;

D.smax = L.smax-1;

D.omin =L.omin;

D.O = L.O;

D.S = L.S;

D.sigma0 = L.sigma0;

for o=1:D.O

[M,N,S] = size(L.octave{o});

D.octave{o} = zeros(M,N,S-1);

for s=1:S-1

D.octave{o}(:,:,s) = L.octave{o}(:,:,s+1) -  L.octave{o}(:,:,s);

end;

end;

===========================================================

function im = appendimages(image1, image2)

% Select the image with the fewest rows and fill in enough empty rows

%   to make it the same height as the other image.

rows1 = size(image1,1);

rows2 = size(image2,1);

if (rows1 < rows2)

image1(rows2,1) = 0;

else

image2(rows1,1) = 0;

end

% Now append both images side-by-side.

im = [image1 image2];

===========================================================

function I = imreadbw(file)

%IMREADBW Summary of this function goes here

%   Detailed explanation goes here

I=im2double(imread(file)) ;

if(size(I,3) > 1)

I = rgb2gray( I ) ;

end

end

===========================================================

function J = smooth(I,s)

%filter

h=fspecial('gaussian',ceil(4*s),s);

%convolution

J=imfilter(I,h);

return;

全文完。