// This unit give several features to work on polygons such as
// found particular points of the polygon.
// Copyright (C) 2005 Pierre Chignac aka peyo
// This program is free software; you can redistribute it and/or 
// modify it under the terms of the GNU General Public License as 
// published by the Free Software Foundation; either version 2 of 
// the License, or (at your option) any later version. This program
// is distributed in the hope that it will be useful, but WITHOUT 
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 
// License for more details. You should have received a copy of the GNU 
// General Public License along with this program; if not, write to the 
// Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.


unit Polygones;

interface

uses
  Classes, SysUtils, Math;

type
  _Ppoint = ^poly_point;

  poly_point = record
    x,y : real;
    next : _Ppoint;
    island : integer;
  end;

  TPolygone = class(TObject)

  public
    id : string;
    count : integer;

  private
    first_item : _Ppoint;
    last_item: _Ppoint;
    maxx, maxy, minx, miny : real;

  public
    constructor Create;
    destructor Destroy;override;
	// ajoute un point au polygone
    procedure AddItem(x,y:real);
	// calcule la surface du ploygone
    function Surface : real;
	// trouvre les coordonées du centre de gravité
    function centre_gravite : poly_point;
	// pareil pour l'isobarycentre
    function isoberycentre : poly_point;
	// et le point isolé
    function point_isole(epsilon : integer) : poly_point;
	// retire tous les points d'un polygone
    procedure Clear;


  private
    function Item(index : integer) : _Ppoint;
    // renvoie l'abscisse du ième point
    function x(index : integer) : real;
    // renvoie l'ordonnée du ième point
    function y(index : integer) : real;
    //inverse la description du polygone par ses points
    procedure Reverse;
    //renvoie 1 si le polynome est décrit dans le sens trigo positif, 0 sinon
    function sens_trigo : boolean;
    //renvoie 1 si p est dans le polynome, 0 sinon
    function inclus(p : poly_point) : boolean;
    //renvoie la distance de p a un polygone;
    function distance(p : poly_point) : real;
    procedure DeleteItem(index : integer);
  end;


implementation


//==========================================================================

//       renvoie la disantance Euclidienne entre deux points A & B

function distance_point_point(xa,ya,xb,yb : real) : real;
begin
  result:=sqrt(sqr(xa-xb)+sqr(ya-yb));
end;

//==========================================================================

function distance_point_segment(xp,yp,xa,ya,xb,yb : real) : real;
var
  tmp, xo, yo, lambda : real;
begin
  Result:=-1;
  //--------------------segment vertical------------------
  if xa=xb then
  begin
    // 1 - le point est sur la droite
    if xp=xa then
    begin
      if ((yp<=max(ya,yb)) and (yp>=min(ya,yb))) then result:=0
        else result:=min(abs(yp-ya),abs(yp-yb));
    exit;
    end
    // 2 - le point n'est pas sur la droite
    else
    begin
      if ((yp<=max(ya,yb)) and (yp>=min(ya,yb))) then result:=abs(xa-xp)
        else result:=min(distance_point_point(xp,yp,xa,ya),distance_point_point(xp,yp,xb,yb));
    exit;
    end;
  end;//seg vert

  //--------------------segment horizontal------------------

  if ya=yb then
  begin
    // 1 - le point est sur la droite
    if yp=ya then
    begin
      if ((xp<=max(xa,xb)) and (xp>=min(xa,xb))) then result:=0
        else result:=min(abs(xp-xa),abs(xp-xb));
    exit;
    end
    // 2 - le point n'est pas sur la droite
    else
    begin
      if ((xp<=max(xa,xb)) and (xp>=min(xa,xb))) then result:=abs(ya-yp)
        else result:=min(distance_point_point(xp,yp,xa,ya),distance_point_point(xp,yp,xb,yb));
    exit;
    end;
  end;//seg hor

  //--------------------segment quelconque------------------

  if xa>xb then
  begin
    tmp:=xa;xa:=xb;xb:=tmp;
    tmp:=ya;ya:=yb;yb:=tmp;
  end;
  //lambda = coefficient directeur de la droite (AB)
  lambda:=(ya-yb)/(xa-xb);
  //  O de coordonnées (xo,yo) est le projeté orthogonal
  //  de P sur la droite (AB)
  xo:=(lambda*xa-ya+1/lambda*xp+yp)/(lambda+1/lambda);
  yo:=-1/lambda*(xo-xp)+yp;
  //si O est sur le segment [AB] alors la distance de P à [AB] est OP
  if ((xo>xa)and(xo<xb)) then result:=distance_point_point(xp,yp,xo,yo)
    //sinon c'est la plus petite distance de P aux extrémités du segment
    else result:=min(distance_point_point(xp,yp,xa,ya),distance_point_point(xp,yp,xb,yb));

  //-----------------------------------------------------------

end;

//==========================================================================

constructor TPolygone.Create;
begin
  inherited create;
  count:=0;
  first_item:=Nil;
  last_item:=Nil;
end;

//==========================================================================

destructor TPolygone.Destroy;
var
  i : integer;
  tmp_item : _Ppoint;
  del_item : _Ppoint;
begin
  if count > 0 then
  begin
    tmp_item := first_item.next;
    del_item:=first_item;
    for i := 1 to count-1 do
    begin
      freemem(del_item);
      del_item:=tmp_item;
      tmp_item:=tmp_item.next;
    end;
    freemem(tmp_item);
    count:=0;
    first_item:=Nil;
    last_item:=Nil;
  end;
  inherited Destroy;
end;

//==========================================================================

function TPolygone.Surface : real;
var
  i : integer;
  mx, my : real;
  tmp_item : _Ppoint;
begin
  tmp_item := first_item;
  mx:=minx;
  my:=miny;
  for i := 0 to count-1 do
  begin
    tmp_item.x:=tmp_item.x-mx+10;
    tmp_item.y:=tmp_item.y-my+10;
    if tmp_item.next<>Nil then tmp_item:=tmp_item.next;
  end;
  tmp_item:=first_item;
  case count of
    1 : Result := -1;
    2 : Result := -2;
    else
    begin
      Result:=0;
      repeat
      begin
        Result:=Result+
                ((tmp_item.x-minx)+(tmp_item.next.x-minx))*
                ((tmp_item.next.y-miny)-(tmp_item.y-miny));
        tmp_item:=tmp_item.next
      end;
      until tmp_item.next=Nil;
      Result:=Result/2;
    end;
  end;
  tmp_item := first_item;
  for i := 0 to count-1 do
  begin
    tmp_item.x:=tmp_item.x+mx-10;
    tmp_item.y:=tmp_item.y+my-10;
    if tmp_item.next<>Nil then tmp_item:=tmp_item.next;
  end;
end;

//==========================================================================

function TPolygone.isoberycentre : poly_point;
var
  i : integer;
begin
  result.x:=0;
  result.y:=0;
  for i := 1 to count-1 do
  begin
    result.x:=result.x+x(i-1);
    result.y:=result.y+y(i-1);
  end;
  result.x:=result.x/(count-1);
  result.y:=result.y/(count-1);
end;

//==========================================================================

function TPolygone.inclus(p : poly_point) : boolean;
var
  i, n_coupes :integer;
  xa, ya, xb, yb: real;
  _A, _B : _Ppoint;
begin
  result:=false;
  _A:=first_item;
  n_coupes:=0;
  for i := 0 to count-2 do
  begin
    _B:=_A.next;
    //mettre A à gauche et B à droite
    if _A.x<_B.x then
      begin
        xa:=_A.x;  ya:=_A.y; xb:=_B.x; yb:=_B.y;
      end
      else
      begin
        xa:=_B.x;  ya:=_B.y; xb:=_A.x; yb:=_A.y;
      end;
    //vérifier l'intersection

      if ((xa=p.x) or (xb=p.x)) then p.x:=p.x+p.x/10000000000;

      if ( (xa<>xb) and (xa<p.x) and (xb>p.x) and ((ya-yb)/(xa-xb)*(p.x-xb)+yb > p.y) ) then
          inc(n_coupes);

    //initialiser le segment suivant
    _A:=_A.next;
  end;
  if n_coupes mod 2 = 0 then result:=false else result:=true;
end;

//==========================================================================

function TPolygone.distance(p : poly_point) : real;
var
  i:integer;
  tmp:real;
  _A, _B : _Ppoint;
begin
  _A:=first_item;
  _B:=_A.next;
  result:=distance_point_segment(p.x,p.y,_A.x,_A.y,_B.x,_B.y);
  for i := 0 to count-2 do
  begin
    _B:=_A.next;
    tmp:=distance_point_segment(p.x,p.y,_A.x,_A.y,_B.x,_B.y);
    if tmp < result then result:=tmp;
    _A:=_A.next;
  end;
end;

//==========================================================================

function TPolygone.point_isole(epsilon : integer) : poly_point;
var
  i, j : integer;
  meilleur_distance : real;
  tmp : real;
  test : poly_point;
  tutu, titi : tstringlist;
begin
// tutu:=tstringlist.create;
//  titi:=tstringlist.create;
//  tutu.loadfromfile('D:/home/Delphi_Projects/polygones/data/ilots_insee/ilots_light/test.mif');
//  titi.loadfromfile('D:/home/Delphi_Projects/polygones/data/ilots_insee/ilots_light/test.mid');
  meilleur_distance:=0;
  for i := 1 to epsilon-1 do
  begin
    test.x:=minx+i*((maxx-minx)/epsilon);
    for j := 1 to epsilon-1 do
    begin
      test.y:=miny+j*((maxy-miny)/epsilon);
      if inclus(test) then
      begin
        tmp:=distance(test);
//        tutu.Add('point '+floattostr(test.x)+' '+floattostr(test.y));
//        titi.add('"'+floattostr(tmp)+'"');
        if tmp>meilleur_distance then
        begin
          meilleur_distance:=tmp;
          result.x:=test.x;
          result.y:=test.y;
        end;
      end;
    end;
  end;
//  tutu.savetofile('D:/home/Delphi_Projects/polygones/data/ilots_insee/ilots_light/test.mif');
//  titi.savetofile('D:/home/Delphi_Projects/polygones/data/ilots_insee/ilots_light/test.mid');
end;

//==========================================================================

function TPolygone.centre_gravite: poly_point;
var
  i : integer;
  p : poly_point;
  mx, my : real;
  tmp_item : _Ppoint;
begin

  tmp_item := first_item;
  p.x:=0;
  p.y:=0;
  mx:=minx;
  my:=miny;
  for i := 0 to count-1 do
  begin
    tmp_item.x:=tmp_item.x-mx+10;
    tmp_item.y:=tmp_item.y-my+10;
    if tmp_item.next<>Nil then tmp_item:=tmp_item.next;
  end;
  for i := 1 to count-1 do
  begin
    p.x:=p.x+(y(i)-y(i-1))*(x(i)*x(i)+x(i)*x(i-1)+x(i-1)*x(i-1));
    p.y:=p.y+(x(i-1)-x(i))*(y(i)*y(i)+y(i)*y(i-1)+y(i-1)*y(i-1));
  end;
  p.x:=p.x/(6*surface);
  p.y:=p.y/(6*surface);
  if (p.x<0) and (p.y<0) then
  begin
    p.x:=p.x*(-1);
    p.y:=p.y*(-1);
  end;
  p.x:=p.x+mx-10;
  p.y:=p.y+my-10;
  tmp_item := first_item;
  for i := 0 to count-1 do
  begin
    tmp_item.x:=tmp_item.x+mx-10;
    tmp_item.y:=tmp_item.y+my-10;
    if tmp_item.next<>Nil then tmp_item:=tmp_item.next;
  end;
  result:=p;
end;

//==========================================================================

function TPolygone.x(index : integer) : real;
var
  i : integer;
  tmp_item : _Ppoint;
begin
  tmp_item := first_item;
  if index > 0 then
    for i := 1 to index do tmp_item:=tmp_item.next;
  result:=tmp_item.x;
end;

//==========================================================================

function TPolygone.y(index : integer) : real;
var
  i : integer;
  tmp_item : _Ppoint;
begin
  tmp_item := first_item;
  if index > 0 then
    for i := 1 to index do tmp_item:=tmp_item.next;
  result:=tmp_item.y;
end;

//==========================================================================

procedure TPolygone.Clear;
var
  i : integer;
begin
  for i := 0 to count-1 do DeleteItem(0);
end;

//==========================================================================

procedure TPolygone.AddItem(x,y:real);
var
  tmp_item : _Ppoint;
begin
  tmp_item:=allocmem(sizeof(poly_point));
  tmp_item.x:=x;
  tmp_item.y:=y;
  tmp_item.next:=Nil;
  if count=0 then
  begin
    minx:=x;miny:=y;maxx:=x;maxy:=y;
    first_item:=tmp_item
  end
  else
  begin
    if x<minx then minx:=x;
    if x>maxx then maxx:=x;
    if y<miny then miny:=y;
    if y>maxy then maxy:=y;
    last_item.next:=tmp_item;
  end;
  last_item:=tmp_item;
  inc(count);
end;

//==========================================================================

procedure TPolygone.DeleteItem(index : integer);
var
  i : integer;
  tmp_item, tmp_item2 : _Ppoint;
begin
  tmp_item:=first_item;
  case count of
  0 : exit;
  1 : begin
        if index =0 then
        begin
          first_item:=nil;
          last_item:=nil;
        end;
      end;//fin case count =1
  2 : begin
        case index of
          0 : first_item:=last_item;
          1 : begin
                tmp_item:=last_item;
                last_item:=first_item;
                first_item.next:=nil;
              end;
        end;
      end;//fin case count = 2
  else
  begin
    case index of
    0 : first_item:=first_item.next;
    1 : begin
          tmp_item:=first_item.next;
          first_item.next:=first_item.next.next;
        end;
    else
    begin
      for i := 0 to index-2 do tmp_item:=tmp_item.next;
      if tmp_item.next.next=nil then
      begin
        last_item:=tmp_item;
        tmp_item:=tmp_item.next;
        last_item.next:=nil;
      end
      else
      begin
        tmp_item2:=tmp_item.next;
        tmp_item.next:=tmp_item.next.next;
        tmp_item:=tmp_item2;
      end;
    end;
    end;
  end;//fin case count else
  end;
  freemem(tmp_item);
  dec(count);
end;

//==========================================================================

function TPolygone.Item(index : integer) : _Ppoint;
var
  i : integer;
  tmp_item : _Ppoint;
begin
  tmp_item := first_item;
  if index > 0 then
    for i := 1 to index do tmp_item:=tmp_item.next;
  result:=tmp_item;
end;

//==========================================================================

function TPolygone.sens_trigo : boolean;
begin
  if surface>0 then Result := True else Result:=False;
end;

//==========================================================================

procedure TPolygone.Reverse;
var
  i : integer;
  tmp_x, tmp_y : real;
begin
  for i := 0 to floor(count/2)-1 do
  begin
    tmp_x:=x(i);
    tmp_y:=y(i);
    Item(i).x:=x(count-i-1);
    Item(i).y:=y(count-i-1);
    Item(count-i-1).x:=tmp_x;
    Item(count-i-1).y:=tmp_y;
  end;
end;

//==========================================================================

end.