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_ns(EMData* circ1, EMData* circ2, vector<int> numr);
00568 
00575         static EMData* Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr);
00576 
00583         static void Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t);
00584 
00591         static EMData* Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr);
00592 
00599         static EMData* Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr);
00600 
00601         static vector<float> Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr );
00602         static void  prb1d(double *b, int npoint, float *pos);
00603 
00604         static void update_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00605         static void sub_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00606 
00607         // helper functions for ali2d_ra
00608         static float ener(EMData* ave, vector<int> numr);
00609 
00610         static float ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot);
00611 
00613         static Dict min_dist_real(EMData* image, const vector<EMData*>& data);
00614 
00616         static Dict min_dist_four(EMData* image, const vector<EMData*>& data);
00617 
00624         static int k_means_cont_table_(int* group1, int* group2, int* stb, long int s1, long int s2, int flag);
00625 
00626         // branch and bound matching algorithm
00627 
00628         
00632         static void initial_prune(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T);
00633 
00637         static bool explore(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T, int partref, int* curintx, int
00638                         size_curintx, int* next, int size_next, int depth);
00639 
00640         
00641         
00642 
00647         static int generatesubmax(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int n_guesses, int LARGEST_CLASS);
00648 
00652         static void search2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* matchlist, int* costlist, int J);
00653         
00654         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*
00655 costlist, int* curbranch);
00656         
00661         static bool sanitycheck(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* output);
00662 
00672         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,
00673         int branchfunc, int LIM);
00674 
00675         
00681         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,
00682         float stmult, int branchfunc, int LIM);
00683         
00684         static int branch_factor_0(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00685         static int branch_factor_2(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int LIM);
00686         static int branch_factor_3(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, int K, int LIM);
00687         static int branch_factor_4(int* costlist, int* matchlist, int J, int T, int nParts, int curlevel, int max_branching, float stmult);
00688         // new code common-lines
00689 
00690         //helper function for the weights calculation by Voronoi to Cml
00691         static vector<double> cml_weights(const vector<float>& cml);
00692 
00694         static vector<int> cml_line_insino(vector<float> Rot, int i_prj, int n_prj);
00695 
00697         static vector<int> cml_line_insino_all(vector<float> Rot, vector<int> seq, int n_prj, int n_lines);
00698 
00700         static vector<double> cml_init_rot(vector<float> Ori);
00701 
00703         static vector<float> cml_update_rot(vector<float> Rot, int iprj, float nph, float th, float nps);
00704 
00706         static vector<double> cml_line_in3d(vector<float> Ori, vector<int> seq, int nprj, int nlines);
00707 
00709         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);
00710         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);
00711 
00713         static double cml_disc(const vector<EMData*>& data, vector<int> com, vector<int> seq, vector<float> weights, int n_lines);
00714 
00720         static void set_line(EMData* img, int posline, EMData* line, int offset, int length);
00721 
00729         static void cml_prepare_line(EMData* sino, EMData* line, int ilf, int ihf, int pos_line, int nblines);
00730 
00731         /* Decimates the image with respect to the image center.
00732          * (i.e) the center of the original image is kept the same
00733          * and then the initial start pixel is calculated with respect to the
00734          * center of the image
00735          * @params(image, x-pixel, y-pixel,z-pixel)
00736          * works for all 3 dimensions
00737         **/
00738         static EMData* decimate(EMData* img, int x_step,int y_step=1,int z_step=1);
00739 
00740         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);
00741 
00742         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");
00743 
00744         static vector<float> histogram(EMData* image, EMData* mask, int nbins = 128, float hmin =0.0f, float hmax = 0.0f );
00745 
00746         static Dict histc(EMData *ref,EMData *img,EMData *mask);
00747 
00748         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);
00749 
00750 
00751         /* 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
00752          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
00753           * sign can be set as +1 or -1 . The unit of frequency ak is 1/Angstrom
00754                   Attention: Envelope function in power spectrum has a form of exp(-b_factor*ak^2)
00755                                           */
00756         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);
00757         static EMData *compress_image_mask(EMData* image, EMData* mask);
00758 
00760         static EMData *reconstitute_image_mask(EMData *image,EMData *mask);
00761 
00762         static vector<float> merge_peaks(vector<float> peak1, vector<float> peak2,float p_size);
00763         static vector<float> pw_extract(vector<float>pw, int n, int iswi,float ps);
00764         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);
00765         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);
00766         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
00767         int *iu, double *s);
00768         static float eval(char * images,EMData * img, vector<int> S,int N, int K,int size);
00769 
00770         /*  VORONOI DIAGRAM */
00771         static vector<double> vrdg(const vector<float>& ph, const vector<float>& th);
00772         static void hsortd(double *theta,double *phi,int *key,int len,int option);
00773         static void voronoidiag(double *theta,double *phi,double* weight,int n);
00774         /*static void angstep(double* thetast,int len);*/
00775         /*static void voronoi(double *phi,double *theta,double *weight,int lenw,int low,int medium,int nt,int last);*/
00776         static void voronoi(double *phi,double *theta,double *weight, int nt);
00777         static void disorder2(double *x,double *y,int *key,int len);
00778         static void ang_to_xyz(double *x,double *y,double *z,int len);
00779         static void flip23(double *x,double *y,double *z,int *key,int k,int len);
00780         struct tmpstruct{
00781                 double theta1,phi1;
00782                 int key1;
00783                 };
00784         static bool cmp1(tmpstruct tmp1,tmpstruct tmp2);
00785         static bool cmp2(tmpstruct tmp1,tmpstruct tmp2);
00786         /**********************************************************/
00787         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00788         static int trmsh3_(int *n0, double *tol, double *x, double *y, double *z__, int *n, int *list, int *lptr,
00789                int *lend, int *lnew, int *indx, int *lcnt, int *near__, int *next, double *dist, int *ier);
00790         static double areav_(int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier);
00791         /**********************************************************/
00792         /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
00793 
00794     /*  Various operation on images */
00795         /* out = img + scalar * img1  */
00796         static EMData* madn_scalar(EMData* img, EMData* img1, float scalar);
00797         /* out = scalar * img  */
00798         static EMData* mult_scalar(EMData* img, float scalar);
00799         /* out = img + img1  */
00800         static EMData* addn_img(EMData* img, EMData* img1);
00801         /* out = img - img1  */
00802         static EMData* subn_img(EMData* img, EMData* img1);
00803         /* out = img * img1  */
00804         static EMData* muln_img(EMData* img, EMData* img1);
00805         /* out = img / img1  */
00806         static EMData* divn_img(EMData* img, EMData* img1);
00807         /* img /= Re(img1) with zero check  */
00808         static EMData* divn_filter(EMData* img, EMData* img1);
00809 
00810         /* img += scalar * img1 */
00811         static void mad_scalar(EMData* img, EMData* img1, float scalar);
00812         /* img *= scalar  */
00813         static void mul_scalar(EMData* img, float scalar);
00814         /* img += img1  */
00815         static void add_img(EMData* img, EMData* img1);
00816         /* img += abs(img1)  */
00817         static void add_img_abs(EMData* img, EMData* img1);
00818         /* img += img1**2  */
00819         static void add_img2(EMData* img, EMData* img1);
00820         /* img -= img1  */
00821         static void sub_img(EMData* img, EMData* img1);
00822         /* img *= img1  */
00823         static void mul_img(EMData* img, EMData* img1);
00824         /* img /= img1  */
00825         static void div_img(EMData* img, EMData* img1);
00826         /* img /= Re(img1) with zero check  */
00827         static void div_filter(EMData* img, EMData* img1);
00828         /* pack absolute values of complex image into  real image with addition of Friedel part  */
00829         static EMData* pack_complex_to_real(EMData* img);
00830 private:
00831         static float ang_n(float peakp, string mode, int maxrin); //this function is used by apmq()
00832 public:
00833 
00839         static vector<float> multiref_polar_ali_2d(EMData* image, const vector< EMData* >& crefim,
00840                 float xrng, float yrng, float step, string mode,
00841                 vector< int >numr, float cnx, float cny);
00842 
00843         /* In this version, we return a list of peaks for all reference images */
00844         static vector<float> multiref_polar_ali_2d_peaklist(EMData* image, const vector< EMData* >& crefim,
00845                 float xrng, float yrng, float step, string mode,
00846                 vector< int >numr, float cnx, float cny);
00847 
00848         /* This is used in ISAC program to assigning particles equally to grops */
00849         static vector<int> assign_groups(std::string matrix_address, int nref, int nima);
00850 
00851         static inline void getvec(float phi, float theta, float& x, float& y, float& z, int option=0) {
00852                 float pi180 = M_PI/180.0f;
00853                 
00854                 if (theta > 180.0f) {
00855                         theta -= 180.0f;
00856                         phi += 180.0f;
00857                 } else if (theta > 90.0f && option == 0) {
00858                         theta = 180.0f - theta;
00859                         phi += 180.0f;
00860                 }
00861 
00862                 phi   *= pi180;
00863                 theta *= pi180;
00864 
00865                 x = sin(theta)*cos(phi);
00866                 y = sin(theta)*sin(phi);
00867                 z = cos(theta);
00868 
00869                 return;
00870         }
00871         
00872         static inline float ang_diff(float v11, float v12, float v13, float v21, float v22, float v23, int& mirror) {
00873                 float v = v11*v21+v12*v22+v13*v23;
00874                 if (v > 1) v = 1;
00875                 if (v < -1) v = -1;
00876                 if (v > 0) { mirror = 1; return acos(v)*180.0f/M_PI; }
00877                 else { mirror = -1; return acos(-v)*180.0f/M_PI; }
00878         }
00879         
00880         /* Find nearest projection angles
00881                 Notice: the input I use is different from python code, which I think is awkward.
00882         */
00883         static int nearest_ang(const vector<float>& vecref, float x, float y, float z);
00884 
00885         /* Assign projection angles to nearest reference projections */
00886         static vector<int> assign_projangles(const vector<float>& projangles, const vector<float>& refangles); 
00887 
00888         /* Assign howmany projection angles to the nearest reference projection */
00889         static vector<int> nearestk_to_refdir(const vector<float>& projangles, const vector<float>& refangles, const int howmany); 
00890 
00891         /* Group projection angles by (phi, theta) */
00892         static vector<int> group_proj_by_phitheta(const vector<float>& projangles, const vector<float>& ref_ang, const int img_per_grp);
00893 
00899         static vector<float> multiref_polar_ali_2d_delta(EMData* image, const vector< EMData* >& crefim,
00900                 float xrng, float yrng, float step, string mode,
00901                 vector< int >numr, float cnx, float cny, float delta_start, float delta);
00902 
00908         static vector<float> multiref_polar_ali_2d_nom(EMData* image, const vector< EMData* >& crefim,
00909                 float xrng, float yrng, float step, string mode,
00910                 vector< int >numr, float cnx, float cny);
00911 
00917         static vector<float> multiref_polar_ali_2d_local(EMData* image, const vector< EMData* >& crefim,
00918                 float xrng, float yrng, float step, float ant, string mode,
00919                 vector< int >numr, float cnx, float cny);
00920         static vector<float> hans(EMData* image, const vector< EMData* >& crefim,
00921                 float xrng, float yrng, float step, float ant, string mode,
00922                 vector< int >numr, float cnx, float cny);
00923 
00930         static vector<float> multiref_polar_ali_helical(EMData* image, const vector< EMData* >& crefim,
00931                 float xrng, float yrng, float step, float psi_max, string mode,
00932                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00933         static vector<float> multiref_polar_ali_helical_local(EMData* image, const vector< EMData* >& crefim,
00934                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00935                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00936         static vector<float> multiref_polar_ali_helical_90(EMData* image, const vector< EMData* >& crefim,
00937                 float xrng, float yrng, float step, float psi_max, string mode,
00938                 vector< int >numr, float cnx, float cny, int ynumber=-1);
00939         static vector<float> multiref_polar_ali_helical_90_local(EMData* image, const vector< EMData* >& crefim,
00940                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00941                 vector< int >numr, float cnx, float cny, int ynumber=-1, float yrnglocal=-1.0);
00942 
00948         static vector<float> multiref_polar_ali_2d_local_psi(EMData* image, const vector< EMData* >& crefim,
00949                 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00950                 vector< int >numr, float cnx, float cny);
00958         static void multiref_peaks_ali2d(EMData* image, EMData* crefim,
00959                 float xrng, float yrng, float step, string mode,
00960                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm);
00961 
00969         static void multiref_peaks_compress_ali2d(EMData* image, EMData* crefim, float xrng, float yrng,
00970              float step, string mode, vector<int>numr, float cnx, float cny, EMData *peaks, EMData *peakm,
00971              EMData *peaks_compress, EMData *peakm_compress);
00972 
00977         static vector<float> ali2d_ccf_list(EMData* image, EMData* crefim, float xrng, float yrng,
00978              float step, string mode, vector<int>numr, float cnx, float cny, double T);
00979      /*
00980         static void multiref_peaks_ali(EMData* image, const vector< EMData* >& crefim,
00981                 float xrng, float yrng, float step, string mode,
00982                 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm,
00983                     int nphi, int ntheta);
00984 */
00985         static vector<float> twoD_fine_ali(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00986 
00987         static vector<float> twoD_fine_ali_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00988 
00989         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);
00990 
00991         static vector<float> twoD_fine_ali_SD(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00992 
00993         static float ccc_images(EMData *, EMData *, EMData *, float , float , float );
00994 
00995         static vector<float> twoD_fine_ali_SD_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00996 
00997         static float ccc_images_G(EMData* image, EMData* refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sx, float sy);
00998 
00999         static EMData* move_points(EMData* img,  float qprob, int ri, int ro);
01000 
01001         static EMData* get_biggest_cluster( EMData* mg );
01002 
01003         //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);
01004         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);
01005         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);
01006         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);
01007 
01008         static inline int mono(int k1, int k2) {
01009 #ifdef _WIN32
01010                 int  mk = _cpp_max(k1,k2);
01011                 return  _cpp_min(k1,k2) + mk*(mk-1)/2;
01012 #else
01013                 int  mk = std::max(k1,k2);
01014                 return  std::min(k1,k2) + mk*(mk-1)/2;
01015 #endif  //_WIN32
01016         }
01017 
01018         static inline int nint180(float arg) {
01019             int res = int(arg + 180.5) - 180;
01020             return res;
01021         }
01022         
01023         static inline double mean(double *x, int n) {
01024                 double s = 0.0;
01025                 for (int i=0; i<n; i++) s+=x[i];
01026                 return s/static_cast<double>(n);
01027         }
01028 
01029         static inline double var(double *x, int n) {
01030                 double s = 0.0;
01031                 double m = mean(x, n);
01032                 for (int i=0; i<n; i++) s += (x[i]-m)*(x[i]-m);
01033                 return s/static_cast<double>(n);
01034         }
01035         
01036         static vector<float> multi_align_error(vector<float> args, vector<float> all_ali_params, int d);
01037         static double multi_align_error_func(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01038         static vector<double> multi_align_error_func2(double* x, vector<float> all_ali_params, int nima, int num_ali, int d);
01039         static void multi_align_error_dfunc(double* x, vector<float> all_ali_params, int nima, int num_ali, double* g, int d);
01040         
01041         static vector<float> cluster_pairwise(EMData* d, int K, float T, float F);
01042         //static vector<float> cluster_equalsize(EMData* d, int m);
01043         static vector<float> cluster_equalsize(EMData* d);
01044         static vector<float> vareas(EMData* d);
01045 
01051         static EMData* get_slice(EMData *vol, int dim, int index);
01052 
01053         static void image_mutation(EMData *img, float mutation_rate);
01054 
01056         static void array_mutation(float* list, int len_list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01057 
01058         static vector<float> list_mutation(vector<float> list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
01059         /*
01060                         To restrict the value to [0, nx)
01061         */
01062         static inline float restrict1(float x, int nx) {
01063                 while ( x < 0.0f )        x += nx;
01064                 while ( x >= (float)(nx) )  x -= nx;
01065                 return x;
01066         }
01067 
01070         static Dict get_transform_params(EMData* image, string xform, string convention);
01071 
01072         static void constrained_helix( vector<EMData*> data, vector<EMData*> fdata, vector<EMData*> refproj, vector<EMData*> rotproj
01073                         , vector<float> dp_dphi_rise_delta, vector<int> nphi_phiwobble_range_ywobble_Dsym_nwx_nwy_nwxc_nwyc
01074                         , bool FindPsi, float psi_max, vector<EMData*> crefim, vector<int> numr, int maxrin, string mode, int cnx, int cny);
01075         static void constrained_helix_test( vector<EMData*> data, vector<EMData*> fdata, vector<EMData*> refproj, vector<EMData*> rotproj
01076                         , vector<float> dp_dphi_rise_delta, vector<int> nphi_phiwobble_range_ywobble_Dsym_nwx_nwy_nwxc_nwyc
01077                         , bool FindPsi, float psi_max, vector<EMData*> crefim, vector<int> numr, int maxrin, string mode, int cnx, int cny);
01078         static Dict predict(float phig, float yg, float dst, float sgn, float ysgn, float dpp, float dphi, bool backpred);
01079         
01080 #endif  //util__sparx_h__

Generated on Tue Jun 11 13:40:45 2013 for EMAN2 by  doxygen 1.3.9.1