00001
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
00032
00033
00034
00035
00036 #include "byteorder.h"
00037 #include "emassert.h"
00038 #include "emdata.h"
00039 #include "util.h"
00040 #include "marchingcubes.h"
00041 #include "randnum.h"
00042
00043 #include <fcntl.h>
00044 #include <iomanip>
00045 #include <sstream>
00046
00047 #include <cstring>
00048
00049 #include <sys/types.h>
00050 #include <gsl/gsl_linalg.h>
00051 #include <algorithm>
00052
00053 #ifndef WIN32
00054 #include <unistd.h>
00055 #include <sys/param.h>
00056 #else
00057 #include <ctime>
00058 #include <io.h>
00059 #define access _access
00060 #define F_OK 00
00061 #endif //WIN32
00062
00063 using namespace std;
00064 using namespace boost;
00065 using namespace EMAN;
00066
00068 int Util::MUTEX_INIT(MUTEX *mutex)
00069 {
00070 #ifdef WIN32
00071 *mutex = CreateMutex(0, FALSE, 0);
00072 return (*mutex==0);
00073 #else
00074 return pthread_mutex_init(mutex, NULL);
00075 #endif
00076 return -1;
00077 }
00078
00079 int Util::MUTEX_LOCK(MUTEX *mutex)
00080 {
00081 #ifdef WIN32
00082 return (WaitForSingleObject(*mutex, INFINITE)==WAIT_FAILED?1:0);
00083 #else
00084 return pthread_mutex_lock(mutex);
00085 #endif
00086 return -1;
00087 }
00088
00089 int Util::MUTEX_UNLOCK(MUTEX *mutex)
00090 {
00091 #ifdef WIN32
00092 return (ReleaseMutex(*mutex)==0);
00093 #else
00094 return pthread_mutex_unlock(mutex);
00095 #endif
00096 return -1;
00097 }
00098
00099
00101 void Util::ap2ri(float *data, size_t n)
00102 {
00103 Assert(n > 0);
00104
00105 if (!data) {
00106 throw NullPointerException("pixel data array");
00107 }
00108
00109 for (size_t i = 0; i < n; i += 2) {
00110 float f = data[i] * sin(data[i + 1]);
00111 data[i] = data[i] * cos(data[i + 1]);
00112 data[i + 1] = f;
00113 }
00114 }
00115
00116 void Util::flip_complex_phase(float *data, size_t n)
00117 {
00118 Assert(n > 0);
00119
00120 if (!data) {
00121 throw NullPointerException("pixel data array");
00122 }
00123
00124 for (size_t i = 0; i < n; i += 2) {
00125 data[i + 1] *= -1;
00126 }
00127 }
00128
00129 void Util::rotate_phase_origin(float *data, size_t nx, size_t ny, size_t nz)
00130 {
00131 if(ny==1 && nz==1) {
00132 return;
00133 }
00134 else if(ny!=1 && nz==1) {
00135 size_t i, j, k, l;
00136 float re;
00137 l=ny/2*nx;
00138 for (i=0; i<ny/2; i++) {
00139 for (j=0; j<nx; j++) {
00140 k=j+i*nx;
00141 re=data[k];
00142 data[k]=data[k+l];
00143 data[k+l]=re;
00144 }
00145 }
00146 }
00147 else {
00148 size_t i, j, k, l, ii, jj;
00149 char * t=(char *)malloc(sizeof(float)*nx);
00150
00151 k=nx*ny*(nz+1)/2;
00152 l=nx*ny*(nz-1)/2;
00153 jj=nx*sizeof(float);
00154 for (j=ii=0; j<nz/2; ++j) {
00155 for (i=0; i<ny; ++i,ii+=nx) {
00156 memcpy(t,data+ii,jj);
00157 if (i<ny/2) {
00158 memcpy(data+ii,data+ii+k,jj);
00159 memcpy(data+ii+k,t,jj);
00160 }
00161 else {
00162 memcpy(data+ii,data+ii+l,jj);
00163 memcpy(data+ii+l,t,jj);
00164 }
00165 }
00166 }
00167 free(t);
00168 }
00169 }
00170
00171 int Util::file_lock_wait(FILE * file)
00172 {
00173 #ifdef WIN32
00174 return 1;
00175 #else
00176
00177 if (!file) {
00178 throw NullPointerException("Tried to lock NULL file");
00179 }
00180
00181 int fdes = fileno(file);
00182
00183 struct flock fl;
00184 fl.l_type = F_WRLCK;
00185 fl.l_whence = SEEK_SET;
00186 fl.l_start = 0;
00187 fl.l_len = 0;
00188 #ifdef WIN32
00189 fl.l_pid = _getpid();
00190 #else
00191 fl.l_pid = getpid();
00192 #endif
00193
00194 #if defined __sgi
00195 fl.l_sysid = getpid();
00196 #endif
00197
00198 int err = 0;
00199 if (fcntl(fdes, F_SETLKW, &fl) == -1) {
00200 LOGERR("file locking error! NFS problem?");
00201
00202 int i = 0;
00203 for (i = 0; i < 5; i++) {
00204 if (fcntl(fdes, F_SETLKW, &fl) != -1) {
00205 break;
00206 }
00207 else {
00208 #ifdef WIN32
00209 Sleep(1000);
00210 #else
00211 sleep(1);
00212 #endif
00213
00214 }
00215 }
00216 if (i == 5) {
00217 LOGERR("Fatal file locking error");
00218 err = 1;
00219 }
00220 }
00221
00222 return err;
00223 #endif
00224 }
00225
00226
00227
00228 bool Util::check_file_by_magic(const void *first_block, const char *magic)
00229 {
00230 if (!first_block || !magic) {
00231 throw NullPointerException("first_block/magic");
00232 }
00233
00234 const char *buf = static_cast < const char *>(first_block);
00235
00236 if (strncmp(buf, magic, strlen(magic)) == 0) {
00237 return true;
00238 }
00239 return false;
00240 }
00241
00242 bool Util::is_file_exist(const string & filename)
00243 {
00244 if (access(filename.c_str(), F_OK) == 0) {
00245 return true;
00246 }
00247 return false;
00248 }
00249
00250
00251 void Util::flip_image(float *data, size_t nx, size_t ny)
00252 {
00253 if (!data) {
00254 throw NullPointerException("image data array");
00255 }
00256 Assert(nx > 0);
00257 Assert(ny > 0);
00258
00259 float *buf = new float[nx];
00260 size_t row_size = nx * sizeof(float);
00261
00262 for (size_t i = 0; i < ny / 2; i++) {
00263 memcpy(buf, &data[i * nx], row_size);
00264 memcpy(&data[i * nx], &data[(ny - 1 - i) * nx], row_size);
00265 memcpy(&data[(ny - 1 - i) * nx], buf, row_size);
00266 }
00267
00268 if( buf )
00269 {
00270 delete[]buf;
00271 buf = 0;
00272 }
00273 }
00274
00275 string Util::str_to_lower(const string& s) {
00276 string ret(s);
00277 std::transform(s.begin(),s.end(),ret.begin(), (int (*)(int) ) std::tolower);
00278 return ret;
00279 }
00280
00281 bool Util::sstrncmp(const char *s1, const char *s2)
00282 {
00283 if (!s1 || !s2) {
00284 throw NullPointerException("Null string");
00285 }
00286
00287 if (strncmp(s1, s2, strlen(s2)) == 0) {
00288 return true;
00289 }
00290
00291 return false;
00292 }
00293
00294 string Util::int2str(int n)
00295 {
00296 char s[32] = {'\0'};
00297 sprintf(s, "%d", n);
00298 return string(s);
00299 }
00300
00301 string Util::get_line_from_string(char **slines)
00302 {
00303 if (!slines || !(*slines)) {
00304 throw NullPointerException("Null string");
00305 }
00306
00307 string result = "";
00308 char *str = *slines;
00309
00310 while (*str != '\n' && *str != '\0') {
00311 result.push_back(*str);
00312 str++;
00313 }
00314 if (*str != '\0') {
00315 str++;
00316 }
00317 *slines = str;
00318
00319 return result;
00320 }
00321
00322 vector<EMData *> Util::svdcmp(const vector<EMData *> &data,int nvec) {
00323 int nimg=data.size();
00324 if (nvec==0) nvec=nimg;
00325 vector<EMData *> ret(nvec);
00326 if (nimg==0) return ret;
00327 int pixels=data[0]->get_xsize()*data[0]->get_ysize()*data[0]->get_zsize();
00328
00329
00330 gsl_vector *work=gsl_vector_alloc(nimg);
00331 gsl_vector *S=gsl_vector_alloc(nimg);
00332 gsl_matrix *A=gsl_matrix_alloc(pixels,nimg);
00333 gsl_matrix *V=gsl_matrix_alloc(nimg,nimg);
00334 gsl_matrix *X=gsl_matrix_alloc(nimg,nimg);
00335
00336 int im,x,y,z,i;
00337 for (im=0; im<nimg; im++) {
00338 for (z=0,i=0; z<data[0]->get_zsize(); z++) {
00339 for (y=0; y<data[0]->get_ysize(); y++) {
00340 for (x=0; x<data[0]->get_xsize(); x++,i++) {
00341 gsl_matrix_set(A,i,im,data[im]->get_value_at(x,y,z));
00342 }
00343 }
00344 }
00345 }
00346
00347
00348 gsl_linalg_SV_decomp_mod (A,X, V, S, work);
00349
00350 for (im=0; im<nvec; im++) {
00351 EMData *a=data[0]->copy_head();
00352 ret[im]=a;
00353 for (z=0,i=0; z<data[0]->get_zsize(); z++) {
00354 for (y=0; y<data[0]->get_ysize(); y++) {
00355 for (x=0; x<data[0]->get_xsize(); x++,i++) {
00356 a->set_value_at(x,y,z,static_cast<float>(gsl_matrix_get(A,i,im)));
00357 }
00358 }
00359 }
00360 }
00361 return ret;
00362 }
00363
00364 bool Util::get_str_float(const char *s, const char *float_var, float *p_val)
00365 {
00366 if (!s || !float_var || !p_val) {
00367 throw NullPointerException("string float");
00368 }
00369 size_t n = strlen(float_var);
00370 if (strncmp(s, float_var, n) == 0) {
00371 *p_val = (float) atof(&s[n]);
00372 return true;
00373 }
00374
00375 return false;
00376 }
00377
00378 bool Util::get_str_float(const char *s, const char *float_var, float *p_v1, float *p_v2)
00379 {
00380 if (!s || !float_var || !p_v1 || !p_v2) {
00381 throw NullPointerException("string float");
00382 }
00383
00384 size_t n = strlen(float_var);
00385 if (strncmp(s, float_var, n) == 0) {
00386 sscanf(&s[n], "%f,%f", p_v1, p_v2);
00387 return true;
00388 }
00389
00390 return false;
00391 }
00392
00393 bool Util::get_str_float(const char *s, const char *float_var,
00394 int *p_v0, float *p_v1, float *p_v2)
00395 {
00396 if (!s || !float_var || !p_v0 || !p_v1 || !p_v2) {
00397 throw NullPointerException("string float");
00398 }
00399
00400 size_t n = strlen(float_var);
00401 *p_v0 = 0;
00402 if (strncmp(s, float_var, n) == 0) {
00403 if (s[n] == '=') {
00404 *p_v0 = 2;
00405 sscanf(&s[n + 1], "%f,%f", p_v1, p_v2);
00406 }
00407 else {
00408 *p_v0 = 1;
00409 }
00410 return true;
00411 }
00412 return false;
00413 }
00414
00415 bool Util::get_str_int(const char *s, const char *int_var, int *p_val)
00416 {
00417 if (!s || !int_var || !p_val) {
00418 throw NullPointerException("string int");
00419 }
00420
00421 size_t n = strlen(int_var);
00422 if (strncmp(s, int_var, n) == 0) {
00423 *p_val = atoi(&s[n]);
00424 return true;
00425 }
00426 return false;
00427 }
00428
00429 bool Util::get_str_int(const char *s, const char *int_var, int *p_v1, int *p_v2)
00430 {
00431 if (!s || !int_var || !p_v1 || !p_v2) {
00432 throw NullPointerException("string int");
00433 }
00434
00435 size_t n = strlen(int_var);
00436 if (strncmp(s, int_var, n) == 0) {
00437 sscanf(&s[n], "%d,%d", p_v1, p_v2);
00438 return true;
00439 }
00440 return false;
00441 }
00442
00443 bool Util::get_str_int(const char *s, const char *int_var, int *p_v0, int *p_v1, int *p_v2)
00444 {
00445 if (!s || !int_var || !p_v0 || !p_v1 || !p_v2) {
00446 throw NullPointerException("string int");
00447 }
00448
00449 size_t n = strlen(int_var);
00450 *p_v0 = 0;
00451 if (strncmp(s, int_var, n) == 0) {
00452 if (s[n] == '=') {
00453 *p_v0 = 2;
00454 sscanf(&s[n + 1], "%d,%d", p_v1, p_v2);
00455 }
00456 else {
00457 *p_v0 = 1;
00458 }
00459 return true;
00460 }
00461 return false;
00462 }
00463
00464 string Util::change_filename_ext(const string & old_filename,
00465 const string & ext)
00466 {
00467 Assert(old_filename != "");
00468 if (ext == "") {
00469 return old_filename;
00470 }
00471
00472 string filename = old_filename;
00473 size_t dot_pos = filename.rfind(".");
00474 if (dot_pos != string::npos) {
00475 filename = filename.substr(0, dot_pos+1);
00476 }
00477 else {
00478 filename = filename + ".";
00479 }
00480 filename = filename + ext;
00481 return filename;
00482 }
00483
00484 string Util::remove_filename_ext(const string& filename)
00485 {
00486 if (filename == "") {
00487 return "";
00488 }
00489
00490 char *buf = new char[filename.size()+1];
00491 strcpy(buf, filename.c_str());
00492 char *old_ext = strrchr(buf, '.');
00493 if (old_ext) {
00494 buf[strlen(buf) - strlen(old_ext)] = '\0';
00495 }
00496 string result = string(buf);
00497 if( buf )
00498 {
00499 delete [] buf;
00500 buf = 0;
00501 }
00502 return result;
00503 }
00504
00505 string Util::sbasename(const string & filename)
00506 {
00507 if (filename == "") {
00508 return "";
00509 }
00510
00511 char s = '/';
00512 #ifdef _WIN32
00513 s = '\\';
00514 #endif
00515 const char * c = strrchr(filename.c_str(), s);
00516 if (!c) {
00517 return filename;
00518 }
00519 else {
00520 c++;
00521 }
00522 return string(c);
00523 }
00524
00525
00526 string Util::get_filename_ext(const string& filename)
00527 {
00528 if (filename == "") {
00529 return "";
00530 }
00531
00532 string result = "";
00533 const char *ext = strrchr(filename.c_str(), '.');
00534 if (ext) {
00535 ext++;
00536 result = string(ext);
00537 }
00538 return result;
00539 }
00540
00541
00542
00543 void Util::calc_least_square_fit(size_t nitems, const float *data_x, const float *data_y,
00544 float *slope, float *intercept, bool ignore_zero,float absmax)
00545 {
00546 Assert(nitems > 0);
00547
00548 if (!data_x || !data_y || !slope || !intercept) {
00549 throw NullPointerException("null float pointer");
00550 }
00551 double sum = 0;
00552 double sum_x = 0;
00553 double sum_y = 0;
00554 double sum_xx = 0;
00555 double sum_xy = 0;
00556
00557 for (size_t i = 0; i < nitems; i++) {
00558 if ((!ignore_zero || (data_x[i] != 0 && data_y[i] != 0))&&(!absmax ||(data_y[i]<absmax && data_y[i]>-absmax))) {
00559 double y = data_y[i];
00560 double x = i;
00561 if (data_x) {
00562 x = data_x[i];
00563 }
00564
00565 sum_x += x;
00566 sum_y += y;
00567 sum_xx += x * x;
00568 sum_xy += x * y;
00569 sum++;
00570 }
00571 }
00572
00573 double div = sum * sum_xx - sum_x * sum_x;
00574 if (div == 0) {
00575 div = 0.0000001f;
00576 }
00577
00578 *intercept = (float) ((sum_xx * sum_y - sum_x * sum_xy) / div);
00579 *slope = (float) ((sum * sum_xy - sum_x * sum_y) / div);
00580 }
00581
00582 Vec3f Util::calc_bilinear_least_square(const vector<float> &p) {
00583 unsigned int i;
00584
00585
00586 double Sx=0,Sy=0,Sxy=0,Sxx=0,Syy=0,Sz=0,Sxz=0,Syz=0,S=0;
00587 for (i=0; i<p.size(); i+=3) {
00588 Sx+=p[i];
00589 Sy+=p[i+1];
00590 Sz+=p[i+2];
00591 Sxx+=p[i]*p[i];
00592 Syy+=p[i+1]*p[i+1];
00593 Sxy+=p[i]*p[i+1];
00594 S+=1.0;
00595 Sxz+=p[i]*p[i+2];
00596 Syz+=p[i+1]*p[i+2];
00597 }
00598 double d=S*Sxy*Sxy - 2*Sx*Sxy*Sy + Sxx*Sy*Sy + Sx*Sx*Syy - S*Sxx*Syy;
00599
00600 Vec3f ret(0,0,0);
00601
00602 ret[0]=static_cast<float>(-((Sxy*Sxz*Sy - Sx*Sxz*Syy + Sx*Sxy*Syz - Sxx*Sy*Syz - Sxy*Sxy*Sz +Sxx*Syy*Sz)/d));
00603 ret[1]=static_cast<float>(-((-Sxz*Sy*Sy + S*Sxz*Syy - S*Sxy*Syz + Sx*Sy*Syz + Sxy*Sy*Sz -Sx*Syy*Sz) /d));
00604 ret[2]=static_cast<float>(-((-S*Sxy*Sxz + Sx*Sxz*Sy - Sx*Sx*Syz + S*Sxx*Syz + Sx*Sxy*Sz -Sxx*Sy*Sz) /d));
00605
00606 return ret;
00607 }
00608
00609 void Util::save_data(const vector < float >&x_array, const vector < float >&y_array,
00610 const string & filename)
00611 {
00612 Assert(x_array.size() > 0);
00613 Assert(y_array.size() > 0);
00614 Assert(filename != "");
00615
00616 if (x_array.size() != y_array.size()) {
00617 LOGERR("array x and array y have different size: %d != %d\n",
00618 x_array.size(), y_array.size());
00619 return;
00620 }
00621
00622 FILE *out = fopen(filename.c_str(), "wb");
00623 if (!out) {
00624 throw FileAccessException(filename);
00625 }
00626
00627 for (size_t i = 0; i < x_array.size(); i++) {
00628 fprintf(out, "%g\t%g\n", x_array[i], y_array[i]);
00629 }
00630 fclose(out);
00631 }
00632
00633 void Util::save_data(float x0, float dx, const vector < float >&y_array,
00634 const string & filename)
00635 {
00636 Assert(dx != 0);
00637 Assert(y_array.size() > 0);
00638 Assert(filename != "");
00639
00640 FILE *out = fopen(filename.c_str(), "wb");
00641 if (!out) {
00642 throw FileAccessException(filename);
00643 }
00644
00645 for (size_t i = 0; i < y_array.size(); i++) {
00646 fprintf(out, "%g\t%g\n", x0 + dx * i, y_array[i]);
00647 }
00648 fclose(out);
00649 }
00650
00651
00652 void Util::save_data(float x0, float dx, float *y_array,
00653 size_t array_size, const string & filename)
00654 {
00655 Assert(dx > 0);
00656 Assert(array_size > 0);
00657 Assert(filename != "");
00658
00659 if (!y_array) {
00660 throw NullPointerException("y array");
00661 }
00662
00663 FILE *out = fopen(filename.c_str(), "wb");
00664 if (!out) {
00665 throw FileAccessException(filename);
00666 }
00667
00668 for (size_t i = 0; i < array_size; i++) {
00669 fprintf(out, "%g\t%g\n", x0 + dx * i, y_array[i]);
00670 }
00671 fclose(out);
00672 }
00673
00674
00675 void Util::sort_mat(float *left, float *right, int *leftPerm, int *rightPerm)
00676
00677 {
00678 float *pivot ; int *pivotPerm;
00679
00680 {
00681 float *pLeft = left; int *pLeftPerm = leftPerm;
00682 float *pRight = right; int *pRightPerm = rightPerm;
00683 float scratch = *left; int scratchPerm = *leftPerm;
00684
00685 while (pLeft < pRight) {
00686 while ((*pRight > scratch) && (pLeft < pRight)) {
00687 pRight--; pRightPerm--;
00688 }
00689 if (pLeft != pRight) {
00690 *pLeft = *pRight; *pLeftPerm = *pRightPerm;
00691 pLeft++; pLeftPerm++;
00692 }
00693 while ((*pLeft < scratch) && (pLeft < pRight)) {
00694 pLeft++; pLeftPerm++;
00695 }
00696 if (pLeft != pRight) {
00697 *pRight = *pLeft; *pRightPerm = *pLeftPerm;
00698 pRight--; pRightPerm--;
00699 }
00700 }
00701 *pLeft = scratch; *pLeftPerm = scratchPerm;
00702 pivot = pLeft; pivotPerm= pLeftPerm;
00703 }
00704 if (left < pivot) {
00705 sort_mat(left, pivot - 1,leftPerm,pivotPerm-1);
00706 }
00707 if (right > pivot) {
00708 sort_mat(pivot + 1, right,pivotPerm+1,rightPerm);
00709 }
00710 }
00711
00712 void Util::set_randnum_seed(unsigned long long seed)
00713 {
00714 Randnum* randnum = Randnum::Instance();
00715 randnum->set_seed(seed);
00716 }
00717
00718 unsigned long long Util::get_randnum_seed()
00719 {
00720 Randnum* randnum = Randnum::Instance();
00721 return randnum->get_seed();
00722 }
00723
00724 int Util::get_irand(int lo, int hi)
00725 {
00726 Randnum* randnum = Randnum::Instance();
00727 return randnum->get_irand(lo, hi);
00728 }
00729
00730 float Util::get_frand(int lo, int hi)
00731 {
00732 return get_frand((float)lo, (float)hi);
00733 }
00734
00735 float Util::get_frand(float lo, float hi)
00736 {
00737 Randnum* randnum = Randnum::Instance();
00738 return randnum->get_frand(lo, hi);
00739 }
00740
00741 float Util::get_frand(double lo, double hi)
00742 {
00743 Randnum* randnum = Randnum::Instance();
00744 return randnum->get_frand(lo, hi);
00745 }
00746
00747 float Util::hypot_fast(int x, int y)
00748 {
00749 static float *mem = (float *)malloc(4*128*128);
00750 static int dim = 0;
00751 x=abs(x);
00752 y=abs(y);
00753
00754 if (x>=dim || y>=dim) {
00755 if (x>2048 || y>2048) return (float)hypot((float)x,(float)y);
00756 dim=x>=dim?x+1:dim;
00757 dim=y>=dim?y+1:dim;
00758 mem=(float*)realloc(mem,4*dim*dim);
00759 for (int y=0; y<dim; y++) {
00760 for (int x=0; x<dim; x++) {
00761 #ifdef _WIN32
00762 mem[x+y*dim]=(float)_hypot((float)x,(float)y);
00763 #else
00764 mem[x+y*dim]=hypot((float)x,(float)y);
00765 #endif
00766 }
00767 }
00768 }
00769
00770 return mem[x+y*dim];
00771 }
00772
00773 short Util::hypot_fast_int(int x, int y)
00774 {
00775 static short *mem = (short *)malloc(2*128*128);
00776 static int dim = 0;
00777 x=abs(x);
00778 y=abs(y);
00779
00780 if (x>=dim || y>=dim) {
00781 if (x>4095 || y>4095) return (short)hypot((float)x,(float)y);
00782 dim=x>=dim?x+1:dim;
00783 dim=y>=dim?y+1:dim;
00784 mem=(short*)realloc(mem,2*dim*dim);
00785 for (int y=0; y<dim; y++) {
00786 for (int x=0; x<dim; x++) {
00787 #ifdef _WIN32
00788 mem[x+y*dim]=(short)Util::round(_hypot((float)x,(float)y));
00789 #else
00790 mem[x+y*dim]=(short)Util::round(hypot((float)x,(float)y));
00791 #endif
00792 }
00793 }
00794 }
00795
00796 return mem[x+y*dim];
00797 }
00798
00799
00800
00801 float Util::fast_exp(const float &f) {
00802 static float *mem = (float *)malloc(sizeof(float)*1000);
00803 static bool needinit=true;
00804
00805 if (needinit) {
00806 needinit=false;
00807 for (int i=0; i<1000; i++) mem[i]=(float)exp(-i/50.0);
00808 }
00809 if (f>0 || f<-19.98) return exp(f);
00810 int g=(int)(-f*50.0+0.5);
00811
00812 return mem[g];
00813 }
00814
00815
00816
00817 float Util::fast_acos(const float &f) {
00818 if (f>=1.0) return 0.0;
00819 if (f<=-1.0) return M_PI;
00820
00821 static float *mem = (float *)malloc(sizeof(float)*2001);
00822 static bool needinit=true;
00823
00824
00825 if (needinit) {
00826 needinit=false;
00827 for (int i=0; i<=2000; i++) mem[i]=(float)acos(i/1000.0-1.0);
00828 }
00829 float f2=f*1000.0f+1000.0f;
00830
00831 int g=(int)(f2+.5);
00832
00833 return mem[g];
00834
00835
00836
00837
00838
00839 }
00840
00841
00842 float Util::get_gauss_rand(float mean, float sigma)
00843 {
00844 Randnum* randnum = Randnum::Instance();
00845 return randnum->get_gauss_rand(mean, sigma);
00846 }
00847
00848 void Util::find_max(const float *data, size_t nitems, float *max_val, int *max_index)
00849 {
00850 Assert(nitems > 0);
00851
00852 if (!data || !max_val || !max_index) {
00853 throw NullPointerException("data/max_val/max_index");
00854 }
00855 float max = -FLT_MAX;
00856 int m = 0;
00857
00858 for (size_t i = 0; i < nitems; i++) {
00859 if (data[i] > max) {
00860 max = data[i];
00861 m = (int)i;
00862 }
00863 }
00864
00865 *max_val = (float)max;
00866
00867 if (max_index) {
00868 *max_index = m;
00869 }
00870 }
00871
00872 void Util::find_min_and_max(const float *data, size_t nitems,
00873 float *max_val, float *min_val,
00874 int *max_index, int *min_index)
00875 {
00876 Assert(nitems > 0);
00877
00878 if (!data || !max_val || !min_val || !max_index || !min_index) {
00879 throw NullPointerException("data/max_val/min_val/max_index/min_index");
00880 }
00881 float max = -FLT_MAX;
00882 float min = FLT_MAX;
00883 int max_i = 0;
00884 int min_i = 0;
00885
00886 for (size_t i = 0; i < nitems; i++) {
00887 if (data[i] > max) {
00888 max = data[i];
00889 max_i = (int)i;
00890 }
00891 if (data[i] < min) {
00892 min = data[i];
00893 min_i = (int)i;
00894 }
00895 }
00896
00897 *max_val = max;
00898 *min_val = min;
00899
00900 if (min_index) {
00901 *min_index = min_i;
00902 }
00903
00904 if (max_index) {
00905 *max_index = max_i;
00906 }
00907
00908 }
00909
00910 Dict Util::get_stats( const vector<float>& data )
00911 {
00912
00913
00914
00915 if (data.size() == 0) EmptyContainerException("Error, attempting to call get stats on an empty container (vector<double>)");
00916
00917 double sum = accumulate(data.begin(), data.end(), 0.0);
00918
00919 double mean = sum / static_cast<double> (data.size());
00920
00921 double std_dev = 0.0, skewness = 0.0, kurtosis = 0.0;
00922
00923 if (data.size() > 1)
00924 {
00925
00926 vector<double> data_mm(data.size());
00927
00928 vector<double> data_mm_sq(data.size());
00929
00930
00931 transform(data.begin(), data.end(), data_mm.begin(), std::bind2nd(std::minus<double>(), mean));
00932
00933
00934 transform(data_mm.begin(), data_mm.end(), data_mm.begin(), data_mm_sq.begin(), std::multiplies<double>());
00935
00936
00937 double square_sum = accumulate(data_mm_sq.begin(), data_mm_sq.end(), 0.0);
00938
00939
00940 std_dev = sqrt(square_sum / static_cast<double>(data.size()-1));
00941 double std_dev_sq = std_dev * std_dev;
00942
00943
00944 double cubic_sum = inner_product(data_mm.begin(), data_mm.end(),data_mm_sq.begin(), 0.0);
00945
00946
00947 double quartic_sum = inner_product(data_mm_sq.begin(), data_mm_sq.end(),data_mm_sq.begin(), 0.0);
00948
00949
00950
00951 skewness = cubic_sum / ((data.size()-1) * std_dev_sq * std_dev );
00952 kurtosis = quartic_sum / ((data.size()-1) * std_dev_sq * std_dev_sq );
00953
00954 }
00955
00956 Dict parms;
00957 parms["mean"] = mean;
00958 parms["std_dev"] = std_dev;
00959 parms["skewness"] = skewness;
00960 parms["kurtosis"] = kurtosis;
00961
00962 return parms;
00963 }
00964
00965
00966 Dict Util::get_stats_cstyle( const vector<float>& data )
00967 {
00968 if (data.size() == 0) EmptyContainerException("Error, attempting to call get stats on an empty container (vector<double>)");
00969
00970 double square_sum = 0.0, sum = 0.0, cube_sum = 0.0, quart_sum = 0.0;
00971 for( vector<float>::const_iterator it = data.begin(); it != data.end(); ++it )
00972 {
00973 double val = *it;
00974 double square = val*val;
00975 quart_sum += square*square;
00976 cube_sum += square*val;
00977 square_sum += square;
00978 sum += val;
00979 }
00980
00981 double mean = sum/(double)data.size();
00982
00983 double std_dev = 0.0, skewness = 0.0, kurtosis = 0.0;
00984
00985 if (data.size() > 1)
00986 {
00987
00988 std_dev = sqrt( (square_sum - mean*sum)/(double)(data.size()-1));
00989
00990 double square_mean = mean*mean;
00991 double cube_mean = mean*square_mean;
00992
00993 double square_std_dev = std_dev*std_dev;
00994
00995
00996
00997 double cubic_sum = cube_sum - 3*square_sum*mean + 3*sum*square_mean - cube_mean*data.size();
00998
00999 skewness = cubic_sum/((data.size()-1)*square_std_dev*std_dev);
01000
01001
01002
01003 double quartic_sum = quart_sum - 4*cube_sum*mean + 6*square_sum*square_mean - 4*sum*cube_mean + square_mean*square_mean*data.size();
01004
01005 kurtosis = quartic_sum /( (data.size()-1)*square_std_dev*square_std_dev);
01006 }
01007
01008 Dict parms;
01009 parms["mean"] = mean;
01010 parms["std_dev"] = std_dev;
01011 parms["skewness"] = skewness;
01012 parms["kurtosis"] = kurtosis;
01013
01014 return parms;
01015 }
01016
01017
01018 int Util::calc_best_fft_size(int low)
01019 {
01020 Assert(low >= 0);
01021
01022
01023 static char *valid = NULL;
01024
01025 if (!valid) {
01026 valid = (char *) calloc(4096, 1);
01027
01028 for (float i2 = 1; i2 < 12.0; i2 += 1.0) {
01029
01030 float f1 = pow((float) 2.0, i2);
01031 for (float i3 = 0; i3 < 8.0; i3 += 1.0) {
01032
01033 float f2 = pow((float) 3.0, i3);
01034 for (float i5 = 0; i5 < 6.0; i5 += 1.0) {
01035
01036 float f3 = pow((float) 5.0, i5);
01037 for (float i7 = 0; i7 < 5.0; i7 += 1.0) {
01038
01039 float f = f1 * f2 * f3 * pow((float) 7.0, i7);
01040 if (f <= 4095.0) {
01041 int n = (int) f;
01042 valid[n] = 1;
01043 }
01044 }
01045 }
01046 }
01047 }
01048 }
01049
01050 for (int i = low; i < 4096; i++) {
01051 if (valid[i]) {
01052 return i;
01053 }
01054 }
01055
01056 LOGERR("Sorry, can only find good fft sizes up to 4096 right now.");
01057
01058 return 1;
01059 }
01060
01061
01062
01063 vector<float> Util::nonconvex(const vector<float>&curve,int first) {
01064
01065 vector<float> ret(curve);
01066 if (first<1) first=1;
01067
01068 int cont=1;
01069 while (cont) {
01070 cont=0;
01071 for (int i=first; i<ret.size()-1; i++) {
01072 float q= (ret[i-1]+ret[i+1])/2.0;
01073 if (ret[i]>q) { ret[i]=q; cont=1; }
01074
01075 }
01076
01077 }
01078
01079 return ret;
01080 }
01081
01082 #ifndef _WIN32
01083 #include <sys/types.h>
01084 #include <sys/socket.h>
01085 #include "byteorder.h"
01086
01087
01088 struct BPKT {
01089 char hdr[4];
01090 int uid;
01091 int len;
01092 int oseq;
01093 int pseq;
01094 unsigned char data[1024]; };
01095
01096 string Util::recv_broadcast(int sock) {
01097
01098
01099
01100
01101 if (ByteOrder::is_host_big_endian()) {
01102 printf("No cache mirroring on Big endian machines yet\n");
01103 return string();
01104 }
01105
01106 BPKT pkt;
01107 string ret;
01108 vector<char> fill;
01109 int obj=-1;
01110 unsigned int i=0;
01111
01112
01113 while (1) {
01114 int l = recv(sock,&pkt,1044,0);
01115 if (l<=0) {
01116 if (obj!=-1) printf("Timeout with incomplete obj %d %d/%d\n",obj,i,(int)fill.size());
01117 return string();
01118 }
01119 if (l<20) {
01120 printf("Bad packet from broadcast");
01121 continue;
01122 }
01123
01124 if (strncmp(pkt.hdr,"EMAN",4)!=0) continue;
01125
01126
01127 if (obj!=pkt.oseq) {
01128 obj=pkt.oseq;
01129 ret.resize(pkt.len);
01130 fill.resize((pkt.len-1)/1024+1);
01131 for (i=0; i<fill.size(); i++) fill[i]=0;
01132 }
01133 if (obj==-1) printf("Something wierd happened. please report\n");
01134
01135
01136 fill[pkt.pseq]=1;
01137 ret.replace(pkt.pseq*1024,l-20,(char *)pkt.data,l-20);
01138
01139
01140 for (i=0; i<fill.size(); i++) {
01141 if (fill[i]!=1) break;
01142 }
01143
01144 fflush(stdout);
01145
01146 if (i==fill.size()) return ret;
01147 }
01148
01149 }
01150 #endif //_WIN32
01151
01152 string Util::get_time_label()
01153 {
01154 time_t t0 = time(0);
01155 struct tm *t = localtime(&t0);
01156 char label[32];
01157 sprintf(label, "%d/%02d/%04d %d:%02d",
01158 t->tm_mon + 1, t->tm_mday, t->tm_year + 1900, t->tm_hour, t->tm_min);
01159 return string(label);
01160 }
01161
01162
01163 void Util::set_log_level(int argc, char *argv[])
01164 {
01165 if (argc > 1 && strncmp(argv[1], "-v", 2) == 0) {
01166 char level_str[32];
01167 strcpy(level_str, argv[1] + 2);
01168 Log::LogLevel log_level = (Log::LogLevel) atoi(level_str);
01169 Log::logger()->set_level(log_level);
01170 }
01171 }
01172
01173 void Util::printMatI3D(MIArray3D& mat, const string str, ostream& out) {
01174
01175
01176 out << "Printing 3D Integer data: " << str << std::endl;
01177 const multi_array_types::size_type* sizes = mat.shape();
01178 int nx = sizes[0], ny = sizes[1], nz = sizes[2];
01179 const multi_array_types::index* indices = mat.index_bases();
01180 int bx = indices[0], by = indices[1], bz = indices[2];
01181 for (int iz = bz; iz < nz+bz; iz++) {
01182 cout << "(z = " << iz << " slice)" << endl;
01183 for (int ix = bx; ix < nx+bx; ix++) {
01184 for (int iy = by; iy < ny+by; iy++) {
01185 cout << setiosflags(ios::fixed) << setw(5)
01186 << mat[ix][iy][iz] << " ";
01187 }
01188 cout << endl;
01189 }
01190 }
01191 }
01192
01193 #include <gsl/gsl_matrix.h>
01194 #include <gsl/gsl_vector.h>
01195 #include <gsl/gsl_linalg.h>
01196
01197 void printmatrix( gsl_matrix* M, const int n, const int m, const string& message = "")
01198 {
01199 cout << message << endl;
01200 for(int i = 0; i < n; ++i ){
01201 for (int j = 0; j < m; ++j ){
01202 cout << gsl_matrix_get(M,i,j) << "\t";
01203 }
01204 cout << endl;
01205 }
01206 }
01207
01208 void printvector( gsl_vector* M, const int n, const string& message = "")
01209 {
01210 cout << message << endl;
01211 for(int i = 0; i < n; ++i ){
01212 cout << gsl_vector_get(M,i) << "\t";
01213 }
01214 cout << endl;
01215 }
01216
01217 float* Util::getBaldwinGridWeights( const int& freq_cutoff, const float& P, const float& r, const float& dfreq, const float& alpha, const float& beta)
01218 {
01219 int i = 0;
01220 int discs = (int)(1+2*freq_cutoff/dfreq);
01221
01222 float* W = new float[discs];
01223
01224 int fc = (int)(2*freq_cutoff + 1);
01225 gsl_matrix* M = gsl_matrix_calloc(fc,fc);
01226
01227 gsl_vector* rhs = gsl_vector_calloc(fc);
01228 cout << i++ << endl;
01229 for(int k = -freq_cutoff; k <= freq_cutoff; ++k){
01230 for(int kp = -freq_cutoff; kp <= freq_cutoff; ++kp){
01231 int kdiff =abs( k-kp);
01232 int evenoddfac = ( kdiff % 2 == 0 ? 1 : -1);
01233
01234 if (kdiff !=0){
01235 float val = sin(M_PI*(float)kdiff*r)/(sin(M_PI*(float)kdiff/(float)P))*(alpha+2.0f*beta*evenoddfac);
01236 gsl_matrix_set(M,int(k+freq_cutoff),int(kp+freq_cutoff),val);
01237 }
01238 }
01239 gsl_matrix_set(M,int(k+freq_cutoff),int(k+freq_cutoff),r*P* (alpha+2*beta));
01240 float val = alpha*sin(M_PI*k*r)/(sin(M_PI*(float)k/(float)P));
01241 if (k!=0) {
01242 gsl_vector_set(rhs,int(k+freq_cutoff),val);
01243 }
01244 }
01245 printmatrix(M,fc,fc,"M");
01246
01247 gsl_vector_set(rhs,int(freq_cutoff),alpha*r*P);
01248 gsl_matrix* V = gsl_matrix_calloc(fc,fc);
01249 gsl_vector* S = gsl_vector_calloc(fc);
01250 gsl_vector* soln = gsl_vector_calloc(fc);
01251 gsl_linalg_SV_decomp(M,V,S,soln);
01252
01253 gsl_linalg_SV_solve(M, V, S, rhs, soln);
01254 printvector(soln,fc,"soln");
01255
01256
01257 int Count=0;
01258 for(float q = (float)(-freq_cutoff); q <= (float)(freq_cutoff); q+= dfreq){
01259 float temp=0;
01260 for(int k = -freq_cutoff; k <= freq_cutoff; ++k){
01261 float dtemp;
01262 if (q!=k) {
01263 dtemp=(1/(float) P)* (float)gsl_vector_get(soln,int(k+freq_cutoff)) * sin(M_PI*(q-k))/sin(M_PI*(q-k)/((float) P));
01264 } else{
01265 dtemp = (1/(float) P)* (float)gsl_vector_get(soln,int(k+freq_cutoff)) * P;
01266 }
01267 temp +=dtemp;
01268 }
01269 W[Count]=temp;
01270 cout << W[Count] << " ";
01271 Count+=1;
01272 }
01273 cout << endl;
01274 return W;
01275 }
01276
01277
01278 void Util::equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane )
01279 {
01280 int x=0,y=1,z=2;
01281 plane[0] = p1[y]*(p2[z]-p3[z])+p2[y]*(p3[z]-p1[z])+p3[y]*(p1[z]-p2[z]);
01282 plane[1] = p1[z]*(p2[x]-p3[x])+p2[z]*(p3[x]-p1[x])+p3[z]*(p1[x]-p2[x]);
01283 plane[2] = p1[x]*(p2[y]-p3[y])+p2[x]*(p3[y]-p1[y])+p3[x]*(p1[y]-p2[y]);
01284 plane[3] = p1[x]*(p2[y]*p3[z]-p3[y]*p2[z])+p2[x]*(p3[y]*p1[z]-p1[y]*p3[z])+p3[x]*(p1[y]*p2[z]-p2[y]*p1[z]);
01285 plane[3] = -plane[3];
01286 }
01287
01288
01289 bool Util::point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3,const Vec2f& point)
01290 {
01291
01292 Vec2f u = p2 - p1;
01293 Vec2f v = p3 - p1;
01294 Vec2f w = point - p1;
01295
01296 float udotu = u.dot(u);
01297 float udotv = u.dot(v);
01298 float udotw = u.dot(w);
01299 float vdotv = v.dot(v);
01300 float vdotw = v.dot(w);
01301
01302 float d = 1.0f/(udotv*udotv - udotu*vdotv);
01303 float s = udotv*vdotw - vdotv*udotw;
01304 s *= d;
01305
01306 float t = udotv*udotw - udotu*vdotw;
01307 t *= d;
01308
01309
01310
01311 if (fabs(s) < Transform::ERR_LIMIT ) s = 0;
01312 if (fabs(t) < Transform::ERR_LIMIT ) t = 0;
01313
01314 if ( fabs((fabs(s)-1.0)) < Transform::ERR_LIMIT ) s = 1;
01315 if ( fabs((fabs(t)-1.0)) < Transform::ERR_LIMIT ) t = 1;
01316
01317
01318
01319
01320 if ( s >= 0 && t >= 0 && (s+t) <= 1 ) return true;
01321 else return false;
01322 }
01323
01324 bool Util::point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point)
01325 {
01326
01327 if (point_is_in_triangle_2d(p1,p2,p4,actual_point)) return true;
01328 else return point_is_in_triangle_2d(p3,p2,p4,actual_point);
01329 }
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358