Unit momenten;

INTERFACE

Uses Generaal;

type
  tensor = record
    xx,yy,xy : double;
  end;

function Lijn_NuldeOrde(figuur : punten) : double;
function Lijn_EersteOrde(figuur : punten) : punt;
function Lijn_TweedeOrde(figuur : punten) : tensor;

IMPLEMENTATION

function stukje(a,b : punt) : double;
begin
  stukje := sqrt(sqr(b.x-a.x) + sqr(b.y-a.y));
end;

function Lijn_NuldeOrde(figuur : punten) : double;
{
  Length of Contour
}
var
  k,L : integer;
  S,ds : double;
  p,q : punt;
begin
  S := 0;
  L := Length(figuur);
  q := figuur[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := figuur[k];
    ds := stukje(p,q);
    S := S + ds;
  end;
  Lijn_NuldeOrde := S;
end;

function Lijn_EersteOrde(figuur : punten) : punt;
{
  Contour as a Curve
}
var
  k,L : integer;
  ds,S : double;
  p,q,c : punt;
begin
  L := Length(figuur);
  S := Lijn_NuldeOrde(figuur);
{
  Central Moments
}
  c.x := 0;
  c.y := 0;
  q := figuur[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := figuur[k];
    ds := stukje(p,q)/(2*S);
    c.x := c.x + ds*(p.x + q.x);
    c.y := c.y + ds*(p.y + q.y);
  end;
  Lijn_EersteOrde := c;
end;

function Lijn_TweedeOrde(figuur : punten) : tensor;
{
  Contour as a Curve
}
var
  k,L : integer;
  ds,S,t : double;
  rand : punten;
  p,q,c : punt;
  h : tensor;
begin
  c := Lijn_EersteOrde(figuur);
  L := Length(figuur);
  SetLength(rand,L);
  for k := 0 to L-1 do
  begin
    rand[k].x := figuur[k].x - c.x;
    rand[k].y := figuur[k].y - c.y;
  end;
  h.xx := 0;
  h.yy := 0;
  h.xy := 0;
  S := Lijn_NuldeOrde(rand);
  q := rand[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rand[k];
    ds := stukje(p,q);
    t := (p.x*p.x + p.x*q.x + q.x*q.x)/3;
    h.xx := h.xx + t*ds/S;
    t := (p.y*p.y + p.y*q.y + q.y*q.y)/3;
    h.yy := h.yy + t*ds/S;
    t := (p.x*p.y + q.x*q.y)/3 + (p.x*q.y + p.y*q.x)/6;
    h.xy := h.xy + t*ds/S;
  end;
  SetLength(rand,0);
  Lijn_TweedeOrde := h;
end;

END.