$ debstats
----- Priority Size(Mo)	Packages
      standard 145	87
      required 48	50
         extra 23	19
      optional 1155	744
     important 31	60
------ Section Size(Mo)	Packages
         games 3.75	4
     otherosfs 8	6
       oldlibs 2.44	10
           net 17	40
           web 46	11
  interpreters 6	4
           x11 157	46
  contrib/misc 0.07	1
          text 56	33
           doc 58	28
         devel 111	42
      libdevel 94	81
          mail 7	9
         sound 10	15
          base 109	56
         admin 27	27
        python 30	23
          math 28	5
         utils 22	38
          libs 218	349
        shells 1.41	2
          perl 28	25
      graphics 60	17
         gnome 232	62
          misc 19	17
       editors 41	8
non-free/utils 0.17	1
-------- Total 1403	960
	

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

.: GNU / Plus de sources pour les développeurs :.