Calcular estadísticas de objetos a partir de los segundos momentos centrales.

Question

Actualmente estoy trabajando en escribir una versión de MATLAB. Accesorios de región función para Octava GNU.He implementado la mayor parte, pero todavía estoy luchando con la implementación de algunas partes.Yo tenía preguntado previamente acerca de segundos momentos centrales de una región.

Esto fue útil en teoría, pero tengo problemas para implementar las sugerencias.Obtengo resultados muy diferentes a los de MATLAB (o al sentido común) y realmente no entiendo por qué.

Considere esta imagen de prueba:

Podemos verlo inclinado a 45 grados del eje X, con ejes menor y mayor de 30 y 100 respectivamente.

Ejecutándolo a través de MATLAB RegionProps La función confirma esto:

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

Mientras tanto, ni siquiera acierto con los ejes.estoy tratando de usar estas fórmulas de Wikipedia.

Mi código hasta ahora es:

momentos_raw.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_central.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);

Y aquí está mi código intentando usarlos.Incluyo comentarios para los valores que obtengo en cada paso:

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

No estoy seguro de qué estoy haciendo mal.Parece que estoy siguiendo esas fórmulas correctamente, pero mis resultados no tienen sentido.No he encontrado ningún error obvio en mis funciones de momento, aunque honestamente, para empezar, mi comprensión de los momentos no es la mejor.

¿Alguien puede ver dónde me estoy desviando?Muchas gracias.

Answer 1

EDITAR:

De acuerdo a Wikipedia:

los valores propios [...] son proporcionala la longitud al cuadrado de los ejes del vector propio.

lo cual se explica por:

axisLength = 4 * sqrt(eigenValue)

---

A continuación se muestra mi versión del código (vectoricé las funciones de momentos):

my_regionprops.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

...y el código para probarlo:

prueba.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

¿Estás seguro del núcleo de tu raw_moments ¿función? Podrías intentar

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

Esto no parece suficiente para causar los problemas que está viendo, pero podría ser al menos una parte.