-
Notifications
You must be signed in to change notification settings - Fork 227
/
compute_grad.m
140 lines (122 loc) · 3.93 KB
/
compute_grad.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
function [grad,gy] = compute_grad(M,options)
% compute_grad - compute the gradient of an image using central differences
%
% grad = compute_grad(M,options);
% [gx,gy] = compute_grad(M,options);
%
% 'options' is a structure:
% - options.h is the sampling step size on both direction (default 1).
% - options.h1 is the sampling step size on X direction (default 1).
% - options.h2 is the sampling step size on Y direction (default 1).
% - options.type is the kind of finite difference.
% type==2 is fwd differences, ie.
% y(i) = (x(i)-x(i-1))/h, with special
% care at boundaries.
% type==1 is forward differences bilinearly interpolated in the
% middle of each pixel (be aware that you have a shift of 1/2 on X and Y for
% the location of the gradient).
% type==1 is backward differences bilinearly interpolated in the
% middle of each pixel (be aware that you have a shift of -1/2 on X and Y for
% the location of the gradient).
%
% Copyright (c) 2004 Gabriel Peyré
% warning, for matlab, X is the second axis !
[gy,gx] = gradient(M);
if nargout==1
grad = zeros([size(M), 2]);
grad(:,:,1) = gx;
grad(:,:,2) = gy;
else
grad=gx;
end
return;
if nargin<2
options.null = 0;
end
if ~isfield(options, 'h1')
options.h1 = 1;
end
h1 = options.h1;
if ~isfield(options, 'h2')
options.h2 = 1;
end
h2 = options.h2;
if isfield(options, 'h')
h1 = options.h;
h2 = options.h;
end
[n,p] = size(M);
if isfield(options, 'type')
type = options.type;
else
type = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% new code, use faster 2D differences
if type==1
% central differences on X
D1 = [M(2:end,:);M(end,:)];
D2 = [M(1,:);M(1:end-1,:)];
grad(:,:,1) = (D1-D2)/(2*h1);
grad(1,:,1) = ( 4*M(2,:) - 3*M(1,:) - M(3,:) )/(2*h1);
grad(end,:,1) = -( 4*M(end-1,:) - 3*M(end,:) - M(end-2,:) )/(2*h1);
% central differences on Y
D1 = [M(:,2:end),M(:,end)];
D2 = [M(:,1),M(:,1:end-1)];
grad(:,:,2) = (D1-D2)/(2*h2);
grad(:,1,2) = ( 4*M(:,2) - 3*M(:,1) - M(:,3) )/(2*h2);
grad(:,end,2) = -( 4*M(:,end-1) - 3*M(:,end) - M(:,end-2) )/(2*h2);
elseif type==2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% accumulate on Y
MM = ( M + [M(:,2:end),M(:,end)] )/2;
% fwd differences on X
D1 = [MM(2:end,:);MM(end,:)];
D2 = MM;
grad(:,:,1) = (D1-D2)/h1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% accumulate on X
MM = ( M + [M(2:end,:);M(end,:)] )/2;
% fwd differences on Y
D1 = [MM(:,2:end),MM(:,end)];
D2 = MM;
grad(:,:,2) = (D1-D2)/h2;
elseif type==3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% accumulate on Y
MM = ( M + [M(:,1),M(:,1:end-1)] )/2;
% fwd differences on X
D1 = MM;
D2 = [MM(1,:);MM(1:end-1,:)];
grad(:,:,1) = (D1-D2)/h;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% accumulate on Y
MM = ( M + [M(1,:);M(1:end-1,:)] )/2;
% fwd differences on Y
D1 = MM;
D2 = [MM(:,1),MM(:,1:end-1)];
grad(:,:,2) = (D1-D2)/h;
else
error('This kind of differences is not supported.');
end
% for 2 arguments return
if nargout==2
gy = grad(:,:,2);
grad = grad(:,:,1);
end
return;
if type~=1
% compute the difference in the center of each square
h = zeros(n,p,2);
for j=1:p-1
h(:,j,1) = ( grad(:,j,1)+grad(:,j+1,1) )/2;
end
for i=1:n-1
h(i,:,2) = ( grad(i,:,2)+grad(i+1,:,2) )/2;
end
if type==2 % fwd differences
grad(1:n-1,1:p-1,:) = h(1:n-1,1:p-1,:);
elseif type==3 % bwd differences
grad(2:n,2:p,:) = h(2:n,2:p,:);
end
end