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
        #pragma warning(disable:4819)
#endif  //_WIN32

#include "sparx/emconstants.h"
#include "exception.h"
#include <vector>
#include <iostream>

#include <string>
using std::string;


#include <boost/multi_array.hpp>
#include <boost/tuple/tuple.hpp>
//#include <boost/math/special_functions/fpclassify.hpp>
#include "vec3.h"

#ifdef WIN32
#include <windows.h>
#define M_PI 3.14159265358979323846f
//#define MAXPATHLEN (MAX_PATH * 4)
#endif

using std::string;
using std::vector;
using std::ostream;
using std::cout;
using std::endl;


namespace EMAN
{
        class EMData;
        class Dict;

        typedef boost::multi_array<int, 3> MIArray3D;

        class Util
        {
                #include "sparx/util_sparx.h"

          public:
                static void ap2ri(float *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 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 int get_randnum_seed();

                static void set_randnum_seed(unsigned long int 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 linear_interpolate(float p1, float p2, float t)
                {
                        return (1-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 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 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 EMData* calc_bessel(const int n, const float& x);

                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 sqrtf(x * x + y * y);
                }

                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 float hypot3(double x, double y, double z)
                {
                        return (float) 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 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_sub_2pi(float x, float y)
                {
                        float r = fmod(fabs(x - y), (float) (2 * 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 int goodf(const float *p_f)
                {
                        //This is the old way to judge a good float, which cause problem on
                        //Fedora Core 64 bit system. Now use isinff() and isnanf() on Linux.
                        //#if defined(_WIN32) || defined(__APPLE__)
                        #if defined(_WIN32)
                        // the first is abnormal zero the second is +-inf or NaN
                        if ((((int *) p_f)[0] & 0x7f800000) == 0 ||
                                (((int *) p_f)[0] & 0x7f800000) == 255) {
                                return 0;
                        }
                        #else
                        using std::fpclassify;
                        int i = fpclassify(*p_f);
                        if ( i == FP_NAN || i == FP_INFINITE) return 0;
                        #endif  //_WIN32 || __APPLE__
                        
                        return 1;
                }
                
                static inline int goodf(const double *p_f)
                {
                        //This is the old way to judge a good float, which cause problem on
                        //Fedora Core 64 bit system. Now use isinff() and isnanf() on Linux.
                        //#if defined(_WIN32) || defined(__APPLE__)
                        #if defined(_WIN32)
                        // the first is abnormal zero the second is +-inf or NaN
                        if ((((int *) p_f)[0] & 0x7f800000) == 0 ||
                                (((int *) p_f)[0] & 0x7f800000) == 255) {
                                return 0;
                        }
                        #else
                        using std::fpclassify;
                        int i = fpclassify(*p_f);
                        if ( i == FP_NAN || i == FP_INFINITE) return 0;
                        #endif  //_WIN32 || __APPLE__
                        
                        return 1;
                }

                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

/* vim: set ts=4 noet nospell: */