唯思有杰

Melkman算法

const epsilon=0.00000001;
type node=record
          x,y:real
          end;
     nodearr=array[1..20000]of node;

type linklist=^dynanode;
     dynanode=record
              dat:integer;
              pre,next:linklist;
              end;
var p0,p1,p2,p3,head,tail:linklist;
var a:nodearr;
    maxy,miny:real;maxi,mini:integer;
    leftset,rightset,list:nodearr;
    k1,k2,k:longint;
    n:longint;
//////////////////////
procedure tail_push(x:longint);
begin
new(p1);
tail^.next:=p1;
p1^.dat:=x;
p1^.pre:=tail;p1^.next:=nil;
tail:=p1;
end;

function tail_pop:longint;
begin
tail_pop:=tail^.dat;
tail:=tail^.pre;
dispose(tail^.next);
tail^.next:=nil;
end;

procedure head_push(x:longint);
begin
new(p2);
head^.pre:=p2;
p2^.dat:=x;
p2^.pre:=nil;p2^.next:=head;
head:=p2;
end;

function head_pop:longint;
begin
head_pop:=head^.dat;
head:=head^.next;
dispose(head^.pre);
head^.pre:=nil;
end;
///////////////////////

function equal(a,b:real):boolean;
begin
if abs(a-b)<epsilon then exit(true) else exit(false);
end;

function cross_product(x1,y1,x2,y2:real):real;
begin
exit(x1*y2-x2*y1);
end;

function left_side(p00,p11,p22:node):boolean;
var ax,ay,bx,by:real;
    product:real;
    p0,p1,p2:node;
begin
p0:=p00;p1:=p11;p2:=p22;
ax:=p1.x-p0.x;ay:=p1.y-p0.y;
bx:=p2.x-p0.x;by:=p2.y-p0.y;
product:=cross_product(ax,ay,bx,by);
if equal(product,0)and( (ax*ax+ay*ay)>=(bx*bx+by*by)) then exit(false);
if equal(product,0)and( (ax*ax+ay*ay)< (bx*bx+by*by)) then exit(true);
if product>0
then exit(true)
else exit(false);
end;

procedure queue_init;
begin
new(p0);
new(p1);
new(p2);
new(p3);
if left_side(list[1],list[2],list[3])
then begin
     head:=p0;
     p0^.dat:=3;p0^.pre:=nil;p0^.next:=p1;
     p1^.dat:=1;p1^.pre:=p0; p1^.next:=p2;
     p2^.dat:=2;p2^.pre:=p1; p2^.next:=p3;
     p3^.dat:=3;p3^.pre:=p2; p3^.next:=nil;
     tail:=p3;
     end
else begin
     head:=p0;
     p0^.dat:=2;p0^.pre:=nil;p0^.next:=p1;
     p1^.dat:=1;p1^.pre:=p0; p1^.next:=p2;
     p2^.dat:=3;p2^.pre:=p1; p2^.next:=p3;
     p3^.dat:=2;p3^.pre:=p2; p3^.next:=nil;
     tail:=p3;
     end;
end;

procedure qsort(var a:nodearr;l,r:longint);
var i,j:longint;
    mid:node;
/////////
procedure swap(var a,b:node);
var t:node;
begin
t:=a;a:=b;b:=t;
end;
//////////
function cmp(aa,bb:node):shortint;
{a&b's relationship:0 for <; 1 for =;  2 for > }
begin
if equal((aa.y),(bb.y))
then begin
     if equal(aa.x,bb.x) then exit(1);
     if aa.x > bb.x then exit(2) else exit(0);
     end;
if aa.y>bb.y then exit(2) else exit(0);
end;

begin
i:=l;j:=r;mid:=a[(i+j)div 2];
repeat
  while cmp(a[i],mid)=0 do inc(i);
  while cmp(mid,a[j])=0 do dec(j);
  if i<=j
  then begin
       swap(a[i],a[j]);
       inc(i);dec(j);
       end;
until i>j;
if l<j then qsort(a,l,j);
if i<r then qsort(a,i,r);
end;

procedure init;
var i:integer;
begin
fillchar(a,sizeof(a),0);
fillchar(leftset,sizeof(leftset),0);
fillchar(rightset,sizeof(rightset),0);
fillchar(list,sizeof(list),0);

maxy:=-maxlongint;miny:=maxlongint;
maxi:=0;mini:=0;

assign(input,'melkman.in');
assign(output,'melkman.out');
reset(input);rewrite(output);

readln(n);{input data}
for i:=1 to n do
  begin
  readln(a[i].x, a[i].y);
  if a[i].y>maxy then begin maxy:=a[i].y;maxi:=i;end;
  if a[i].y<miny then begin miny:=a[i].y;mini:=i;end;
  end;
{sort}
k1:=0;k2:=0;
for i:=1 to n do
  begin
  if left_side(a[mini],a[maxi],a[i])
  then begin inc(k1);leftset[k1]:=a[i];end
  else begin inc(k2);rightset[k2]:=a[i];end;
  end;
;
qsort(leftset,1,k1);
qsort(rightset,1,k2);
list:=rightset;
k:=k2;
for i:=k1 downto 1 do
  begin
  inc(k);
  list[k]:=leftset[i];
  end;{now, list saves sorted nodes}
end;

procedure melkman;
var i:longint;
    null:longint;
begin
queue_init;
i:=4;
while i<=n do
  begin
  if (left_side(list[(tail^.pre)^.dat ], list[ tail^.dat ] ,list[i])=true)
     and(left_side(list[(head^.next)^.dat] ,list[ head^.dat ] , list[i])=false)
  then begin inc(i);continue;end;
  while (left_side(list[(tail^.pre)^.dat ], list[ tail^.dat ] ,list[i] )=false) do null:=tail_pop;
  tail_push(i);
  while (left_side(list[(head^.next)^.dat] ,list[ head^.dat ] , list[i])=true)  do null:=head_pop;
  head_push(i);
  inc(i);
  end;
end;

procedure print;
begin
p1:=head;
repeat
writeln(p1^.dat,':(',list[p1^.dat].x:0:1,',',list[p1^.dat].y:0:1,')');
p1:=p1^.next;
until p1^.next=nil;
readln;readln;
end;

begin
init;
melkman;
print;
close(input);close(output);
end.
 

posted on 2010-08-07 17:53 陈斯杰 阅读(713) 评论(0)  编辑 收藏 引用 所属分类: 计算几何