00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
00061
00062
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
00069 for(i=n; i < argc; i++)
00070 sl += strlen(argv[i]);
00071
00072
00073 if(sl) {
00074 newstr = new char[sl + 8 + argc - n];
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
00085 return newstr;
00086 }
00087
00088
00089
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
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
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
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
00156
00157
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
00173 if((namestrt = strrchr(full,'/')) != NULL && strlen(namestrt) > 0) {
00174 namestrt++;
00175 } else {
00176 namestrt = full;
00177 }
00178
00179
00180 *name = stringdup(namestrt);
00181
00182
00183 pathlen = strlen(full) - strlen(*name);
00184 *path = new char[pathlen + 1];
00185 strncpy(*path,full,pathlen);
00186 (*path)[pathlen] = '\0';
00187 }
00188
00189
00190 char *str_tokenize(const char *newcmd, int *argc, char *argv[]) {
00191 char *cmd;
00192 const char *cmdstart;
00193 cmdstart = newcmd;
00194
00195
00196
00197
00198 while (cmdstart != NULL &&
00199 (*cmdstart == ' ' ||
00200 *cmdstart == ',' ||
00201 *cmdstart == ';' ||
00202 *cmdstart == '\t' ||
00203 *cmdstart == '\n')) {
00204 cmdstart++;
00205 }
00206
00207 cmd = stringdup(cmdstart);
00208 *argc = 0;
00209
00210
00211 argv[*argc] = strtok(cmd, " ,;\t\n");
00212
00213
00214 while(argv[*argc] != NULL) {
00215
00216 if(argv[*argc][0] == '#') {
00217 break;
00218 } else {
00219 (*argc)++;
00220 }
00221
00222
00223 argv[*argc] = strtok(NULL," ,;\t\n");
00224 }
00225
00226 return (*argc > 0 ? argv[0] : (char *) NULL);
00227 }
00228
00229
00230
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) {
00273 if (errno != EINTR)
00274 printf("select() error while attempting to read text input.\n");
00275 return FALSE;
00276 } else if (ret == 0) {
00277 return FALSE;
00278 }
00279 return TRUE;
00280 #endif
00281 }
00282
00283
00284
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
00315
00316
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
00326
00327 float * vec_normalize(float *vect) {
00328 float len2 = vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2];
00329
00330
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
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
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
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
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
00380
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
00387 if ((s1_length >= 1.0f ) || (s1_length <= -1.0f))
00388 return 1;
00389 else
00390 return 0;
00391 }
00392
00393
00394
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
00405
00406
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
00421
00422 return 57.2958f * (float) atan2(psin, pcos);
00423 }
00424
00425
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);
00451
00452 char *tmpfilebuf = new char[1024];
00453
00454
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
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;
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; ;
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
00588 long pgsz = sysconf(_SC_PAGESIZE);
00589 long physpgs = sysconf(_SC_PHYS_PAGES);
00590 return ((pgsz / 1024) * physpgs) / 1024;
00591 #elif defined(__APPLE__)
00592
00593
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;
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;
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; ;
00632 val += strtol(pos, (char **)NULL, 10);
00633 }
00634 pos=strstr(meminfobuf,"Buffers:");
00635 if (pos != NULL) {
00636 pos += 8; ;
00637 val += strtol(pos, (char **)NULL, 10);
00638 }
00639 pos=strstr(meminfobuf,"Cached:");
00640 if (pos != NULL) {
00641 pos += 8; ;
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
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
00663
00664
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;
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;
00689 }
00690
00691
00692
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;
00701 }
00702
00703
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
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
00734
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;
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
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
00783
00784
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;
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++;
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;
00827 } else {
00828 activelist[testpt]=0;
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;
00859 attempts++;
00860 } else if (popul==(N+1)) {
00861 converged=1;
00862 break;
00863 }
00864 if (attempts>10) {
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