EMAN2
util.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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> // using accumulate, inner_product, transform
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) {    //1D, do nothing
00132                 return;
00133         }
00134         else if(ny!=1 && nz==1) {       //2D, rotate vertically by ny/2
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 {  //3D, in the y,z plane, swaps quadrants I,III and II,IV, this is the 'rotation' in y and z
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         // Allocate the working space
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         // This calculates the SVD
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 // various sums used in the final solution
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 // Adapted by PRB from a macro definition posted on SourceForge by evildave
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);           // We won't cache anything bigger than 4096^2
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);           // We won't cache anything bigger than 4096^2
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 // Uses a precached table to return a good approximate to exp(x)
00800 // if outside the cached range, uses regular exp
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 // Uses a precached table to return a good approximate to acos(x)
00816 // tolerates values outside the -1 - 1 domain by clamping to PI,0
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 // This version interpolates, but is slower
00836 /*int g=(int)f2;
00837 f2-=g;
00838 return mem[g+1]*f2+mem[g]*(1.0-f2);*/
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         // Note that this is a heavy STL approach using generic algorithms - some memory could be saved
00913         // using plain c style code, as in get_stats_cstyle below
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                 // read mm is "minus_mean"
00926                 vector<double> data_mm(data.size());
00927                 // read ts as "then squared"
00928                 vector<double> data_mm_sq(data.size());
00929 
00930                 // Subtract the mean from the data and store it in data_mm
00931                 transform(data.begin(), data.end(), data_mm.begin(), std::bind2nd(std::minus<double>(), mean));
00932 
00933                 // Get the square of the data minus the mean and store it in data_mm_sq
00934                 transform(data_mm.begin(), data_mm.end(), data_mm.begin(), data_mm_sq.begin(), std::multiplies<double>());
00935 
00936                 // Get the sum of the squares for the calculation of the standard deviation
00937                 double square_sum = accumulate(data_mm_sq.begin(), data_mm_sq.end(), 0.0);
00938 
00939                 //Calculate teh standard deviation
00940                 std_dev = sqrt(square_sum / static_cast<double>(data.size()-1));
00941                 double std_dev_sq = std_dev * std_dev;
00942 
00943                 // The numerator for the skewness fraction, as defined in http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
00944                 double cubic_sum = inner_product(data_mm.begin(), data_mm.end(),data_mm_sq.begin(), 0.0);
00945 
00946                 // The numerator for the kurtosis fraction, as defined in http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
00947                 double quartic_sum = inner_product(data_mm_sq.begin(), data_mm_sq.end(),data_mm_sq.begin(), 0.0);
00948 
00949                 // Finalize the calculation of the skewness and kurtosis, as defined in
00950                 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
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                 // The standard deviation is calculated here
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                 // This is the numerator of the skewness fraction, if you expand the brackets, as defined in
00996                 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
00997                 double cubic_sum = cube_sum - 3*square_sum*mean + 3*sum*square_mean - cube_mean*data.size();
00998                 // Complete the skewness fraction
00999                 skewness = cubic_sum/((data.size()-1)*square_std_dev*std_dev);
01000 
01001                 // This is the numerator of the kurtosis fraction, if you expand the brackets, as defined in
01002                 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
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                 // Complete the kurtosis fraction
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         //array containing valid sizes <1024 for speed
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 // This takes a 1-D curve and makes it non-convex, by iteratively constraining each point to be
01062 // lower than the mean of the surrounding 2 values
01063 vector<float> Util::nonconvex(const vector<float>&curve,int first) {
01064         
01065         vector<float> ret(curve);
01066         if (first<1) first=1;           // we need one point at each end as an anchor
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 //                      printf("%1.2f (%1.2f) ",ret[i],q);
01075                 }
01076 //              printf("\n");
01077         }
01078         
01079         return ret;
01080 }
01081 
01082 vector<float> Util::windowdot(const vector<float>&curveA,const vector<float>&curveB,int window,int normA) {
01083         if (window<=0) { printf("Warning, %d window in windowdot\n",window); window=1; }
01084         
01085         int ws=window*2+1;
01086         vector<float> ret(curveA.size(),0.0f);
01087         vector<float> suba(ws,0.0f);
01088         vector<float> subb(ws,0.0f);
01089         
01090 
01091         for (int i=window; i<curveA.size()-window; i++) {
01092                 double asig=0,bsig=0,amean=0,bmean=0;
01093                 
01094                 // extract subcurves and compute stats
01095                 for (int j=0,k=i-window; k<=i+window; j++,k++) {
01096                         suba[j]=curveA[k];
01097                         subb[j]=curveB[k];
01098                         amean+=suba[j];
01099                         bmean+=subb[j];
01100                         asig +=suba[j]*suba[j];
01101                         bsig +=subb[j]*subb[j];
01102                 }
01103                 amean/=(double)ws;
01104                 bmean/=(double)ws;
01105                 asig=sqrt(asig/(double)ws-amean*amean);
01106                 bsig=sqrt(bsig/(double)ws-bmean*bmean);
01107 //              printf("%lf %lf %lf %lf\n",amean,asig,bmean,bsig);
01108                 
01109                 double dot=0;
01110                 // normalize vector(s) & compute dot
01111                 for (int j=0; j<ws; j++) {
01112                         subb[j]=(subb[j]-bmean)/bsig;
01113                         if (normA) suba[j]=(suba[j]-amean)/asig;
01114                         
01115                         dot+=suba[j]*subb[j];
01116                 }
01117                 ret[i]=dot/(float)ws;
01118         }
01119         return ret;
01120 }
01121 
01122 
01123 #ifndef _WIN32
01124 #include <sys/types.h>
01125 #include <sys/socket.h>
01126 #include "byteorder.h"
01127 // This is the structure of a broadcast packet
01128 // Transmitter assumes little endian
01129 struct BPKT {
01130         char hdr[4];
01131         int uid;        // user id on transmitter
01132         int len;        // length of total object
01133         int oseq;       // object id
01134         int pseq;       // packet id within object
01135         unsigned char data[1024]; };
01136 
01137 string Util::recv_broadcast(int sock) {
01138 //      struct sockaddr_in sadr = { AF_INET, 9989, INADDR_ANY};
01139 //      int sock=socket(AF_INET,SOCK_DGRAM,0);
01140 //      if (bind(sock,&sadr,sizeof(sockaddr_in))) return string();
01141 
01142         if (ByteOrder::is_host_big_endian()) {
01143                 printf("No cache mirroring on Big endian machines yet\n");
01144                 return string();        // FIXME: no support for big endian hosts
01145         }
01146 
01147         BPKT pkt;
01148         string ret;
01149         vector<char> fill;
01150         int obj=-1;
01151         unsigned int i=0;
01152 //      printf ("Listening\n");
01153 
01154         while (1) {
01155                 int l = recv(sock,&pkt,1044,0);
01156                 if (l<=0) {
01157                         if (obj!=-1) printf("Timeout with incomplete obj %d  %d/%d\n",obj,i,(int)fill.size());
01158                         return string();                // probably a timeout
01159                 }
01160                 if (l<20) {
01161                         printf("Bad packet from broadcast");
01162                         continue;
01163                 }
01164 
01165                 if (strncmp(pkt.hdr,"EMAN",4)!=0) continue;
01166 
01167                 // new object coming in
01168                 if (obj!=pkt.oseq) {
01169                         obj=pkt.oseq;
01170                         ret.resize(pkt.len);
01171                         fill.resize((pkt.len-1)/1024+1);
01172                         for (i=0; i<fill.size(); i++) fill[i]=0;
01173                 }
01174                 if (obj==-1) printf("Something wierd happened. please report\n");
01175 
01176                 // copy the packet into the output buffer
01177                 fill[pkt.pseq]=1;
01178                 ret.replace(pkt.pseq*1024,l-20,(char *)pkt.data,l-20);
01179 
01180                 // see if we got everything
01181                 for (i=0; i<fill.size(); i++) {
01182                         if (fill[i]!=1) break;
01183                 }
01184 //              printf("\t\t\tObj %d  %d/%d      \r",obj,i,(int)fill.size());
01185                 fflush(stdout);
01186 
01187                 if (i==fill.size()) return ret;         // Yea !  We got a good packet
01188         }
01189 
01190 }
01191 #endif  //_WIN32
01192 
01193 string Util::get_time_label()
01194 {
01195         time_t t0 = time(0);
01196         struct tm *t = localtime(&t0);
01197         char label[32];
01198         sprintf(label, "%d/%02d/%04d %d:%02d",
01199                         t->tm_mon + 1, t->tm_mday, t->tm_year + 1900, t->tm_hour, t->tm_min);
01200         return string(label);
01201 }
01202 
01203 
01204 void Util::set_log_level(int argc, char *argv[])
01205 {
01206         if (argc > 1 && strncmp(argv[1], "-v", 2) == 0) {
01207                 char level_str[32];
01208                 strcpy(level_str, argv[1] + 2);
01209                 Log::LogLevel log_level = (Log::LogLevel) atoi(level_str);
01210                 Log::logger()->set_level(log_level);
01211         }
01212 }
01213 
01214 void Util::printMatI3D(MIArray3D& mat, const string str, ostream& out) {
01215         // Note: Don't need to check if 3-D because 3D is part of
01216         //       the MIArray3D typedef.
01217         out << "Printing 3D Integer data: " << str << std::endl;
01218         const multi_array_types::size_type* sizes = mat.shape();
01219         int nx = sizes[0], ny = sizes[1], nz = sizes[2];
01220         const multi_array_types::index* indices = mat.index_bases();
01221         int bx = indices[0], by = indices[1], bz = indices[2];
01222         for (int iz = bz; iz < nz+bz; iz++) {
01223                 cout << "(z = " << iz << " slice)" << endl;
01224                 for (int ix = bx; ix < nx+bx; ix++) {
01225                         for (int iy = by; iy < ny+by; iy++) {
01226                                 cout << setiosflags(ios::fixed) << setw(5)
01227                                          << mat[ix][iy][iz] << "  ";
01228                         }
01229                         cout << endl;
01230                 }
01231         }
01232 }
01233 
01234 #include <gsl/gsl_matrix.h>
01235 #include <gsl/gsl_vector.h>
01236 #include <gsl/gsl_linalg.h>
01237 
01238 void printmatrix( gsl_matrix* M, const int n, const int m, const string& message = "")
01239 {
01240         cout << message << endl;
01241         for(int i = 0; i < n; ++i ){
01242                 for (int j = 0; j < m; ++j ){
01243                         cout << gsl_matrix_get(M,i,j) << "\t";
01244                 }
01245                 cout << endl;
01246         }
01247 }
01248 
01249 void printvector( gsl_vector* M, const int n, const string& message = "")
01250 {
01251         cout << message << endl;
01252         for(int i = 0; i < n; ++i ){
01253                 cout << gsl_vector_get(M,i) << "\t";
01254         }
01255         cout << endl;
01256 }
01257 
01258 float* Util::getBaldwinGridWeights( const int& freq_cutoff, const float& P, const float& r, const float& dfreq, const float& alpha, const float& beta)
01259 {
01260         int i = 0;
01261         int discs = (int)(1+2*freq_cutoff/dfreq);
01262 
01263         float*  W = new float[discs];
01264 
01265         int fc = (int)(2*freq_cutoff + 1);
01266         gsl_matrix* M = gsl_matrix_calloc(fc,fc);
01267 
01268         gsl_vector* rhs = gsl_vector_calloc(fc);
01269         cout << i++ << endl;
01270         for(int k = -freq_cutoff; k <= freq_cutoff; ++k){
01271                 for(int kp = -freq_cutoff; kp <= freq_cutoff; ++kp){
01272                         int kdiff =abs( k-kp);
01273                         int evenoddfac = ( kdiff % 2 == 0 ? 1 : -1);
01274 
01275                         if (kdiff !=0){
01276                                 float val = sin(M_PI*(float)kdiff*r)/(sin(M_PI*(float)kdiff/(float)P))*(alpha+2.0f*beta*evenoddfac);
01277                                 gsl_matrix_set(M,int(k+freq_cutoff),int(kp+freq_cutoff),val);
01278                         }
01279                 }
01280                 gsl_matrix_set(M,int(k+freq_cutoff),int(k+freq_cutoff),r*P* (alpha+2*beta));
01281                 float val = alpha*sin(M_PI*k*r)/(sin(M_PI*(float)k/(float)P));
01282                 if (k!=0) {
01283                         gsl_vector_set(rhs,int(k+freq_cutoff),val);
01284                 }
01285         }
01286         printmatrix(M,fc,fc,"M");
01287 
01288         gsl_vector_set(rhs,int(freq_cutoff),alpha*r*P);
01289         gsl_matrix* V = gsl_matrix_calloc(fc,fc);
01290         gsl_vector* S = gsl_vector_calloc(fc);
01291         gsl_vector* soln = gsl_vector_calloc(fc);
01292         gsl_linalg_SV_decomp(M,V,S,soln);
01293 
01294         gsl_linalg_SV_solve(M, V, S, rhs, soln); // soln now runs from -freq_cutoff to + freq_cutoff
01295         printvector(soln,fc,"soln");
01296 
01297         // we want to solve for W, which ranges from -freq_cutoff to +freq_cutoff in steps of dfreq                            2
01298         int Count=0;
01299         for(float q = (float)(-freq_cutoff); q <= (float)(freq_cutoff); q+= dfreq){
01300                 float temp=0;
01301                 for(int k = -freq_cutoff; k <= freq_cutoff; ++k){
01302                         float dtemp;
01303                         if (q!=k) {
01304                                 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));
01305                         } else{
01306                                 dtemp = (1/(float) P)* (float)gsl_vector_get(soln,int(k+freq_cutoff))  * P;
01307                         }
01308                         temp +=dtemp;
01309                 }
01310                 W[Count]=temp;
01311                 cout << W[Count] << " ";
01312                 Count+=1;
01313         }
01314         cout << endl;
01315         return W;
01316 }
01317 
01318 
01319 void Util::equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane )
01320 {
01321         int x=0,y=1,z=2;
01322         plane[0] = p1[y]*(p2[z]-p3[z])+p2[y]*(p3[z]-p1[z])+p3[y]*(p1[z]-p2[z]);
01323         plane[1] = p1[z]*(p2[x]-p3[x])+p2[z]*(p3[x]-p1[x])+p3[z]*(p1[x]-p2[x]);
01324         plane[2] = p1[x]*(p2[y]-p3[y])+p2[x]*(p3[y]-p1[y])+p3[x]*(p1[y]-p2[y]);
01325         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]);
01326         plane[3] = -plane[3];
01327 }
01328 
01329 
01330 bool Util::point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3,const Vec2f& point)
01331 {
01332 
01333         Vec2f u = p2 - p1;
01334         Vec2f v = p3 - p1;
01335         Vec2f w = point - p1;
01336 
01337         float udotu = u.dot(u);
01338         float udotv = u.dot(v);
01339         float udotw = u.dot(w);
01340         float vdotv = v.dot(v);
01341         float vdotw = v.dot(w);
01342 
01343         float d = 1.0f/(udotv*udotv - udotu*vdotv);
01344         float s = udotv*vdotw - vdotv*udotw;
01345         s *= d;
01346 
01347         float t = udotv*udotw - udotu*vdotw;
01348         t *= d;
01349 
01350         // We've done a few multiplications, so detect when there are tiny residuals that may throw off the final
01351         // decision
01352         if (fabs(s) < Transform::ERR_LIMIT ) s = 0;
01353         if (fabs(t) < Transform::ERR_LIMIT ) t = 0;
01354 
01355         if ( fabs((fabs(s)-1.0)) < Transform::ERR_LIMIT ) s = 1;
01356         if ( fabs((fabs(t)-1.0)) < Transform::ERR_LIMIT ) t = 1;
01357 
01358 //      cout << "s and t are " << s << " " << t << endl;
01359 
01360         // The final decision, if this is true then we've hit the jackpot
01361         if ( s >= 0 && t >= 0 && (s+t) <= 1 ) return true;
01362         else return false;
01363 }
01364 
01365 bool Util::point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point)
01366 {
01367 
01368         if (point_is_in_triangle_2d(p1,p2,p4,actual_point)) return true;
01369         else return point_is_in_triangle_2d(p3,p2,p4,actual_point);
01370 }
01371 
01372 /*
01373 Dict Util::get_isosurface(EMData * image, float surface_palue, bool smooth)
01374 {
01375         if (image->get_ndim() != 3) {
01376                 throw ImageDimensionException("3D only");
01377         }
01378 
01379         MarchingCubes * mc = new MarchingCubes(image, smooth);
01380         mc->set_surface_palue(surface_palue);
01381 
01382         Dict d;
01383         if(smooth) {
01384                 d.put("points", *(mc->get_points()));
01385                 d.put("faces", *(mc->get_faces()));
01386                 d.put("normals", *(mc->get_normalsSm()));
01387         }
01388         else {
01389                 d.put("points", *(mc->get_points()));
01390                 d.put("faces", *(mc->get_faces()));
01391                 d.put("normals", *(mc->get_normals()));
01392         }
01393 
01394         delete mc;
01395         mc = 0;
01396 
01397         return d;
01398 }
01399 */