Поиск ближайших соседей радиальных сегментов

Question

Во-первых, не бойтесь, судя по этому вопросу;)

Я пытаюсь реализовать дескриптор формы в 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 

Для меня это звучит как каскадная сортировка, сначала сортировка по восходящему расстоянию, а затем по восходящему индексу, но ближайшие соседи, которые я нахожу, не соответствуют бумаге.

В качестве примера я покажу вам, как ближайших соседей ПОЛУЧИТЬ для сегмента б {4,1}, это один с пометкой «EX» на рисунке 1d)

я получаю следующее список ближайших соседей для Ь {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

Answer 1

В документе говорится о «соседях по региону»; интерпретация того, что это «ближайшие соседи» в смысле евклидова расстояния, неверна. Это просто районы, которые являются соседями определенной области, а метод их поиска тривиален:

Области имеют 2 координаты: (c, s), где c обозначает, какой концентрический круг они являются частью, от 1 в центре до C на краю, а s обозначает, какой сектор они являются частью, начиная с 1, начиная с угла 0 ° до S, заканчивающегося под углом 360 °.

Каждая область, чьи координаты c и s отличаются не более 1 от координат области, является соседней областью (номера сегментов обертываются от S до 1.) В зависимости от местоположения региона существует 3 случая: (как показано на фиг. 1d)

• область является средняя область (отмечены MI), например, область б (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)
  Существует только один соседний круг и 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), например, область б (3,1)
  Существует только один круг, и соседние два соседних секторов, так что 6 регионов в общей сложности:
  каждый регион в круг 2 или 3 и сектор 8, 1 или 2:
  b(2,8), b(2,1), b(2,2), b(3,8), b(3,1), b(3,2)

Answer 2

чтобы добавить немного MATLAB к great answer из @ 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,[])]; 

1. Первый определяет an inline if для использования в анонимных функциях, чтобы иметь возможность определять функции в отдельных файлах (это необходимо для случая c). Это имеет синтаксис iif(condition1,result1,condition2,result2,...), где каждое условие проверяется друг за другом, и возвращается результат, соответствующий первому условию, дающему true.

2. Второй определяет радиальные соседних индексов с, если: возвращение [с-1, с + 1], если ни одна из них не является неопределенным (что привело бы к массива связанным нарушений)

3. третий определяет периодическую индексацию для угловых секторов, такую, что для S=4nsfun(2)==[1, 3] и nsfun(4)==[3, 1].

4. Я просто добавил пример, когда для данной действуют пара c,sneighbs(c,s) вернет подмассив b(1:C,1:S), которые являются соседями: сначала вверх-вниз/соседи влево-вправо, затем (до) четыре угловых.

Answer 3

После того, как проводить слишком много времени на это, я понял это, уловка, чтобы вычислить соседние области элемента Ь {C, S} был:

• принимать все сегменты для соседние кольца с, в то числе с самого по себе, т.е. c-1, c, c+1, если это внутренний самым круг только тогда c,c+1 если это внешнее большинство принимает c-1, c

• Вычислить расстояние евклидового из б {C, S} для все выбранной центроиды s с предыдущего шага, включая сам пункт b {c, s}

• Отсортируйте расстояния в порядке возрастания и возьмите первые N сегментов.

Дескриптор работает достаточно хорошо, однако это вращательное инвариантность зависит от оси с самой высокой плотностью в некоторых случаях это является очень чувствительным к шуму/окклюзии, т.е. дескриптор может составлять +/- 1 вращение выкл.

Вот полная реализация дескрипторов CBSM в MATLAB, не оптимизированы в любом случае (я слышал MATLAB ненавидит за-петель):

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); 

Выход:

Только для пинков, некоторые страшные картины происходило при кодировании его: