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 #ifndef eman__util_h__
00037 #define eman__util_h__ 1
00038
00039 #ifdef _WIN32
00040 #pragma warning(disable:4819)
00041 #include <cfloat>
00042 #endif //_WIN32
00043
00044 #include "sparx/emconstants.h"
00045 #include "exception.h"
00046 #include <vector>
00047 #include <iostream>
00048
00049 #include <string>
00050 using std::string;
00051
00052
00053 #include <boost/multi_array.hpp>
00054 #include <boost/tuple/tuple.hpp>
00055
00056 #include "vec3.h"
00057
00058 #ifdef WIN32
00059 #include <windows.h>
00060 #include <process.h>
00061 #define M_PI 3.14159265358979323846f
00062 #define MUTEX HANDLE
00063
00064 #else
00065 #include <pthread.h>
00066 #define MUTEX pthread_mutex_t
00067 #endif
00068
00069 using std::string;
00070 using std::vector;
00071 using std::ostream;
00072 using std::cout;
00073 using std::endl;
00074
00075
00076 namespace EMAN
00077 {
00078 class EMData;
00079 class Dict;
00080
00081 typedef boost::multi_array<int, 3> MIArray3D;
00082
00086 class Util
00087 {
00089 #include "sparx/util_sparx.h"
00090
00091
00092 public:
00093
00094
00095 static int MUTEX_INIT(MUTEX *mutex);
00096 static int MUTEX_LOCK(MUTEX *mutex);
00097 static int MUTEX_UNLOCK(MUTEX *mutex);
00098
00104 static void ap2ri(float *data, size_t n);
00105
00110 static void flip_complex_phase(float *data, size_t n);
00111
00119 static void rotate_phase_origin(float *data, size_t nx, size_t ny, size_t nz);
00120
00128 static int file_lock_wait(FILE * file);
00129
00130
00136 static bool check_file_by_magic(const void *first_block, const char *magic);
00137
00141 static bool is_file_exist(const string & filename);
00142
00148 static void flip_image(float *data, size_t nx, size_t ny);
00149
00155 static vector<EMData *> svdcmp(const vector<EMData *> &data,int nvec);
00156
00157
00162 static string str_to_lower(const string& s);
00163
00172 static bool sstrncmp(const char *s1, const char *s2);
00173
00178 static string int2str(int n);
00179
00188 static string get_line_from_string(char **str);
00189
00201 static bool get_str_float(const char *str, const char *float_var, float *p_val);
00202
00203
00216 static bool get_str_float(const char *str, const char *float_var,
00217 float *p_v1, float *p_v2);
00218
00237 static bool get_str_float(const char *str, const char *float_var,
00238 int *p_nvalues, float *p_v1, float *p_v2);
00239
00251 static bool get_str_int(const char *str, const char *int_var, int *p_val);
00252
00265 static bool get_str_int(const char *str, const char *int_var,
00266 int *p_v1, int *p_v2);
00267
00286 static bool get_str_int(const char *str, const char *int_var,
00287 int *p_nvalues, int *p_v1, int *p_v2);
00288
00299 static string change_filename_ext(const string& old_filename,
00300 const string & new_ext);
00301
00308 static string remove_filename_ext(const string& filename);
00309
00315 static string get_filename_ext(const string& filename);
00316
00323 static string sbasename(const string & filename);
00324
00337 static void calc_least_square_fit(size_t nitems, const float *data_x,
00338 const float *data_y, float *p_slope,
00339 float *p_intercept, bool ignore_zero,float absmax=0);
00340
00347 static Vec3f calc_bilinear_least_square(const vector<float> &points);
00348
00359 static void save_data(const vector < float >&x_array,
00360 const vector < float >&y_array,
00361 const string & filename);
00362
00371 static void save_data(float x0, float dx,
00372 const vector < float >&y_array,
00373 const string & filename);
00374
00384 static void save_data(float x0, float dx, float *y_array,
00385 size_t array_size, const string & filename);
00386
00387
00396 static void sort_mat(float *left, float *right, int *leftPerm,
00397 int *rightPerm);
00398
00399
00403 static unsigned long long get_randnum_seed();
00404
00408 static void set_randnum_seed(unsigned long long seed);
00409
00415 static int get_irand(int low, int high);
00416
00422 static float get_frand(int low, int high);
00423
00429 static float get_frand(float low, float high);
00430
00436 static float get_frand(double low, double high);
00437
00444 static float get_gauss_rand(float mean, float sigma);
00445
00450 static inline int round(float x)
00451 {
00452 if (x < 0) {
00453 return (int) (x - 0.5f);
00454 }
00455 return (int) (x + 0.5f);
00456 }
00457
00462 static inline int round(double x)
00463 {
00464 if (x < 0) {
00465 return (int) (x - 0.5);
00466 }
00467 return (int) (x + 0.5);
00468 }
00469
00476 static inline float linear_interpolate(float p1, float p2, float t)
00477 {
00478 return (1-t) * p1 + t * p2;
00479 }
00480
00490 static inline float bilinear_interpolate(float p1, float p2, float p3,
00491 float p4, float t, float u)
00492 {
00493 return (1-t) * (1-u) * p1 + t * (1-u) * p2 + (1-t) * u * p3 + t * u * p4;
00494 }
00495
00511 static inline float trilinear_interpolate(float p1, float p2, float p3,
00512 float p4, float p5, float p6,
00513 float p7, float p8, float t,
00514 float u, float v)
00515 {
00516 return ((1 - t) * (1 - u) * (1 - v) * p1 + t * (1 - u) * (1 - v) * p2
00517 + (1 - t) * u * (1 - v) * p3 + t * u * (1 - v) * p4
00518 + (1 - t) * (1 - u) * v * p5 + t * (1 - u) * v * p6
00519 + (1 - t) * u * v * p7 + t * u * v * p8);
00520 }
00521
00528 static void find_max(const float *data, size_t nitems,
00529 float *p_max_val, int *p_max_index = 0);
00530
00541 static void find_min_and_max(const float *data, size_t nitems,
00542 float *p_max_val, float *p_min_val,
00543 int *p_max_index = 0, int *p_min_index = 0);
00544
00545
00550 static Dict get_stats( const vector<float>& data );
00551
00555 static Dict get_stats_cstyle( const vector<float>& data );
00556
00563 static int calc_best_fft_size(int low);
00564
00572 static vector<float> nonconvex(const vector<float>&curve,int first=3);
00573
00574 static EMData* calc_bessel(const int n, const float& x);
00575
00580 static inline int square(int n)
00581 {
00582 return (n * n);
00583 }
00584
00589 static inline float square(float x)
00590 {
00591 return (x * x);
00592 }
00593
00598 static inline float square(double x)
00599 {
00600 return (float)(x * x);
00601 }
00602
00608 static inline float square_sum(float x, float y)
00609 {
00610 return (float)(x * x + y * y);
00611 }
00612
00618 static inline float hypot2(float x, float y)
00619 {
00620
00621 #ifdef _WIN32
00622 return (float) _hypot(x, y);
00623 #else
00624 return (float) hypot(x, y);
00625 #endif //_WIN32
00626
00627 }
00628
00635 static inline int hypot3sq(int x, int y, int z)
00636 {
00637 return ((x * x + y * y + z * z));
00638 }
00639
00646 static inline float hypot3sq(float x, float y, float z)
00647 {
00648 return (x * x + y * y + z * z);
00649 }
00650
00657 static inline float hypot3(int x, int y, int z)
00658 {
00659 return sqrtf((float)(x * x + y * y + z * z));
00660 }
00661
00668 static inline float hypot3(float x, float y, float z)
00669 {
00670 return sqrtf(x * x + y * y + z * z);
00671 }
00672
00679 static inline float hypot3(double x, double y, double z)
00680 {
00681 return (float) sqrt(x * x + y * y + z * z);
00682 }
00683
00689 static float hypot_fast(int x, int y);
00690
00696 static short hypot_fast_int(int x, int y);
00697
00704 static inline int fast_floor(float x)
00705 {
00706 if (x < 0) {
00707 return ((int) x - 1);
00708 }
00709 return (int) x;
00710 }
00716 static float fast_exp(const float &f) ;
00717
00723 static float fast_acos(const float &f) ;
00724
00733 static inline float agauss(float a, float dx, float dy, float dz, float d)
00734 {
00735 return (a * exp(-(dx * dx + dy * dy + dz * dz) / d));
00736 }
00737
00743 static inline int get_min(int f1, int f2)
00744 {
00745 return (f1 < f2 ? f1 : f2);
00746 }
00747
00754 static inline int get_min(int f1, int f2, int f3)
00755 {
00756 if (f1 <= f2 && f1 <= f3) {
00757 return f1;
00758 }
00759 if (f2 <= f1 && f2 <= f3) {
00760 return f2;
00761 }
00762 return f3;
00763 }
00764
00770 static inline float get_min(float f1, float f2)
00771 {
00772 return (f1 < f2 ? f1 : f2);
00773 }
00774
00781 static inline float get_min(float f1, float f2, float f3)
00782 {
00783 if (f1 <= f2 && f1 <= f3) {
00784 return f1;
00785 }
00786 if (f2 <= f1 && f2 <= f3) {
00787 return f2;
00788 }
00789 return f3;
00790 }
00791
00792
00800 static inline float get_min(float f1, float f2, float f3, float f4)
00801 {
00802 float m = f1;
00803 if (f2 < m) {
00804 m = f2;
00805 }
00806 if (f3 < m) {
00807 m = f3;
00808 }
00809 if (f4 < m) {
00810 m = f4;
00811 }
00812 return m;
00813 }
00814
00820 static inline float get_max(float f1, float f2)
00821 {
00822 return (f1 < f2 ? f2 : f1);
00823 }
00824
00831 static inline float get_max(float f1, float f2, float f3)
00832 {
00833 if (f1 >= f2 && f1 >= f3) {
00834 return f1;
00835 }
00836 if (f2 >= f1 && f2 >= f3) {
00837 return f2;
00838 }
00839 return f3;
00840 }
00841
00849 static inline float get_max(float f1, float f2, float f3, float f4)
00850 {
00851 float m = f1;
00852 if (f2 > m) {
00853 m = f2;
00854 }
00855 if (f3 > m) {
00856 m = f3;
00857 }
00858 if (f4 > m) {
00859 m = f4;
00860 }
00861 return m;
00862 }
00863
00871 static inline float angle_sub_2pi(float x, float y)
00872 {
00873 float r = fmod(fabs(x - y), (float) (2 * M_PI));
00874 if (r > M_PI) {
00875 r = (float) (2.0 * M_PI - r);
00876 }
00877
00878 return r;
00879 }
00880
00888 static inline float angle_sub_pi(float x, float y)
00889 {
00890 float r = fmod(fabs(x - y), (float) M_PI);
00891 if (r > M_PI / 2.0) {
00892 r = (float)(M_PI - r);
00893 }
00894 return r;
00895 }
00896
00903 static inline float angle_err_ri(float r1,float i1,float r2,float i2) {
00904 if ((r1==0 && i1==0) || (r2==0 && i2==0)) return 0;
00905
00906
00907 return fast_acos((r1*r2+i1*i2)/(float)(hypot(r1,i1)*hypot(r2,i2)));
00908 }
00909
00916 static inline int goodf(const float *p_f)
00917 {
00918
00919
00920
00921 if (((*((int *)p_f)>>23)&255)==255) return 0;
00922
00923 return 1;
00924 }
00925
00926 static inline int goodf(const double *p_f)
00927 {
00928
00929
00930
00931 if (((*((long long *)p_f)>>52)&2047)==2047) return 0;
00932
00933 return 1;
00934 }
00935
00936 #ifndef _WIN32
00937 static string recv_broadcast(int port);
00938 #endif //_WIN32
00939
00943 static string get_time_label();
00944
00951 static void set_log_level(int argc, char *argv[]);
00952
00963 static inline float eman_copysign(float a, float b)
00964 {
00965 #ifndef WIN32
00966 return copysign(a, b);
00967 #else
00968 int flip = -1;
00969 if ((a <= 0 && b <= 0) || (a > 0 && b > 0)) {
00970 flip = 1;
00971 }
00972 return a * flip;
00973 #endif
00974 }
00975
00988 static inline float eman_erfc(float x)
00989 {
00990 #ifndef WIN32
00991 return (float)erfc(x);
00992 #else
00993 static double a[] = { -1.26551223, 1.00002368,
00994 0.37409196, 0.09678418,
00995 -0.18628806, 0.27886807,
00996 -1.13520398, 1.48851587,
00997 -0.82215223, 0.17087277
00998 };
00999
01000 double result = 1;
01001 double z = fabs(x);
01002 if (z > 0) {
01003 double t = 1 / (1 + 0.5 * z);
01004 double f1 = t * (a[4] + t * (a[5] + t * (a[6] +
01005 t * (a[7] + t * (a[8] + t * a[9])))));
01006 result = t * exp((-z * z) + a[0] + t * (a[1] + t * (a[2] + t * (a[3] + f1))));
01007
01008 if (x < 0) {
01009 result = 2 - result;
01010 }
01011 }
01012 return (float)result;
01013 #endif
01014 }
01015
01027 static void equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane );
01028
01029
01030
01045 static bool point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point);
01046
01058 static bool point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& actual_point);
01059
01066 static void printMatI3D(MIArray3D& mat,
01067 const string str = string(""),
01068 ostream& out = std::cout);
01075 template<class T> static inline T sgn(T& val) {
01076 return (val > 0) ? T(+1) : T(-1);
01077 }
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01102 static float* getBaldwinGridWeights( const int& freq_cutoff, const float& P, const float& r, const float& dfreq = 1, const float& alpha=0.5, const float& beta = 0.2);
01103
01110 static bool IsPower2(int x) {
01111 return ( (x>1) && (x & (x-1))==0 );
01112 }
01113
01114
01115 static void apply_precision(float& value, const float& precision) {
01116 float c = ceilf(value);
01117 float f = (float)fast_floor(value);
01118 if (fabs(value - c) < precision) value = c;
01119 else if (fabs(value - f) < precision) value = f;
01120 }
01121 };
01122 }
01123
01124 #endif
01125
01126