Calculando estatísticas de objetos a partir dos segundos momentos centrais

Question

Atualmente estou trabalhando na escrita de uma versão do MATLAB RegionProps função para Oitava GNU.Já implementei a maior parte, mas ainda estou lutando com a implementação de algumas partes.Eu tive perguntado anteriormente sobre a segundos momentos centrais de uma região.

Isso foi útil teoricamente, mas estou tendo problemas para implementar as sugestões.Obtenho resultados totalmente diferentes dos do MATLAB (ou do bom senso) e realmente não entendo o porquê.

Considere esta imagem de teste:

Podemos vê-lo inclinado a 45 graus em relação ao eixo X, com eixos menor e maior de 30 e 100, respectivamente.

Executando-o através do MATLAB RegionProps função confirma isso:

MajorAxisLength: 101.3362
MinorAxisLength: 32.2961
Eccentricity: 0.9479
Orientation: -44.9480

Entretanto, nem acerto os eixos.estou tentando usar essas fórmulas da Wikipédia.

Meu código até agora é:

raw_moments.m:

function outmom = raw_moments(im,i,j)

  total = 0;
  total = int32(total);
  im = int32(im);

  [height,width] = size(im);

  for x = 1:width;
     for y = 1:height;
        amount = (x ** i) * (y ** j) * im(y,x);
        total = total + amount;
     end;
  end;

  outmom = total;

momentos_centrais.m:

function cmom = central_moments(im,p,q);

  total = 0;
  total = double(total);
  im = int32(im);

  rawm00 = raw_moments(im,0,0);

  xbar = double(raw_moments(im,1,0)) / double(rawm00);
  ybar = double(raw_moments(im,0,1)) / double(rawm00);

  [height,width] = size(im);

  for x = 1:width;
    for y = 1:height;
      amount = ((x - xbar) ** p) *  ((y - ybar) ** q) * double(im(y,x));
      total = total + double(amount);
    end;
  end;

  cmom = double(total);

E aqui está meu código tentando usá-los.Incluo comentários para os valores que recebo em cada etapa:

inim = logical(imread('135deg100by30ell.png'));

cm00 = central_moments(inim,0,0);          % 2567

up20 = central_moments(inim,2,0) / cm00;   % 353.94
up02 = central_moments(inim,0,2) / cm00;   % 352.89
up11 = central_moments(inim,1,1) / cm00;   % 288.31

covmat = [up20, up11; up11, up02];
%[ 353.94  288.31
%  288.31  352.89 ]

eigvals = eig(covmat);          % [65.106 641.730]

minoraxislength = eigvals(1);   % 65.106
majoraxislength = eigvals(2);   % 641.730

Não tenho certeza do que estou fazendo de errado.Parece que estou seguindo essas fórmulas corretamente, mas meus resultados são absurdos.Não encontrei nenhum erro óbvio em minhas funções de momento, embora, honestamente, minha compreensão de momentos não seja a melhor para começar.

Alguém pode ver onde estou me extraviando?Muito obrigado.

Answer 1

EDITAR:

De acordo com Wikipédia:

os autovalores [...] são proporcionalao quadrado do comprimento dos eixos dos autovetores.

que é explicado por:

axisLength = 4 * sqrt(eigenValue)

---

Abaixo é mostrada minha versão do código (vetorizei as funções dos momentos):

minha_regiãoprops.m

function props = my_regionprops(im)
    cm00 = central_moments(im, 0, 0);
    up20 = central_moments(im, 2, 0) / cm00;
    up02 = central_moments(im, 0, 2) / cm00;
    up11 = central_moments(im, 1, 1) / cm00;

    covMat = [up20 up11 ; up11 up02];
    [V,D] = eig( covMat );
    [D,order] = sort(diag(D), 'descend');        %# sort cols high to low
    V = V(:,order);

    %# D(1) = (up20+up02)/2 + sqrt(4*up11^2 + (up20-up02)^2)/2;
    %# D(2) = (up20+up02)/2 - sqrt(4*up11^2 + (up20-up02)^2)/2;

    props = struct();
    props.MajorAxisLength = 4*sqrt(D(1));
    props.MinorAxisLength = 4*sqrt(D(2));
    props.Eccentricity = sqrt(1 - D(2)/D(1));
    %# props.Orientation = -atan(V(2,1)/V(1,1)) * (180/pi);      %# sign?
    props.Orientation = -atan(2*up11/(up20-up02))/2 * (180/pi);
end

function cmom = central_moments(im,i,j)
    rawm00 = raw_moments(im,0,0);
    centroids = [raw_moments(im,1,0)/rawm00 , raw_moments(im,0,1)/rawm00];
    cmom = sum(sum( (([1:size(im,1)]-centroids(2))'.^j * ...
                     ([1:size(im,2)]-centroids(1)).^i) .* im ));
end

function outmom = raw_moments(im,i,j)
    outmom = sum(sum( ((1:size(im,1))'.^j * (1:size(im,2)).^i) .* im ));
end

...e o código para testá-lo:

teste.m

I = imread('135deg100by30ell.png');
I = logical(I);

>> p = regionprops(I, {'Eccentricity' 'MajorAxisLength' 'MinorAxisLength' 'Orientation'})
p = 
    MajorAxisLength: 101.34
    MinorAxisLength: 32.296
       Eccentricity: 0.94785
        Orientation: -44.948

>> props = my_regionprops(I)
props = 
    MajorAxisLength: 101.33
    MinorAxisLength: 32.275
       Eccentricity: 0.94792
        Orientation: -44.948

%# these values are by hand only ;)
subplot(121), imshow(I), imdistline(gca, [17 88],[9 82]);
subplot(122), imshow(I), imdistline(gca, [43 67],[59 37]);

Answer 2

Você tem certeza sobre o núcleo do seu raw_moments função? Você pode tentar

amount = ((x-1) ** i) * ((y-1) ** j) * im(y,x);

Isso não parece o suficiente para causar os problemas que você está vendo, mas pode ser pelo menos uma parte.