Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

utilities.C

Go to the documentation of this file.
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 

Generated on Tue Dec 3 02:46:06 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002