-
Notifications
You must be signed in to change notification settings - Fork 0
/
sp_find_sift_grid.m
134 lines (110 loc) · 4.09 KB
/
sp_find_sift_grid.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
function sift_arr = sp_find_sift_grid(I, grid_x, grid_y, patch_size, sigma_edge)
% parameters
num_angles = 8;
num_bins = 4;
num_samples = num_bins * num_bins;
alpha = 9;
if nargin < 5
sigma_edge = 1;
end
angle_step = 2 * pi / num_angles;
angles = 0:angle_step:2*pi;
angles(num_angles+1) = []; % bin centers
[hgt wid] = size(I);
num_patches = numel(grid_x);
sift_arr = zeros(num_patches, num_samples * num_angles);
[G_X,G_Y]=gen_dgauss(sigma_edge);
I_X = filter2(G_X, I, 'same'); % vertical edges
I_Y = filter2(G_Y, I, 'same'); % horizontal edges
I_mag = sqrt(I_X.^2 + I_Y.^2); % gradient magnitude
I_theta = atan2(I_Y,I_X);
I_theta(find(isnan(I_theta))) = 0; % necessary????
% make default grid of samples (centered at zero, width 2)
interval = 2/num_bins:2/num_bins:2;
interval = interval - (1/num_bins + 1);
[sample_x sample_y] = meshgrid(interval, interval);
sample_x = reshape(sample_x, [1 num_samples]);
sample_y = reshape(sample_y, [1 num_samples]);
% make orientation images
I_orientation = zeros(hgt, wid, num_angles);
% for each histogram angle
for a=1:num_angles
% compute each orientation channel
tmp = cos(I_theta - angles(a)).^alpha;
tmp = tmp .* (tmp > 0);
% weight by magnitude
I_orientation(:,:,a) = tmp .* I_mag;
end
% for all patches
for i=1:num_patches
r = patch_size/2;
cx = grid_x(i) + r - 0.5;
cy = grid_y(i) + r - 0.5;
% find coordinates of sample points (bin centers)
sample_x_t = sample_x * r + cx;
sample_y_t = sample_y * r + cy;
sample_res = sample_y_t(2) - sample_y_t(1);
% find window of pixels that contributes to this descriptor
x_lo = grid_x(i);
x_hi = grid_x(i) + patch_size - 1;
y_lo = grid_y(i);
y_hi = grid_y(i) + patch_size - 1;
% find coordinates of pixels
[sample_px, sample_py] = meshgrid(x_lo:x_hi,y_lo:y_hi);
num_pix = numel(sample_px);
sample_px = reshape(sample_px, [num_pix 1]);
sample_py = reshape(sample_py, [num_pix 1]);
% find (horiz, vert) distance between each pixel and each grid sample
dist_px = abs(repmat(sample_px, [1 num_samples]) - repmat(sample_x_t, [num_pix 1]));
dist_py = abs(repmat(sample_py, [1 num_samples]) - repmat(sample_y_t, [num_pix 1]));
% find weight of contribution of each pixel to each bin
weights_x = dist_px/sample_res;
weights_x = (1 - weights_x) .* (weights_x <= 1);
weights_y = dist_py/sample_res;
weights_y = (1 - weights_y) .* (weights_y <= 1);
weights = weights_x .* weights_y;
% % make sure that the weights for each pixel sum to one?
% tmp = sum(weights,2);
% tmp = tmp + (tmp == 0);
% weights = weights ./ repmat(tmp, [1 num_samples]);
% make sift descriptor
curr_sift = zeros(num_angles, num_samples);
for a = 1:num_angles
tmp = reshape(I_orientation(y_lo:y_hi,x_lo:x_hi,a),[num_pix 1]);
tmp = repmat(tmp, [1 num_samples]);
curr_sift(a,:) = sum(tmp .* weights);
end
sift_arr(i,:) = reshape(curr_sift, [1 num_samples * num_angles]);
% % visualization
% if sigma_edge >= 3
% subplot(1,2,1);
% rescale_and_imshow(I(y_lo:y_hi,x_lo:x_hi) .* reshape(sum(weights,2), [y_hi-y_lo+1,x_hi-x_lo+1]));
% subplot(1,2,2);
% rescale_and_imshow(curr_sift);
% pause;
% end
end
function G=gen_gauss(sigma)
if all(size(sigma)==[1, 1])
% isotropic gaussian
f_wid = 4 * ceil(sigma) + 1;
G = fspecial('gaussian', f_wid, sigma);
% G = normpdf(-f_wid:f_wid,0,sigma);
% G = G' * G;
else
% anisotropic gaussian
f_wid_x = 2 * ceil(sigma(1)) + 1;
f_wid_y = 2 * ceil(sigma(2)) + 1;
G_x = normpdf(-f_wid_x:f_wid_x,0,sigma(1));
G_y = normpdf(-f_wid_y:f_wid_y,0,sigma(2));
G = G_y' * G_x;
end
function [GX,GY]=gen_dgauss(sigma)
% laplacian of size sigma
%f_wid = 4 * floor(sigma);
%G = normpdf(-f_wid:f_wid,0,sigma);
%G = G' * G;
G = gen_gauss(sigma);
[GX,GY] = gradient(G);
GX = GX * 2 ./ sum(sum(abs(GX)));
GY = GY * 2 ./ sum(sum(abs(GY)));