EMAN2
util_sparx.h
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
00007  * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
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 
00040 #ifndef util__sparx_h__
00041 #define util__sparx_h__
00042 
00043 public:
00044 
00045 static int coveig(int n, float *covmat, float *eigval, float *eigvec);
00046 
00048 /* Functions WTF, WTM and BPCQ accept as first parameter images in real and Fourier space. Output image is always in real space. */
00049 static Dict coveig_for_py(int ncov, const vector<float>& covmatpy);
00050 
00051 static void WTF(EMData* PROJ,vector<float> SS,float SNR,int K);
00052 
00053 static void WTM(EMData* PROJ, vector<float> SS,int DIAMETER,int NUMP);
00054 
00055 static Dict CANG(float PHI, float THETA, float PSI);
00056 
00057 static void BPCQ(EMData* B, EMData *CUBE, const int radius);
00058 
00059 static vector<float> infomask(EMData* Vol, EMData* mask, bool);
00060 
00061 static vector<float> helixshiftali(vector<EMData*> ctx, vector<vector<float> > pcoords, int nsegms, float maxincline, int kang, int search_rng, int nxc);
00062 
00063 static vector<float> curhelixshiftali(vector<EMData*> ctx, vector<vector<float> > pcoords, int nsegms,  int search_rng, int nx, int ny);
00064 
00065 static void colreverse(float* beg, float* end, int nx);
00066 
00067 static void slicereverse(float* beg, float* end, int nx,int ny);
00068 
00103 static void cyclicshift(EMData* image, Dict params);
00104 
00105 static Dict im_diff(EMData* V1, EMData* V2, EMData* mask=0);
00106 
00119 static EMData* TwoDTestFunc(int Size, float p, float q,  float a, float b,
00120                    int flag=0, float alphaDeg=0); //PRB
00121 
00122 
00134 static void spline_mat(float *x, float *y, int n,  float *xq, float *yq, int m); //PRB
00135 
00148 static void spline(float *x, float *y, int n, float yp1, float ypn, float *y2);
00149   // PRB
00161 static void splint( float *xa, float *ya, float *y2a, int n,
00162                                      float *xq, float *yq, int m);
00163 
00164 
00175 static void Radialize(int *PermMatTr,  float * kValsSorted,
00176             float *weightofkvalsSorted, int Size, int *SizeReturned);
00177 
00178 
00179 
00180 
00181 class sincBlackman
00182 {
00183         protected:
00184                 int M; 
00185                 float fc; 
00186                 int ntable;
00187                 vector<float> sBtable;
00188                 virtual void build_sBtable(); 
00189                 float fltb;
00190         public:
00191                 sincBlackman(int M_, float fc_, int ntable_ = 1999);
00192                 virtual ~sincBlackman() {};
00193 
00194                 inline  float sBwin_tab(float x) const {
00195                         float xt;
00196                         if(x<0.0f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00197                         return sBtable[ (int) xt];
00198                 }
00200                 int get_sB_size() const { return M; }
00201 };
00202 
00203 
00204 
00228 class KaiserBessel
00229 {
00230         protected:
00231                 float alpha, v, r; 
00232                 int N; 
00233                 int K; 
00234                 float vtable; 
00235                 int ntable;
00236                 vector<float> i0table;
00237                 float dtable; 
00238                 float alphar; 
00239                 float fac; 
00240                 float vadjust;
00241                 float facadj; 
00242                 virtual void build_I0table(); 
00243                 float fltb;
00244         public:
00245                 KaiserBessel(float alpha_, int K, float r_,
00246                                      float v_, int N_, float vtable_=0.f,
00247                                          int ntable_ = 5999);
00248                 virtual ~KaiserBessel() {};
00250                 float I0table_maxerror();
00251                 vector<float> dump_table() {
00252                         return i0table;
00253                 }
00255                 virtual float sinhwin(float x) const;
00257                 virtual float i0win(float x) const;
00259                 inline float i0win_tab(float x) const {
00260                 /*float absx = fabs(x);
00261                         int loc = int(round(absx*fltb));
00262                         return i0table[loc];*/
00263                         float xt;
00264                         if(x<0.f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00265                         return i0table[ (int) xt];
00266                         /*return i0table[ (int) (fabs(x)*fltb+0.5f)];
00267                                 if (absx > vtable) return 0.f;
00268                                 float loc = absx/dtable;
00269                                 return i0table[int(loc + 0.5f)]; */
00270                 }
00272                 int get_window_size() const { return K; }
00274                 class kbsinh_win {
00275                         KaiserBessel& kb;
00276                         public:
00277                         kbsinh_win(KaiserBessel& kb_) : kb(kb_) {}
00278                         float operator()(float x) const {
00279                                 return kb.sinhwin(x);
00280                         }
00281                         int get_window_size() const {return kb.get_window_size();}
00282                 };
00284                 kbsinh_win get_kbsinh_win() {
00285                         return kbsinh_win(*this);
00286                 }
00288                 class kbi0_win {
00289                         KaiserBessel& kb;
00290                         public:
00291                         kbi0_win(KaiserBessel& kb_) : kb(kb_) {}
00292                         float operator()(float x) const {
00293                                 return kb.i0win(x);
00294                         }
00295                         int get_window_size() const {return kb.get_window_size();}
00296                 };
00298                 kbi0_win get_kbi0_win() {
00299                         return kbi0_win(*this);
00300                 }
00301 };
00302 
00303 class FakeKaiserBessel : public KaiserBessel {
00304         public:
00305                 FakeKaiserBessel(float alpha, int K, float r_,
00306                                          float v_, int N_, float vtable_=0.f,
00307                                                  int ntable_ = 5999)
00308         : KaiserBessel(alpha, K, r_, v_, N_, vtable_, ntable_) {
00309                         build_I0table();
00310                 }
00311                 float sinhwin(float x) const;
00312                 float i0win(float x) const;
00313                 void build_I0table();
00314 };
00315 
00327                 static vector<float>
00328                 even_angles(float delta, float t1=0, float t2=90, float p1=0, float p2=359.999);
00329 
00330 
00383                 static float quadri(float x, float y, int nx, int ny, float* image);
00384 
00439                 static float quadri_background(float x, float y, int nx, int ny, float* image, int xnew, int ynew);
00440 
00441                 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
00442                 /* Commented by Zhengfan Yang on 04/20/07
00443                 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
00444                 I made the following changes to get_pixel_conv():
00445                 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
00446                 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
00447                 3. Unfold the 'for' loop
00448                 4. Reduce the usage of multiplications through some bracketing (from 98 times to 57 times in 2D case, from 1029 times to 400 times in 3D case)
00449 
00450                 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
00451                 size, say N=5, you can easily modify it by referring my code.
00452                 */
00453                 static float get_pixel_conv_new(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb);
00454 
00455                 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
00456                 /* Commented by Zhengfan Yang on 04/20/07
00457                 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
00458                 I made the following changes to get_pixel_conv():
00459                 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
00460                 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
00461                 3. Unfold the 'for' loop
00462                 4. Reduce the usage of multiplications through some bracketing (from 98 times to 57 times in 2D case, from 1029 times to 400 times in 3D case)
00463 
00464                 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
00465                 size, say N=5, you can easily modify it by referring my code.
00466                 */
00467         static float get_pixel_conv_new_background(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb, int xnew, int ynew);
00468 
00469                 static std::complex<float> extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb);
00470 
00471                 /*static float quadris(float x, float y, int nx, int ny, float* image);*/
00472                 static float bilinear(float xold, float yold, int nsam, int nrow, float* xim);
00473 
00474 
00483                 static float triquad(float r, float s, float t, float* fdata);
00484 
00492                 class Gaussian {
00493                         float sigma;
00494                         float rttwopisigma;
00495                         float twosigma2;
00496                         public:
00497                         Gaussian(float sigma_ = 1.0) : sigma(sigma_) {
00498                                 rttwopisigma = sqrtf(static_cast<float>(twopi)*sigma);
00499                                 twosigma2 = 2*sigma*sigma;
00500                         }
00501                         inline float operator()(float x) const {
00502                                 return exp(-x*x/(twosigma2))/rttwopisigma;
00503                         }
00504                 };
00505         /*static void alrq(float *xim,  int nsam , int nrow , int *numr,
00506                              float *circ, int lcirc, int nring, char mode);*/
00507         static EMData* Polar2D(EMData* image, vector<int> numr, string mode);
00508         static EMData* Polar2Dm(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode);
00509         /*static void alrq_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
00510                             int  *numr, float *circ, int lcirc, int  nring, char  mode);*/
00511         static void alrl_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
00512                             int  *numr, float *circ, int lcirc, int  nring, char  mode);
00513         /*static void alrq_msi(EMData* image,float cns2, float cnr2,
00514                            int  *numr, float *circ, int lcirc, int  nring, char  mode, Util::KaiserBessel&
00515                                                kb);*/
00516         static EMData* Polar2Dmi(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode, Util::KaiserBessel& kb);
00517 
00518         static void  fftr_q(float  *xcmplx, int nv);
00519         static void  fftr_d(double *xcmplx, int nv);
00520         static void  fftc_q(float  *br, float  *bi, int ln, int ks);
00521         static void  fftc_d(double *br, double *bi, int ln, int ks);
00522 
00524         static void  Frngs(EMData* circ, vector<int> numr);
00525         static void  Normalize_ring(EMData* ring, const vector<int>& numr);
00526 
00528         static void  Frngs_inv(EMData* circ, vector<int> numr);
00529 
00531         static void  Applyws(EMData* circ, vector<int> numr, vector<float> wr);
00532 
00533         /*
00534                 A little note about different Crosrng:
00535                 Basically, they all do cross-correlation function to two images in polar coordinates
00536                 Crosrng_e is the original one
00537                 Crosrng_ew is the one that you could apply weights to different rings
00538                 Crosrng_ms assumes the user already apply weights to circ1, it also returns both
00539                            straight and mirrored positions simultaneously.
00540                 Crosrng_msg differs from the previous ones in that it returns the cross-correlation
00541                             function entirely instead of the peak value and position, thus makes it
00542                             possible to use the gridding method to determine the peak position
00543                 Crosrng_msg_s is same as Crosrng_msg except that it only checks straight position
00544                 Crosrng_msg_m is same as Crosrng_msg except that it only checks mirrored position
00545           */
00546         static Dict Crosrng_e(EMData* circ1, EMData* circ2, vector<int> numr, int neg);
00547         static Dict Crosrng_ew(EMData* circ1, EMData* circ2, vector<int> numr, vector<float> w, int neg);
00548 
00549         static Dict Crosrng_ms(EMData* circ1, EMData* circ2, vector<int> numr);
00550         static Dict Crosrng_ms_delta(EMData* circ1, EMData* circ2, vector<int> numr, float delta_start, float delta);
00551 
00557         static Dict Crosrng_sm_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, int flag, float psimax);
00558 
00564         static Dict Crosrng_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, float psimax);
00565  
00572         static Dict Crosrng_ns(EMData* circ1, EMData* circ2, vector<int> numr);
00573 
00580         static EMData* Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr);
00581 
00588         static void Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t);
00589 
00596         static EMData* Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr);
00597 
00604         static EMData* Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr);
00605 
00606         static vector<float> Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr );
00607         static void  prb1d(double *b, int npoint, float *pos);
00608 
00609         static void update_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00610         static void sub_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00611 
00612         // helper functions for ali2d_ra
00613         static float ener(EMData* ave, vector<int> numr);
00614 
00615         static float ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot);
00616 
00618         static Dict min_dist_real(EMData* image, const vector<EMData*>& data);
00619 
00621         static Dict min_dist_four(EMData* image, const vector<EMData*>& data);
00622 
00629         static int k_means_cont_table_(int* group1, int* group2, int* stb, long int s1, long int s2, int flag);
00630 
00631         // branch and bound matching algorithm
00632 
00633         
00637         static void initial_prune(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T);
00638 
00642         static bool explore(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T, int partref, int* curintx, int
00643                         size_curintx, int* next, int size_next, int depth);
00644 
00645         
00646         
00647 
00652         static int generatesubmax(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int n_guesses, int LARGEST_CLASS);
00653 
00657         static void search2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* matchlist, int* costlist, int J);
00658         
00659         static void explore2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* curintx, int size_curintx, int* next, int size_next, int depth, int J, int* matchlist, int*
00660 costlist, int* curbranch);
00661         
00666         static bool sanitycheck(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* output);
00667 
00677         static vector<int> bb_enumerateMPI_(int* argParts, int* dimClasses, int nParts, int K, int T, int n_guesses, int LARGEST_CLASS, int J, int max_branching, float stmult,
00678         int branchfunc, int LIM);
00679 
00680         
00686         static int* branchMPI(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int curlevel,int n_guesses, int LARGEST_CLASS, int J, int max_branching,
00687         float stmult, int branchfunc, int LIM);
00688         
00689         static int branch_factor_0(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00690         static int branch_factor_2(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00691         static int branch_factor_3(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int K, int LIM);
00692         static int branch_factor_4(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, float stmult);
00693         // new code common-lines
00694 
00695         //helper function for the weights calculation by Voronoi to Cml
00696         static vector<double> cml_weights(const vector<float>& cml);
00697 
00699         static vector<int> cml_line_insino(vector<float> Rot, int i_prj, int n_prj);
00700 
00702         static vector<int> cml_line_insino_all(vector<float> Rot, vector<int> seq, int n_prj, int n_lines);
00703 
00705         static vector<double> cml_init_rot(vector<float> Ori);
00706 
00708         static vector<float> cml_update_rot(vector<float> Rot, int iprj, float nph, float th, float nps);
00709 
00711         static vector<double> cml_line_in3d(vector<float> Ori, vector<int> seq, int nprj, int nlines);
00712 
00714         static vector<double> cml_spin_psi(const vector<EMData*>& data, vector<int> com, vector<float> weights, int iprj, vector<int> iw, int n_psi, int d_psi, int n_prj);
00715         static vector<double> cml_spin_psi_now(const vector<EMData*>& data, vector<int> com, int iprj, vector<int> iw, int n_psi, int d_psi, int n_prj);
00716 
00718         static double cml_disc(const vector<EMData*>& data, vector<int> com, vector<int> seq, vector<float> weights, int n_lines);
00719 
00725         static void set_line(EMData* img, int posline, EMData* line, int offset, int length);
00726 
00734         static void cml_prepare_line(EMData* sino, EMData* line, int ilf, int ihf, int pos_line, int nblines);
00735 
00736         /* Decimates the image with respect to the image center.
00737          * (i.e) the center of the original image is kept the same
00738          * and then the initial start pixel is calculated with respect to the
00739          * center of the image
00740          * @params(image, x-pixel, y-pixel,z-pixel)
00741          * works for all 3 dimensions
00742         **/
00743         static EMData* decimate(EMData* img, int x_step,int y_step=1,int z_step=1);
00744 
00745         static EMData* window(EMData* img,int new_nx ,int new_ny=1, int new_nz=1, int x_offset=0, int y_offset=0, int z_offset=0);
00746 
00747         static EMData* pad(EMData* img, int new_nx, int new_ny=1, int new_nz=1, int x_offset=0, int y_offset=0, int z_offset=0, const char *params="average");
00748 
00749         static vector<float> histogram(EMData* image, EMData* mask, int nbins = 128, float hmin =0.0f, float hmax = 0.0f );
00750 
00751         static Dict histc(EMData *ref,EMData *img,EMData *mask);
00752 
00753         static float hist_comp_freq(float PA,float PB,size_t size_img, int hist_len, EMData *img, vector<float> ref_freq_hist, EMData *mask, float ref_h_diff, float ref_h_min);
00754 
00755 
00756         /* The unit in the ctf function: dz: Angstrom, cs: CM  Ps: Angstrom, Voltage: Kv,dza: Angstrom, azz: degree wgh: None unit. b_factor: Angstrom^2
00757          The CTF function takes form of   *sin(-quadpi*(dz*lambda*ak^2-cs*lambda^3*ak^4/2.)-wgh)*exp(-b_factor*ak^2)*sign
00758           * sign can be set as +1 or -1 . The unit of frequency ak is 1/Angstrom
00759                   Attention: Envelope function in power spectrum has a form of exp(-b_factor*ak^2)
00760                                           */
00761         static float   tf(float dzz, float ak, float voltage = 300.0f, float cs = 2.0f, float wgh = 0.1f, float b_factor = 0.0f, float sign = -1.0f);
00762         static EMData *compress_image_mask(EMData* image, EMData* mask);
00763 
00765         static EMData *reconstitute_image_mask(EMData *image,EMData *mask);
00766 
00767         static vector<float> merge_peaks(vector<float> peak1, vector<float> peak2,float p_size);
00768         static vector<float> pw_extract(vector<float>pw, int n, int iswi,float ps);
00769         static vector<float> call_cl1(long int *k,long int *n, float *ps, long int *iswi, float *pw, float *q2, double *q, double *x, double *res, double *cu, double *s, long int *iu);
00770         static vector<float> lsfit(long int *ks,long int *n, long int *klm2d, long int *iswi, float *q1,double *q, double *x, double *res, double *cu, double *s,long int *iu);
00771         static void cl1(long int *k, long int *l, long int *m, long int *n, long int *klm2d,double *q, double *x, double *res, double *cu, long
00772         int *iu, double *s);
00773         static float eval(char * images,EMData * img, vector<int> S,int N, int K,int size);
00774 
00775         /*  VORONOI DIAGRAM */
00776         static vector<double> vrdg(const vector<float>& ph, const vector<float>& th);
00777         static void hsortd(double *theta,double *phi,int *key,int len,int option);
00778         static void voronoidiag(double *theta,double *phi,double* weight,int n);
00779         /*static void angstep(double* thetast,int len);*/
00780         /*static void voronoi(double *phi,double *theta,double *weight,int lenw,int low,int medium,int nt,int last);*/
00781         static void voronoi(double *phi,double *theta,double *weight, int nt);
00782         static void disorder2(double *x,double *y,int *key,int len);
00783         static void ang_to_xyz(double *x,double *y,double *z,int len);
00784         static void flip23(double *x,double *y,double *z,int *key,int k,int len);
00785         struct tmpstruct{
00786                 double theta1,phi1;
00787                 int key1;
00788                 };
00789         static bool cmp1(tmpstruct tmp1,tmpstruct tmp2);
00790         static bool cmp2(tmpstruct tmp1,tmpstruct tmp2);
00791         /**********************************************************/
00792         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00793         static int trmsh3_(int *n0, double *tol, double *x, double *y, double *z__, int *n, int *list, int *lptr,
00794                int *lend, int *lnew, int *indx, int *lcnt, int *near__, int *next, double *dist, int *ier);
00795         static double areav_(int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier);
00796         /**********************************************************/
00797         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00798 
00799     /*  Various operation on images */
00800         /* out = img + scalar * img1  */
00801         static EMData* madn_scalar(EMData* img, EMData* img1, float scalar);
00802         /* out = scalar * img  */
00803         static EMData* mult_scalar(EMData* img, float scalar);
00804         /* out = img + img1  */
00805         static EMData* addn_img(EMData* img, EMData* img1);
00806         /* out = img - img1  */
00807         static EMData* subn_img(EMData* img, EMData* img1);
00808         /* out = img * img1  */
00809         static EMData* muln_img(EMData* img, EMData* img1);
00810         /* out = img / img1  */
00811         static EMData* divn_img(EMData* img, EMData* img1);
00812         /* out = |img|^2  */
00813         static EMData* squaren_img(EMData* img);
00814         /* img /= Re(img1) with zero check  */
00815         static EMData* divn_filter(EMData* img, EMData* img1);
00816 
00817         /* img += scalar * img1 */
00818         static void mad_scalar(EMData* img, EMData* img1, float scalar);
00819         /* img *= scalar  */
00820         static void mul_scalar(EMData* img, float scalar);
00821         /* img += img1  */
00822         static void add_img(EMData* img, EMData* img1);
00823         /* img += abs(img1)  */
00824         static void add_img_abs(EMData* img, EMData* img1);
00825         /* img += img1**2  */
00826         static void add_img2(EMData* img, EMData* img1);
00827         /* img -= img1  */
00828         static void sub_img(EMData* img, EMData* img1);
00829         /* img *= img1  */
00830         static void mul_img(EMData* img, EMData* img1);
00831         /* img /= img1  */
00832         static void div_img(EMData* img, EMData* img1);
00833         /* img = |img|^2  */
00834         static void square_img(EMData* img);
00835         /* img /= Re(img1) with zero check  */
00836         static void div_filter(EMData* img, EMData* img1);
00837 
00838         //utility for sxlocres
00839         static void set_freq(EMData* freqvol, EMData* temp, EMData* mask, float cutoff, float freq);
00840 
00841 
00842         /* pack absolute values of complex image into  real image with addition of Friedel part  */
00843         static EMData* pack_complex_to_real(EMData* img);
00844 private:
00845         static float ang_n(float peakp, string mode, int maxrin); //this function is used by apmq()
00846 public:
00847 
00853         static vector<float> multiref_polar_ali_2d(EMData* image, const vector< EMData* >& crefim,
00854                 float xrng, float yrng, float step, string mode,
00855                 vector< int >numr, float cnx, float cny);
00856 
00857         /* In this version, we return a list of peaks for all reference images */
00858         static vector<float> multiref_polar_ali_2d_peaklist(EMData* image, const vector< EMData* >& crefim,
00859                 float xrng, float yrng, float step, string mode,
00860                 vector< int >numr, float cnx, float cny);
00861 
00862         /* In this version, we return a list of peaks for local subset of reference images */
00863         static vector<float> multiref_polar_ali_2d_peaklist_local(EMData* image, const vector< EMData* >& crefim,
00864                 float xrng, float yrng, float step, float ant, string mode,
00865                 vector< int >numr, float cnx, float cny);
00866 
00867         /* This is used in ISAC program to assigning particles equally to grops */
00868         static vector<int> assign_groups(std::string matrix_address, int nref, int nima);
00869 
00870         static inline void getvec(float phi, float theta, float& x, float& y, float& z, int option=0) {
00871                 float pi180 = M_PI/180.0f;
00872                 
00873                 if (theta > 180.0f) {
00874                         theta -= 180.0f;
00875                         phi += 180.0f;
00876                 } else if (theta > 90.0f && option == 0) {
00877                         theta = 180.0f - theta;
00878                         phi += 180.0f;
00879                 }
00880 
00881                 phi   *= pi180;
00882                 theta *= pi180;
00883 
00884                 x = sin(theta)*cos(phi);
00885                 y = sin(theta)*sin(phi);
00886                 z = cos(theta);
00887 
00888                 return;
00889         }
00890         
00891         static inline float ang_diff(float v11, float v12, float v13, float v21, float v22, float v23, int& mirror) {
00892                 float v = v11*v21+v12*v22+v13*v23;
00893                 if (v > 1) v = 1;
00894                 if (v < -1) v = -1;
00895                 if (v > 0) { mirror = 1; return acos(v)*180.0f/M_PI; }
00896                 else { mirror = -1; return acos(-v)*180.0f/M_PI; }
00897         }
00898         
00899         /* Find nearest projection angles
00900                 Notice: the input I use is different from python code, which I think is awkward.
00901         */
00902         static int nearest_ang(const vector<float>& vecref, float x, float y, float z);
00903 
00904         /* Assign projection angles to nearest reference projections */
00905         static vector<int> assign_projangles(const vector<float>& projangles, const vector<float>& refangles); 
00906 
00907         /* Assign howmany projection angles to the nearest reference projection */
00908         static vector<int> nearestk_to_refdir(const vector<float>& projangles, const vector<float>& refangles, const int howmany); 
00909 
00910         /* Group projection angles by (phi, theta) */
00911         static vector<int> group_proj_by_phitheta(const vector<float>& projangles, const vector<float>& ref_ang, const int img_per_grp);
00912 
00918         static vector<float> multiref_polar_ali_2d_delta(EMData* image, const vector< EMData* >& crefim,
00919                 float xrng, float yrng, float step, string mode,
00920                 vector< int >numr, float cnx, float cny, float delta_start, float delta);
00921 
00927         static vector<float> multiref_polar_ali_2d_nom(EMData* image, const vector< EMData* >& crefim,
00928                 float xrng, float yrng, float step, string mode,
00929                 vector< int >numr, float cnx, float cny);
00930 
00936         static vector<float> multiref_polar_ali_2d_local(EMData* image, const vector< EMData* >& crefim,
00937                 float xrng, float yrng, float step, float ant, string mode,
00938                 vector< int >numr, float cnx, float cny);
00939         /* Returns first match with peak greater than previousmax or the best match in whole space (when there are no peaks > previousmax).
00940          * The reference rings are checked in random order.
00941          * */
00942         static vector<float> shc(EMData* image, const vector< EMData* >& crefim,
00943                 float xrng, float yrng, float step, float ant, string mode,
00944                 vector< int >numr, float cnx, float cny, string sym);
00945         static vector<float> shc_multipeaks(EMData* image, const vector< EMData* >& crefim,
00946                     float xrng, float yrng, float step, float ant, string mode,
00947                     vector<int>numr, float cnx, float cny, int max_peaks_count);
00948 
00955         static vector<float> multiref_polar_ali_helical(EMData* image, const vector< EMData* >& crefim,
00956                 float xrng, float yrng, float step, float psi_max, string mode,
00957                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00958         static vector<float> multiref_polar_ali_helical_local(EMData* image, const vector< EMData* >& crefim,
00959                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00960                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00961         static vector<float> multiref_polar_ali_helical_90(EMData* image, const vector< EMData* >& crefim,
00962                 float xrng, float yrng, float step, float psi_max, string mode,
00963                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00964         static vector<float> multiref_polar_ali_helical_90_local(EMData* image, const vector< EMData* >& crefim,
00965                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00966                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00968         static vector<float> multiref_polar_ali_helicon_local(EMData* image, const vector< EMData* >& crefim,
00969                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00970                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00971         static vector<float> multiref_polar_ali_helicon_90_local(EMData* image, const vector< EMData* >& crefim,
00972                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00973                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00974 
00980         static vector<float> multiref_polar_ali_2d_local_psi(EMData* image, const vector< EMData* >& crefim,
00981                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00982                 vector< int >numr, float cnx, float cny);
00990         static void multiref_peaks_ali2d(EMData* image, EMData* crefim,
00991                 float xrng, float yrng, float step, string mode,
00992                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm);
00993 
01001         static void multiref_peaks_compress_ali2d(EMData* image, EMData* crefim, float xrng, float yrng,
01002              float step, string mode, vector<int>numr, float cnx, float cny, EMData *peaks, EMData *peakm,
01003              EMData *peaks_compress, EMData *peakm_compress);
01004 
01009         static vector<float> ali2d_ccf_list(EMData* image, EMData* crefim, float xrng, float yrng,
01010              float step, string mode, vector<int>numr, float cnx, float cny, double T);
01011      /*
01012         static void multiref_peaks_ali(EMData* image, const vector< EMData* >& crefim,
01013                 float xrng, float yrng, float step, string mode,
01014                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm,
01015                     int nphi, int ntheta);
01016 */
01017         static vector<float> twoD_fine_ali(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
01018 
01019         static vector<float> twoD_fine_ali_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
01020 
01021         static vector<float> twoD_to_3D_ali(EMData* volft, Util::KaiserBessel& kb, EMData *refim, EMData* mask, float phi, float theta, float psi, float sxs, float sxy);
01022 
01023         static vector<float> twoD_fine_ali_SD(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
01024 
01025         static float ccc_images(EMData *, EMData *, EMData *, float , float , float );
01026 
01027         static vector<float> twoD_fine_ali_SD_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
01028 
01029         static float ccc_images_G(EMData* image, EMData* refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sx, float sy);
01030 
01031         static float local_inner_product(EMData* image1, EMData* image2, int lx, int ly, int lz, int w);
01032 
01033         static EMData* move_points(EMData* img,  float qprob, int ri, int ro);
01034 
01035         static EMData* get_biggest_cluster( EMData* mg );
01036 
01037         //static EMData* ctf_img(int nx, int ny, int nz, float dz, float ps, float voltage=300.0f,float cs=2.0f,float wgh=0.1f,float b_factor=0.0f,float dza=0.0f,float azz=0.0f,float sign=-1.0f);
01038         static EMData* ctf_img(int nx, int ny, int nz, float dz, float ps, float voltage,float cs,float wgh,float b_factor,float dza,float azz,float sign);
01039         static EMData* ctf_rimg(int nx, int ny, int nz, float dz, float ps, float voltage,float cs,float wgh,float b_factor,float dza,float azz,float sign);
01040         static EMData* ctf2_rimg(int nx, int ny, int nz, float dz, float ps, float voltage,float cs,float wgh,float b_factor,float dza,float azz,float sign);
01041 
01042         static inline int mono(int k1, int k2) {
01043 #ifdef _WIN32
01044                 int  mk = _cpp_max(k1,k2);
01045                 return  _cpp_min(k1,k2) + mk*(mk-1)/2;
01046 #else
01047                 int  mk = std::max(k1,k2);
01048                 return  std::min(k1,k2) + mk*(mk-1)/2;
01049 #endif  //_WIN32
01050         }
01051 
01052         static inline int nint180(float arg) {
01053             int res = int(arg + 180.5) - 180;
01054             return res;
01055         }
01056         
01057         static inline double mean(double *x, int n) {
01058                 double s = 0.0;
01059                 for (int i=0; i<n; i++) s+=x[i];
01060                 return s/static_cast<double>(n);
01061         }
01062 
01063         static inline double var(double *x, int n) {
01064                 double s = 0.0;
01065                 double m = mean(x, n);
01066                 for (int i=0; i<n; i++) s += (x[i]-m)*(x[i]-m);
01067                 return s/static_cast<double>(n);
01068         }
01069         
01070         static vector<float> multi_align_error(vector<float> args, vector<float> all_ali_params, int d);
01071         static double multi_align_error_func(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01072         static vector<double> multi_align_error_func2(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01073         static void multi_align_error_dfunc(double* x, vector<float> all_ali_params, int nima, int num_ali, double* g, int d);
01074         
01075         static vector<float> cluster_pairwise(EMData* d, int K, float T, float F);
01076         //static vector<float> cluster_equalsize(EMData* d, int m);
01077         static vector<float> cluster_equalsize(EMData* d);
01078         static vector<float> vareas(EMData* d);
01079 
01085         static EMData* get_slice(EMData *vol, int dim, int index);
01086 
01087         static void image_mutation(EMData *img, float mutation_rate);
01088 
01090         static void array_mutation(float* list, int len_list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01091 
01092         static vector<float> list_mutation(vector<float> list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01093         /*
01094                         To restrict the value to [0, nx)
01095         */
01096         static inline float restrict1(float x, int nx) {
01097                 while ( x < 0.0f )        x += nx;
01098                 while ( x >= (float)(nx) )  x -= nx;
01099                 return x;
01100         }
01101 
01104         static Dict get_transform_params(EMData* image, string xform, string convention);
01105 
01106         static void constrained_helix_exhaustive( vector<EMData*> data, vector<EMData*> fdata, vector<EMData*> refproj, vector<EMData*> rotproj
01107                         , vector<float> dp_dphi_rise_delta, vector<int> nphi_phiwobble_range_ywobble_Dsym_nwx_nwy_nwxc_nwyc
01108                         , bool FindPsi, float psi_max, vector<EMData*> crefim, vector<int> numr, int maxrin, string mode, int cnx, int cny);
01109         
01110         static std::vector<float> diff_between_matrix_of_3D_parameters_angles( std::vector<float> all_params, std::vector<float> rotations );
01111 
01115         static std::vector<int> max_clique(std::vector<int> edges);
01116 
01117 #endif  //util__sparx_h__