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 }

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 }

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 }