Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 void colreverse(float* beg, float* end, int nx);
00062 
00063 static void slicereverse(float* beg, float* end, int nx,int ny);
00064 
00099 static void cyclicshift(EMData* image, Dict params);
00100 
00101 static Dict im_diff(EMData* V1, EMData* V2, EMData* mask=0);
00102 
00115 static EMData* TwoDTestFunc(int Size, float p, float q,  float a, float b,
00116                    int flag=0, float alphaDeg=0); //PRB
00117 
00118 
00130 static void spline_mat(float *x, float *y, int n,  float *xq, float *yq, int m); //PRB
00131 
00144 static void spline(float *x, float *y, int n, float yp1, float ypn, float *y2);
00145   // PRB
00157 static void splint( float *xa, float *ya, float *y2a, int n,
00158                                      float *xq, float *yq, int m);
00159 
00160 
00171 static void Radialize(int *PermMatTr,  float * kValsSorted,
00172             float *weightofkvalsSorted, int Size, int *SizeReturned);
00173 
00174 
00175 
00176 class sincBlackman
00177 {
00178         protected:
00179                 int M; 
00180                 float fc; 
00181                 int ntable;
00182                 vector<float> sBtable;
00183                 virtual void build_sBtable(); 
00184                 float fltb;
00185         public:
00186                 sincBlackman(int M_, float fc_, int ntable_ = 1999);
00187                 virtual ~sincBlackman() {};
00188 
00189                 inline  float sBwin_tab(float x) const {
00190                         float xt;
00191                         if(x<0.0f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00192                         return sBtable[ (int) xt];
00193                 }
00195                 int get_sB_size() const { return M; }
00196 };
00197 
00198 
00199 
00223 class KaiserBessel
00224 {
00225         protected:
00226                 float alpha, v, r; 
00227                 int N; 
00228                 int K; 
00229                 float vtable; 
00230                 int ntable;
00231                 vector<float> i0table;
00232                 float dtable; 
00233                 float alphar; 
00234                 float fac; 
00235                 float vadjust;
00236                 float facadj; 
00237                 virtual void build_I0table(); 
00238                 float fltb;
00239         public:
00240                 KaiserBessel(float alpha_, int K, float r_,
00241                                      float v_, int N_, float vtable_=0.f,
00242                                          int ntable_ = 5999);
00243                 virtual ~KaiserBessel() {};
00245                 float I0table_maxerror();
00246                 vector<float> dump_table() {
00247                         return i0table;
00248                 }
00250                 virtual float sinhwin(float x) const;
00252                 virtual float i0win(float x) const;
00254                 inline float i0win_tab(float x) const {
00255                 /*float absx = fabs(x);
00256                         int loc = int(round(absx*fltb));
00257                         return i0table[loc];*/
00258                         float xt;
00259                         if(x<0.f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00260                         return i0table[ (int) xt];
00261                         /*return i0table[ (int) (fabs(x)*fltb+0.5f)];
00262                                 if (absx > vtable) return 0.f;
00263                                 float loc = absx/dtable;
00264                                 return i0table[int(loc + 0.5f)]; */
00265                 }
00267                 int get_window_size() const { return K; }
00269                 class kbsinh_win {
00270                         KaiserBessel& kb;
00271                         public:
00272                         kbsinh_win(KaiserBessel& kb_) : kb(kb_) {}
00273                         float operator()(float x) const {
00274                                 return kb.sinhwin(x);
00275                         }
00276                         int get_window_size() const {return kb.get_window_size();}
00277                 };
00279                 kbsinh_win get_kbsinh_win() {
00280                         return kbsinh_win(*this);
00281                 }
00283                 class kbi0_win {
00284                         KaiserBessel& kb;
00285                         public:
00286                         kbi0_win(KaiserBessel& kb_) : kb(kb_) {}
00287                         float operator()(float x) const {
00288                                 return kb.i0win(x);
00289                         }
00290                         int get_window_size() const {return kb.get_window_size();}
00291                 };
00293                 kbi0_win get_kbi0_win() {
00294                         return kbi0_win(*this);
00295                 }
00296 };
00297 
00298 class FakeKaiserBessel : public KaiserBessel {
00299         public:
00300                 FakeKaiserBessel(float alpha, int K, float r_,
00301                                          float v_, int N_, float vtable_=0.f,
00302                                                  int ntable_ = 5999)
00303         : KaiserBessel(alpha, K, r_, v_, N_, vtable_, ntable_) {
00304                         build_I0table();
00305                 }
00306                 float sinhwin(float x) const;
00307                 float i0win(float x) const;
00308                 void build_I0table();
00309 };
00310 
00322                 static vector<float>
00323                 even_angles(float delta, float t1=0, float t2=90, float p1=0, float p2=359.999);
00324 
00325 
00378                 static float quadri(float x, float y, int nx, int ny, float* image);
00379 
00434                 static float quadri_background(float x, float y, int nx, int ny, float* image, int xnew, int ynew);
00435 
00436                 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
00437                 /* Commented by Zhengfan Yang on 04/20/07
00438                 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
00439                 I made the following changes to get_pixel_conv():
00440                 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
00441                 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
00442                 3. Unfold the 'for' loop
00443                 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)
00444 
00445                 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
00446                 size, say N=5, you can easily modify it by referring my code.
00447                 */
00448                 static float get_pixel_conv_new(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb);
00449 
00450                 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
00451                 /* Commented by Zhengfan Yang on 04/20/07
00452                 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
00453                 I made the following changes to get_pixel_conv():
00454                 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
00455                 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
00456                 3. Unfold the 'for' loop
00457                 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)
00458 
00459                 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
00460                 size, say N=5, you can easily modify it by referring my code.
00461                 */
00462         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);
00463 
00464                 static std::complex<float> extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb);
00465 
00466                 /*static float quadris(float x, float y, int nx, int ny, float* image);*/
00467                 static float bilinear(float xold, float yold, int nsam, int nrow, float* xim);
00468 
00469 
00478                 static float triquad(float r, float s, float t, float* fdata);
00479 
00487                 class Gaussian {
00488                         float sigma;
00489                         float rttwopisigma;
00490                         float twosigma2;
00491                         public:
00492                         Gaussian(float sigma_ = 1.0) : sigma(sigma_) {
00493                                 rttwopisigma = sqrtf(static_cast<float>(twopi)*sigma);
00494                                 twosigma2 = 2*sigma*sigma;
00495                         }
00496                         inline float operator()(float x) const {
00497                                 return exp(-x*x/(twosigma2))/rttwopisigma;
00498                         }
00499                 };
00500         /*static void alrq(float *xim,  int nsam , int nrow , int *numr,
00501                              float *circ, int lcirc, int nring, char mode);*/
00502         static EMData* Polar2D(EMData* image, vector<int> numr, string mode);
00503         static EMData* Polar2Dm(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode);
00504         /*static void alrq_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
00505                             int  *numr, float *circ, int lcirc, int  nring, char  mode);*/
00506         static void alrl_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
00507                             int  *numr, float *circ, int lcirc, int  nring, char  mode);
00508         /*static void alrq_msi(EMData* image,float cns2, float cnr2,
00509                            int  *numr, float *circ, int lcirc, int  nring, char  mode, Util::KaiserBessel&
00510                                                kb);*/
00511         static EMData* Polar2Dmi(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode, Util::KaiserBessel& kb);
00512 
00513         static void  fftr_q(float  *xcmplx, int nv);
00514         static void  fftr_d(double *xcmplx, int nv);
00515         static void  fftc_q(float  *br, float  *bi, int ln, int ks);
00516         static void  fftc_d(double *br, double *bi, int ln, int ks);
00517 
00519         static void  Frngs(EMData* circ, vector<int> numr);
00520         static void  Normalize_ring(EMData* ring, const vector<int>& numr);
00521 
00523         static void  Frngs_inv(EMData* circ, vector<int> numr);
00524 
00526         static void  Applyws(EMData* circ, vector<int> numr, vector<float> wr);
00527 
00528         /*
00529                 A little notes about different Crosrng:
00530                 Basically, they all do cross-correlation function to two images in polar coordinates
00531                 Crosrng_e is the original one
00532                 Crosrng_ew is the one that you could apply weights to different rings
00533                 Crosrng_ms assumes the user already apply weights to circ1, it also returns both
00534                            straight and mirrored positions simultaneously.
00535                 Crosrng_msg differs from the previous ones in that it returns the cross-correlation
00536                             function entirely instead of the peak value and position, thus makes it
00537                             possible to use the gridding method to determine the peak position
00538                 Crosrng_msg_s is same as Crosrng_msg except that it only checks straight position
00539                 Crosrng_msg_m is same as Crosrng_msg except that it only checks mirrored position
00540           */
00541         static Dict Crosrng_e(EMData* circ1, EMData* circ2, vector<int> numr, int neg);
00542         static Dict Crosrng_ew(EMData* circ1, EMData* circ2, vector<int> numr, vector<float> w, int neg);
00543 
00544         static Dict Crosrng_ms(EMData* circ1, EMData* circ2, vector<int> numr);
00545         static Dict Crosrng_ms_delta(EMData* circ1, EMData* circ2, vector<int> numr, float delta_start, float delta);
00546 
00552         static Dict Crosrng_sm_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, int flag, float psimax);
00553 
00559         static Dict Crosrng_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, float psimax);
00560  
00567         static Dict Crosrng_psi_0_180_no_mirror(EMData* circ1, EMData* circ2, vector<int> numr, float psi_max);
00568         static Dict Crosrng_ns(EMData* circ1, EMData* circ2, vector<int> numr);
00569 
00576         static EMData* Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr);
00577 
00584         static void Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t);
00585 
00592         static EMData* Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr);
00593 
00600         static EMData* Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr);
00601 
00602         static vector<float> Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr );
00603         static void  prb1d(double *b, int npoint, float *pos);
00604 
00605         static void update_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00606         static void sub_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00607 
00608         // helper functions for ali2d_ra
00609         static float ener(EMData* ave, vector<int> numr);
00610 
00611         static float ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot);
00612 
00614         static Dict min_dist_real(EMData* image, const vector<EMData*>& data);
00615 
00617         static Dict min_dist_four(EMData* image, const vector<EMData*>& data);
00618 
00625         static int k_means_cont_table_(int* group1, int* group2, int* stb, long int s1, long int s2, int flag);
00626 
00627         // branch and bound matching algorithm
00628 
00629         
00633         static void initial_prune(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T);
00634 
00638         static bool explore(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T, int partref, int* curintx, int
00639                         size_curintx, int* next, int size_next, int depth);
00640 
00641         
00642         
00643 
00648         static int generatesubmax(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int n_guesses, int LARGEST_CLASS);
00649 
00653         static void search2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* matchlist, int* costlist, int J);
00654         
00655         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*
00656 costlist, int* curbranch);
00657         
00662         static bool sanitycheck(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* output);
00663 
00673         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,
00674         int branchfunc, int LIM);
00675 
00676         
00682         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,
00683         float stmult, int branchfunc, int LIM);
00684         
00685         static int branch_factor_0(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00686         static int branch_factor_2(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00687         static int branch_factor_3(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int K, int LIM);
00688         static int branch_factor_4(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, float stmult);
00689         // new code common-lines
00690 
00691         //helper function for the weights calculation by Voronoi to Cml
00692         static vector<double> cml_weights(const vector<float>& cml);
00693 
00695         static vector<int> cml_line_insino(vector<float> Rot, int i_prj, int n_prj);
00696 
00698         static vector<int> cml_line_insino_all(vector<float> Rot, vector<int> seq, int n_prj, int n_lines);
00699 
00701         static vector<double> cml_init_rot(vector<float> Ori);
00702 
00704         static vector<float> cml_update_rot(vector<float> Rot, int iprj, float nph, float th, float nps);
00705 
00707         static vector<double> cml_line_in3d(vector<float> Ori, vector<int> seq, int nprj, int nlines);
00708 
00710         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);
00711         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);
00712 
00714         static double cml_disc(const vector<EMData*>& data, vector<int> com, vector<int> seq, vector<float> weights, int n_lines);
00715 
00721         static void set_line(EMData* img, int posline, EMData* line, int offset, int length);
00722 
00730         static void cml_prepare_line(EMData* sino, EMData* line, int ilf, int ihf, int pos_line, int nblines);
00731 
00732         /* Decimates the image with respect to the image center.
00733          * (i.e) the center of the original image is kept the same
00734          * and then the initial start pixel is calculated with respect to the
00735          * center of the image
00736          * @params(image, x-pixel, y-pixel,z-pixel)
00737          * works for all 3 dimensions
00738         **/
00739         static EMData* decimate(EMData* img, int x_step,int y_step=1,int z_step=1);
00740 
00741         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);
00742 
00743         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");
00744 
00745         static vector<float> histogram(EMData* image, EMData* mask, int nbins = 128, float hmin =0.0f, float hmax = 0.0f );
00746 
00747         static Dict histc(EMData *ref,EMData *img,EMData *mask);
00748 
00749         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);
00750 
00751 
00752         /* 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
00753          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
00754           * sign can be set as +1 or -1 . The unit of frequency ak is 1/Angstrom
00755                   Attention: Envelope function in power spectrum has a form of exp(-b_factor*ak^2)
00756                                           */
00757         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);
00758         static EMData *compress_image_mask(EMData* image, EMData* mask);
00759 
00761         static EMData *reconstitute_image_mask(EMData *image,EMData *mask);
00762 
00763         static vector<float> merge_peaks(vector<float> peak1, vector<float> peak2,float p_size);
00764         static vector<float> pw_extract(vector<float>pw, int n, int iswi,float ps);
00765         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);
00766         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);
00767         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
00768         int *iu, double *s);
00769         static float eval(char * images,EMData * img, vector<int> S,int N, int K,int size);
00770 
00771         /*  VORONOI DIAGRAM */
00772         static vector<double> vrdg(const vector<float>& ph, const vector<float>& th);
00773         static void hsortd(double *theta,double *phi,int *key,int len,int option);
00774         static void voronoidiag(double *theta,double *phi,double* weight,int n);
00775         /*static void angstep(double* thetast,int len);*/
00776         /*static void voronoi(double *phi,double *theta,double *weight,int lenw,int low,int medium,int nt,int last);*/
00777         static void voronoi(double *phi,double *theta,double *weight, int nt);
00778         static void disorder2(double *x,double *y,int *key,int len);
00779         static void ang_to_xyz(double *x,double *y,double *z,int len);
00780         static void flip23(double *x,double *y,double *z,int *key,int k,int len);
00781         struct tmpstruct{
00782                 double theta1,phi1;
00783                 int key1;
00784                 };
00785         static bool cmp1(tmpstruct tmp1,tmpstruct tmp2);
00786         static bool cmp2(tmpstruct tmp1,tmpstruct tmp2);
00787         /**********************************************************/
00788         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00789         static int trmsh3_(int *n0, double *tol, double *x, double *y, double *z__, int *n, int *list, int *lptr,
00790                int *lend, int *lnew, int *indx, int *lcnt, int *near__, int *next, double *dist, int *ier);
00791         static double areav_(int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier);
00792         /**********************************************************/
00793         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00794 
00795     /*  Various operation on images */
00796         /* out = img + scalar * img1  */
00797         static EMData* madn_scalar(EMData* img, EMData* img1, float scalar);
00798         /* out = scalar * img  */
00799         static EMData* mult_scalar(EMData* img, float scalar);
00800         /* out = img + img1  */
00801         static EMData* addn_img(EMData* img, EMData* img1);
00802         /* out = img - img1  */
00803         static EMData* subn_img(EMData* img, EMData* img1);
00804         /* out = img * img1  */
00805         static EMData* muln_img(EMData* img, EMData* img1);
00806         /* out = img / img1  */
00807         static EMData* divn_img(EMData* img, EMData* img1);
00808         /* img /= Re(img1) with zero check  */
00809         static EMData* divn_filter(EMData* img, EMData* img1);
00810 
00811         /* img += scalar * img1 */
00812         static void mad_scalar(EMData* img, EMData* img1, float scalar);
00813         /* img *= scalar  */
00814         static void mul_scalar(EMData* img, float scalar);
00815         /* img += img1  */
00816         static void add_img(EMData* img, EMData* img1);
00817         /* img += abs(img1)  */
00818         static void add_img_abs(EMData* img, EMData* img1);
00819         /* img += img1**2  */
00820         static void add_img2(EMData* img, EMData* img1);
00821         /* img -= img1  */
00822         static void sub_img(EMData* img, EMData* img1);
00823         /* img *= img1  */
00824         static void mul_img(EMData* img, EMData* img1);
00825         /* img /= img1  */
00826         static void div_img(EMData* img, EMData* img1);
00827         /* img /= Re(img1) with zero check  */
00828         static void div_filter(EMData* img, EMData* img1);
00829         /* pack absolute values of complex image into  real image with addition of Friedel part  */
00830         static EMData* pack_complex_to_real(EMData* img);
00831 private:
00832         static float ang_n(float peakp, string mode, int maxrin); //this function is used by apmq()
00833 public:
00834 
00840         static vector<float> multiref_polar_ali_2d(EMData* image, const vector< EMData* >& crefim,
00841                 float xrng, float yrng, float step, string mode,
00842                 vector< int >numr, float cnx, float cny);
00843 
00844         /* In this version, we return a list of peaks for all reference images */
00845         static vector<float> multiref_polar_ali_2d_peaklist(EMData* image, const vector< EMData* >& crefim,
00846                 float xrng, float yrng, float step, string mode,
00847                 vector< int >numr, float cnx, float cny);
00848 
00849         /* This is used in ISAC program to assigning particles equally to grops */
00850         static vector<int> assign_groups(const vector<float>& d, int nref, int nima);
00851 
00852         static inline void getvec(float phi, float theta, float& x, float& y, float& z, int option=0) {
00853                 float pi180 = M_PI/180.0f;
00854                 
00855                 if (theta > 180.0f) {
00856                         theta -= 180.0f;
00857                         phi += 180.0f;
00858                 } else if (theta > 90.0f && option == 0) {
00859                         theta = 180.0f - theta;
00860                         phi += 180.0f;
00861                 }
00862 
00863                 phi   *= pi180;
00864                 theta *= pi180;
00865 
00866                 x = sin(theta)*cos(phi);
00867                 y = sin(theta)*sin(phi);
00868                 z = cos(theta);
00869 
00870                 return;
00871         }
00872         
00873         static inline float ang_diff(float v11, float v12, float v13, float v21, float v22, float v23, int& mirror) {
00874                 float v = v11*v21+v12*v22+v13*v23;
00875                 if (v > 1) v = 1;
00876                 if (v < -1) v = -1;
00877                 if (v > 0) { mirror = 1; return acos(v)*180.0f/M_PI; }
00878                 else { mirror = -1; return acos(-v)*180.0f/M_PI; }
00879         }
00880         
00881         /* Find nearest projection angles
00882                 Notice: the input I use is different from python code, which I think is awkward.
00883         */
00884         static int nearest_ang(const vector<float>& vecref, float x, float y, float z);
00885 
00886         /* Assign projection angles to nearest reference projections */
00887         static vector<int> assign_projangles(const vector<float>& projangles, const vector<float>& refangles); 
00888 
00889         /* Assign howmany projection angles to the nearest reference projection */
00890         static vector<int> nearestk_to_refdir(const vector<float>& projangles, const vector<float>& refangles, const int howmany); 
00891 
00892         /* Group projection angles by (phi, theta) */
00893         static vector<int> group_proj_by_phitheta(const vector<float>& projangles, const vector<float>& ref_ang, const int img_per_grp);
00894 
00900         static vector<float> multiref_polar_ali_2d_delta(EMData* image, const vector< EMData* >& crefim,
00901                 float xrng, float yrng, float step, string mode,
00902                 vector< int >numr, float cnx, float cny, float delta_start, float delta);
00903 
00909         static vector<float> multiref_polar_ali_2d_nom(EMData* image, const vector< EMData* >& crefim,
00910                 float xrng, float yrng, float step, string mode,
00911                 vector< int >numr, float cnx, float cny);
00912 
00918         static vector<float> multiref_polar_ali_2d_local(EMData* image, const vector< EMData* >& crefim,
00919                 float xrng, float yrng, float step, float ant, string mode,
00920                 vector< int >numr, float cnx, float cny);
00921 
00928         static vector<float> multiref_polar_ali_helical(EMData* image, const vector< EMData* >& crefim,
00929                 float xrng, float yrng, float step, float psi_max, string mode,
00930                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00931         static vector<float> multiref_polar_ali_helical_local(EMData* image, const vector< EMData* >& crefim,
00932                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00933                 vector< int >numr, float cnx, float cny, int ynumber=-1, bool mirror_only=false, float yrnglocal=-1.0, bool CONS=false);
00934         static vector<float> multiref_polar_ali_helical_90(EMData* image, const vector< EMData* >& crefim,
00935                 float xrng, float yrng, float step, float psi_max, string mode,
00936                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00937         static vector<float> multiref_polar_ali_helical_90_local(EMData* image, const vector< EMData* >& crefim,
00938                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00939                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00940 
00946         static vector<float> multiref_polar_ali_2d_local_psi(EMData* image, const vector< EMData* >& crefim,
00947                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00948                 vector< int >numr, float cnx, float cny);
00949 
00957         static void multiref_peaks_ali2d(EMData* image, EMData* crefim,
00958                 float xrng, float yrng, float step, string mode,
00959                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm);
00960 
00968         static void multiref_peaks_compress_ali2d(EMData* image, EMData* crefim, float xrng, float yrng,
00969              float step, string mode, vector<int>numr, float cnx, float cny, EMData *peaks, EMData *peakm,
00970              EMData *peaks_compress, EMData *peakm_compress);
00971 
00976         static vector<float> ali2d_ccf_list(EMData* image, EMData* crefim, float xrng, float yrng,
00977              float step, string mode, vector<int>numr, float cnx, float cny, double T);
00978      /*
00979         static void multiref_peaks_ali(EMData* image, const vector< EMData* >& crefim,
00980                 float xrng, float yrng, float step, string mode,
00981                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm,
00982                     int nphi, int ntheta);
00983 */
00984         static vector<float> twoD_fine_ali(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00985 
00986         static vector<float> twoD_fine_ali_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00987 
00988         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);
00989 
00990         static vector<float> twoD_fine_ali_SD(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00991 
00992         static float ccc_images(EMData *, EMData *, EMData *, float , float , float );
00993 
00994         static vector<float> twoD_fine_ali_SD_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00995 
00996         static float ccc_images_G(EMData* image, EMData* refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sx, float sy);
00997 
00998         static EMData* move_points(EMData* img,  float qprob, int ri, int ro);
00999 
01000         static EMData* get_biggest_cluster( EMData* mg );
01001 
01002         //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);
01003         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);
01004 
01005         static inline int mono(int k1, int k2) {
01006 #ifdef _WIN32
01007                 int  mk = _cpp_max(k1,k2);
01008                 return  _cpp_min(k1,k2) + mk*(mk-1)/2;
01009 #else
01010                 int  mk = std::max(k1,k2);
01011                 return  std::min(k1,k2) + mk*(mk-1)/2;
01012 #endif  //_WIN32
01013         }
01014 
01015         static inline int nint180(float arg) {
01016             int res = int(arg + 180.5) - 180;
01017             return res;
01018         }
01019         
01020         static inline double mean(double *x, int n) {
01021                 double s = 0.0;
01022                 for (int i=0; i<n; i++) s+=x[i];
01023                 return s/static_cast<double>(n);
01024         }
01025 
01026         static inline double var(double *x, int n) {
01027                 double s = 0.0;
01028                 double m = mean(x, n);
01029                 for (int i=0; i<n; i++) s += (x[i]-m)*(x[i]-m);
01030                 return s/static_cast<double>(n);
01031         }
01032         
01033         static vector<float> multi_align_error(vector<float> args, vector<float> all_ali_params, int d);
01034         static double multi_align_error_func(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01035         static vector<double> multi_align_error_func2(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01036         static void multi_align_error_dfunc(double* x, vector<float> all_ali_params, int nima, int num_ali, double* g, int d);
01037         
01038         static vector<float> cluster_pairwise(EMData* d, int K, float T, float F);
01039         //static vector<float> cluster_equalsize(EMData* d, int m);
01040         static vector<float> cluster_equalsize(EMData* d);
01041         static vector<float> vareas(EMData* d);
01042 
01048         static EMData* get_slice(EMData *vol, int dim, int index);
01049 
01050         static void image_mutation(EMData *img, float mutation_rate);
01051 
01053         static void array_mutation(float* list, int len_list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01054 
01055         static vector<float> list_mutation(vector<float> list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01056         /*
01057                         To restrict the value to [0, nx)
01058         */
01059         static inline float restrict1(float x, int nx) {
01060                 while ( x < 0.0f )        x += nx;
01061                 while ( x >= (float)(nx) )  x -= nx;
01062                 return x;
01063         }
01064 
01065 #endif  //util__sparx_h__

Generated on Fri Aug 10 16:35:18 2012 for EMAN2 by  doxygen 1.3.9.1