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

Generated on Thu May 3 10:06:29 2012 for EMAN2 by  doxygen 1.4.7