Calcolo statistiche oggetto della seconda centrale momenti

Question

Attualmente sto lavorando sulla scrittura di una versione di MATLAB RegionProps funzione per GNU Octave.Ho implementato, ma io sono ancora alle prese con l'attuazione di alcune parti.Ho avuto chiesto in precedenza circa il secondo momenti centrali di una regione.

Questo è stato utile in teoria, ma sto avendo problemi di implementazione delle proposte.I risultati selvaggiamente diverso da MATLAB (o di senso comune per quella materia) e davvero non capisco perché.

Considerare questa immagine di test:

Possiamo vedere inclinato a 45 gradi dall'asse X, con piccoli e grandi assi di 30 e 100, rispettivamente.

Esecuzione attraverso MATLAB s RegionProps funzione di conferma:

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

Nel frattempo, anche non ottenere gli assi a destra.Sto cercando di utilizzare queste formule da Wikipedia.

Il mio codice finora è:

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;

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

Ed ecco il mio codice di tentare di utilizzare questi.Io sono commenti per i valori ho ad ogni passo:

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

Io non sono sicuro di quello che sto facendo di sbagliato.Mi sembra di essere seguito le formule correttamente, ma i miei risultati sono sciocchezze.Non ho trovato errori evidenti nel mio momento di funzioni, anche se onestamente la mia comprensione di momenti non è il massimo per iniziare con.

Chiunque può vedere dove sto andando fuori strada?Vi ringrazio molto.

Answer 1

EDIT:

Secondo Wikipedia:

il eignevalues [...] sono proporzionale al quadrato della lunghezza dell'autovettore assi.

che viene spiegata da:

axisLength = 4 * sqrt(eigenValue)

---

Qui di seguito è mostrata la mia versione del codice (ho vettorializzare i momenti funzioni):

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

...e il codice per il test:

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

Sei sicuro circa il nucleo della vostra funzione raw_moments? Si potrebbe provare

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

Questa non sembra sufficiente a causare i problemi che stai vedendo, ma potrebbe essere almeno una parte.