utilities.C

00001 /***************************************************************************
00002 *cr 
00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 
00004 *cr University of Illinois 
00005 *cr All Rights Reserved 
00006 *cr 
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: utilities.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.177 $ $Date: 2021年09月23日 15:20:20 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 *
00019 * General utility routines and definitions.
00020 *
00021 ***************************************************************************/
00022 
00023 #include <string.h>
00024 #include <ctype.h>
00025 #include <math.h>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 
00029 #if defined(_MSC_VER)
00030 #include <windows.h>
00031 #include <conio.h>
00032 #else
00033 #include <unistd.h>
00034 #include <sys/time.h>
00035 #include <errno.h>
00036 
00037 #if defined(ARCH_AIX4)
00038 #include <strings.h>
00039 #endif
00040 
00041 #if defined(__irix)
00042 #include <bstring.h>
00043 #endif
00044 
00045 #if defined(__hpux)
00046 #include <time.h>
00047 #endif // HPUX
00048 #endif // _MSC_VER
00049 
00050 #if defined(AIXUSEPERFSTAT)
00051 #include <libperfstat.h>
00052 #endif
00053 
00054 #if defined(__APPLE__)
00055 #include <sys/sysctl.h>
00056 #endif
00057 
00058 #include "utilities.h"
00059 
00060 // given an argc, argv pair, take all the arguments from the Nth one on
00061 // and combine them into a single string with spaces separating words. This
00062 // allocates space for the string, which must be freed by the user.
00063 char *combine_arguments(int argc, const char **argv, int n) {
00064 char *newstr = NULL;
00065 
00066 if(argc > 0 && n < argc && n >= 0) {
00067 int i, sl = 0;
00068 // find out the length of the words we must combine
00069 for(i=n; i < argc; i++)
00070 sl += strlen(argv[i]);
00071 
00072 // combine the words together
00073 if(sl) {
00074 newstr = new char[sl + 8 + argc - n]; // extra buffer added
00075 *newstr = '0円';
00076 for(i=n; i < argc; i++) {
00077 if(i != n)
00078 strcat(newstr," ");
00079 strcat(newstr, argv[i]);
00080 }
00081 }
00082 }
00083 
00084 // return the string, or NULL if a problem occurred
00085 return newstr;
00086 }
00087 
00088 
00089 // duplicate a string using c++ new call
00090 char *stringdup(const char *s) {
00091 char *rs;
00092 
00093 if(!s)
00094 return NULL;
00095 
00096 rs = new char[strlen(s) + 1];
00097 strcpy(rs,s);
00098 
00099 return rs;
00100 }
00101 
00102 
00103 // convert a string to upper case
00104 char *stringtoupper(char *s) {
00105 if (s != NULL) {
00106 int i;
00107 int sz = strlen(s);
00108 for(i=0; i<sz; i++)
00109 s[i] = toupper(s[i]);
00110 }
00111 
00112 return s;
00113 }
00114 
00115 void stripslashes(char *str) {
00116 while (strlen(str) > 0 && str[strlen(str) - 1] == '/') {
00117 str[strlen(str) - 1] = '0円';
00118 }
00119 }
00120 
00121 // do upper-case comparison
00122 int strupcmp(const char *a, const char *b) {
00123 char *ua, *ub;
00124 int retval;
00125 
00126 ua = stringtoupper(stringdup(a));
00127 ub = stringtoupper(stringdup(b));
00128 
00129 retval = strcmp(ua,ub);
00130 
00131 delete [] ub;
00132 delete [] ua;
00133 
00134 return retval;
00135 }
00136 
00137 
00138 // do upper-case comparison, up to n characters
00139 int strupncmp(const char *a, const char *b, int n) {
00140 #if defined(ARCH_AIX3) || defined(ARCH_AIX4) || defined(_MSC_VER)
00141 while (n-- > 0) {
00142 if (toupper(*a) != toupper(*b)) {
00143 return toupper(*b) - toupper(*a);
00144 }
00145 if (*a == 0) return 0;
00146 a++; b++;
00147 }
00148 return 0;
00149 #else
00150 return strncasecmp(a, b, n);
00151 #endif
00152 }
00153 
00154 
00155 // break a file name up into path + name, returning both in the specified
00156 // character pointers. This creates storage for the new strings
00157 // by allocating space for them.
00158 void breakup_filename(const char *full, char **path, char **name) {
00159 const char *namestrt;
00160 int pathlen;
00161 
00162 if(full == NULL) {
00163 *path = *name = NULL;
00164 return;
00165 } else if (strlen(full) == 0) {
00166 *path = new char[1];
00167 *name = new char[1];
00168 (*path)[0] = (*name)[0] = '0円';
00169 return;
00170 }
00171 
00172 // find start of final file name
00173 if((namestrt = strrchr(full,'/')) != NULL && strlen(namestrt) > 0) {
00174 namestrt++;
00175 } else {
00176 namestrt = full;
00177 }
00178 
00179 // make a copy of the name
00180 *name = stringdup(namestrt);
00181 
00182 // make a copy of the path
00183 pathlen = strlen(full) - strlen(*name);
00184 *path = new char[pathlen + 1];
00185 strncpy(*path,full,pathlen);
00186 (*path)[pathlen] = '0円';
00187 } 
00188 
00189 // break a configuration line up into tokens.
00190 char *str_tokenize(const char *newcmd, int *argc, char *argv[]) {
00191 char *cmd; 
00192 const char *cmdstart;
00193 cmdstart = newcmd;
00194 
00195 // guarantee that the command string we return begins on the first
00196 // character returned by strtok(), otherwise the subsequent delete[]
00197 // calls will reference invalid memory blocks
00198 while (cmdstart != NULL &&
00199 (*cmdstart == ' ' ||
00200 *cmdstart == ',' ||
00201 *cmdstart == ';' ||
00202 *cmdstart == '\t' ||
00203 *cmdstart == '\n')) {
00204 cmdstart++; // advance pointer to first command character
00205 } 
00206 
00207 cmd = stringdup(cmdstart);
00208 *argc = 0;
00209 
00210 // initialize tokenizing calls
00211 argv[*argc] = strtok(cmd, " ,;\t\n");
00212 
00213 // loop through words until end-of-string, or comment character, found
00214 while(argv[*argc] != NULL) {
00215 // see if the token starts with '#'
00216 if(argv[*argc][0] == '#') {
00217 break; // don't process any further tokens
00218 } else {
00219 (*argc)++; // another token in list
00220 }
00221 
00222 // scan for next token
00223 argv[*argc] = strtok(NULL," ,;\t\n");
00224 }
00225 
00226 return (*argc > 0 ? argv[0] : (char *) NULL);
00227 }
00228 
00229 
00230 // get the time of day from the system clock, and store it (in seconds)
00231 double time_of_day(void) {
00232 #if defined(_MSC_VER)
00233 double t;
00234 
00235 t = GetTickCount(); 
00236 t = t / 1000.0;
00237 
00238 return t;
00239 #else
00240 struct timeval tm;
00241 struct timezone tz;
00242 
00243 gettimeofday(&tm, &tz);
00244 return((double)(tm.tv_sec) + (double)(tm.tv_usec)/1000000.0);
00245 #endif
00246 }
00247 
00248 
00249 int vmd_check_stdin(void) {
00250 #if defined(_MSC_VER)
00251 if (_kbhit() != 0)
00252 return TRUE;
00253 else
00254 return FALSE;
00255 #else
00256 fd_set readvec;
00257 struct timeval timeout;
00258 int ret, stdin_fd;
00259 
00260 timeout.tv_sec = 0;
00261 timeout.tv_usec = 0;
00262 stdin_fd = 0;
00263 FD_ZERO(&readvec);
00264 FD_SET(stdin_fd, &readvec);
00265 
00266 #if !defined(ARCH_AIX3)
00267 ret = select(16, &readvec, NULL, NULL, &timeout);
00268 #else
00269 ret = select(16, (int *)(&readvec), NULL, NULL, &timeout);
00270 #endif
00271 
00272 if (ret == -1) { // got an error
00273 if (errno != EINTR) // XXX: this is probably too lowlevel to be converted to Inform.h
00274 printf("select() error while attempting to read text input.\n");
00275 return FALSE;
00276 } else if (ret == 0) {
00277 return FALSE; // select timed out
00278 }
00279 return TRUE;
00280 #endif
00281 }
00282 
00283 
00284 // return the username of the currently logged-on user
00285 char *vmd_username(void) {
00286 #if defined(_MSC_VER)
00287 char username[1024];
00288 unsigned long size = 1023;
00289 
00290 if (GetUserName((char *) &username, &size)) {
00291 return stringdup(username);
00292 }
00293 else { 
00294 return stringdup("Windows User");
00295 }
00296 #else
00297 #if defined(ARCH_FREEBSD) || defined(ARCH_FREEBSDAMD64) || defined(__APPLE__) || defined(__linux)
00298 return stringdup(getlogin());
00299 #else
00300 return stringdup(cuserid(NULL));
00301 #endif 
00302 #endif
00303 }
00304 
00305 int vmd_getuid(void) {
00306 #if defined(_MSC_VER)
00307 return 0;
00308 #else
00309 return getuid(); 
00310 #endif
00311 }
00312 
00313 
00314 // take three 3-vectors and compute x2 cross x3; with the results
00315 // in x1. x1 must point to different memory than x2 or x3
00316 // This returns a pointer to x1
00317 float * cross_prod(float *x1, const float *x2, const float *x3)
00318 {
00319 x1[0] = x2[1]*x3[2] - x3[1]*x2[2];
00320 x1[1] = -x2[0]*x3[2] + x3[0]*x2[2];
00321 x1[2] = x2[0]*x3[1] - x3[0]*x2[1];
00322 return x1;
00323 }
00324 
00325 // normalize a vector, and return a pointer to it
00326 // Warning: it changes the value of the vector!!
00327 float * vec_normalize(float *vect) {
00328 float len2 = vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2];
00329 
00330 // prevent division by zero
00331 if (len2 > 0) {
00332 float rescale = 1.0f / sqrtf(len2);
00333 vect[0] *= rescale;
00334 vect[1] *= rescale;
00335 vect[2] *= rescale;
00336 }
00337 
00338 return vect;
00339 }
00340 
00341 
00342 // find and return the norm of a 3-vector
00343 float norm(const float *vect) {
00344 return sqrtf(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2]);
00345 }
00346 
00347 
00348 // determine if a triangle is degenerate or not
00349 int tri_degenerate(const float * v0, const float * v1, const float * v2) {
00350 float s1[3], s2[3], s1_length, s2_length;
00351 
00352 /*
00353 various rendering packages have amusingly different ideas about what
00354 constitutes a degenerate triangle. -1 and 1 work well. numbers
00355 below 0.999 and -0.999 show up in OpenGL
00356 numbers as low as 0.98 have worked in POVRay with certain models while
00357 numbers as high as 0.999999 have produced massive holes in other
00358 models
00359 -matt 11/13/96
00360 */
00361 
00362 /**************************************************************/
00363 /* turn the triangle into 2 normalized vectors. */
00364 /* If the dot product is 1 or -1 then */
00365 /* the triangle is degenerate */
00366 /**************************************************************/
00367 s1[0] = v0[0] - v1[0];
00368 s1[1] = v0[1] - v1[1];
00369 s1[2] = v0[2] - v1[2];
00370 
00371 s2[0] = v0[0] - v2[0];
00372 s2[1] = v0[1] - v2[1];
00373 s2[2] = v0[2] - v2[2];
00374 
00375 s1_length = sqrtf(s1[0]*s1[0] + s1[1]*s1[1] + s1[2]*s1[2]);
00376 s2_length = sqrtf(s2[0]*s2[0] + s2[1]*s2[1] + s2[2]*s2[2]);
00377 
00378 /**************************************************************/
00379 /* invert to avoid divides: */
00380 /* 1.0/v1_length * 1.0/v2_length */
00381 /**************************************************************/
00382 
00383 s2_length = 1.0f / (s1_length*s2_length);
00384 s1_length = s2_length * (s1[0]*s2[0] + s1[1]*s2[1] + s1[2]*s2[2]);
00385 
00386 // and add it to the list if it's not degenerate
00387 if ((s1_length >= 1.0f ) || (s1_length <= -1.0f)) 
00388 return 1;
00389 else
00390 return 0;
00391 }
00392 
00393 
00394 // compute the angle (in degrees 0 to 180 ) between two vectors a & b
00395 float angle(const float *a, const float *b) {
00396 float ab[3];
00397 cross_prod(ab, a, b);
00398 float psin = sqrtf(dot_prod(ab, ab));
00399 float pcos = dot_prod(a, b);
00400 return 57.2958f * (float) atan2(psin, pcos);
00401 }
00402 
00403 
00404 // Compute the dihedral angle for the given atoms, returning a value between
00405 // -180 and 180.
00406 // faster, cleaner implementation based on atan2
00407 float dihedral(const float *a1,const float *a2,const float *a3,const float *a4)
00408 {
00409 float r1[3], r2[3], r3[3], n1[3], n2[3];
00410 vec_sub(r1, a2, a1);
00411 vec_sub(r2, a3, a2);
00412 vec_sub(r3, a4, a3);
00413 
00414 cross_prod(n1, r1, r2);
00415 cross_prod(n2, r2, r3);
00416 
00417 float psin = dot_prod(n1, r3) * sqrtf(dot_prod(r2, r2));
00418 float pcos = dot_prod(n1, n2);
00419 
00420 // atan2f would be faster, but we'll have to workaround the lack
00421 // of existence on some platforms.
00422 return 57.2958f * (float) atan2(psin, pcos);
00423 }
00424 
00425 // compute the distance between points a & b
00426 float distance(const float *a, const float *b) {
00427 return sqrtf(distance2(a,b));
00428 }
00429 
00430 char *vmd_tempfile(const char *s) {
00431 char *envtxt, *TempDir;
00432 
00433 if((envtxt = getenv("VMDTMPDIR")) != NULL) {
00434 TempDir = stringdup(envtxt);
00435 } else {
00436 #if defined(_MSC_VER)
00437 if ((envtxt = getenv("TMP")) != NULL) {
00438 TempDir = stringdup(envtxt);
00439 }
00440 else if ((envtxt = getenv("TEMP")) != NULL) {
00441 TempDir = stringdup(envtxt);
00442 }
00443 else {
00444 TempDir = stringdup("c:\\\\");
00445 }
00446 #else
00447 TempDir = stringdup("/tmp");
00448 #endif
00449 }
00450 stripslashes(TempDir); // strip out ending '/' chars.
00451 
00452 char *tmpfilebuf = new char[1024];
00453 
00454 // copy in temp string
00455 strcpy(tmpfilebuf, TempDir);
00456 
00457 #if defined(_MSC_VER)
00458 strcat(tmpfilebuf, "\\");
00459 strncat(tmpfilebuf, s, 1022 - strlen(TempDir));
00460 #else
00461 strcat(tmpfilebuf, "/");
00462 strncat(tmpfilebuf, s, 1022 - strlen(TempDir));
00463 #endif
00464 
00465 tmpfilebuf[1023] = '0円';
00466 
00467 delete [] TempDir;
00468 
00469 // return converted string
00470 return tmpfilebuf;
00471 }
00472 
00473 
00474 int vmd_delete_file(const char * path) {
00475 #if defined(_MSC_VER)
00476 if (DeleteFile(path) == 0) 
00477 return -1;
00478 else 
00479 return 0; 
00480 #else
00481 return unlink(path);
00482 #endif
00483 }
00484 
00485 void vmd_sleep(int secs) {
00486 #if defined(_MSC_VER)
00487 Sleep(secs * 1000);
00488 #else 
00489 sleep(secs);
00490 #endif
00491 }
00492 
00493 void vmd_msleep(int msecs) {
00494 #if defined(_MSC_VER)
00495 Sleep(msecs);
00496 #else 
00497 struct timeval timeout;
00498 timeout.tv_sec = 0;
00499 timeout.tv_usec = 1000 * msecs;
00500 select(0, NULL, NULL, NULL, &timeout);
00501 #endif // _MSC_VER
00502 }
00503 
00504 int vmd_system(const char* cmd) {
00505 return system(cmd);
00506 }
00507 
00508 
00512 long vmd_random(void) {
00513 #ifdef _MSC_VER
00514 return rand();
00515 #else
00516 return random();
00517 #endif
00518 }
00519 
00520 void vmd_srandom(unsigned int seed) {
00521 #ifdef _MSC_VER
00522 srand(seed);
00523 #else
00524 srandom(seed);
00525 #endif
00526 }
00527 
00530 float vmd_random_gaussian() {
00531 static bool cache = false;
00532 static float cached_value;
00533 const float RAND_FACTOR = 2.f/float(VMD_RAND_MAX);
00534 float r, s, w;
00535 
00536 if (cache) {
00537 cache = false;
00538 return cached_value;
00539 }
00540 do {
00541 r = RAND_FACTOR*vmd_random()-1.f; 
00542 s = RAND_FACTOR*vmd_random()-1.f;
00543 w = r*r+s*s;
00544 } while (w >= 1.f);
00545 w = sqrtf(-2.f*logf(w)/w);
00546 cached_value = s * w;
00547 cache = true;
00548 return (r*w);
00549 }
00550 
00551 
00554 long vmd_get_total_physmem_mb(void) {
00555 #if defined(_MSC_VER)
00556 MEMORYSTATUS memstat;
00557 GlobalMemoryStatus(&memstat);
00558 if (memstat.dwLength != sizeof(memstat))
00559 return -1; /* memstat result is wrong size! */
00560 return memstat.dwTotalPhys/(1024 * 1024);
00561 #elif defined(__linux)
00562 FILE *fp;
00563 char meminfobuf[1024], *pos;
00564 size_t len;
00565 
00566 fp = fopen("/proc/meminfo", "r");
00567 if (fp != NULL) {
00568 len = fread(meminfobuf,1,1024, fp);
00569 meminfobuf[1023] = 0;
00570 fclose(fp);
00571 if (len > 0) {
00572 pos=strstr(meminfobuf,"MemTotal:");
00573 if (pos == NULL) 
00574 return -1;
00575 pos += 9; /* skip tag */;
00576 return strtol(pos, (char **)NULL, 10)/1024L;
00577 }
00578 } 
00579 return -1;
00580 #elif defined(AIXUSEPERFSTAT) && defined(_AIX)
00581 perfstat_memory_total_t minfo;
00582 perfstat_memory_total(NULL, &minfo, sizeof(perfstat_memory_total_t), 1);
00583 return minfo.real_total*(4096/1024)/1024;
00584 #elif defined(_AIX)
00585 return (sysconf(_SC_AIX_REALMEM) / 1024);
00586 #elif defined(_SC_PAGESIZE) && defined(_SC_PHYS_PAGES)
00587 /* SysV Unix */
00588 long pgsz = sysconf(_SC_PAGESIZE);
00589 long physpgs = sysconf(_SC_PHYS_PAGES);
00590 return ((pgsz / 1024) * physpgs) / 1024;
00591 #elif defined(__APPLE__)
00592 /* MacOS X uses BSD sysctl */
00593 /* use hw.memsize, as it's a 64-bit value */
00594 int rc;
00595 uint64_t membytes;
00596 size_t len = sizeof(membytes);
00597 if (sysctlbyname("hw.memsize", &membytes, &len, NULL, 0)) 
00598 return -1;
00599 return (membytes / (1024*1024));
00600 #else
00601 return -1; /* unrecognized system, no method to get this info */
00602 #endif
00603 }
00604 
00605 
00606 
00609 long vmd_get_avail_physmem_mb(void) {
00610 #if defined(_MSC_VER)
00611 MEMORYSTATUS memstat;
00612 GlobalMemoryStatus(&memstat);
00613 if (memstat.dwLength != sizeof(memstat))
00614 return -1; /* memstat result is wrong size! */ 
00615 return memstat.dwAvailPhys / (1024 * 1024);
00616 #elif defined(__linux)
00617 FILE *fp;
00618 char meminfobuf[1024], *pos;
00619 size_t len;
00620 long val;
00621 
00622 fp = fopen("/proc/meminfo", "r");
00623 if (fp != NULL) {
00624 len = fread(meminfobuf,1,1024, fp);
00625 meminfobuf[1023] = 0;
00626 fclose(fp);
00627 if (len > 0) {
00628 val = 0L;
00629 pos=strstr(meminfobuf,"MemFree:");
00630 if (pos != NULL) {
00631 pos += 8; /* skip tag */;
00632 val += strtol(pos, (char **)NULL, 10);
00633 }
00634 pos=strstr(meminfobuf,"Buffers:");
00635 if (pos != NULL) {
00636 pos += 8; /* skip tag */;
00637 val += strtol(pos, (char **)NULL, 10);
00638 }
00639 pos=strstr(meminfobuf,"Cached:");
00640 if (pos != NULL) {
00641 pos += 8; /* skip tag */;
00642 val += strtol(pos, (char **)NULL, 10);
00643 }
00644 return val/1024L;
00645 } else {
00646 return -1;
00647 }
00648 } else {
00649 return -1;
00650 }
00651 #elif defined(AIXUSEPERFSTAT) && defined(_AIX)
00652 perfstat_memory_total_t minfo;
00653 perfstat_memory_total(NULL, &minfo, sizeof(perfstat_memory_total_t), 1);
00654 return minfo.real_free*(4096/1024)/1024;
00655 #elif defined(_SC_PAGESIZE) && defined(_SC_AVPHYS_PAGES)
00656 /* SysV Unix */
00657 long pgsz = sysconf(_SC_PAGESIZE);
00658 long avphyspgs = sysconf(_SC_AVPHYS_PAGES);
00659 return ((pgsz / 1024) * avphyspgs) / 1024;
00660 #elif defined(__APPLE__)
00661 #if 0
00662 /* BSD sysctl */
00663 /* hw.usermem isn't really the amount of free memory, it's */
00664 /* really more a measure of the non-kernel memory */
00665 int rc;
00666 int membytes;
00667 size_t len = sizeof(membytes);
00668 if (sysctlbyname("hw.usermem", &membytes, &len, NULL, 0)) 
00669 return -1;
00670 return (membytes / (1024*1024));
00671 #else
00672 return -1;
00673 #endif
00674 #else
00675 return -1; /* unrecognized system, no method to get this info */
00676 #endif
00677 }
00678 
00679 
00681 long vmd_get_avail_physmem_percent(void) {
00682 double total, avail;
00683 total = (double) vmd_get_total_physmem_mb();
00684 avail = (double) vmd_get_avail_physmem_mb();
00685 if (total > 0.0 && avail >= 0.0)
00686 return (long) (avail / (total / 100.0));
00687 
00688 return -1; /* return an error */
00689 }
00690 
00691 
00692 // returns minimum distance for Poisson disk sampler
00693 float correction(int nrays) {
00694 float N=(float)nrays;
00695 float eightPi=VMD_PIF*8.0f;
00696 float denom=(float)( N*(3.0f*sqrtf(3)) );
00697 float ans=sqrtf(eightPi/denom);
00698 float minD=sqrtf((ans*ans) + powf(((VMD_PIF/6.0f)*ans),2) );
00699 
00700 return 1.1275f*minD; //with padding
00701 }
00702 
00703 // print poisson points to .xyz file
00704 void print_xyz(float* population, int nrays) {
00705 float x,y,z;
00706 FILE* fp=fopen("validate.xyz","w");
00707 fprintf(fp,"%d\n",nrays);
00708 for (int i=0; i<=nrays; i++) {
00709 x=cosf(population[i*2+0])*sinf(population[i*2+1]);
00710 y=sinf(population[i*2+0])*sinf(population[i*2+1]);
00711 z=cosf(population[i*2+1]);
00712 fprintf(fp,"C %2.6f %2.6f %2.6f\n",x,y,z);
00713 }
00714 fclose(fp);
00715 return;
00716 }
00717 
00718 // compute arc distance between two points on unit sphere
00719 float arcdistance(float lambda1, float lambda2, float phi1, float phi2) {
00720 float sl1,cl1,sl2,cl2;
00721 float sp1,cp1,sp2,cp2;
00722 
00723 sincosf(lambda1, &sl1, &cl1);
00724 sincosf(phi1, &sp1, &cp1);
00725 sincosf(lambda2, &sl2, &cl2);
00726 sincosf(phi2, &sp2, &cp2);
00727 
00728 float cos_Ang=(float)( ((cl1*sp1)*(cl2*sp2)) + ((sl1*sp1)*(sl2*sp2)) + (cp1*cp2) );
00729 
00730 return acosf(cos_Ang);
00731 }
00732 
00733 // generate K candidates, return the farthest from all previously
00734 // committed samples
00735 int k_candidates(int k, int nrays, int idx, int testpt, float minD, float* candidates, float* population) {
00736 static const float RAND_MAX_INV=1.0f/float(VMD_RAND_MAX);
00737 int bestIdx=-1;
00738 int count=0; //score for each test point
00739 float bestDist=0.0f;
00740 float lambda1=population[testpt*2+0];
00741 float phi1=population[testpt*2+1];
00742 float currDist;
00743 
00744 float minLambda=VMD_TWOPIF*minD;
00745 float minPhi=VMD_PIF*minD;
00746 
00747 float dp, dl;
00748 for (int i=0; i<k; i++) {
00749 dl=(float)(lambda1-(minLambda+((float)(RAND_MAX_INV*vmd_random()))*minLambda));
00750 dp=(float)(phi1-(minPhi+((float)(RAND_MAX_INV*vmd_random()))*minPhi));
00751 
00752 candidates[i*2+0]=(float)(lambda1+(dl));
00753 candidates[i*2+1]=(float)(phi1+(dp));
00754 }
00755 
00756 for (int j=0; j<k; j++) {
00757 currDist=0.0f;
00758 count=0;
00759 for(int jj=0; jj<idx; jj++) {
00760 float dist=arcdistance(population[jj*2+0],candidates[j*2+0],
00761 population[jj*2+1],candidates[j*2+1]);
00762 
00763 if ((dist-minD)>0.00000001f) {
00764 currDist+=dist;
00765 count++;
00766 }
00767 }
00768 
00769 //pick best candidate
00770 if (count==idx) {
00771 if (currDist>bestDist) {
00772 bestIdx=j;
00773 bestDist=currDist;
00774 count=0;
00775 }
00776 count=0;
00777 }
00778 }
00779 return bestIdx;
00780 }
00781 
00782 // Poisson disk sampler -- fills 2D array with N 
00783 // lambda phi pairs which embed in the unit sphere when 
00784 // converted to cartesian coordinates
00785 int poisson_sample_on_sphere(float *population, int N, int k, int verbose) {
00786 int converged=0; 
00787 int result=0; 
00788 int popul=0;
00789 int fc=0; 
00790 int testpt=0; 
00791 int attempts=0;
00792 int tr=0;
00793 
00794 float minD=correction(N);
00795 static const float RAND_MAX_INV=1.0f/float(VMD_RAND_MAX);
00796 
00797 float *candidates=NULL;
00798 int *activelist=NULL;
00799 
00800 candidates=new float[k*2];
00801 activelist=new int[N];
00802 for (int ii=0; ii<k; ii++) {
00803 candidates[ii*2+0]=0.0f;
00804 candidates[ii*2+1]=0.0f;
00805 }
00806 for (int jj=0; jj<N; jj++) {
00807 population[jj*2+0]=0.0f;
00808 population[jj*2+1]=0.0f;
00809 activelist[jj]=1; //set all to active
00810 }
00811 int numactive=N;
00812 
00813 vmd_srandom(512346);
00814 population[0]=(float)((RAND_MAX_INV*vmd_random())*VMD_TWOPI);
00815 population[1]=(float)((RAND_MAX_INV*vmd_random())*VMD_PI);
00816 popul++; //for consistency -- popul incremented after each addition
00817 
00818 while (converged==0) {
00819 while (activelist[testpt]==1) {
00820 result=k_candidates(k,N,popul,testpt,minD,candidates,population);
00821 if (result!=-1) {
00822 population[popul*2+0]=candidates[result*2+0];
00823 population[popul*2+1]=candidates[result*2+1];
00824 popul++;
00825 if (popul==(N+1)) { converged=1; break; }
00826 testpt=vmd_random()%popul; //pick new test pt
00827 } else {
00828 activelist[testpt]=0; //if a new pt can't be added around testpt, deactivate
00829 numactive--;
00830 break;
00831 }
00832 if (popul==(N+1)) {
00833 converged=1;
00834 break;
00835 }
00836 }
00837 
00838 if (popul==(N+1)) {
00839 converged=1;
00840 break;
00841 } else if (fc==N*N) {
00842 if (verbose) printf("Poisson sampler failed at population=%d\n",popul);
00843 break;
00844 } else if (numactive>=2) {
00845 while (activelist[testpt]!=1) {
00846 testpt=vmd_random()%popul;
00847 tr++;
00848 if (tr>N*1000) { break; }
00849 }
00850 } else {
00851 if (verbose) printf("Error in Poisson disk sampling procedure.\n");
00852 }
00853 
00854 if (numactive<=2 && popul!=(N+1)) {
00855 for (int f=0; f<N; f++)
00856 activelist[f]=1;
00857 
00858 popul=1; fc=0; numactive=N; testpt=0; //failed to converge, reset
00859 attempts++;
00860 } else if (popul==(N+1)) {
00861 converged=1;
00862 break;
00863 }
00864 if (attempts>10) { //hard exit after 10 attempts
00865 if (verbose) printf("Poisson sampler failed to converge (population=%d)\n",popul);
00866 break;
00867 }
00868 }
00869 
00870 delete [] candidates;
00871 delete [] activelist;
00872 
00873 if (popul==(N+1) && converged==1) {
00874 if (verbose) printf("Poisson sampler done. (%d/10 attempts)\n",attempts+1);
00875 return 1;
00876 }
00877 
00878 if (converged==0) {
00879 if (verbose) printf("Poisson sampler failed -- exiting.");
00880 return -1;
00881 }
00882 
00883 return 0;
00884 }
00885