EMAN2
util.h
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 #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 //#include <boost/math/special_functions/fpclassify.hpp>
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 //#define MAXPATHLEN (MAX_PATH * 4)
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                 //Functions do Cross-Platform Mutex 
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                 //static void file_unlock(FILE * file);
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 
00584                 static vector<float> windowdot(const vector<float>&curveA,const vector<float>&curveB,int window,int normA);
00585 
00586                 
00587                 static EMData* calc_bessel(const int n, const float& x);
00588 
00590                 inline double angle3(double x1,double y1, double z1,double x2, double y2, double z2) {
00591                         return acos((x1*x2+y1*y2+z1*z2)/(hypot3(x1,y1,z1)*hypot3(x2,y2,z2)));
00592                 }
00593                 
00598                 static inline int square(int n)
00599                 {
00600                         return (n * n);
00601                 }
00602 
00607                 static inline float square(float x)
00608                 {
00609                         return (x * x);
00610                 }
00611 
00616                 static inline float square(double x)
00617                 {
00618                         return (float)(x * x);
00619                 }
00620 
00626                 static inline float square_sum(float x, float y)
00627                 {
00628                         return (float)(x * x + y * y);
00629                 }
00630 
00636                 static inline float hypot2(float x, float y)
00637                 {
00638 //                      return sqrtf(x * x + y * y);
00639 #ifdef  _WIN32
00640                         return (float) _hypot(x, y);
00641 #else
00642                         return (float) hypot(x, y);
00643 #endif  //_WIN32
00644                         
00645                 }
00646 
00653                 static inline int hypot3sq(int x, int y, int z)
00654                 {
00655                         return ((x * x + y * y + z * z));
00656                 }
00657 
00664                 static inline float hypot3sq(float x, float y, float z)
00665                 {
00666                         return (x * x + y * y + z * z);
00667                 }
00668 
00675                 static inline float hypot3(int x, int y, int z)
00676                 {
00677                         return sqrtf((float)(x * x + y * y + z * z));
00678                 }
00679 
00686                 static inline float hypot3(float x, float y, float z)
00687                 {
00688                         return sqrtf(x * x + y * y + z * z);
00689                 }
00690 
00697                 static inline double hypot3(double x, double y, double z)
00698                 {
00699                         return (double) sqrt(x * x + y * y + z * z);
00700                 }
00701 
00707                 static float hypot_fast(int x, int y);
00708 
00714                 static short hypot_fast_int(int x, int y);
00715 
00722                 static inline int fast_floor(float x)
00723                 {
00724                         if (x < 0) {
00725                                 return ((int) x - 1);
00726                         }
00727                         return (int) x;
00728                 }
00734                 static float fast_exp(const float &f) ;
00735 
00741                 static float fast_acos(const float &f) ;
00742 
00751                 static inline float agauss(float a, float dx, float dy, float dz, float d)
00752                 {
00753                         return (a * exp(-(dx * dx + dy * dy + dz * dz) / d));
00754                 }
00755 
00761                 static inline int get_min(int f1, int f2)
00762                 {
00763                         return (f1 < f2 ? f1 : f2);
00764                 }
00765 
00772                 static inline int get_min(int f1, int f2, int f3)
00773                 {
00774                         if (f1 <= f2 && f1 <= f3) {
00775                                 return f1;
00776                         }
00777                         if (f2 <= f1 && f2 <= f3) {
00778                                 return f2;
00779                         }
00780                         return f3;
00781                 }
00782 
00788                 static inline float get_min(float f1, float f2)
00789                 {
00790                         return (f1 < f2 ? f1 : f2);
00791                 }
00792 
00799                 static inline float get_min(float f1, float f2, float f3)
00800                 {
00801                         if (f1 <= f2 && f1 <= f3) {
00802                                 return f1;
00803                         }
00804                         if (f2 <= f1 && f2 <= f3) {
00805                                 return f2;
00806                         }
00807                         return f3;
00808                 }
00809 
00810 
00818                 static inline float get_min(float f1, float f2, float f3, float f4)
00819                 {
00820                         float m = f1;
00821                         if (f2 < m) {
00822                                 m = f2;
00823                         }
00824                         if (f3 < m) {
00825                                 m = f3;
00826                         }
00827                         if (f4 < m) {
00828                                 m = f4;
00829                         }
00830                         return m;
00831                 }
00832 
00838                 static inline float get_max(float f1, float f2)
00839                 {
00840                         return (f1 < f2 ? f2 : f1);
00841                 }
00842 
00849                 static inline float get_max(float f1, float f2, float f3)
00850                 {
00851                         if (f1 >= f2 && f1 >= f3) {
00852                                 return f1;
00853                         }
00854                         if (f2 >= f1 && f2 >= f3) {
00855                                 return f2;
00856                         }
00857                         return f3;
00858                 }
00859 
00867                 static inline float get_max(float f1, float f2, float f3, float f4)
00868                 {
00869                         float m = f1;
00870                         if (f2 > m) {
00871                                 m = f2;
00872                         }
00873                         if (f3 > m) {
00874                                 m = f3;
00875                         }
00876                         if (f4 > m) {
00877                                 m = f4;
00878                         }
00879                         return m;
00880                 }
00881 
00889                 static inline float angle_sub_2pi(float x, float y)
00890                 {
00891                         float r = fmod(fabs(x - y), (float) (2 * M_PI));
00892                         if (r > M_PI) {
00893                                 r = (float) (2.0 * M_PI - r);
00894                         }
00895 
00896                         return r;
00897                 }
00898 
00906                 static inline float angle_sub_pi(float x, float y)
00907                 {
00908                         float r = fmod(fabs(x - y), (float) M_PI);
00909                         if (r > M_PI / 2.0) {
00910                                 r = (float)(M_PI - r);
00911                         }
00912                         return r;
00913                 }
00914 
00921                 static inline float angle_err_ri(float r1,float i1,float r2,float i2) {
00922                         if ((r1==0 && i1==0) || (r2==0 && i2==0)) return 0;
00923 //                      printf("%f\t%f\t%f\n",r1*r2+i1*i2,hypot(r1,i1),hypot(r2,i2));
00924 //                      return acos((r1*r2+i1*i2)/(hypot(r1,i1)*hypot(r2,i2)));
00925                         return fast_acos((r1*r2+i1*i2)/(float)(hypot(r1,i1)*hypot(r2,i2)));             // fast_acos also tolerates values very slightly above 1
00926                 }
00927 
00934                 static inline int goodf(const float *p_f)
00935                 {
00936                 // The old code used 'fpclassify' which was demonstrably broken on many platforms, 
00937                 // causing many irritating problems
00938 
00939                         if (((*((int *)p_f)>>23)&255)==255) return 0;
00940 
00941                         return 1;
00942                 }
00943 
00944                 static inline int goodf(const double *p_f)
00945                 {
00946                 // The old code used 'fpclassify' which was demonstrably broken on many platforms, 
00947                 // causing many irritating problems
00948 
00949                         if (((*((long long *)p_f)>>52)&2047)==2047) return 0;
00950 
00951                         return 1;
00952                 }
00953 
00954 #ifndef _WIN32
00955                 static string recv_broadcast(int port);
00956 #endif  //_WIN32
00957 
00961                 static string get_time_label();
00962 
00969                 static void set_log_level(int argc, char *argv[]);
00970 
00981                 static inline float eman_copysign(float a, float b)
00982                 {
00983 #ifndef WIN32
00984                         return copysign(a, b);
00985 #else
00986                         int flip = -1;
00987                         if ((a <= 0 && b <= 0) || (a > 0 && b > 0)) {
00988                                 flip = 1;
00989                         }
00990                         return a * flip;
00991 #endif
00992                 }
00993 
01006                 static inline float eman_erfc(float x)
01007                 {
01008 #ifndef WIN32
01009                         return (float)erfc(x);
01010 #else
01011                         static double a[] = { -1.26551223, 1.00002368,
01012                                                                   0.37409196, 0.09678418,
01013                                                                   -0.18628806, 0.27886807,
01014                                                                   -1.13520398, 1.48851587,
01015                                                                   -0.82215223, 0.17087277
01016                         };
01017 
01018                         double result = 1;
01019                         double z = fabs(x);
01020                         if (z > 0) {
01021                                 double t = 1 / (1 + 0.5 * z);
01022                                 double f1 = t * (a[4] + t * (a[5] + t * (a[6] +
01023                                                         t * (a[7] + t * (a[8] + t * a[9])))));
01024                                 result = t * exp((-z * z) + a[0] + t * (a[1] + t * (a[2] + t * (a[3] + f1))));
01025 
01026                                 if (x < 0) {
01027                                         result = 2 - result;
01028                                 }
01029                         }
01030                         return (float)result;
01031 #endif
01032                 }
01033 
01045                 static void equation_of_plane(const Vec3f& p1, const Vec3f& p2, const Vec3f& p3, float * plane );
01046 
01047 
01048 
01063                 static bool point_is_in_convex_polygon_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& p4,const Vec2f& actual_point);
01064 
01076                 static bool point_is_in_triangle_2d(const Vec2f& p1, const Vec2f& p2, const Vec2f& p3, const Vec2f& actual_point);
01077 
01084                 static void printMatI3D(MIArray3D& mat,
01085                                                                 const string str = string(""),
01086                                                                 ostream& out = std::cout);
01093                 template<class T> static inline T sgn(T& val) {
01094                         return (val > 0) ? T(+1) : T(-1);
01095                 }
01096 
01097 //              /** Get the isosurface value for 3D image.
01098 //               *
01099 //               * @param[in] image 3D image data
01100 //               * @param[in] surface_value threshhold for isosurface valuse
01101 //               * @param[in] smooth boolean to specify whether the smooth value needed
01102 //               *
01103 //               * @return Dict to wrap the points(float array), and triangles(int array)
01104 //               *
01105 //               * @exception ImageDimensionException 3D image only
01106 //               * */
01107 //              static Dict get_isosurface(EMData * image, float surface_value, bool smooth);
01108 
01120                 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);
01121 
01128                 static bool IsPower2(int x) {
01129                         return ( (x>1) && (x & (x-1))==0 );
01130                 }
01131 
01132 
01133                 static void apply_precision(float& value, const float& precision) {
01134                         float c = ceilf(value);
01135                         float f = (float)fast_floor(value);
01136                         if (fabs(value - c) < precision) value = c;
01137                         else if (fabs(value - f) < precision) value = f;
01138                 }
01139         };
01140 }
01141 
01142 #endif
01143 
01144 /* vim: set ts=4 noet nospell: */