util.h

Go to the documentation of this file.

/*
 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
 * Copyright (c) 2000-2006 Baylor College of Medicine
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * */
 
#ifndef eman__util_h__
#define eman__util_h__ 1
 
#ifdef WIN32
#include <windows.h>
        #include <process.h>
        #define M_PI 3.14159265358979323846f
        #define MUTEX HANDLE
//#define MAXPATHLEN (MAX_PATH * 4)
#else
#include <pthread.h>
#define MUTEX pthread_mutex_t
#endif
 
#include <vector>
#include <iostream>
 
#include <string>
using std::string;
 
#include "sparx/emconstants.h"
#include "exception.h"
 
#include <boost/multi_array.hpp>
#include "vec3.h"
 
using std::string;
using std::vector;
using std::ostream;
using std::cout;
using std::endl;
 
const int good_fft_sizes[] = { 16, 24, 32, 36, 40, 44, 48, 52, 56, 60, 64, 72, 84, 96, 100, 104, 112, 120, 128, 132, 140, 168, 180, 192, 196, 
        208, 216, 220, 224, 240, 256, 260, 288, 300, 320, 352, 360, 384, 416, 440, 448, 480, 512, 540, 560, 576, 588, 600, 630, 640, 648, 672, 
        686, 700, 720, 750, 756, 768, 784, 800, 810, 840, 864, 882, 896, 900, 960, 972, 980, 1000, 1008, 1024, 1050, 1080, 1120, 1134, 1152, 
        1176, 1200, 1250, 1260, 1280, 1296, 1344, 1350, 1372, 1400, 1440, 1458, 1470, 1500, 1512, 1536, 1568, 1600, 1620, 1680, 1728, 1750, 
        1764, 1792, 1800, 1890, 1920, 1944, 1960, 2000, 2016, 2048, 2058, 2100, 2160, 2240, 2250, 2268, 2304, 2352, 2400, 2430, 2450, 2500, 
        2520, 2560, 2592, 2646, 2688, 2700, 2744, 2800, 2880, 2916, 2940, 3000, 3024, 3072, 3136, 3150, 3200, 3240, 3360, 3402, 3430, 3456, 
        3500, 3528, 3584, 3600, 3750, 3780, 3840, 3888, 3920, 4000, 4032, 4050, 4096, 4116, 4200, 4320, 4374, 4410, 4480, 4500, 4536, 4608, 
        4704, 4800, 4802, 4860, 4900, 5000, 5040, 5120, 5184, 5250, 5292, 5376, 5400, 5488, 5600, 5670, 5760, 5832, 5880, 6000, 6048, 6144, 
        6174, 6250, 6272, 6300, 6400, 6480, 6720, 6750, 6804, 6860, 6912, 7000, 7056, 7168, 7200, 7290, 7350, 7500, 7560, 7680, 7776, 7840, 
        7938, 8000, 8064, 8100, 8192, 8232, 8400, 8640, 8748, 8750, 8820, 8960, 9000, 9072, 9216, 9408, 9450, 9600, 9604, 9720, 9800, 10000, 
        10080, 10206, 10240, 10290, 10368, 10500, 10584, 10752, 10800, 10976, 11200, 11250, 11340, 11520, 11664, 11760, 12000, 12096, 
        12150, 12250, 12288, 12348, 12500, 12544, 12600, 12800, 12960, 13122, 13230, 13440, 13500, 13608, 13720, 13824, 14000, 14112, 
        14336, 14400, 14406, 14580, 14700, 15000, 15120, 15360, 15552, 15680, 15750, 15876, 16000, 16128, 16200, 16384 };
 
 
namespace EMAN
{
        class EMData;
        class Dict;
 
        typedef boost::multi_array<int, 3> MIArray3D;
 
        class Util
        {
                #include "sparx/util_sparx.h"
                #include "sphire/util_sphire.h"
                
          public:
 
                //Functions do Cross-Platform Mutex 
                static int MUTEX_INIT(MUTEX *mutex);
                static int MUTEX_LOCK(MUTEX *mutex);
                static int MUTEX_UNLOCK(MUTEX *mutex);
 
                static inline bool is_nan(const float number)
                {
                        return std::isnan(number);
                }
 
                static void ap2ri(float *data, size_t n);
                static void ap2ri(double *data, size_t n);
 
                static void flip_complex_phase(float *data, size_t n);
 
                static void rotate_phase_origin(float *data, size_t nx, size_t ny, size_t nz);
 
                static int file_lock_wait(FILE * file);
                //static void file_unlock(FILE * file);
 
                static bool check_file_by_magic(const void *first_block, const char *magic);
 
                static bool is_file_exist(const string & filename);
 
                static void flip_image(float *data, size_t nx, size_t ny);
 
                static vector<EMData *> svdcmp(const vector<EMData *> &data,int nvec);
 
                static string str_to_lower(const string& s);
 
                static void replace_non_ascii(char *str, int max_size, char repl_char='?');
 
                static bool sstrncmp(const char *s1, const char *s2);
 
                static string int2str(int n);
 
                static string get_line_from_string(char **str);
 
                static bool get_str_float(const char *str, const char *float_var, float *p_val);
 
                static bool get_str_float(const char *str, const char *float_var,
                                                                  float *p_v1, float *p_v2);
 
                static bool get_str_float(const char *str, const char *float_var,
                                                                  int *p_nvalues, float *p_v1, float *p_v2);
 
                static bool get_str_int(const char *str, const char *int_var, int *p_val);
 
                static bool get_str_int(const char *str, const char *int_var,
                                                                int *p_v1, int *p_v2);
 
                static bool get_str_int(const char *str, const char *int_var,
                                                                int *p_nvalues, int *p_v1, int *p_v2);
 
                static string change_filename_ext(const string& old_filename,
                                                                                  const string & new_ext);
 
                static string remove_filename_ext(const string& filename);
 
                static string get_filename_ext(const string& filename);
 
                static string sbasename(const string & filename);
 
                static void calc_least_square_fit(size_t nitems, const float *data_x,
                                                                                  const float *data_y, float *p_slope,
                                                                                  float *p_intercept, bool ignore_zero,float absmax=0);
 
                static Vec3f calc_bilinear_least_square(const vector<float> &points);
 
                static void save_data(const vector < float >&x_array,
                                                          const vector < float >&y_array,
                                                          const string & filename);
 
                static void save_data(float x0, float dx,
                                                          const vector < float >&y_array,
                                                          const string & filename);
 
                static void save_data(float x0, float dx, float *y_array,
                                                          size_t array_size, const string & filename);
 
                static void sort_mat(float *left, float *right, int *leftPerm,
                                                         int *rightPerm);
 
                static unsigned long long get_randnum_seed();
 
                static void set_randnum_seed(unsigned long long seed);
 
                static int get_irand(int low, int high);
 
                static float get_frand(int low, int high);
 
                static float get_frand(float low, float high);
 
                static float get_frand(double low, double high);
 
                static float get_gauss_rand(float mean, float sigma);
 
                static inline int round(float x)
                {
                        if (x < 0) {
                                return (int) (x - 0.5f);
                        }
                        return (int) (x + 0.5f);
                }
 
           static inline int round(double x)
                {
                        if (x < 0) {
                                return (int) (x - 0.5);
                        }
                        return (int) (x + 0.5);
                }
 
                static inline float log_2(float x)
                {
                        return (float) (log(x) / log(2.0f));
                }
 
                static inline double log_2(double x)
                {
                        return (double) (log(x) / log(2.0));
                }
 
                static inline float log_2(int i)
                {
                        return (float) (log((float)i) / log(2.0f));
                }
 
                static inline float linear_interpolate(float p1, float p2, float t)
                {
                        return (1-t) * p1 + t  * p2;
                }
                
                static inline std::complex<float> linear_interpolate_cmplx(std::complex<float> p1, std::complex<float> p2, float t)
                {
                        return (1.0f-t) * p1 + t * p2;
                }
 
                static inline float bilinear_interpolate(float p1, float p2, float p3,
                                                                                                 float p4, float t, float u)
                {
                        return (1-t) * (1-u) * p1 + t * (1-u) * p2 + (1-t) * u * p3 + t * u * p4;
                }
 
                static inline std::complex<float> bilinear_interpolate_cmplx(std::complex<float> p1, std::complex<float> p2, std::complex<float> p3,
                         std::complex<float> p4, float t, float u)
                {
                        return (1.0f-t) * (1.0f-u) * p1 + t * (1.0f-u) * p2 + (1.0f-t) * u * p3 + t * u * p4;
                }
 
                static inline std::complex<float> gauss_interpolate_cmplx(std::complex<float> p1, std::complex<float> p2, std::complex<float> p3,
                         std::complex<float> p4, float t, float u)
                {
                        float c4 = Util::fast_gauss_B2(1.0f-t)*Util::fast_gauss_B2(1.0f-u);
                        float c3 = Util::fast_gauss_B2(t)*Util::fast_gauss_B2(1.0f-u);
                        float c2 = Util::fast_gauss_B2(1.0f-t)*Util::fast_gauss_B2(u);
                        float c1 = Util::fast_gauss_B2(t)*Util::fast_gauss_B2(u);
                        return (c1 * p1 + c2 * p2 + c3 * p3 + c4 * p4)/(c1+c2+c3+c4);
                }
 
                static inline std::complex<float> gauss3_interpolate_cmplx(std::complex<float> *p, float t, float u)
                {
                        std::complex<float> sm=0;
                        float wt=0;
                        for (int y=0,i=0; y<3; y++) {
                                for (int x=0; x<3; x++,i++) {
                                        float c=Util::fast_gauss_B2(t-x)*Util::fast_gauss_B2(u-y);
                                        sm+=c*p[i];
                                        wt+=c;
                                }
                        }
                        return sm/wt;
                }
 
                
                static inline float trilinear_interpolate(float p1, float p2, float p3,
                                                          float p4, float p5, float p6,
                                                         float p7, float p8, float t,
                                                                  float u, float v)
                {
                        return ((1 - t) * (1 - u) * (1 - v) * p1 + t * (1 - u) * (1 - v) * p2
                                        + (1 - t) * u * (1 - v) * p3 + t * u * (1 - v) * p4
                                        + (1 - t) * (1 - u) * v * p5 + t * (1 - u) * v * p6
                                        + (1 - t) * u * v * p7 + t * u * v * p8);
                }
 
                static inline std::complex<float> trilinear_interpolate_complex(std::complex<float> p1, std::complex<float> p2, std::complex<float> p3,
                                                          std::complex<float> p4, std::complex<float> p5, std::complex<float> p6,
                                                         std::complex<float> p7, std::complex<float> p8, float t,
                                                                  float u, float v)
                {
                        return ((1 - t) * (1 - u) * (1 - v) * p1 + t * (1 - u) * (1 - v) * p2
                                        + (1 - t) * u * (1 - v) * p3 + t * u * (1 - v) * p4
                                        + (1 - t) * (1 - u) * v * p5 + t * (1 - u) * v * p6
                                        + (1 - t) * u * v * p7 + t * u * v * p8);
                }
 
                
                static void find_max(const float *data, size_t nitems,
                                                         float *p_max_val, int *p_max_index = 0);
 
                static void find_min_and_max(const float *data, size_t nitems,
                                                                         float *p_max_val, float *p_min_val,
                                                                         int *p_max_index = 0, int *p_min_index = 0);
 
                static Dict get_stats( const vector<float>& data );
 
                static Dict get_stats_cstyle( const vector<float>& data );
 
                static int calc_best_fft_size(int low);
 
                static vector<float> nonconvex(const vector<float>&curve,int first=3);
 
                static vector<float> windowdot(const vector<float>&curveA,const vector<float>&curveB,int window,int normA);
 
                
                static EMData* calc_bessel(const int n, const float& x);
 
                static inline double angle3(double x1,double y1, double z1,double x2, double y2, double z2) {
                        return fast_acos((x1*x2+y1*y2+z1*z2)/(hypot3(x1,y1,z1)*hypot3(x2,y2,z2))); //fast_acos tolerates roundoff error slightly outside -1 to 1 range
                }
 
                static inline float angle3(float x1,float y1, float z1,float x2, float y2, float z2) {
                        return fast_acos((x1*x2+y1*y2+z1*z2)/(hypot3(x1,y1,z1)*hypot3(x2,y2,z2)));
                }
 
                static inline int square(int n)
                {
                        return (n * n);
                }
 
                static inline float square(float x)
                {
                        return (x * x);
                }
 
                static inline float square(double x)
                {
                        return (float)(x * x);
                }
 
                static inline float square_sum(float x, float y)
                {
                        return (float)(x * x + y * y);
                }
 
                static inline float hypot2(float x, float y)
                {
                        return (float) hypot(x, y);
                }
                
                static inline int hypot2sq(int x, int y)
                {
                        return (x * x + y * y);
                }
                
                static inline float hypot2sq(float x, float y)
                {
                        return (float)(x * x + y * y);
                }
 
                static inline int hypot3sq(int x, int y, int z)
                {
                        return ((x * x + y * y + z * z));
                }
 
                static inline float hypot3sq(float x, float y, float z)
                {
                        return (x * x + y * y + z * z);
                }
 
                static inline float hypot3(int x, int y, int z)
                {
                        return sqrtf((float)(x * x + y * y + z * z));
                }
 
                static inline float hypot3(float x, float y, float z)
                {
                        return sqrtf(x * x + y * y + z * z);
                }
 
                static inline double hypot3(double x, double y, double z)
                {
                        return (double) sqrt(x * x + y * y + z * z);
                }
 
                static float hypot_fast(int x, int y);
 
                static short hypot_fast_int(int x, int y);
 
                static inline int fast_floor(float x)
                {
                        if (x < 0) {
                                return ((int) x - 1);
                        }
                        return (int) x;
                }
 
                static float fast_exp(const float &f) ;
 
                static float fast_acos(const float &f) ;
 
                static float fast_gauss_B2(const float &f) ;
 
                
                static inline float agauss(float a, float dx, float dy, float dz, float d)
                {
                        return (a * exp(-(dx * dx + dy * dy + dz * dz) / d));
                }
 
                static inline int get_min(int f1, int f2)
                {
                        return (f1 < f2 ? f1 : f2);
                }
 
                static inline int get_min(int f1, int f2, int f3)
                {
                        if (f1 <= f2 && f1 <= f3) {
                                return f1;
                        }
                        if (f2 <= f1 && f2 <= f3) {
                                return f2;
                        }
                        return f3;
                }
 
                static inline float get_min(float f1, float f2)
                {
                        return (f1 < f2 ? f1 : f2);
                }
 
                static inline float get_min(float f1, float f2, float f3)
                {
                        if (f1 <= f2 && f1 <= f3) {
                                return f1;
                        }
                        if (f2 <= f1 && f2 <= f3) {
                                return f2;
                        }
                        return f3;
                }
 
                static inline float get_min(float f1, float f2, float f3, float f4)
                {
                        float m = f1;
                        if (f2 < m) {
                                m = f2;
                        }
                        if (f3 < m) {
                                m = f3;
                        }
                        if (f4 < m) {
                                m = f4;
                        }
                        return m;
                }
 
                static inline float get_max(float f1, float f2)
                {
                        return (f1 < f2 ? f2 : f1);
                }
 
                static inline float get_max(float f1, float f2, float f3)
                {
                        if (f1 >= f2 && f1 >= f3) {
                                return f1;
                        }
                        if (f2 >= f1 && f2 >= f3) {
                                return f2;
                        }
                        return f3;
                }
 
                static inline float get_max(float f1, float f2, float f3, float f4)
                {
                        float m = f1;
                        if (f2 > m) {
                                m = f2;
                        }
                        if (f3 > m) {
                                m = f3;
                        }
                        if (f4 > m) {
                                m = f4;
                        }
                        return m;
                }
                
                static inline float angle_norm_2pi(float in)
                {
                        float m = fmod((float)in, (float)(2.0*M_PI));
                        
                        return m<0?m+2.0*M_PI:m;
                }
 
                static inline float angle_norm_pi(float in)
                {
                        float m = fmod((float)in, (float)(2.0*M_PI));
                        if (m<-M_PI) m+=2.0*M_PI;
                        return m>M_PI?m-2.0*M_PI:m;
                }
 
 
                static inline float angle_sub_2pi(float x, float y)
                {
                        float r = fmod(fabs(x - y), (float) (2.0 * M_PI));
                        if (r<0) r+=2.0*M_PI;
                        if (r > M_PI) {
                                r = (float) (2.0 * M_PI - r);
                        }
 
                        return r;
                }
 
                static inline float angle_sub_pi(float x, float y)
                {
                        float r = fmod(fabs(x - y), (float) M_PI);
                        if (r > M_PI / 2.0) {
                                r = (float)(M_PI - r);
                        }
                        return r;
                }
 
                static inline float angle_err_ri(float r1,float i1,float r2,float i2) {
                        if ((r1==0 && i1==0) || (r2==0 && i2==0)) return 0;
//                      printf("%f\t%f\t%f\n",r1*r2+i1*i2,hypot(r1,i1),hypot(r2,i2));
//                      return acos((r1*r2+i1*i2)/(hypot(r1,i1)*hypot(r2,i2)));
                        return fast_acos((r1*r2+i1*i2)/(float)(hypot(r1,i1)*hypot(r2,i2)));             // fast_acos also tolerates values very slightly above 1
                }
 
                static inline int goodf(const float *p_f)
                {
                // The old code used 'fpclassify' which was demonstrably broken on many platforms, 
                // causing many irritating problems
 
                        if (((*((int *)p_f)>>23)&255)==255) return 0;
 
                        return 1;
                }
 
                static inline int goodf(const double *p_f)
                {
                // The old code used 'fpclassify' which was demonstrably broken on many platforms, 
                // causing many irritating problems
 
                        if (((*((long long *)p_f)>>52)&2047)==2047) return 0;
 
                        return 1;
                }
 
#ifndef _WIN32
                static string recv_broadcast(int port);
#endif  //_WIN32
 
                static string get_time_label();
 
                static void set_log_level(int argc, char *argv[]);
 
                static inline float eman_copysign(float a, float b)
                {
#ifndef WIN32
                        return copysign(a, b);
#else
                        int flip = -1;
                        if ((a <= 0 && b <= 0) || (a > 0 && b > 0)) {
                                flip = 1;
                        }
                        return a * flip;
#endif
                }
 
                static inline float eman_erfc(float x)
                {
#ifndef WIN32
                        return (float)erfc(x);
#else
                        static double a[] = { -1.26551223, 1.00002368,
                                                                  0.37409196, 0.09678418,
                                                                  -0.18628806, 0.27886807,
                                                                  -1.13520398, 1.48851587,
                                                                  -0.82215223, 0.17087277
                        };
 
                        double result = 1;
                        double z = fabs(x);
                        if (z > 0) {
                                double t = 1 / (1 + 0.5 * z);
                                double f1 = t * (a[4] + t * (a[5] + t * (a[6] +
                                                        t * (a[7] + t * (a[8] + t * a[9])))));
                                result = t * exp((-z * z) + a[0] + t * (a[1] + t * (a[2] + t * (a[3] + f1))));
 
                                if (x < 0) {
                                        result = 2 - result;
                                }
                        }
                        return (float)result;
#endif
                }
 
                static void equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane );
 
                static bool point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point);
 
                static bool point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& actual_point);
 
                static void printMatI3D(MIArray3D& mat,
                                                                const string str = string(""),
                                                                ostream& out = std::cout);
                template<class T> static inline T sgn(T& val) {
                        return (val > 0) ? T(+1) : T(-1);
                }
 
//              /** Get the isosurface value for 3D image.
//               *
//               * @param[in] image 3D image data
//               * @param[in] surface_value threshhold for isosurface valuse
//               * @param[in] smooth boolean to specify whether the smooth value needed
//               *
//               * @return Dict to wrap the points(float array), and triangles(int array)
//               *
//               * @exception ImageDimensionException 3D image only
//               * */
//              static Dict get_isosurface(EMData * image, float surface_value, bool smooth);
 
                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);
 
                static bool IsPower2(int x) {
                        return ( (x>1) && (x & (x-1))==0 );
                }
 
 
                static void apply_precision(float& value, const float& precision) {
                        float c = ceilf(value);
                        float f = (float)fast_floor(value);
                        if (fabs(value - c) < precision) value = c;
                        else if (fabs(value - f) < precision) value = f;
                }
        };
}
 
#endif