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 #define M_PI 3.14159265358979323846f
00061
00062 #endif
00063
00064 using std::string;
00065 using std::vector;
00066 using std::ostream;
00067 using std::cout;
00068 using std::endl;
00069
00070
00071 namespace EMAN
00072 {
00073 class EMData;
00074 class Dict;
00075
00076 typedef boost::multi_array<int, 3> MIArray3D;
00077
00081 class Util
00082 {
00084 #include "sparx/util_sparx.h"
00085
00086 public:
00092 static void ap2ri(float *data, size_t n);
00093
00098 static void flip_complex_phase(float *data, size_t n);
00099
00107 static void rotate_phase_origin(float *data, size_t nx, size_t ny, size_t nz);
00108
00116 static int file_lock_wait(FILE * file);
00117
00118
00124 static bool check_file_by_magic(const void *first_block, const char *magic);
00125
00129 static bool is_file_exist(const string & filename);
00130
00136 static void flip_image(float *data, size_t nx, size_t ny);
00137
00143 static vector<EMData *> svdcmp(const vector<EMData *> &data,int nvec);
00144
00145
00150 static string str_to_lower(const string& s);
00151
00160 static bool sstrncmp(const char *s1, const char *s2);
00161
00166 static string int2str(int n);
00167
00176 static string get_line_from_string(char **str);
00177
00189 static bool get_str_float(const char *str, const char *float_var, float *p_val);
00190
00191
00204 static bool get_str_float(const char *str, const char *float_var,
00205 float *p_v1, float *p_v2);
00206
00225 static bool get_str_float(const char *str, const char *float_var,
00226 int *p_nvalues, float *p_v1, float *p_v2);
00227
00239 static bool get_str_int(const char *str, const char *int_var, int *p_val);
00240
00253 static bool get_str_int(const char *str, const char *int_var,
00254 int *p_v1, int *p_v2);
00255
00274 static bool get_str_int(const char *str, const char *int_var,
00275 int *p_nvalues, int *p_v1, int *p_v2);
00276
00287 static string change_filename_ext(const string& old_filename,
00288 const string & new_ext);
00289
00296 static string remove_filename_ext(const string& filename);
00297
00303 static string get_filename_ext(const string& filename);
00304
00311 static string sbasename(const string & filename);
00312
00325 static void calc_least_square_fit(size_t nitems, const float *data_x,
00326 const float *data_y, float *p_slope,
00327 float *p_intercept, bool ignore_zero,float absmax=0);
00328
00335 static Vec3f calc_bilinear_least_square(const vector<float> &points);
00336
00347 static void save_data(const vector < float >&x_array,
00348 const vector < float >&y_array,
00349 const string & filename);
00350
00359 static void save_data(float x0, float dx,
00360 const vector < float >&y_array,
00361 const string & filename);
00362
00372 static void save_data(float x0, float dx, float *y_array,
00373 size_t array_size, const string & filename);
00374
00375
00384 static void sort_mat(float *left, float *right, int *leftPerm,
00385 int *rightPerm);
00386
00387
00391 static unsigned long int get_randnum_seed();
00392
00396 static void set_randnum_seed(unsigned long int seed);
00397
00403 static int get_irand(int low, int high);
00404
00410 static float get_frand(int low, int high);
00411
00417 static float get_frand(float low, float high);
00418
00424 static float get_frand(double low, double high);
00425
00432 static float get_gauss_rand(float mean, float sigma);
00433
00438 static inline int round(float x)
00439 {
00440 if (x < 0) {
00441 return (int) (x - 0.5f);
00442 }
00443 return (int) (x + 0.5f);
00444 }
00445
00450 static inline int round(double x)
00451 {
00452 if (x < 0) {
00453 return (int) (x - 0.5);
00454 }
00455 return (int) (x + 0.5);
00456 }
00457
00464 static inline float linear_interpolate(float p1, float p2, float t)
00465 {
00466 return (1-t) * p1 + t * p2;
00467 }
00468
00478 static inline float bilinear_interpolate(float p1, float p2, float p3,
00479 float p4, float t, float u)
00480 {
00481 return (1-t) * (1-u) * p1 + t * (1-u) * p2 + (1-t) * u * p3 + t * u * p4;
00482 }
00483
00499 static inline float trilinear_interpolate(float p1, float p2, float p3,
00500 float p4, float p5, float p6,
00501 float p7, float p8, float t,
00502 float u, float v)
00503 {
00504 return ((1 - t) * (1 - u) * (1 - v) * p1 + t * (1 - u) * (1 - v) * p2
00505 + (1 - t) * u * (1 - v) * p3 + t * u * (1 - v) * p4
00506 + (1 - t) * (1 - u) * v * p5 + t * (1 - u) * v * p6
00507 + (1 - t) * u * v * p7 + t * u * v * p8);
00508 }
00509
00516 static void find_max(const float *data, size_t nitems,
00517 float *p_max_val, int *p_max_index = 0);
00518
00529 static void find_min_and_max(const float *data, size_t nitems,
00530 float *p_max_val, float *p_min_val,
00531 int *p_max_index = 0, int *p_min_index = 0);
00532
00533
00538 static Dict get_stats( const vector<float>& data );
00539
00543 static Dict get_stats_cstyle( const vector<float>& data );
00544
00551 static int calc_best_fft_size(int low);
00552
00553
00554 static EMData* calc_bessel(const int n, const float& x);
00555
00560 static inline int square(int n)
00561 {
00562 return (n * n);
00563 }
00564
00569 static inline float square(float x)
00570 {
00571 return (x * x);
00572 }
00573
00578 static inline float square(double x)
00579 {
00580 return (float)(x * x);
00581 }
00582
00588 static inline float square_sum(float x, float y)
00589 {
00590 return (float)(x * x + y * y);
00591 }
00592
00598 static inline float hypot2(float x, float y)
00599 {
00600 return sqrtf(x * x + y * y);
00601 }
00602
00609 static inline int hypot3sq(int x, int y, int z)
00610 {
00611 return ((x * x + y * y + z * z));
00612 }
00613
00620 static inline float hypot3sq(float x, float y, float z)
00621 {
00622 return (x * x + y * y + z * z);
00623 }
00624
00631 static inline float hypot3(int x, int y, int z)
00632 {
00633 return sqrtf((float)(x * x + y * y + z * z));
00634 }
00635
00642 static inline float hypot3(float x, float y, float z)
00643 {
00644 return sqrtf(x * x + y * y + z * z);
00645 }
00646
00653 static inline float hypot3(double x, double y, double z)
00654 {
00655 return (float) sqrt(x * x + y * y + z * z);
00656 }
00657
00663 static float hypot_fast(int x, int y);
00664
00670 static short hypot_fast_int(int x, int y);
00671
00678 static inline int fast_floor(float x)
00679 {
00680 if (x < 0) {
00681 return ((int) x - 1);
00682 }
00683 return (int) x;
00684 }
00690 static float fast_exp(const float &f) ;
00691
00697 static float fast_acos(const float &f) ;
00698
00707 static inline float agauss(float a, float dx, float dy, float dz, float d)
00708 {
00709 return (a * exp(-(dx * dx + dy * dy + dz * dz) / d));
00710 }
00711
00717 static inline int get_min(int f1, int f2)
00718 {
00719 return (f1 < f2 ? f1 : f2);
00720 }
00721
00728 static inline int get_min(int f1, int f2, int f3)
00729 {
00730 if (f1 <= f2 && f1 <= f3) {
00731 return f1;
00732 }
00733 if (f2 <= f1 && f2 <= f3) {
00734 return f2;
00735 }
00736 return f3;
00737 }
00738
00744 static inline float get_min(float f1, float f2)
00745 {
00746 return (f1 < f2 ? f1 : f2);
00747 }
00748
00755 static inline float get_min(float f1, float f2, float f3)
00756 {
00757 if (f1 <= f2 && f1 <= f3) {
00758 return f1;
00759 }
00760 if (f2 <= f1 && f2 <= f3) {
00761 return f2;
00762 }
00763 return f3;
00764 }
00765
00766
00774 static inline float get_min(float f1, float f2, float f3, float f4)
00775 {
00776 float m = f1;
00777 if (f2 < m) {
00778 m = f2;
00779 }
00780 if (f3 < m) {
00781 m = f3;
00782 }
00783 if (f4 < m) {
00784 m = f4;
00785 }
00786 return m;
00787 }
00788
00794 static inline float get_max(float f1, float f2)
00795 {
00796 return (f1 < f2 ? f2 : f1);
00797 }
00798
00805 static inline float get_max(float f1, float f2, float f3)
00806 {
00807 if (f1 >= f2 && f1 >= f3) {
00808 return f1;
00809 }
00810 if (f2 >= f1 && f2 >= f3) {
00811 return f2;
00812 }
00813 return f3;
00814 }
00815
00823 static inline float get_max(float f1, float f2, float f3, float f4)
00824 {
00825 float m = f1;
00826 if (f2 > m) {
00827 m = f2;
00828 }
00829 if (f3 > m) {
00830 m = f3;
00831 }
00832 if (f4 > m) {
00833 m = f4;
00834 }
00835 return m;
00836 }
00837
00845 static inline float angle_sub_2pi(float x, float y)
00846 {
00847 float r = fmod(fabs(x - y), (float) (2 * M_PI));
00848 if (r > M_PI) {
00849 r = (float) (2.0 * M_PI - r);
00850 }
00851
00852 return r;
00853 }
00854
00862 static inline float angle_sub_pi(float x, float y)
00863 {
00864 float r = fmod(fabs(x - y), (float) M_PI);
00865 if (r > M_PI / 2.0) {
00866 r = (float)(M_PI - r);
00867 }
00868 return r;
00869 }
00870
00877 static inline float angle_err_ri(float r1,float i1,float r2,float i2) {
00878 if ((r1==0 && i1==0) || (r2==0 && i2==0)) return 0;
00879
00880
00881 return fast_acos((r1*r2+i1*i2)/(float)(hypot(r1,i1)*hypot(r2,i2)));
00882 }
00883
00890 static inline int goodf(const float *p_f)
00891 {
00892 #ifdef _WIN32
00893
00894 if ((((int *) p_f)[0] & 0x7f800000) == 0 ||
00895 (((int *) p_f)[0] & 0x7f800000) == 255) {
00896 return 0;
00897 }
00898 #else //_WIN32
00899 using std::fpclassify;
00900 int i = fpclassify(*p_f);
00901 if ( i == FP_NAN || i == FP_INFINITE) return 0;
00902 #endif //_WIN32
00903
00904 return 1;
00905 }
00906
00907 static inline int goodf(const double *p_f)
00908 {
00909 #ifdef _WIN32
00910 int i = _fpclass(*p_f);
00911 if ( i == _FPCLASS_SNAN || i == _FPCLASS_QNAN || i == _FPCLASS_NINF || i == _FPCLASS_PINF) return 0;
00912 #else //_WIN32
00913 using std::fpclassify;
00914 int i = fpclassify(*p_f);
00915 if ( i == FP_NAN || i == FP_INFINITE) return 0;
00916 #endif //_WIN32
00917
00918 return 1;
00919 }
00920
00921 #ifndef _WIN32
00922 static string recv_broadcast(int port);
00923 #endif //_WIN32
00924
00928 static string get_time_label();
00929
00936 static void set_log_level(int argc, char *argv[]);
00937
00948 static inline float eman_copysign(float a, float b)
00949 {
00950 #ifndef WIN32
00951 return copysign(a, b);
00952 #else
00953 int flip = -1;
00954 if ((a <= 0 && b <= 0) || (a > 0 && b > 0)) {
00955 flip = 1;
00956 }
00957 return a * flip;
00958 #endif
00959 }
00960
00973 static inline float eman_erfc(float x)
00974 {
00975 #ifndef WIN32
00976 return (float)erfc(x);
00977 #else
00978 static double a[] = { -1.26551223, 1.00002368,
00979 0.37409196, 0.09678418,
00980 -0.18628806, 0.27886807,
00981 -1.13520398, 1.48851587,
00982 -0.82215223, 0.17087277
00983 };
00984
00985 double result = 1;
00986 double z = fabs(x);
00987 if (z > 0) {
00988 double t = 1 / (1 + 0.5 * z);
00989 double f1 = t * (a[4] + t * (a[5] + t * (a[6] +
00990 t * (a[7] + t * (a[8] + t * a[9])))));
00991 result = t * exp((-z * z) + a[0] + t * (a[1] + t * (a[2] + t * (a[3] + f1))));
00992
00993 if (x < 0) {
00994 result = 2 - result;
00995 }
00996 }
00997 return (float)result;
00998 #endif
00999 }
01000
01012 static void equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane );
01013
01014
01015
01030 static bool point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point);
01031
01043 static bool point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& actual_point);
01044
01051 static void printMatI3D(MIArray3D& mat,
01052 const string str = string(""),
01053 ostream& out = std::cout);
01060 template<class T> static inline T sgn(T& val) {
01061 return (val > 0) ? T(+1) : T(-1);
01062 }
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01087 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);
01088
01095 static bool IsPower2(int x) {
01096 return ( (x>1) && (x & (x-1))==0 );
01097 }
01098
01099
01100 static void apply_precision(float& value, const float& precision) {
01101 float c = ceilf(value);
01102 float f = (float)fast_floor(value);
01103 if (fabs(value - c) < precision) value = c;
01104 else if (fabs(value - f) < precision) value = f;
01105 }
01106 };
01107 }
01108
01109 #endif
01110
01111