Поиск локальных максимумов для 3d-матрицы

Question

Я пытаюсь реализовать преобразование Хафа в Matlab, чтобы найти круги в picture.

В аккумуляторной матрице глобальный максимум равен 105 по адресу A(32,31,24). Поэтому я могу это получить: max circle

Проблема в том, как я могу найти локальные максимумы, чтобы найти остальные круги?

Я написал это, чтобы найти A(i,j,k), который больше, чем 26 соседних точек (26-Connected воксельных окрестностей):

[i j k]=find(A~=0) ; 
f=0; 
for s=1:size(i) 
    if(i(s)~=100&&j(s)~=100&&k(s)~=141&&i(s)~=1&&j(s)~=1&&k(s)~=1) 
    if  (A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s)-1)) 
     f=f+1 
     i(s) 
     j(s) 
     k(s) 
    end 
    end 
end 

Почему я никогда не получил их правильно i, j, k и f всегда 0? Думаю, я должен хотя бы найти (32,31,24) и f=1?

Может ли кто-нибудь мне помочь? Большое вам спасибо!

Полный код здесь:

im=imread('C:\Users\dell\Desktop\tp-complet\Hough\four.png'); 

sigma = 0.3; 
gausFilter = fspecial('gaussian',[5 5],sigma); 
sobelFilter=fspecial('sobel'); 
img=imfilter(im,gausFilter,'replicate'); 

ims=edge(img,'sobel'); 

rmax=size(im,1); 
cmax=size(im,2); 
radmax=round(sqrt(rmax^2+cmax^2)); 

for i=1:rmax 
    for j=1:cmax 
     for k=1:radmax 
      A(i,j,k)=0; 
     end 
    end 
end 

[r c]=find(ims==1); 
length=size(r); 

for k=1:length 
    for l=1:rmax 
    for m=1:cmax 
     if((l~=r(k))&&(m~=c(k))) 
      x=sqrt((l-r(k))^2+(m-c(k))^2); 
      x=round(x); 
      A(l,m,x)=A(l,m,x)+1; 
     end 
    end 
    end 
end 

[i j k]=find(A~=0) ; 
f=0; 
for s=1:size(i) 
    if(i(s)~=100&&j(s)~=100&&k(s)~=141&&i(s)~=1&&j(s)~=1&&k(s)~=1) 
    if  (A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s))&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s))&&A(i(s),j(s),k(s))>A(i(s),j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s)+1)&&A(i(s),j(s),k(s))>A(i(s),j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s),k(s)-1)&&A(i(s),j(s),k(s))>A(i(s),j(s)-1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s),j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)-1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)-1,j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)+1,k(s)-1)&&A(i(s),j(s),k(s))>A(i(s)+1,j(s)-1,k(s)-1)) 
     f=f+1 
     i(s) 
     j(s) 
     k(s) 
    end 
    end 
end 

Answer 1

не имеют набор инструментов для обработки изображений, поэтому не может легко запустить свой код. Однако, быстрый обзор показывает мне, что вы неправильно используете функцию find. Если у вас есть 3 выхода (как у вас есть), find не возвращает 3D-индексы, а скорее строку, столбец и вектор значений. Это означает, что k - это просто гигантский вектор 1, поэтому ваш for цикл if k(s)~=1 заявление никогда не будет удовлетворен. Вы должны сделать что-то вроде [i,j,k]=ind2sub(size(A),find(A~=0)), если вам нужны 3D-индексы.

Кроме того, FYI, это:

for i=1:rmax 
    for j=1:cmax 
     for k=1:radmax 
      A(i,j,k)=0; 
     end 
    end 
end 

Может быть заменен A=zeros(rmax,cmax,radmax);.