util_sparx.h

Go to the documentation of this file.

/*
 * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
 * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

#ifndef util__sparx_h__
#define util__sparx_h__

public:

static int coveig(int n, float *covmat, float *eigval, float *eigvec);

static Dict coveig_for_py(int ncov, const vector<float>& covmatpy);

static Dict ExpMinus4YSqr(float ymax,int nsamples);

static void WTM(EMData* PROJ, vector<float> SS,int DIAMETER,int NUMP);

static void WTF(EMData* PROJ,vector<float> SS,float SNR,int K,vector<float> exptable);

static Dict CANG(float PHI, float THETA, float PSI);

static void BPCQ(EMData* B, EMData *CUBE,vector<float> DM);

static vector<float> infomask(EMData* Vol, EMData* mask, bool);

static void colreverse(float* beg, float* end, int nx);

static void slicereverse(float* beg, float* end, int nx,int ny);

static void cyclicshift(EMData* image, Dict params);

static Dict im_diff(EMData* V1, EMData* V2, EMData* mask=0);

static EMData* TwoDTestFunc(int Size, float p, float q,  float a, float b,
                   int flag=0, float alphaDeg=0); //PRB


static void spline_mat(float *x, float *y, int n,  float *xq, float *yq, int m); //PRB

static void spline(float *x, float *y, int n, float yp1, float ypn, float *y2);
  // PRB
static void splint( float *xa, float *ya, float *y2a, int n,
                                     float *xq, float *yq, int m);


static void Radialize(int *PermMatTr,  float * kValsSorted,
            float *weightofkvalsSorted, int Size, int *SizeReturned);



class sincBlackman
{
        protected:
                int M; 
                float fc; 
                int ntable;
                vector<float> sBtable;
                virtual void build_sBtable(); 
                float fltb;
        public:
                sincBlackman(int M_, float fc_, int ntable_ = 1999);
                virtual ~sincBlackman() {};

                inline  float sBwin_tab(float x) const {
                        float xt;
                        if(x<0.0f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
                        return sBtable[ (int) xt];
                }
                int get_sB_size() const { return M; }
};



class KaiserBessel
{
        protected:
                float alpha, v, r; 
                int N; 
                int K; 
                float vtable; 
                int ntable;
                vector<float> i0table;
                float dtable; 
                float alphar; 
                float fac; 
                float vadjust;
                float facadj; 
                virtual void build_I0table(); 
                float fltb;
        public:
                KaiserBessel(float alpha_, int K, float r_,
                                     float v_, int N_, float vtable_=0.f,
                                         int ntable_ = 5999);
                virtual ~KaiserBessel() {};
                float I0table_maxerror();
                vector<float> dump_table() {
                        return i0table;
                }
                virtual float sinhwin(float x) const;
                virtual float i0win(float x) const;
                inline float i0win_tab(float x) const {
                /*float absx = fabs(x);
                        int loc = int(round(absx*fltb));
                        return i0table[loc];*/
                        float xt;
                        if(x<0.f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
                        return i0table[ (int) xt];
                        /*return i0table[ (int) (fabs(x)*fltb+0.5f)];
                                if (absx > vtable) return 0.f;
                                float loc = absx/dtable;
                                return i0table[int(loc + 0.5f)]; */
                }
                int get_window_size() const { return K; }
                class kbsinh_win {
                        KaiserBessel& kb;
                        public:
                        kbsinh_win(KaiserBessel& kb_) : kb(kb_) {}
                        float operator()(float x) const {
                                return kb.sinhwin(x);
                        }
                        int get_window_size() const {return kb.get_window_size();}
                };
                kbsinh_win get_kbsinh_win() {
                        return kbsinh_win(*this);
                }
                class kbi0_win {
                        KaiserBessel& kb;
                        public:
                        kbi0_win(KaiserBessel& kb_) : kb(kb_) {}
                        float operator()(float x) const {
                                return kb.i0win(x);
                        }
                        int get_window_size() const {return kb.get_window_size();}
                };
                kbi0_win get_kbi0_win() {
                        return kbi0_win(*this);
                }
};

class FakeKaiserBessel : public KaiserBessel {
        public:
                FakeKaiserBessel(float alpha, int K, float r_,
                                         float v_, int N_, float vtable_=0.f,
                                                 int ntable_ = 5999)
        : KaiserBessel(alpha, K, r_, v_, N_, vtable_, ntable_) {
                        build_I0table();
                }
                float sinhwin(float x) const;
                float i0win(float x) const;
                void build_I0table();
};

                static vector<float>
                even_angles(float delta, float t1=0, float t2=90, float p1=0, float p2=359.999);


                static float quadri(float x, float y, int nx, int ny, float* image);

                static float quadri_background(float x, float y, int nx, int ny, float* image, int xnew, int ynew);

                static float get_pixel_conv_new(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb);

        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);

                static std::complex<float> extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb);

                /*static float quadris(float x, float y, int nx, int ny, float* image);*/
                static float bilinear(float xold, float yold, int nsam, int nrow, float* xim);


                static float triquad(float r, float s, float t, float* fdata);

                class Gaussian {
                        float sigma;
                        float rttwopisigma;
                        float twosigma2;
                        public:
                        Gaussian(float sigma_ = 1.0) : sigma(sigma_) {
                                rttwopisigma = sqrtf(static_cast<float>(twopi)*sigma);
                                twosigma2 = 2*sigma*sigma;
                        }
                        inline float operator()(float x) const {
                                return exp(-x*x/(twosigma2))/rttwopisigma;
                        }
                };
        /*static void alrq(float *xim,  int nsam , int nrow , int *numr,
                             float *circ, int lcirc, int nring, char mode);*/
        static EMData* Polar2D(EMData* image, vector<int> numr, string mode);
        static EMData* Polar2Dm(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode);
        /*static void alrq_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
                            int  *numr, float *circ, int lcirc, int  nring, char  mode);*/
        static void alrl_ms(float *xim, int    nsam, int  nrow, float cns2, float cnr2,
                            int  *numr, float *circ, int lcirc, int  nring, char  mode);
        /*static void alrq_msi(EMData* image,float cns2, float cnr2,
                           int  *numr, float *circ, int lcirc, int  nring, char  mode, Util::KaiserBessel&
                                               kb);*/
        static EMData* Polar2Dmi(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode, Util::KaiserBessel& kb);

        static void  fftr_q(float  *xcmplx, int nv);
        static void  fftr_d(double *xcmplx, int nv);
        static void  fftc_q(float  *br, float  *bi, int ln, int ks);
        static void  fftc_d(double *br, double *bi, int ln, int ks);
        static void  Frngs(EMData* circ, vector<int> numr);
        static void  Normalize_ring(EMData* ring, const vector<int>& numr);
        static void  Frngs_inv(EMData* circ, vector<int> numr);

        /*
                A little notes about different Crosrng:
                Basically, they all do cross-correlation function to two images in polar coordinates
                Crosrng_e is the original one
                Crosrng_ew is the one that you could apply weights to different rings
                Crosrng_ms assumes the user already apply weights to circ1, it also returns both
                           straight and mirrored positions simultaneously.
                Crosrng_msg differs from the previous ones in that it returns the cross-correlation
                            function entirely instead of the peak value and position, thus makes it
                            possible to use the gridding method to determine the peak position
                Crosrng_msg_s is same as Crosrng_msg except that it only checks straight position
                Crosrng_msg_m is same as Crosrng_msg except that it only checks mirrored position
          */
        static Dict Crosrng_e(EMData* circ1, EMData* circ2, vector<int> numr, int neg);
        static Dict Crosrng_ew(EMData* circ1, EMData* circ2, vector<int> numr, vector<float> w, int neg);

        static Dict Crosrng_ms(EMData* circ1, EMData* circ2, vector<int> numr);
        static Dict Crosrng_sm_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, int flag);
        static Dict Crosrng_ns(EMData* circ1, EMData* circ2, vector<int> numr);

        static EMData* Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr);
        static void Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t);
        static EMData* Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr);
        static EMData* Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr);

        static vector<float> Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr );
        static void  prb1d(double *b, int npoint, float *pos);

        static void update_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
        static void sub_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);

        static float ener(EMData* ave, vector<int> numr);

        static float ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot);

        // k-means helper
        static Dict min_dist_real(EMData* image, const vector<EMData*>& data);
        static Dict min_dist_four(EMData* image, const vector<EMData*>& data);
        static int k_means_cont_table_(int* grp1, int* grp2, int* stb, long int s1, long int s2, int flag);

        // new code common-lines
        static vector<double> cml_weights(const vector<float>& cml);
        static vector<int> cml_line_insino(vector<float> Rot, int i_prj, int n_prj);
        static vector<int> cml_line_insino_all(vector<float> Rot, vector<int> seq, int n_prj, int n_lines);
        static vector<double> cml_init_rot(vector<float> Ori);
        static vector<float> cml_update_rot(vector<float> Rot, int iprj, float nph, float th, float nps);
        static vector<double> cml_line_in3d(vector<float> Ori, vector<int> seq, int nprj, int nlines);
        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);
        static double cml_disc(const vector<EMData*>& data, vector<int> com, vector<int> seq, vector<float> weights, int n_lines);
        static void set_line(EMData* img, int posline, EMData* line, int offset, int length);
        static void cml_prepare_line(EMData* sino, EMData* line, int ilf, int ihf, int pos_line, int nblines);
        /* Decimates the image with respect to the image center.
         * (i.e) the center of the original image is kept the same
         * and then the initial start pixel is calculated with respect to the
         * center of the image
         * @params(image, x-pixel, y-pixel,z-pixel)
         * works for all 3 dimensions
        **/
        static EMData* decimate(EMData* img, int x_step,int y_step=1,int z_step=1);

        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);

        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, char *params="average");

        static vector<float> histogram(EMData* image, EMData* mask, int nbins = 128, float hmin =0.0f, float hmax = 0.0f );

        static Dict histc(EMData *ref,EMData *img,EMData *mask);

        static float hist_comp_freq(float PA,float PB,int size_img, int hist_len, EMData *img, vector<float> ref_freq_hist, EMData *mask, float ref_h_diff, float ref_h_min);


        /* 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
         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
          * sign can be set as +1 or -1 . The unit of frequency ak is 1/Angstrom
                  Attention: Envelope function in power spectrum has a form of exp(-b_factor*ak^2)
                                          */
        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);
        static EMData *compress_image_mask(EMData* image, EMData* mask);
        static EMData *reconstitute_image_mask(EMData *image,EMData *mask);
        static vector<float> merge_peaks(vector<float> peak1, vector<float> peak2,float p_size);
        static vector<float> pw_extract(vector<float>pw, int n, int iswi,float ps);
        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);
        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);
        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
        int *iu, double *s);
        static float eval(char * images,EMData * img, vector<int> S,int N, int K,int size);

        /*  VORONOI DIAGRAM */
        static vector<double> vrdg(const vector<float>& ph, const vector<float>& th);
        static void hsortd(double *theta,double *phi,int *key,int len,int option);
        static void voronoidiag(double *theta,double *phi,double* weight,int n);
        /*static void angstep(double* thetast,int len);*/
        /*static void voronoi(double *phi,double *theta,double *weight,int lenw,int low,int medium,int nt,int last);*/
        static void voronoi(double *phi,double *theta,double *weight, int nt);
        static void disorder2(double *x,double *y,int *key,int len);
        static void ang_to_xyz(double *x,double *y,double *z,int len);
        static void flip23(double *x,double *y,double *z,int *key,int k,int len);
        struct tmpstruct{
                double theta1,phi1;
                int key1;
                };
        static bool cmp1(tmpstruct tmp1,tmpstruct tmp2);
        static bool cmp2(tmpstruct tmp1,tmpstruct tmp2);
        /**********************************************************/
        /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/
        static int trmsh3_(int *n0, double *tol, double *x, double *y, double *z__, int *n, int *list, int *lptr,
               int *lend, int *lnew, int *indx, int *lcnt, int *near__, int *next, double *dist, int *ier);
        static double areav_(int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier);
        /**********************************************************/
        /* ######### STRIDPACK USED COMMANDS FOR VORONOI #########################*/

    /*  Various operation on images */
        /* out = img + scalar * img1  */
        static EMData* madn_scalar(EMData* img, EMData* img1, float scalar);
        /* out = scalar * img  */
        static EMData* mult_scalar(EMData* img, float scalar);
        /* out = img + img1  */
        static EMData* addn_img(EMData* img, EMData* img1);
        /* out = img - img1  */
        static EMData* subn_img(EMData* img, EMData* img1);
        /* out = img * img1  */
        static EMData* muln_img(EMData* img, EMData* img1);
        /* out = img / img1  */
        static EMData* divn_img(EMData* img, EMData* img1);
        /* img /= Re(img1) with zero check  */
        static EMData* divn_filter(EMData* img, EMData* img1);

        /* img += scalar * img1 */
        static void mad_scalar(EMData* img, EMData* img1, float scalar);
        /* img *= scalar  */
        static void mul_scalar(EMData* img, float scalar);
        /* img += img1  */
        static void add_img(EMData* img, EMData* img1);
        /* img += img1**2  */
        static void add_img2(EMData* img, EMData* img1);
        /* img -= img1  */
        static void sub_img(EMData* img, EMData* img1);
        /* img *= img1  */
        static void mul_img(EMData* img, EMData* img1);
        /* img /= img1  */
        static void div_img(EMData* img, EMData* img1);
        /* img /= Re(img1) with zero check  */
        static void div_filter(EMData* img, EMData* img1);
        /* pack absolute values of complex image into  real image with addition of Friedel part  */
        static EMData* pack_complex_to_real(EMData* img);
private:
        static float ang_n(float peakp, string mode, int maxrin); //this function is used by apmq()
public:
        static vector<float> multiref_polar_ali_2d(EMData* image, const vector< EMData* >& crefim,
                float xrng, float yrng, float step, string mode,
                vector< int >numr, float cnx, float cny);
        static vector<float> multiref_polar_ali_2d_nom(EMData* image, const vector< EMData* >& crefim,
                float xrng, float yrng, float step, string mode,
                vector< int >numr, float cnx, float cny);
        static vector<float> multiref_polar_ali_2d_local(EMData* image, const vector< EMData* >& crefim,
                float xrng, float yrng, float step, float ant, string mode,
                vector< int >numr, float cnx, float cny);
        static vector<float> multiref_polar_ali_2d_local_psi(EMData* image, const vector< EMData* >& crefim,
                float xrng, float yrng, float step, float ant, float psi_max, string mode,
                vector< int >numr, float cnx, float cny);

        static void multiref_peaks_ali2d(EMData* image, EMData* crefim,
                float xrng, float yrng, float step, string mode,
                vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm);

        static void multiref_peaks_compress_ali2d(EMData* image, EMData* crefim, float xrng, float yrng,
             float step, string mode, vector<int>numr, float cnx, float cny, EMData *peaks, EMData *peakm,
             EMData *peaks_compress, EMData *peakm_compress);

        static vector<float> ali2d_ccf_list(EMData* image, EMData* crefim, float xrng, float yrng,
             float step, string mode, vector<int>numr, float cnx, float cny, double T);
     /*
        static void multiref_peaks_ali(EMData* image, const vector< EMData* >& crefim,
                float xrng, float yrng, float step, string mode,
                vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm,
                    int nphi, int ntheta);
*/
        static vector<float> twoD_fine_ali(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);

        static vector<float> twoD_fine_ali_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);

        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);

        static vector<float> twoD_fine_ali_SD(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);

        static float ccc_images(EMData *, EMData *, EMData *, float , float , float );

        static vector<float> twoD_fine_ali_SD_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);

        static float ccc_images_G(EMData* image, EMData* refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sx, float sy);

        static EMData* move_points(EMData* img,  float qprob, int ri, int ro);

        static EMData* get_biggest_cluster( EMData* mg );

        //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);
        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);

        static inline int mono(int k1, int k2) {
#ifdef _WIN32
                int  mk = _cpp_max(k1,k2);
                return  _cpp_min(k1,k2) + mk*(mk-1)/2;
#else
                int  mk = std::max(k1,k2);
                return  std::min(k1,k2) + mk*(mk-1)/2;
#endif  //_WIN32
        }

        static inline int nint180(float arg) {
            int res = int(arg + 180.5) - 180;
            return res;
        }

        static vector<float> cluster_pairwise(EMData* d, int K, float T, float F);
        //static vector<float> cluster_equalsize(EMData* d, int m);
        static vector<float> cluster_equalsize(EMData* d);
        static vector<float> vareas(EMData* d);
        static EMData* get_slice(EMData *vol, int dim, int index);
        static void image_mutation(EMData *img, float mutation_rate);
        /*
                        To restrict the value to [0, nx)
        */
        static inline float restrict1(float x, int nx) {
                while ( x < 0.0f )        x += nx;
                while ( x >= (float)(nx) )  x -= nx;
                return x;
        }

#endif  //util__sparx_h__