00001 #include "cs.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
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
00091
00092
00093
00094
00095
00096
00097
00098
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
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
00129
00130
00131
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
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
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++ ){
00164
00165
00166
00167
00168
00169
00170
00171
00172
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
00182
00183
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
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
00237
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
00256
00257
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) ;
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) ;
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) ;
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) ;
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) ;
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) ;
00280 ww[rr] += a_rr_pr_cc = ATR(ss,rr_pr,rr_cc,pr_cc) ;
00281
00282
00283 #if 0
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 }
00333
00334
00335
00336 *wt = ww ; return 1 ;
00337 }