凯鲁嘎吉 用书写铭记日常，最迷人的不在远方

# Gaussian field consensus论文解读及MATLAB实现

## 一、Introduction

An image pair and its putative correspondences. Blue and red lines represent inliers and outliers respectively.

## 五、Code

results.m

function criterion=results()
%  GFC;
% 读取图像及对应关系数据
% church
% ImgName1 = './TestData/church1.jpg' ;
% ImgName2 = './TestData/church2.jpg' ;
% eye
ImgName1 = './TestData/14_a.jpg' ;
ImgName2 = './TestData/14_b.jpg' ;

t0=cputime;
[precision, recall, accuracy, f1]=demo_GFC(I1, I2, X, Y, CorrectIndex);
run_time=cputime-t0;
criterion=[precision; recall; accuracy; f1; run_time];


demo_GFC.m

function [precision, recall, accuracy, f1]=demo_GFC(I1, I2, X, Y, CorrectIndex)
label=CorrectIndex;
delta =0.5;
[index] = GFC_match(X,Y,delta);
[precision, recall, accuracy, f1]  = evaluatePR(label, index, size(X,1));

% Plot results
%原始结果
% N=size(X);
% temp=1:N;
% temp=temp';
% plot_matches(I1, I2, X, Y, temp, label);
%实验结果
plot_matches(I1, I2, X, Y, index, label);


function A = adjacency(DATA, TYPE, PARAM, DISTANCEFUNCTION);

% Compute the adjacency graph of the data set DATA
%
% A = adjacency(DATA, TYPE, PARAM, DISTANCEFUNCTION);
%
% DATA - NxK matrix. Data points are rows.
% TYPE - string 'nn' or string 'epsballs'.
% PARAM - integer if TYPE='nn', real number if TYPE='epsballs'.
% DISTANCEFUNCTION - function mapping a (DxM) and a (D x N) matrix
%                    to an M x N distance matrix (D:dimensionality)
% Returns: A, sparse symmetric NxN matrix of distances between the
%
% Example:
%
%   A contains the adjacency matrix for the data
%   set X. For each point, the distances to 6 adjacent points are
%   stored. N
%
% Note: the adjacency relation is symmetrized, i.e. if
% point a is adjacent to point b, then point b is also considered to be
%
%
% Author:
%
% Mikhail Belkin
% misha@math.uchicago.edu
%
% Modified by: Vikas Sindhwani
% June 2004

if (nargin < 3) | (strcmp(TYPE,'nn') & strcmp(TYPE,'epsballs')) | ~isreal(PARAM)

disp(sprintf('ERROR: Too few arguments given or incorrect arguments.\n'));
disp(sprintf('USAGE:\n A = laplacian(DATA, TYPE, PARAM)'));
disp(sprintf('DATA - the data matrix. Data points are rows.'));
disp(sprintf('Nearest neigbors: TYPE =''nn''    PARAM = number of nearest neigbors'));
disp(sprintf('Epsilon balls: TYPE =''epsballs''    PARAM = redius of the ball\n'));
return;
end

n = size(DATA,1);
%disp (sprintf ('DATA: %d points in %d dimensional space.',n,size (DATA,2)));

switch TYPE
case {'nn'}
% disp(sprintf('Creating the adjacency matrix. Nearest neighbors, N=%d.', PARAM));
case{'eps', 'epsballs'}
%disp(sprintf('Creating the adjacency matrix. Epsilon balls, eps=%f.', PARAM));
end;

A = sparse(n,n);
step = 100;

if (strcmp(TYPE,'nn'))
for i1=1:step:n
i2 = i1+step-1;
if (i2> n)
i2=n;
end;
XX= DATA(i1:i2,:);
dt = feval(DISTANCEFUNCTION, XX',DATA');
[Z,I] = sort ( dt,2);

for i=i1:i2
if ( mod(i, 500) ==0)
%disp(sprintf('%d points processed.', i));
end;
for j=2:PARAM+1
A(i,I(i-i1+1,j))= Z(i-i1+1,j);
A(I(i-i1+1,j),i)= Z(i-i1+1,j);
end;
end

end;

% epsilon balls
else
for i1=1:step:n
i2 = i1+step-1;
if (i2> n)
i2=n;
end;

XX= DATA(i1:i2,:);
dt = feval(DISTANCEFUNCTION, XX',DATA');
[Z,I] = sort ( dt,2 );

for i=i1:i2
%  if ( mod(i, 500) ==0) disp(sprintf('%d points processed.', i)); end;
j=2;
while ( (Z(i-i1+1,j) < PARAM))
j = j+1;
jj = I(i-i1+1,j);
A(i,jj)= Z(i-i1+1,j);
A(jj,i)= Z(i-i1+1,j);
end;
end
end;
end;


con_K.m

function K=con_K(x,y,beta)
if nargin<3
error('Error! Not enough input parameters.');
end
ks=-2 * beta^2;
[n, d]=size(x);
[m, d]=size(y);
K=repmat(x,[1 1 m])-permute(repmat(y,[1 1 n]),[3 2 1]);
K=squeeze(sum(K.^2,2));
K=K/ks;
K=exp(K);


costfun_GFC.m

function [E, G] = costfun_GFC(param, X, Y,  U, beta,lambda,sigma2)
[N, D] = size(X);
M = size(U, 2);
Alpha = reshape(param, [M D]);
options=ml_options('Kernel','rbf', 'KernelParam', beta,'NN',5);
options.GraphWeights= 'heat';
options.GraphWeightParam=sqrt(sigma2);
L=laplacian(X,'nn',options);
E=lambda * trace(Alpha'*U'*L*U*Alpha);
V = Y-(X+ U*Alpha);
a = -2 / N / (2*sigma2)^(D/2);
F = exp(-sum(V.^2, 2) / (2*sigma2));
E = E + a * sum(F);
G =  -a * U' * ( V .* repmat(F, [1 D]) / sigma2 ) + 2*lambda * U'* L * U *Alpha;


euclidean.m

function d = euclidean(a,b,df)
% EUCLIDEAN - computes Euclidean distance matrix
%
% E = euclidean(A,B)
%
%    A - (DxM) matrix
%    B - (DxN) matrix
%    df = 1, force diagonals to be zero; 0 (default), do not force
%
% Returns:
%    E - (MxN) Euclidean distances between vectors in A and B
%
%
% Description :
%    This fully vectorized (VERY FAST!) m-file computes the
%    Euclidean distance between two vectors by:
%
%                 ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
%
% Example :
%    A = rand(400,100); B = rand(400,200);
%    d = distance(A,B);

% Author   : Roland Bunschoten
%            University of Amsterdam
%            Intelligent Autonomous Systems (IAS) group
%            Kruislaan 403  1098 SJ Amsterdam
%            tel.(+31)20-5257524
%            bunschot@wins.uva.nl
% Last Rev : Wed Oct 20 08:58:08 MET DST 1999
% Tested   : PC Matlab v5.2 and Solaris Matlab v5.3

% Copyright notice: You are free to modify, extend and distribute
%    this code granted that the author of the original code is
%    mentioned as the original author of the code.

% Fixed by JBT (3/18/00) to work for 1-dimensional vectors
% and to warn for imaginary numbers.  Also ensures that
% output is all real, and allows the option of forcing diagonals to
% be zero.

if (nargin < 2)
error('Not enough input arguments');
end

if (nargin < 3)
df = 0;    % by default, do not force 0 on the diagonal
end

if (size(a,1) ~= size(b,1))
error('A and B should be of same dimensionality');
end

if ~(isreal(a)*isreal(b))
disp('Warning: running distance.m with imaginary numbers.  Results may be off.');
end

if (size(a,1) == 1)
a = [a; zeros(1,size(a,2))];
b = [b; zeros(1,size(b,2))];
end

aa=sum(a.*a); bb=sum(b.*b); ab=a'*b;
d = sqrt(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);

% make sure result is all real
d = real(d);

% force 0 on the diagonal?
if (df==1)
d = d.*(1-eye(size(d)));
end


evaluatePR.m

function [precision, recall, accuracy, f1] = evaluatePR(CorrectIndex, Index, siz)
tmp=zeros(1, siz);
tmp(Index) = 1;
tmp(CorrectIndex) = tmp(CorrectIndex)+1;
GFCCorrect = find(tmp == 2);
NumCorrectIndex = length(CorrectIndex);
NumGFCIndex = length(Index);
NumGFCCorrect = length(GFCCorrect);
% corrRate = NumCorrectIndex/siz;
precision = NumGFCCorrect/NumGFCIndex;
TP=NumGFCCorrect;
FP=NumGFCIndex-TP;
recall = NumGFCCorrect/NumCorrectIndex;
FN=NumCorrectIndex-TP;
TN=siz-NumGFCIndex-FN;
accuracy=(TP+TN)/(TP+TN+FP+FN);
f1=(2*precision*recall)/(precision+recall);


GFC.m

function [idt, V , param] = GFC(X, Y, delta)
[N1, D] = size(X);
eta=0.98;
history.x = [ ];
history.fval = [ ];
beta =1.5;
lambda = 3;
sigma2=power(det(X'*X/N1), 1/(2^D));
M=round(sqrt(N1));
x0 = zeros(M*D, 1);
options = optimset( 'display','iter');
options = optimset(options, 'outputfcn',@GFCoutfun);
options = optimset( options, 'LargeScale','off');
options = optimset(options, 'MaxFunEvals', 100);
options = optimset(options, 'MaxIter', 100);
U=get_U(X,beta,M);
param = fminunc(@(x)costfun_GFC(x, X, Y,  U, beta, lambda, sigma2), x0, options);
Alpha = reshape(param, [M D]);
V=X+U*Alpha;
Pb =  exp(-sum((Y-V).^2, 2) / (sigma2)) ;
idt = find(Pb > delta);
function stop = GFCoutfun(x,optimValues,state,varargin)
stop = false;
switch state
case 'init'
case 'iter'
history.fval = [history.fval; optimValues.fval];
history.x = [history.x; reshape(x,1,length(x))];
Alpha = reshape(x, [M D]);
V=X+U*Alpha;
sigma2=sigma2 * eta;
case 'done'
otherwise
end
end
end

function U=get_U(X,beta,M)
tmp_X = unique(X, 'rows');
idx = randperm(size(tmp_X,1));
idx = idx(1:min(M,size(tmp_X,1)));
ctrl_pts=tmp_X(idx,:);
U = con_K(X, ctrl_pts, beta);
end


GFC_match.m

function [idt]=GFC_match(X,Y,delta)
Ni=2;
Xk=X;
k=1;
s=1;
while s
X2=Xk;Y2=Y;
normal.xm=0; normal.ym=0;
normal.xscale=1; normal.yscale=1;
[nX, nY, normal]=norm2s(X2,Y2);
[idt, trans] = GFC(nX,nY,delta);
trans=(trans)*normal.yscale+repmat(normal.ym,size(Y2,1),1);
Xk = trans;
if k==Ni
s=0;
else
k=k+1;
end
end
end


laplacian.m

function L = laplacian(DATA, TYPE, options)

% Calculate the graph laplacian of the adjacency graph of data set DATA.
%
% L = laplacian(DATA, TYPE, PARAM)
%
% DATA - NxK matrix. Data points are rows.
% TYPE - string 'nn' or string 'epsballs'
% options - Data structure containing the following fields
% NN - integer if TYPE='nn' (number of nearest neighbors),
%       or size of 'epsballs'
%
% DISTANCEFUNCTION - distance function used to make the graph
% WEIGHTTYPPE='binary' | 'distance' | 'heat'
% WEIGHTPARAM= width for heat kernel
% NORMALIZE= 0 | 1 whether to return normalized graph laplacian or not
%
% Returns: L, sparse symmetric NxN matrix
%
% Author:
%
% Mikhail Belkin
% misha@math.uchicago.edu
%
% Modified by: Vikas Sindhwani (vikass@cs.uchicago.edu)
% June 2004

% disp('Computing Graph Laplacian.');

NN=options.NN;
DISTANCEFUNCTION=options.GraphDistanceFunction;
WEIGHTTYPE=options.GraphWeights;
WEIGHTPARAM=options.GraphWeightParam;
NORMALIZE=options.GraphNormalize;

% calculate the adjacency matrix for DATA
A = adjacency(DATA, TYPE, NN, DISTANCEFUNCTION);

W = A;

% disassemble the sparse matrix
[A_i, A_j, A_v] = find(A);

switch WEIGHTTYPE

case 'distance'
for i = 1: size(A_i)
W(A_i(i), A_j(i)) = A_v(i);
end;

case 'binary'
disp('Laplacian : Using Binary weights ');
for i = 1: size(A_i)
W(A_i(i), A_j(i)) = 1;
end;

case 'heat'
%     disp(['Laplacian : Using Heat Kernel sigma : ' num2str(WEIGHTPARAM)]);
t=WEIGHTPARAM;
for i = 1: size(A_i)
W(A_i(i), A_j(i)) = exp(-A_v(i)^2/(2*t*t));
end;

otherwise
error('Unknown Weighttype');
end

D = sum(W(:,:),2);

if NORMALIZE==0
L = spdiags(D,0,speye(size(W,1)))-W;
else % normalized laplacian
D=diag(sqrt(1./D));
L=eye(size(W,1))-D*W*D;
end


ml_options.m

function options = ml_options(varargin)

% ML_OPTIONS - Generate/alter options structure for training classifiers
% ----------------------------------------------------------------------------------------%
%   options = ml_options('PARAM1',VALUE1,'PARAM2',VALUE2,...)
%
%   Creates an options structure "options" in which the named parameters
%   have the specified values.  Any unspecified parameters are set to
%   default values specified below.
%   options = ml_options (with no input arguments) creates an options structure
%   "options" where all the fields are set to default values specified below.
%
%   Example:
%          options=ml_options('Kernel','rbf','KernelParam',0.5,'NN',6);
%
%   "options" structure is as follows:
%
%            Field Name:  Description                                         : default
% -------------------------------------------------------------------------------------
%              'Kernel':  'linear' | 'rbf' | 'poly'                           : 'linear'
%         'KernelParam':  --    | sigma | degree                              :  1
%                  'NN':  number of nearest neighbor                          :  6
%'GraphDistanceFuncion':  distance function for graph: 'euclidean' | 'cosine' : 'euclidean'
%        'GraphWeights':  'binary' | 'distance' | 'heat'                      : 'binary'
%    'GraphWeightParam':  e.g For heat kernel, width to use                   :  1
%      'GraphNormalize':  Use normalized Graph laplacian (1) or not (0)       :  1
%          'ClassEdges':  Disconnect Edges across classes:yes(1) no (0)       :  0
%             'gamma_A':  RKHS norm regularization parameter (Ambient)        :  1
%             'gamma_I':  Manifold regularization parameter  (Intrinsic)      :  1
% -------------------------------------------------------------------------------------
%
% Acknowledgement: Adapted from Anton Schwaighofer's software:
%               http://www.cis.tugraz.at/igi/aschwaig/software.html
%
% Author: Vikas Sindhwani (vikass@cs.uchicago.edu)
% June 2004
% ----------------------------------------------------------------------------------------%

% options default values
options = struct('Kernel','linear', ...
'KernelParam',1, ...
'NN',6,...
'GraphDistanceFunction','euclidean',...
'GraphWeights', 'binary', ...
'GraphWeightParam',1, ...
'GraphNormalize',1, ...
'ClassEdges',0,...
'gamma_A',1.0,...
'gamma_I',1.0);

numberargs = nargin;

Names = fieldnames(options);
[m,n] = size(Names);
names = lower(Names);

i = 1;
while i <= numberargs
arg = varargin{i};
if isstr(arg)
break;
end
if ~isempty(arg)
if ~isa(arg,'struct')
error(sprintf('Expected argument %d to be a string parameter name or an options structure.', i));
end
for j = 1:m
if any(strcmp(fieldnames(arg),Names{j,:}))
val = getfield(arg, Names{j,:});
else
val = [];
end
if ~isempty(val)
[valid, errmsg] = checkfield(Names{j,:},val);
if valid
options = setfield(options, Names{j,:},val);
else
error(errmsg);
end
end
end
end
i = i + 1;
end

% A finite state machine to parse name-value pairs.
if rem(numberargs-i+1,2) ~= 0
error('Arguments must occur in name-value pairs.');
end
expectval = 0;
while i <= numberargs
arg = varargin{i};
if ~expectval
if ~isstr(arg)
error(sprintf('Expected argument %d to be a string parameter name.', i));
end
lowArg = lower(arg);
j = strmatch(lowArg,names);
if isempty(j)
error(sprintf('Unrecognized parameter name ''%s''.', arg));
elseif length(j) > 1
% Check for any exact matches (in case any names are subsets of others)
k = strmatch(lowArg,names,'exact');
if length(k) == 1
j = k;
else
msg = sprintf('Ambiguous parameter name ''%s'' ', arg);
msg = [msg '(' Names{j(1),:}];
for k = j(2:length(j))'
msg = [msg ', ' Names{k,:}];
end
msg = sprintf('%s).', msg);
error(msg);
end
end
expectval = 1;
else
[valid, errmsg] = checkfield(Names{j,:}, arg);
if valid
options = setfield(options, Names{j,:}, arg);
else
error(errmsg);
end
expectval = 0;
end
i = i + 1;
end

function [valid, errmsg] = checkfield(field,value)
% CHECKFIELD Check validity of structure field contents.
%   [VALID, MSG] = CHECKFIELD('field',V) checks the contents of the specified
%   value V to be valid for the field 'field'.
%

valid = 1;
errmsg = '';
if isempty(value)
return
end
isFloat = length(value==1) & isa(value, 'double');
isPositive = isFloat & (value>=0);
isString = isa(value, 'char');
range = [];
requireInt = 0;
switch field
case 'NN'
requireInt = 1;
range=[1 Inf];
case 'GraphNormalize'
requireInt = 1;
range=[0 1];
case 'ClassEdges'
requireInt = 1;
range=[0 1];
case {'Kernel', 'GraphWeights', 'GraphDistanceFunction'}
if ~isString,
valid = 0;
errmsg = sprintf('Invalid value for %s parameter: Must be a string', field);
end
case {'gamma_A', 'gamma_I','GraphWeightParam'}
range = [0 Inf];
case 'KernelParam'
valid = 1;
otherwise
valid = 0;
error('Unknown field name for Options structure.')
end

if ~isempty(range),
if (value<range(1)) | (value>range(2)),
valid = 0;
errmsg = sprintf('Invalid value for %s parameter: Must be scalar in the range [%g..%g]', ...
field, range(1), range(2));
end
end

if requireInt & ((value-round(value))~=0),
valid = 0;
errmsg = sprintf('Invalid value for %s parameter: Must be integer', ...
field);
end


norm2s.m

function  [X, Y, normal] =norm2s(x,y)
x = double(x);
y = double(y);

n=size(x,1);
m=size(y,1);

normal.xm=mean(x);
normal.ym=mean(y);

x=x-repmat(normal.xm,n,1);
y=y-repmat(normal.ym,m,1);

normal.xscale=sqrt(sum(sum(x.^2,2))/n);
normal.yscale=sqrt(sum(sum(y.^2,2))/m);

X=x/normal.xscale;
Y=y/normal.yscale;


plot_matches.m

function plot_matches(I1, I2, X, Y, VFCIndex, CorrectIndex)
%   PLOT_MATCHES(I1, I2, X, Y, VFCINDEX, CORRECTINDEX)
%   only plots the ture positives with blue lines, false positives with red
%   lines, and false negatives with green lines. For visibility, it plots at
%   most NUMPLOT (Default value is 50) matches proportionately.
%
% Input:
%   I1, I2: Tow input images.
%
%   X, Y: Coordinates of intrest points in I1, I2 respectively.
%
%   VFCIndex: Indexes preserved by VFC.
%
%   CorrectIndex: Ground truth indexes.
%

% Authors: Jiayi Ma (jyma2010@gmail.com)
% Date:    04/17/2012

% Define the maximum number of matches to plot

%blue = true positive+true negative, green = false negative, red = false positive
NumPlot = 50;

TruePos = intersect(VFCIndex, CorrectIndex);%Ture positive
FalsePos = setdiff(VFCIndex, CorrectIndex); %False positive
FalseNeg = setdiff(CorrectIndex, VFCIndex); %False negative

NumPos = length(TruePos)+length(FalsePos)+length(FalseNeg);
if NumPos > NumPlot
t_p = length(TruePos)/NumPos;
n1 = round(t_p*NumPlot);
f_p = length(FalsePos)/NumPos;
n2 = ceil(f_p*NumPlot);
f_n = length(FalseNeg)/NumPos;
n3 = ceil(f_n*NumPlot);
else
n1 = length(TruePos);
n2 = length(FalsePos);
n3 = length(FalseNeg);
end

per = randperm(length(TruePos));
TruePos = TruePos(per(1:n1));
per = randperm(length(FalsePos));
FalsePos = FalsePos(per(1:n2));
per = randperm(length(FalseNeg));
FalseNeg = FalseNeg(per(1:n3));

interval = 20;
WhiteInterval = 255*ones(size(I1,1), interval, 3);
figure;imagesc(cat(2, I1, WhiteInterval, I2)) ;
hold on ;
line([X(FalsePos,1)'; Y(FalsePos,1)'+size(I1,2)+interval], [X(FalsePos,2)' ;  Y(FalsePos,2)'],'linewidth', 1, 'color', 'r') ;
line([X(FalseNeg,1)'; Y(FalseNeg,1)'+size(I1,2)+interval], [X(FalseNeg,2)' ;  Y(FalseNeg,2)'],'linewidth', 1, 'color', 'g') ;
line([X(TruePos,1)'; Y(TruePos,1)'+size(I1,2)+interval], [X(TruePos,2)' ;  Y(TruePos,2)'],'linewidth', 1, 'color', 'b') ;
axis equal ;axis off  ;
drawnow;


criterion =

0.9259
1.0000
0.9604
0.9615
0.8125


## 六、Reference

posted on 2019-07-09 22:03  凯鲁嘎吉  阅读(672)  评论(0编辑  收藏