00001 #include "cs.h"
00002 
00003 /*----------------------------------------------------
00004   Compute the convex hull of a bunch of 3-vectors
00005   Inputs:
00006     npt = number of vectors
00007     xyz = array of coordinates of 3-vectors;
00008           the i-th vector is stored in
00009             xyz[3*i] xyz[3*i+1] xyz[3*i+2]
00010   Output:
00011     *ijk = pointer to malloc()-ed array of triangles;
00012            the j-th triangle is stored in
00013              ijk[3*j] ijk[3*j+1] ijk[3*j+2]
00014            where the integer index i refers to the
00015            i-th 3-vector input
00016 
00017   Return value is the number of triangles.  If this
00018   is zero, something bad happened.
00019 
00020   Example:
00021     int ntri , *tri , nvec ;
00022     float vec[something] ;
00023     ntri = qhull_wrap( nvec , vec , &tri ) ;
00024 
00025   This function just executes the Geometry Center
00026   program qhull to compute the result.  This program
00027   should be in the user's path, or this function
00028   will fail (return 0).
00029 ------------------------------------------------------*/
00030 
00031 int qhull_wrap( int npt , float * xyz , int ** ijk )
00032 {
00033    int ii,jj , nfac , *fac ;
00034    int fd ; FILE *fp ;
00035    char qbuf[128] ;
00036 
00037 #ifndef DONT_USE_MKSTEMP
00038    char fname[] = "/tmp/afniXXXXXX" ;
00039 #else
00040    char *fname ;
00041 #endif
00042 
00043    if( npt < 3 || xyz == NULL || ijk == NULL ){
00044       fprintf(stderr,"qhull_wrap: bad inputs\n") ;
00045       return 0 ;
00046    }
00047 
00048 #ifndef DONT_USE_MKSTEMP
00049    fd = mkstemp( fname ) ;
00050    if( fd == -1 ){ fprintf(stderr,"qhull_wrap: mkstemp fails\n"); return 0; }
00051    fp = fdopen( fd , "w" ) ;
00052    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: fdopen fails\n"); close(fd); return 0; }
00053 #else
00054    fname = tmpnam(NULL) ;
00055    if( fname == NULL ){ fprintf(stderr,"qhull_wrap: tmpnam fails\n"); return 0; }
00056    fp = fopen( fname , "w" ) ;
00057    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: fopen fails\n"); return 0; }
00058 #endif
00059 
00060    fprintf(fp,"3\n%d\n",npt) ;
00061    for( ii=0 ; ii < npt ; ii++ )
00062       fprintf(fp,"%g %g %g\n",xyz[3*ii],xyz[3*ii+1],xyz[3*ii+2]) ;
00063 
00064    fclose(fp) ;
00065 
00066    sprintf(qbuf,"qhull -i -Pp < %s",fname) ;
00067    fp = popen( qbuf , "r" ) ;
00068    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: popen fails\n"); remove(fname); return 0; }
00069 
00070    jj = fscanf(fp,"%d",&nfac) ;
00071    if( jj != 1 || nfac < 1 ){ fprintf(stderr,"qhull_wrap: 1st fscanf fails\n"); pclose(fp); remove(fname); return 0; }
00072 
00073    fac = (int *) malloc( sizeof(int)*3*nfac ) ;
00074    if( fac == NULL ){ fprintf(stderr,"qhull_wrap: malloc fails\n"); pclose(fp); remove(fname); return 0; }
00075 
00076    for( ii=0 ; ii < nfac ; ii++ ){
00077       jj = fscanf(fp,"%d %d %d",fac+(3*ii),fac+(3*ii+1),fac+(3*ii+2)) ;
00078       if( jj < 3 ){
00079          fprintf(stderr,"qhull_wrap: fscanf fails at ii=%d\n",ii) ;
00080          pclose(fp); remove(fname); free(fac); return 0;
00081       }
00082    }
00083 
00084    pclose(fp); remove(fname);
00085 
00086    *ijk = fac ; return nfac ;
00087 }
00088 
00089 /*------------------------------------------------------------------
00090    Compute the Voronoi areas for a collection of points on the
00091    surface of the unit sphere:
00092      npt    = number of points
00093      pol[i] = polar angle (co-latitude) of i-th point (radians)
00094      azi[i] = azimuthal angle (longitude) of i-th point (radians)
00095      *wt    = malloc()-ed array holding the output; the area attached
00096               to the i-th input point is in wt[i]
00097    Return value is 0 if an error transpired, or 1 if all went OK.
00098    If return is 0, then *wt is not changed.
00099 --------------------------------------------------------------------*/
00100 
00101 int sphere_voronoi_angles( int npt , float *pol , float *azi , float **wt )
00102 {
00103    float *xyz ;
00104    double cp,sp,ca,sa ;
00105    int ii ;
00106 
00107    if( npt < 3 || pol == NULL || azi == NULL || wt == NULL ){
00108       fprintf(stderr,"sphere_voronoi_angles: bad inputs\n"); return 0;
00109    }
00110 
00111    /* make 3-vectors of the points on the sphere */
00112 
00113    xyz = (float *) malloc( sizeof(float) * (3*npt) ) ;
00114 
00115    for( ii=0 ; ii < npt ; ii++ ){
00116       cp = cos(pol[ii]) ; sp = sin(pol[ii]) ;
00117       ca = cos(azi[ii]) ; sa = sin(azi[ii]) ;
00118       xyz[3*ii]   = sp * ca ;
00119       xyz[3*ii+1] = sp * sa ;
00120       xyz[3*ii+2] = cp ;
00121    }
00122 
00123    ii = sphere_voronoi_vectors( npt , xyz , wt ) ;
00124    free(xyz) ; return ii ;
00125 }
00126 
00127 /*------------------------------------------------------------------------
00128    Same as above, but points on surface of sphere are specified by
00129    unit 3-vectors; i-th vector is in xyz[3*ii], xyz[3*ii+1], xyz[3*ii+2]
00130    N.B.: if the 3-vectors are NOT unit vectors, the results of this
00131          routine are erroneous!
00132 --------------------------------------------------------------------------*/
00133 
00134 int sphere_voronoi_vectors( int npt , float *xyz , float **wt )
00135 {
00136    int ntri , *tri , ii,jj , pp,qq,rr ;
00137    float *ww ;
00138    double cp,sp,ca,sa ;
00139    double xpq,ypq,zpq , xpr,ypr,zpr , xqr,yqr,zqr , xcc,ycc,zcc ;
00140    double xpp,ypp,zpp , xqq,yqq,zqq , xrr,yrr,zrr , xnn,ynn,znn ;
00141    double pp_pq , pp_pr , pp_cc ,
00142           qq_pq , qq_qr , qq_cc ,
00143           rr_qr , rr_cc , rr_pr ,
00144           pq_cc , qr_cc , pr_cc , ss ;
00145    double a_pp_pq_cc , a_pp_pr_cc ,
00146           a_qq_pq_cc , a_qq_qr_cc ,
00147           a_rr_qr_cc , a_rr_pr_cc  ;
00148 
00149    if( npt < 3 || xyz == NULL || wt == NULL ){
00150       fprintf(stderr,"sphere_voronoi: bad inputs\n"); return 0;
00151    }
00152 
00153    /* compute convex hull triangular facets */
00154 
00155    ntri = qhull_wrap( npt , xyz , &tri ) ;
00156    if( ntri == 0 ){ fprintf(stderr,"sphere_voronoi: qhull_wrap fails\n"); free(xyz); return 0; }
00157 
00158    /* initialize output */
00159 
00160    ww = (float *) malloc( sizeof(float) * npt ) ;
00161    for( ii=0 ; ii < npt ; ii++ ) ww[ii] = 0.0 ;
00162 
00163    for( jj=0 ; jj < ntri ; jj++ ){  /* loop over triangles */
00164 
00165 
00166       /* triangle vertices pp,qq,rr      * pp
00167                                         / \
00168                                        /   \
00169                                    pq *     * pr
00170                                      /   *   \
00171                                     /   cc    \
00172                                 qq *-----------* rr   */
00173 
00174       pp  = tri[3*jj  ] ;
00175       xpp = xyz[3*pp  ] ; ypp = xyz[3*pp+1] ; zpp = xyz[3*pp+2] ;
00176       qq  = tri[3*jj+1] ;
00177       xqq = xyz[3*qq  ] ; yqq = xyz[3*qq+1] ; zqq = xyz[3*qq+2] ;
00178       rr  = tri[3*jj+2] ;
00179       xrr = xyz[3*rr  ] ; yrr = xyz[3*rr+1] ; zrr = xyz[3*rr+2] ;
00180 
00181       /* midpoints pq,pr,qr, and centroid cc */
00182       /*** Q: should centroid be replaced by normal?
00183               what about orientation, largeness?     ***/
00184 
00185       xpq = 0.5*(xpp+xqq) ; ypq = 0.5*(ypp+yqq) ; zpq = 0.5*(zpp+zqq) ;
00186       xpr = 0.5*(xpp+xrr) ; ypr = 0.5*(ypp+yrr) ; zpr = 0.5*(zpp+zrr) ;
00187       xqr = 0.5*(xqq+xrr) ; yqr = 0.5*(yqq+yrr) ; zqr = 0.5*(zqq+zrr) ;
00188 
00189       xcc = 0.3333333*(xpp+xqq+xrr) ;
00190       ycc = 0.3333333*(ypp+yqq+yrr) ;
00191       zcc = 0.3333333*(zpp+zqq+zrr) ;
00192 
00193 #undef SCL
00194 #define SCL(a,b,c) 1.0/sqrt(a*a+b*b+c*c)
00195 
00196 #undef USE_NORMAL
00197 #ifdef USE_NORMAL
00198 # define XCROSS(a,b,c,x,y,z) ((b)*(z)-(c)*(y))
00199 # define YCROSS(a,b,c,x,y,z) ((c)*(x)-(a)*(z))
00200 # define ZCROSS(a,b,c,x,y,z) ((a)*(y)-(b)*(x))
00201       { double apq=xpp-xqq , bpq=ypp-yqq , cpq=zpp-zqq ,
00202                aqr=xqq-xrr , bqr=yqq-yrr , cqr=zqq-zrr  ;
00203 
00204         xnn = XCROSS(apq,bpq,cpq,aqr,bqr,cqr) ,
00205         ynn = YCROSS(apq,bpq,cpq,aqr,bqr,cqr) ,
00206         znn = ZCROSS(apq,bpq,cpq,aqr,bqr,cqr)  ;
00207 
00208         cp = SCL(xnn,ynn,znn) ; xnn *= cp ; ynn *= cp ; znn *= cp ;
00209         if( xnn*xcc + ynn*ycc + znn*zcc < 0 ){
00210            xnn = -xnn ; ynn = -ynn ; znn = -znn ;
00211         }
00212       }
00213 
00214 # define xVV xnn
00215 # define yVV ynn
00216 # define zVV znn
00217 
00218 #else
00219 
00220 # define xVV xcc
00221 # define yVV ycc
00222 # define zVV zcc
00223 
00224 #endif
00225 
00226       /* project pq,pr,qr,cc to sphere (nn is already on sphere) */
00227 
00228       cp = SCL(xpq,ypq,zpq) ; xpq *= cp ; ypq *= cp ; zpq *= cp ;
00229       cp = SCL(xpr,ypr,zpr) ; xpr *= cp ; ypr *= cp ; zpr *= cp ;
00230       cp = SCL(xqr,yqr,zqr) ; xqr *= cp ; yqr *= cp ; zqr *= cp ;
00231       cp = SCL(xcc,ycc,zcc) ; xcc *= cp ; ycc *= cp ; zcc *= cp ;
00232 
00233 #undef  ANG
00234 #define ANG(u1,u2,u3,v1,v2,v3) acos(u1*v1+u2*v2+u3*v3)
00235 
00236       /* compute distance on surface between points:
00237          aa_bb = between points aa and bb, from the picture above */
00238 
00239       pp_pq = ANG( xpp,ypp,zpp , xpq,ypq,zpq ) ;
00240       pp_cc = ANG( xpp,ypp,zpp , xVV,yVV,zVV ) ;
00241       pp_pr = ANG( xpp,ypp,zpp , xpr,ypr,zpr ) ;
00242 
00243       qq_pq = ANG( xqq,yqq,zqq , xpq,ypq,zpq ) ;
00244       qq_qr = ANG( xqq,yqq,zqq , xqr,yqr,zqr ) ;
00245       qq_cc = ANG( xqq,yqq,zqq , xVV,yVV,zVV ) ;
00246 
00247       rr_qr = ANG( xrr,yrr,zrr , xqr,yqr,zqr ) ;
00248       rr_pr = ANG( xrr,yrr,zrr , xpr,ypr,zpr ) ;
00249       rr_cc = ANG( xrr,yrr,zrr , xVV,yVV,zVV ) ;
00250 
00251       pq_cc = ANG( xpq,ypq,zpq , xVV,yVV,zVV ) ;
00252       qr_cc = ANG( xqr,yqr,zqr , xVV,yVV,zVV ) ;
00253       pr_cc = ANG( xpr,ypr,zpr , xVV,yVV,zVV ) ;
00254 
00255       /* for each vertex,
00256          compute areas of 2 subtriangles it touches,
00257          add these into the area weight for that vertex */
00258 
00259 #undef  SS
00260 #define SS(a,b,c) ((a+b+c)/2)
00261 #undef  ATR
00262 #define ATR(s,a,b,c) (4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2))))
00263 
00264       ss      = SS(pp_pq,pp_cc,pq_cc) ;     /* subtriangle pp,pq,cc */
00265       ww[pp] += a_pp_pq_cc = ATR(ss,pp_pq,pp_cc,pq_cc) ;
00266 
00267       ss      = SS(pp_pr,pp_cc,pr_cc) ;     /* subtriangle pp,pr,cc */
00268       ww[pp] += a_pp_pr_cc = ATR(ss,pp_pr,pp_cc,pr_cc) ;
00269 
00270       ss      = SS(qq_pq,qq_cc,pq_cc) ;     /* subtriangle qq,pq,cc */
00271       ww[qq] += a_qq_pq_cc = ATR(ss,qq_pq,qq_cc,pq_cc) ;
00272 
00273       ss      = SS(qq_qr,qq_cc,qr_cc) ;     /* subtriangle qq,qr,cc */
00274       ww[qq] += a_qq_qr_cc = ATR(ss,qq_qr,qq_cc,qr_cc) ;
00275 
00276       ss      = SS(rr_qr,rr_cc,qr_cc) ;     /* subtriangle rr,qr,cc */
00277       ww[rr] += a_rr_qr_cc = ATR(ss,rr_qr,rr_cc,qr_cc) ;
00278 
00279       ss      = SS(rr_pr,rr_cc,pr_cc) ;     /* subtriangle rr,pr,cc */
00280       ww[rr] += a_rr_pr_cc = ATR(ss,rr_pr,rr_cc,pr_cc) ;
00281 
00282 
00283 #if 0          /* debugging printouts */
00284 # undef  DDD
00285 # define DDD(x,y,z,a,b,c) sqrt((x-a)*(x-a)+(y-b)*(y-b)+(z-c)*(z-c))
00286 
00287       cp = DDD(xpp,ypp,zpp,xqq,yqq,zqq) ;
00288       sp = DDD(xpp,ypp,zpp,xrr,yrr,zrr) ;
00289       ca = DDD(xqq,yqq,zqq,xrr,yrr,zrr) ;
00290       sa = (cp+sp+ca)/2 ;
00291       ss = sqrt(sa*(sa-cp)*(sa-sp)*(sa-ca)) ;
00292 
00293       fprintf(stderr,"triangle %d: pp=%d qq=%d rr=%d  AREA=%6.3f PLANAR=%6.3f\n"
00294                      "  xpp=%6.3f ypp=%6.3f zpp=%6.3f\n"
00295                      "  xqq=%6.3f yqq=%6.3f zqq=%6.3f\n"
00296                      "  xrr=%6.3f yrr=%6.3f zrr=%6.3f\n"
00297                      "  xpq=%6.3f ypq=%6.3f zpq=%6.3f\n"
00298                      "  xqr=%6.3f yqr=%6.3f zqr=%6.3f\n"
00299                      "  xpr=%6.3f ypr=%6.3f zpr=%6.3f\n"
00300                      "  xcc=%6.3f ycc=%6.3f zcc=%6.3f\n"
00301 #ifdef USE_NORMAL
00302                      "  xnn=%6.3f ynn=%6.3f znn=%6.3f\n"
00303 #endif
00304                      "  pp_pq=%6.3f pp_pr=%6.3f pp_cc=%6.3f\n"
00305                      "  qq_pq=%6.3f qq_qr=%6.3f qq_cc=%6.3f\n"
00306                      "  rr_qr=%6.3f rr_cc=%6.3f rr_pr=%6.3f\n"
00307                      "  pq_cc=%6.3f qr_cc=%6.3f pr_cc=%6.3f\n"
00308                      "  a_pp_pq_cc=%6.3f a_pp_pr_cc=%6.3f\n"
00309                      "  a_qq_pq_cc=%6.3f a_qq_qr_cc=%6.3f\n"
00310                      "  a_rr_qr_cc=%6.3f a_rr_pr_cc=%6.3f\n" ,
00311                jj,pp,qq,rr ,
00312                a_pp_pq_cc+a_pp_pr_cc+a_qq_pq_cc+a_qq_qr_cc+a_rr_qr_cc+a_rr_pr_cc , ss ,
00313                xpp, ypp, zpp,
00314                xqq, yqq, zqq,
00315                xrr, yrr, zrr,
00316                xpq, ypq, zpq,
00317                xqr, yqr, zqr,
00318                xpr, ypr, zpr,
00319                xcc, ycc, zcc,
00320 #ifdef USE_NORMAL
00321                xnn, ynn, znn,
00322 #endif
00323                pp_pq, pp_pr, pp_cc,
00324                qq_pq, qq_qr, qq_cc,
00325                rr_qr, rr_cc, rr_pr,
00326                pq_cc, qr_cc, pr_cc,
00327                a_pp_pq_cc, a_pp_pr_cc,
00328                a_qq_pq_cc, a_qq_qr_cc,
00329                a_rr_qr_cc, a_rr_pr_cc  ) ;
00330 #endif
00331 
00332    } /* end of loop over triangles */
00333 
00334    /* exit stage left */
00335 
00336    *wt = ww ; return 1 ;
00337 }