/********************** CONVEX HULL  ***************

Author: PD Dr. Siegfried Gabler
	ZUMA
	Centre for Survey Research and Methodology
	Mannheim
	Germany
	phone int-49-621-1246-281
	fax   int-49-621-1246-100
	e-mail gabler@zuma-mannheim.de
Date: Tue, 18 Feb 1997

Below there is a procedure CHUELLE(x) which computes the vertices of a
cloud of points in the plane. CHUELLEG(x,i) is  for  graphics . An
example is included. Sorry that parts are written in German. The
procedures are part of my program on "Correspondence Analysis".
The code is provided without performance guarantees
for public, non-commercial use.
*****************************************************/

/*
An example:
*/
library pgraph;
graphset;
_plctrl=-1;
_psymsiz=1;

_ctype=6;   /* 2th column of _pline */
_ccolor=12; /* 8th column of _pline */
_cthick=5;  /* 9th column of _pline */

x=rndn(1000,2);
i=400�600;
call chuelleg(x,i);
xy(x�.,1!,x�.,2!);


/*
---------------------------------------------------------------------
*/

/* CHUELLE.SRC 25.8.1993/12.12.95  Siegfried Gabler

Purpose: Computes the vertices of  n  points

Format:  i= chuelle(x)

Input:   x -- nx2 matrix; rows are the coordinates of the points

Output:  i -- Kx1 vector; number of rows of x being vertices

Remark:  Jarvi's procedure;

References: Cormen/Leiserson/Rivest (1990). Introduction to Algorithms
            MIT Press S.905
----------------------------------------------------------------------
-- */ proc chuelle(x); local
ind,n,s,imin,xmin1,xmin2,xmin3,xmin,xmax,xmax1,xmax2,xmax3,
rmin,rmax,j,jmin,jmax,j1,j2,emin,emax,e1;

n=rows(x);
if x==meanc(x)';
    retp(1);
else;
s=seqa(1,1,n);

ind=sortmc(x~s,2�1);
x=ind�.,1 2!;
ind=ind�.,3!;
imin=1;

j=1;
j1=imin;
j2=imin;
do until j==0;                         /* j:Zeilennummer von x */
    xmin=x�j1,.!;
    xmax=x�j2,.!;
    xmin1=x-xmin;
    xmax1=x-xmax;
    xmin2=delif(xmin1~s, (sumc( (xmin1.==0)' ).>1));
    xmax2=delif(xmax1~s, (sumc( (xmax1.==0)' ).>1));
    xmin3=trimr(xmin2',0,1)';
    xmax3=trimr(xmax2',0,1)';
    emin=(xmin3�.,2!.<0); /* untere Quadranten */
    emax=(xmax3�.,2!.>0); /* obere Quadranten */
    rmin= arccos( xmin3�.,1!./sqrt(sumc(xmin3'^2)) );
    rmax= arccos( -xmax3�.,1!./sqrt(sumc(xmax3'^2)) );
    rmin=rmin+emin.*(2*pi-2*rmin);
    rmax=rmax+emax.*(2*pi-2*rmax);

    jmin=minindc(rmin);
    e1=(rmin�jmin!.==rmin);   /* Falls mehrere Punkte auf Kanten
    liegen */ if sumc(e1)>1;
        jmin=maxindc( sumc( ((xmin3.*e1)^2 )' ) );
    endif;

    jmax=maxindc(rmax);
    e1=(rmax�jmax!.==rmax);   /* Falls mehrere Punkte auf Kanten
    liegen */ if sumc(e1)>1;
        jmax=maxindc( sumc( ((xmax3.*e1)^2 )' ) );
    endif;

    j1=xmin2�jmin,3!;
    j2=xmax2�jmax,3!;
    if (j1==j2);
        j=0;
        imin=imin�j1;
    elseif (j1/=j2) and (j1==imin�1!);
        j=0;
    elseif (j1/=j2) and (j2==imin�rows(imin)!);
        j=0;
    else;
        if (sumc(j1.==imin)>0) and (sumc(j2.==imin)==0);
            imin=j2�imin;
        elseif (sumc(j1.==imin)==0) and (sumc(j2.==imin)>0);
            imin=imin�j1;
        elseif (sumc(j1.==imin)>0) and (sumc(j2.==imin)>0);
            j=0;
        else;
            imin=j2�imin�j1;
        endif;
    endif;

endo;

retp(ind�imin!);
endif;
endp;
/*
----------------------------------------------------------------------
-- */

/* CHUELLEG.SRC 25.8.1993

Purpose: Connects the vertices of point clouds by lines;
         lines are saved in the global variable _PLINE

Format:  chuelleg(x,i)

Input:   x -- nx2 Matrix; rows are the coordinates of the points
         i -- 0 oder Kx1 vector;
              partition of the K clouds of points in the sequel of x
             0  means: only one cloud of points (SUMC(i)=n)

----------------------------------------------------------------------
----- */

proc(0)=chuelleg(x,n);
local nc,i,x1,imin1,i1,o1,z1,pl;

if n==0;             /* Nur eine konvexe Huelle */
    nc=0�rows(x);
else;
    if sumc(n)==rows(x);
        nc=0�cumsumc(n);
    else;
        "In Prozedur CHUELLEG(Arg1,Arg2)";
        "ist Zahl der Zeilen von Arg1 ungleich SUMC(Arg2)";
        end;
    endif;
endif;

i=1;
do until i>rows(n);
    x1=x�nc�i!+1:nc�i+1!,.!;
    imin1=chuelle(x1);
    i1=rotater(imin1',-1)';
    o1=ones(rows(imin1),1);
    z1=_cthick*ones(rows(imin1),1);
    pl=o1~_ctype*o1~x1�imin1,1!~x1�imin1,2!~x1�i1,1!~x1�i1,2!~z1;
    pl=pl~_ccolor*o1~shiftr(z1',1,0)'; if _pline==0;
        _pline=pl;
    else;
        _pline=_pline�pl;
    endif;
    i=i+1;
endo;

endp;
/*
----------------------------------------------------------------------
*/

