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