寻找径向线段的最近邻居
Finding nearest neighbours of radial segments
首先，不要被这个问题的外表吓到 ;)
我正在尝试在 matlab 中实现一个称为圆形模糊形状模型的形状描述符，其中一部分是为每个径向段获取最近邻居列表，如图 1d)
我在 MATLAB 中进行了直接而简单的实现，但我卡在了算法的第 5 步和第 6 步，主要是因为我无法理解定义：
Xb{c,s} = {b1, ..., b{c*s}} as the sorted set of the elements in B*
so that d(b*{c,s}, bi*) <= d(b*{c,s}, bj*), i<j
对我来说这听起来像是级联排序，先按距离升序排序，然后按索引升序排序，但我找到的最近邻居不是根据论文。
作为示例，我向您展示了我为段 b{4,1} 获得的最近邻居，这是图 1d)
中标记为 "EX" 的那个
我得到以下 b{4,1} 的最近邻居列表：b{3,2}, b{3,1}, b{3,8}, b{2,1}, b{2,8}
根据论文正确的是：b{4,2}, b{4,8}, b{3,2}, b{3,1}, b{3,8}
然而，我的点实际上是最接近所选线段的一组（按欧氏距离测量）！距离 b{4,1} <=> b{2,1} 小于 b{4,1} <=> b{4,2} 或 b{4,1} <=> b{4,8}...
这是我的（丑陋但直接的）MATLAB 代码：
width = 734;
height = 734;
assert(width == height, 'Image must be square in size!');
% Radius of the correlogram
R = width;
% Number of circles in correlogram
C = 4;
% Number of sections in correlogram
S = 8;
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 360/S;
% set of bins for the circular description of I
B = zeros(C, S);
% centroid coordinates for bins
B_star = zeros(C,S,2);
% calculate centroids of bins
for c=1:C
for s=1:S
alpha = deg2rad(max(s-1, 0)*g + g/2);
r = d*max((c-1),0) + d/2;
B_star(c,s,1) = r*cos(alpha);
B_star(c,s,2) = r*sin(alpha);
end
end
% create sorted list of bin numbers which fullfill
% d(b{c,s}*, bi*) <= d(b{c,s}, bj*) where i<j
% B_star_dists is a simple square distance matrix for getting
% the distance between two centroids c_i,s_i and c_j,s_j
B_star_dists = zeros(C*S, C*S);
for i=1:C*S
[c_i, s_i] = ind2sub([C,S], i);
% x,y centroid coordinates for point i
b_star_i = [B_star(c_i, s_i, 1), B_star(c_i, s_i, 2)];
for j=1:C*S
[c_j, s_j] = ind2sub([C,S], j);
% x,y centroid coordinates for point j
b_star_j = [B_star(c_j, s_j, 1), B_star(c_j, s_j, 2)];
% store the euclidean distance between these two centroids
% in the distance matrix.
B_star_dists(i,j) = norm(b_star_i - b_star_j);
end
end
% calculate nearest neighbour "centroids" for each centroid
% B_NN is a cell array, B{idx} gives an array of indexes to the
% nearest neighbour centroids.
B_NN = cell(C*S, 1);
for i=1:C*S
[c_i, s_i] = ind2sub([C,S], i);
% get a (C*S)x2 matrix of all distances, the first column are the array
% indexes and the second column are the distances e.g
% 1 d1
% 2 d2
% .. ..
% CS d{c,s}
dists = [transpose(1:C*S), B_star_dists(:, i)];
% sort ascending by the distances first (e.g second column) then
% sort ascending by the array index (e.g first column)
dists = sortrows(dists, [2,1]);
% middle section has nine neighbours, set as default
neighbour_count = 9;
if c_i == 1
% inner region has S+3 neighbours
neighbour_count = S+3;
elseif c_i == C
% outer most ring has 6 neighbours
neighbour_count = 6;
end
B_NN{i} = dists(1:neighbour_count,1);
end
% FROM HERE ON JUST VISUALIZATION CODE
figure(1);
hold on;
for c=1:C
% plot circles
r = c*d;
plot(r*cos(0:pi/50:2*pi), r*sin(0:pi/50:2*pi), 'k:');
end
for s=1:S
% plot lines
line_len = C*d;
alpha = deg2rad(s*g);
start_pt = [0, 0];
end_pt = start_pt + line_len.*[cos(alpha), sin(alpha)];
plot([start_pt(1), end_pt(1)], [start_pt(2), end_pt(2)], 'k-');
end
for c=1:C
% plot centroids of segments
for s=1:S
segment_centroid = B_star(c,s, :);
plot(segment_centroid(1), segment_centroid(2), '.k');
end
end
% plot some nearest neighbours
% list of [C;S]
plot_nn = [4;1];
for i = 1:size(plot_nn,2)
start_c = plot_nn(1,i);
start_s = plot_nn(2,i);
start_pt = [B_star(start_c, start_s,1), B_star(start_c, start_s,2)];
start_idx = sub2ind([C, S], start_c, start_s);
plot(start_pt(1), start_pt(2), 'xb');
nn_idx_list = B_NN{start_idx};
for j = 1:length(nn_idx_list)
nn_idx = nn_idx_list(j);
[nn_c, nn_s] = ind2sub([C, S], nn_idx);
nn_pt = [B_star(nn_c, nn_s,1), B_star(nn_c, nn_s,2)];
plot(nn_pt(1), nn_pt(2), 'xr');
end
end
可以找到完整的论文here
论文讲的是"region neighbours"；这些是欧几里得距离感中的 "nearest neighbours" 的解释是不正确的。它们只是与某个区域相邻的区域，找到它们的方法很简单：
这些区域有 2 个坐标：(c,s)，其中 c 表示它们属于哪个同心圆，从中心的 1 到边缘的 C，s 表示它们属于哪个扇区，从 1 开始于角度 0° 到 S 结束于角度 360°。
每个区域的c和s坐标与区域坐标最多相差1，是一个相邻区域（段号从S到1。）根据区域的位置，有3种情况：(如图1d)
该区域是中间区域（标有MI）e.g.区域 b(2,4)
有2个相邻的圆圈和2个相邻的扇区，所以总共有9个区域：
圆圈 1、2 或 3 和扇区 3、4 或 5 中的每个区域：
b(1,3), b(2,3), b(3,3), b(1,4), b(2,4), b(3,4), b(1,5), b(2,5), b(3,5)
该区域是内部区域（标记为IN）例如区域 b(1,8)
只有1个相邻圆圈和2个相邻扇区，但所有内部区域都是邻居，所以总共S + 3个区域：
圆圈 2 和扇区 7、8 或 1 中的每个区域：
b(2,7), b(2,8), b(2,1)
以及内圈中的每个区域：
b(1,1), b(1,2), b(1,3), b(1,4), b(1,5), b(1,6), b(1,7), b(1,8)
该区域是外部区域（标记为EX）例如区域 b(3,1)
只有一个相邻的圆和2个相邻的扇区，所以总共有6个区域：
圆圈 2 或 3 和扇区 8、1 或 2 中的每个区域：
b(2,8), b(2,1), b(2,2), b(3,8), b(3,1), b(3,2)
要向 @m69 的 添加一点 matlab，您可以像这样自动执行邻居索引：
%assume C and S defined according to the answer of @m69
iif=@(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
ncfun=@(c) iif(c==1,c+1,c==C,c-1,true,[c-1, c+1]);
nsfun=@(s)mod([s-1, s+1]-1,S)+1;
neighbs=@(c,s) [b(c,nsfun(s)), b(ncfun(c),s)', reshape(b(ncfun(c),nsfun(s)),1,[])];
第一个定义了 用于匿名函数，以避免在单独的文件中定义函数（这对于 c 情况是需要的）。这具有语法 iif(condition1,result1,condition2,result2,...)，其中每个条件一个接一个地被测试，并且对应于产生 true 的第一个条件的结果是 returned.
第二个使用 if 定义径向邻居索引：return [c-1, c+1] 除非其中任何一个未定义（这将导致数组边界违规)
第三个为 angular 扇区定义了一个周期性索引，例如 S=4 nsfun(2)==[1, 3] 和 nsfun(4)==[3, 1].
我刚刚添加了一个示例，其中对于给定的有效一对 c,s neighbs(c,s) 将 return b(1:C,1:S) 的子数组作为邻居：首先是 up-down/left-right 个邻居，然后是（最多）四个角落的邻居。
在这上面花了太多时间之后我想通了，计算元素 b{c,s} 相邻区域的技巧是：
取 c 的相邻环的所有段，包括 c 本身，即 c-1, c, c+1，如果它是最内圈，则仅 c,c+1 如果它是最外圈c-1, c
计算b{c,s}到上一步选中的所有质心的欧氏距离，包括b{c,s}本身
将距离升序排列，取前N段
描述符工作得很好，但是它的旋转不变性取决于密度最高的轴，在某些情况下这对 noise/occlusion 非常敏感，即描述符可能偏离 +/- 1 旋转。
这是 MATLAB 中完整的 CBSM 描述符实现，没有以任何方式优化（我听说 MATLAB 讨厌 for 循环）：
csbm.m
function [descriptor] = csbm(I, R, C, S)
% Input:
% ------
% I = The image, binary, must be square in sice
% R = The radius of the correlogram, could be half the image width
% C = Number of circles
% S = Number of segments per circle
% Output:
% -------
% A C*S-by-1 row vector, the descriptor
[width, height] = size(I);
assert(width == height, 'Image must be square in size!');
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 2*pi/S;
% centroid coordinates for bins
B_star = zeros(C*S,2);
% initialize the descriptor
descriptor = zeros(C*S,1);
% calculate centroids of bins
for c=1:C
for s=1:S
alpha = max(s-1, 0)*g + g/2;
r = d*max((c-1),0) + d/2;
ind = (c-1)*S+s; %sub2ind(cs_size, c,s);
B_star(ind, :) = [r*cos(alpha), r*sin(alpha)];
end
end
% calculate nearest neighbor regions
B_NN = cell(C*S,1);
for c=1:C
min_circle = max(1,c-1);
max_circle = min(C,c+1);
start_sel = (min_circle-1)*S+1;
end_sel = (max_circle)*S;
centroids = B_star(start_sel:end_sel, :);
for s=1:S
current_ind = (c-1)*S+s;
centroid = B_star(current_ind, :);
num_centroids = length(centroids);
distances = sqrt(sum((centroids-repmat(centroid, num_centroids,1)).^2,2));
distances = horzcat(distances, transpose((start_sel:end_sel)));
distances = sortrows(distances, [1,2]);
neighbour_count = 9;
if c == 1
% inner region has S+3 neighbours
neighbour_count = S+3;
elseif c == C
% outer most ring has 6 neighbours
neighbour_count = 6;
end
B_NN{current_ind} = distances(1:neighbour_count, 2);
end
end
for x=1:width
x_centered = x-width/2;
for y=1:width
if I(x,y) == 0
continue;
end
y_centered = y-width/2;
% bin the image point
r = sqrt(x_centered^2 + y_centered^2);
a = wrapTo2Pi(atan2(y_centered, x_centered));
% get bin
c_pixel = max(1, floor(r/d));
s_pixel = max(1, floor(a/g));
if c_pixel > C
continue;
end
ind_pixel = (c_pixel-1)*S+s_pixel;
pt_pixel = [x_centered, y_centered];
% get neighbours of this bin
neighbours = B_NN{ind_pixel};
% calculate distance to neighbours
nn_dists = sqrt(sum((B_star(neighbours, :) - repmat(pt_pixel, length(neighbours), 1)).^2,2));
D = sum(1./nn_dists);
% update the probabilities vector
descriptor(neighbours) = descriptor(neighbours) + 1./(nn_dists.*D);
end
end
% normalize the vector v
descriptor = descriptor./sum(descriptor);
% make it rotationally invariant
G = zeros(S/2, 2*C);
for s=1:S/2
for c=1:C
G(s,c) = descriptor((c-1)*S+s);
G(s,c+C) = descriptor((c-1)*S+s+S/2);
end
end
[~, max_G_idx] = max(sum(G,2));
L_G = 0;
R_G = 0;
for j=1:C
for k=1:S
if k > (max_G_idx) && k < (max_G_idx+S/2)
L_G = L_G + descriptor((j-1)*S+k);
elseif k ~= max_G_idx && k ~= (max_G_idx+S/2)
R_G = R_G + descriptor((j-1)*S+k);
end
end
end
if L_G > R_G
% B is rotated k=i+S/2-1 positions to the left:
fprintf('rotate %d to the left\n', max_G_idx+S/2-1);
rotate_by = -(max_G_idx+S/2-1);
else
% B is rotated k=i-1 positions to the right:
fprintf('rotate %d to the right\n', max_G_idx-1);
rotate_by = -(max_G_idx-1);
end
% segments are grouped by circle
% so for every circle we get all segments and circular shift them
for c=1:C
range_sel = ((c-1)*S+1):(c*S);
segments = descriptor(range_sel);
descriptor(range_sel) = circshift(segments, [rotate_by,0]);
end
end
plot_descriptor.m
function plot_descriptor(R,C,S, descriptor)
% Input:
% ------
% R Radius for the correlogram in pixels, can be arbitrary
% C Number of circles
% S Number of segments per circle
% descriptor The C*S-by-1 descriptor vector
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 2*pi/S;
% full image
[x,y] = meshgrid(-R:R);
[theta, rho] = cart2pol(x,y);
theta = wrapTo2Pi(theta);
brightness = zeros(size(rho));
min_val = min(descriptor);
max_val = max(descriptor);
scale_fact = 1/(max_val-min_val);
for c=1:C
rInner = (c-1)*d;
rOuter = c*d;
for s=1:S
idx = (c-1)*S+s;
minTheta = (s-1)*g;
maxTheta = (s*g);
matching_theta = find(theta > minTheta & theta < maxTheta);
matching_rho = find(rho > rInner & rho < rOuter);
matching_idx = intersect(matching_theta, matching_rho);
intensity = descriptor(idx)*scale_fact;
brightness(matching_idx) = intensity;
end
end
figure; imshow(mat2gray(brightness));
如何运行:
I = imread('bat-18.gif');
I = transpose(im2bw(I, graythresh(I)));
I = edge(I);
[width, height] = size(I);
R = width/2;
C = 24;
S = 24;
descriptor = csbm(I, R, C, S);
plot_descriptor(R,C,S,descriptor);
输出：
为了好玩，编码时出现了一些可怕的画面：
首先，不要被这个问题的外表吓到 ;)
我正在尝试在 matlab 中实现一个称为圆形模糊形状模型的形状描述符，其中一部分是为每个径向段获取最近邻居列表，如图 1d)
我在 MATLAB 中进行了直接而简单的实现，但我卡在了算法的第 5 步和第 6 步，主要是因为我无法理解定义：
Xb{c,s} = {b1, ..., b{c*s}} as the sorted set of the elements in B*
so that d(b*{c,s}, bi*) <= d(b*{c,s}, bj*), i<j
对我来说这听起来像是级联排序，先按距离升序排序，然后按索引升序排序，但我找到的最近邻居不是根据论文。
作为示例，我向您展示了我为段 b{4,1} 获得的最近邻居，这是图 1d)
中标记为 "EX" 的那个我得到以下 b{4,1} 的最近邻居列表：b{3,2}, b{3,1}, b{3,8}, b{2,1}, b{2,8}
根据论文正确的是：b{4,2}, b{4,8}, b{3,2}, b{3,1}, b{3,8}
然而，我的点实际上是最接近所选线段的一组（按欧氏距离测量）！距离 b{4,1} <=> b{2,1} 小于 b{4,1} <=> b{4,2} 或 b{4,1} <=> b{4,8}...
这是我的（丑陋但直接的）MATLAB 代码：
width = 734;
height = 734;
assert(width == height, 'Image must be square in size!');
% Radius of the correlogram
R = width;
% Number of circles in correlogram
C = 4;
% Number of sections in correlogram
S = 8;
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 360/S;
% set of bins for the circular description of I
B = zeros(C, S);
% centroid coordinates for bins
B_star = zeros(C,S,2);
% calculate centroids of bins
for c=1:C
for s=1:S
alpha = deg2rad(max(s-1, 0)*g + g/2);
r = d*max((c-1),0) + d/2;
B_star(c,s,1) = r*cos(alpha);
B_star(c,s,2) = r*sin(alpha);
end
end
% create sorted list of bin numbers which fullfill
% d(b{c,s}*, bi*) <= d(b{c,s}, bj*) where i<j
% B_star_dists is a simple square distance matrix for getting
% the distance between two centroids c_i,s_i and c_j,s_j
B_star_dists = zeros(C*S, C*S);
for i=1:C*S
[c_i, s_i] = ind2sub([C,S], i);
% x,y centroid coordinates for point i
b_star_i = [B_star(c_i, s_i, 1), B_star(c_i, s_i, 2)];
for j=1:C*S
[c_j, s_j] = ind2sub([C,S], j);
% x,y centroid coordinates for point j
b_star_j = [B_star(c_j, s_j, 1), B_star(c_j, s_j, 2)];
% store the euclidean distance between these two centroids
% in the distance matrix.
B_star_dists(i,j) = norm(b_star_i - b_star_j);
end
end
% calculate nearest neighbour "centroids" for each centroid
% B_NN is a cell array, B{idx} gives an array of indexes to the
% nearest neighbour centroids.
B_NN = cell(C*S, 1);
for i=1:C*S
[c_i, s_i] = ind2sub([C,S], i);
% get a (C*S)x2 matrix of all distances, the first column are the array
% indexes and the second column are the distances e.g
% 1 d1
% 2 d2
% .. ..
% CS d{c,s}
dists = [transpose(1:C*S), B_star_dists(:, i)];
% sort ascending by the distances first (e.g second column) then
% sort ascending by the array index (e.g first column)
dists = sortrows(dists, [2,1]);
% middle section has nine neighbours, set as default
neighbour_count = 9;
if c_i == 1
% inner region has S+3 neighbours
neighbour_count = S+3;
elseif c_i == C
% outer most ring has 6 neighbours
neighbour_count = 6;
end
B_NN{i} = dists(1:neighbour_count,1);
end
% FROM HERE ON JUST VISUALIZATION CODE
figure(1);
hold on;
for c=1:C
% plot circles
r = c*d;
plot(r*cos(0:pi/50:2*pi), r*sin(0:pi/50:2*pi), 'k:');
end
for s=1:S
% plot lines
line_len = C*d;
alpha = deg2rad(s*g);
start_pt = [0, 0];
end_pt = start_pt + line_len.*[cos(alpha), sin(alpha)];
plot([start_pt(1), end_pt(1)], [start_pt(2), end_pt(2)], 'k-');
end
for c=1:C
% plot centroids of segments
for s=1:S
segment_centroid = B_star(c,s, :);
plot(segment_centroid(1), segment_centroid(2), '.k');
end
end
% plot some nearest neighbours
% list of [C;S]
plot_nn = [4;1];
for i = 1:size(plot_nn,2)
start_c = plot_nn(1,i);
start_s = plot_nn(2,i);
start_pt = [B_star(start_c, start_s,1), B_star(start_c, start_s,2)];
start_idx = sub2ind([C, S], start_c, start_s);
plot(start_pt(1), start_pt(2), 'xb');
nn_idx_list = B_NN{start_idx};
for j = 1:length(nn_idx_list)
nn_idx = nn_idx_list(j);
[nn_c, nn_s] = ind2sub([C, S], nn_idx);
nn_pt = [B_star(nn_c, nn_s,1), B_star(nn_c, nn_s,2)];
plot(nn_pt(1), nn_pt(2), 'xr');
end
end
可以找到完整的论文here
论文讲的是"region neighbours"；这些是欧几里得距离感中的 "nearest neighbours" 的解释是不正确的。它们只是与某个区域相邻的区域，找到它们的方法很简单：
这些区域有 2 个坐标：(c,s)，其中 c 表示它们属于哪个同心圆，从中心的 1 到边缘的 C，s 表示它们属于哪个扇区，从 1 开始于角度 0° 到 S 结束于角度 360°。
每个区域的c和s坐标与区域坐标最多相差1，是一个相邻区域（段号从S到1。）根据区域的位置，有3种情况：(如图1d)
该区域是中间区域（标有MI）e.g.区域 b(2,4)
有2个相邻的圆圈和2个相邻的扇区，所以总共有9个区域：
圆圈 1、2 或 3 和扇区 3、4 或 5 中的每个区域：
b(1,3), b(2,3), b(3,3), b(1,4), b(2,4), b(3,4), b(1,5), b(2,5), b(3,5)该区域是内部区域（标记为IN）例如区域 b(1,8)
只有1个相邻圆圈和2个相邻扇区，但所有内部区域都是邻居，所以总共S + 3个区域：
圆圈 2 和扇区 7、8 或 1 中的每个区域：
b(2,7), b(2,8), b(2,1)
以及内圈中的每个区域：
b(1,1), b(1,2), b(1,3), b(1,4), b(1,5), b(1,6), b(1,7), b(1,8)该区域是外部区域（标记为EX）例如区域 b(3,1)
只有一个相邻的圆和2个相邻的扇区，所以总共有6个区域：
圆圈 2 或 3 和扇区 8、1 或 2 中的每个区域：
b(2,8), b(2,1), b(2,2), b(3,8), b(3,1), b(3,2)
要向 @m69 的
%assume C and S defined according to the answer of @m69
iif=@(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
ncfun=@(c) iif(c==1,c+1,c==C,c-1,true,[c-1, c+1]);
nsfun=@(s)mod([s-1, s+1]-1,S)+1;
neighbs=@(c,s) [b(c,nsfun(s)), b(ncfun(c),s)', reshape(b(ncfun(c),nsfun(s)),1,[])];
第一个定义了
c情况是需要的）。这具有语法iif(condition1,result1,condition2,result2,...)，其中每个条件一个接一个地被测试，并且对应于产生true的第一个条件的结果是 returned.第二个使用 if 定义径向邻居索引：return [c-1, c+1] 除非其中任何一个未定义（这将导致数组边界违规)
第三个为 angular 扇区定义了一个周期性索引，例如
S=4nsfun(2)==[1, 3]和nsfun(4)==[3, 1].我刚刚添加了一个示例，其中对于给定的有效一对
c,sneighbs(c,s)将 returnb(1:C,1:S)的子数组作为邻居：首先是 up-down/left-right 个邻居，然后是（最多）四个角落的邻居。
在这上面花了太多时间之后我想通了，计算元素 b{c,s} 相邻区域的技巧是：
取 c 的相邻环的所有段，包括 c 本身，即
c-1, c, c+1，如果它是最内圈，则仅c,c+1如果它是最外圈c-1, c计算b{c,s}到上一步选中的所有质心的欧氏距离，包括b{c,s}本身
将距离升序排列，取前N段
描述符工作得很好，但是它的旋转不变性取决于密度最高的轴，在某些情况下这对 noise/occlusion 非常敏感，即描述符可能偏离 +/- 1 旋转。
这是 MATLAB 中完整的 CBSM 描述符实现，没有以任何方式优化（我听说 MATLAB 讨厌 for 循环）：
csbm.m
function [descriptor] = csbm(I, R, C, S)
% Input:
% ------
% I = The image, binary, must be square in sice
% R = The radius of the correlogram, could be half the image width
% C = Number of circles
% S = Number of segments per circle
% Output:
% -------
% A C*S-by-1 row vector, the descriptor
[width, height] = size(I);
assert(width == height, 'Image must be square in size!');
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 2*pi/S;
% centroid coordinates for bins
B_star = zeros(C*S,2);
% initialize the descriptor
descriptor = zeros(C*S,1);
% calculate centroids of bins
for c=1:C
for s=1:S
alpha = max(s-1, 0)*g + g/2;
r = d*max((c-1),0) + d/2;
ind = (c-1)*S+s; %sub2ind(cs_size, c,s);
B_star(ind, :) = [r*cos(alpha), r*sin(alpha)];
end
end
% calculate nearest neighbor regions
B_NN = cell(C*S,1);
for c=1:C
min_circle = max(1,c-1);
max_circle = min(C,c+1);
start_sel = (min_circle-1)*S+1;
end_sel = (max_circle)*S;
centroids = B_star(start_sel:end_sel, :);
for s=1:S
current_ind = (c-1)*S+s;
centroid = B_star(current_ind, :);
num_centroids = length(centroids);
distances = sqrt(sum((centroids-repmat(centroid, num_centroids,1)).^2,2));
distances = horzcat(distances, transpose((start_sel:end_sel)));
distances = sortrows(distances, [1,2]);
neighbour_count = 9;
if c == 1
% inner region has S+3 neighbours
neighbour_count = S+3;
elseif c == C
% outer most ring has 6 neighbours
neighbour_count = 6;
end
B_NN{current_ind} = distances(1:neighbour_count, 2);
end
end
for x=1:width
x_centered = x-width/2;
for y=1:width
if I(x,y) == 0
continue;
end
y_centered = y-width/2;
% bin the image point
r = sqrt(x_centered^2 + y_centered^2);
a = wrapTo2Pi(atan2(y_centered, x_centered));
% get bin
c_pixel = max(1, floor(r/d));
s_pixel = max(1, floor(a/g));
if c_pixel > C
continue;
end
ind_pixel = (c_pixel-1)*S+s_pixel;
pt_pixel = [x_centered, y_centered];
% get neighbours of this bin
neighbours = B_NN{ind_pixel};
% calculate distance to neighbours
nn_dists = sqrt(sum((B_star(neighbours, :) - repmat(pt_pixel, length(neighbours), 1)).^2,2));
D = sum(1./nn_dists);
% update the probabilities vector
descriptor(neighbours) = descriptor(neighbours) + 1./(nn_dists.*D);
end
end
% normalize the vector v
descriptor = descriptor./sum(descriptor);
% make it rotationally invariant
G = zeros(S/2, 2*C);
for s=1:S/2
for c=1:C
G(s,c) = descriptor((c-1)*S+s);
G(s,c+C) = descriptor((c-1)*S+s+S/2);
end
end
[~, max_G_idx] = max(sum(G,2));
L_G = 0;
R_G = 0;
for j=1:C
for k=1:S
if k > (max_G_idx) && k < (max_G_idx+S/2)
L_G = L_G + descriptor((j-1)*S+k);
elseif k ~= max_G_idx && k ~= (max_G_idx+S/2)
R_G = R_G + descriptor((j-1)*S+k);
end
end
end
if L_G > R_G
% B is rotated k=i+S/2-1 positions to the left:
fprintf('rotate %d to the left\n', max_G_idx+S/2-1);
rotate_by = -(max_G_idx+S/2-1);
else
% B is rotated k=i-1 positions to the right:
fprintf('rotate %d to the right\n', max_G_idx-1);
rotate_by = -(max_G_idx-1);
end
% segments are grouped by circle
% so for every circle we get all segments and circular shift them
for c=1:C
range_sel = ((c-1)*S+1):(c*S);
segments = descriptor(range_sel);
descriptor(range_sel) = circshift(segments, [rotate_by,0]);
end
end
plot_descriptor.m
function plot_descriptor(R,C,S, descriptor)
% Input:
% ------
% R Radius for the correlogram in pixels, can be arbitrary
% C Number of circles
% S Number of segments per circle
% descriptor The C*S-by-1 descriptor vector
% "width" of ring segments
d = R/C;
% angle of one segment in degrees
g = 2*pi/S;
% full image
[x,y] = meshgrid(-R:R);
[theta, rho] = cart2pol(x,y);
theta = wrapTo2Pi(theta);
brightness = zeros(size(rho));
min_val = min(descriptor);
max_val = max(descriptor);
scale_fact = 1/(max_val-min_val);
for c=1:C
rInner = (c-1)*d;
rOuter = c*d;
for s=1:S
idx = (c-1)*S+s;
minTheta = (s-1)*g;
maxTheta = (s*g);
matching_theta = find(theta > minTheta & theta < maxTheta);
matching_rho = find(rho > rInner & rho < rOuter);
matching_idx = intersect(matching_theta, matching_rho);
intensity = descriptor(idx)*scale_fact;
brightness(matching_idx) = intensity;
end
end
figure; imshow(mat2gray(brightness));
如何运行:
I = imread('bat-18.gif');
I = transpose(im2bw(I, graythresh(I)));
I = edge(I);
[width, height] = size(I);
R = width/2;
C = 24;
S = 24;
descriptor = csbm(I, R, C, S);
plot_descriptor(R,C,S,descriptor);
输出：
为了好玩，编码时出现了一些可怕的画面：