00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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 Dict ExpMinus4YSqr(float ymax,int nsamples);
00051
00052 static void WTM(EMData* PROJ, vector<float> SS,int DIAMETER,int NUMP);
00053
00054 static void WTF(EMData* PROJ,vector<float> SS,float SNR,int K,vector<float> exptable);
00055
00056 static Dict CANG(float PHI, float THETA, float PSI);
00057
00058 static void BPCQ(EMData* B, EMData *CUBE,vector<float> DM);
00059
00060 static vector<float> infomask(EMData* Vol, EMData* mask, bool);
00061
00062 static void colreverse(float* beg, float* end, int nx);
00063
00064 static void slicereverse(float* beg, float* end, int nx,int ny);
00065
00100 static void cyclicshift(EMData* image, Dict params);
00101
00102 static Dict im_diff(EMData* V1, EMData* V2, EMData* mask=0);
00103
00116 static EMData* TwoDTestFunc(int Size, float p, float q, float a, float b,
00117 int flag=0, float alphaDeg=0);
00118
00119
00131 static void spline_mat(float *x, float *y, int n, float *xq, float *yq, int m);
00132
00145 static void spline(float *x, float *y, int n, float yp1, float ypn, float *y2);
00146
00158 static void splint( float *xa, float *ya, float *y2a, int n,
00159 float *xq, float *yq, int m);
00160
00161
00172 static void Radialize(int *PermMatTr, float * kValsSorted,
00173 float *weightofkvalsSorted, int Size, int *SizeReturned);
00174
00175
00176
00177 class sincBlackman
00178 {
00179 protected:
00180 int M;
00181 float fc;
00182 int ntable;
00183 vector<float> sBtable;
00184 virtual void build_sBtable();
00185 float fltb;
00186 public:
00187 sincBlackman(int M_, float fc_, int ntable_ = 1999);
00188 virtual ~sincBlackman() {};
00189
00190 inline float sBwin_tab(float x) const {
00191 float xt;
00192 if(x<0.0f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00193 return sBtable[ (int) xt];
00194 }
00196 int get_sB_size() const { return M; }
00197 };
00198
00199
00200
00224 class KaiserBessel
00225 {
00226 protected:
00227 float alpha, v, r;
00228 int N;
00229 int K;
00230 float vtable;
00231 int ntable;
00232 vector<float> i0table;
00233 float dtable;
00234 float alphar;
00235 float fac;
00236 float vadjust;
00237 float facadj;
00238 virtual void build_I0table();
00239 float fltb;
00240 public:
00241 KaiserBessel(float alpha_, int K, float r_,
00242 float v_, int N_, float vtable_=0.f,
00243 int ntable_ = 5999);
00244 virtual ~KaiserBessel() {};
00246 float I0table_maxerror();
00247 vector<float> dump_table() {
00248 return i0table;
00249 }
00251 virtual float sinhwin(float x) const;
00253 virtual float i0win(float x) const;
00255 inline float i0win_tab(float x) const {
00256
00257
00258
00259 float xt;
00260 if(x<0.f) xt = -x*fltb+0.5f; else xt = x*fltb+0.5f;
00261 return i0table[ (int) xt];
00262
00263
00264
00265
00266 }
00268 int get_window_size() const { return K; }
00270 class kbsinh_win {
00271 KaiserBessel& kb;
00272 public:
00273 kbsinh_win(KaiserBessel& kb_) : kb(kb_) {}
00274 float operator()(float x) const {
00275 return kb.sinhwin(x);
00276 }
00277 int get_window_size() const {return kb.get_window_size();}
00278 };
00280 kbsinh_win get_kbsinh_win() {
00281 return kbsinh_win(*this);
00282 }
00284 class kbi0_win {
00285 KaiserBessel& kb;
00286 public:
00287 kbi0_win(KaiserBessel& kb_) : kb(kb_) {}
00288 float operator()(float x) const {
00289 return kb.i0win(x);
00290 }
00291 int get_window_size() const {return kb.get_window_size();}
00292 };
00294 kbi0_win get_kbi0_win() {
00295 return kbi0_win(*this);
00296 }
00297 };
00298
00299 class FakeKaiserBessel : public KaiserBessel {
00300 public:
00301 FakeKaiserBessel(float alpha, int K, float r_,
00302 float v_, int N_, float vtable_=0.f,
00303 int ntable_ = 5999)
00304 : KaiserBessel(alpha, K, r_, v_, N_, vtable_, ntable_) {
00305 build_I0table();
00306 }
00307 float sinhwin(float x) const;
00308 float i0win(float x) const;
00309 void build_I0table();
00310 };
00311
00323 static vector<float>
00324 even_angles(float delta, float t1=0, float t2=90, float p1=0, float p2=359.999);
00325
00326
00379 static float quadri(float x, float y, int nx, int ny, float* image);
00380
00435 static float quadri_background(float x, float y, int nx, int ny, float* image, int xnew, int ynew);
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 static float get_pixel_conv_new(int nx, int ny, int nz, float delx, float dely, float delz, float* data, Util::KaiserBessel& kb);
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463 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);
00464
00465 static std::complex<float> extractpoint2(int nx, int ny, float nuxnew, float nuynew, EMData *fimage, Util::KaiserBessel& kb);
00466
00467
00468 static float bilinear(float xold, float yold, int nsam, int nrow, float* xim);
00469
00470
00479 static float triquad(float r, float s, float t, float* fdata);
00480
00488 class Gaussian {
00489 float sigma;
00490 float rttwopisigma;
00491 float twosigma2;
00492 public:
00493 Gaussian(float sigma_ = 1.0) : sigma(sigma_) {
00494 rttwopisigma = sqrtf(static_cast<float>(twopi)*sigma);
00495 twosigma2 = 2*sigma*sigma;
00496 }
00497 inline float operator()(float x) const {
00498 return exp(-x*x/(twosigma2))/rttwopisigma;
00499 }
00500 };
00501
00502
00503 static EMData* Polar2D(EMData* image, vector<int> numr, string mode);
00504 static EMData* Polar2Dm(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode);
00505
00506
00507 static void alrl_ms(float *xim, int nsam, int nrow, float cns2, float cnr2,
00508 int *numr, float *circ, int lcirc, int nring, char mode);
00509
00510
00511
00512 static EMData* Polar2Dmi(EMData* image, float cns2, float cnr2, vector<int> numr, string cmode, Util::KaiserBessel& kb);
00513
00514 static void fftr_q(float *xcmplx, int nv);
00515 static void fftr_d(double *xcmplx, int nv);
00516 static void fftc_q(float *br, float *bi, int ln, int ks);
00517 static void fftc_d(double *br, double *bi, int ln, int ks);
00518
00520 static void Frngs(EMData* circ, vector<int> numr);
00521 static void Normalize_ring(EMData* ring, const vector<int>& numr);
00522
00524 static void Frngs_inv(EMData* circ, vector<int> numr);
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 static Dict Crosrng_e(EMData* circ1, EMData* circ2, vector<int> numr, int neg);
00540 static Dict Crosrng_ew(EMData* circ1, EMData* circ2, vector<int> numr, vector<float> w, int neg);
00541
00542 static Dict Crosrng_ms(EMData* circ1, EMData* circ2, vector<int> numr);
00543 static Dict Crosrng_ms_delta(EMData* circ1, EMData* circ2, vector<int> numr, float delta_start, float delta);
00544
00550 static Dict Crosrng_sm_psi(EMData* circ1, EMData* circ2, vector<int> numr, float psi, int flag);
00551
00558 static Dict Crosrng_psi_0_180(EMData* circ1, EMData* circ2, vector<int> numr, float psi_max);
00559 static Dict Crosrng_ns(EMData* circ1, EMData* circ2, vector<int> numr);
00560
00567 static EMData* Crosrng_msg(EMData* circ1, EMData* circ2, vector<int> numr);
00568
00575 static void Crosrng_msg_vec(EMData* circ1, EMData* circ2, vector<int> numr, float *q, float *t);
00576
00583 static EMData* Crosrng_msg_s(EMData* circ1, EMData* circ2, vector<int> numr);
00584
00591 static EMData* Crosrng_msg_m(EMData* circ1, EMData* circ2, vector<int> numr);
00592
00593 static vector<float> Crosrng_msg_vec_p(EMData* circ1, EMData* circ2, vector<int> numr );
00594 static void prb1d(double *b, int npoint, float *pos);
00595
00596 static void update_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00597 static void sub_fav(EMData* ave,EMData* dat, float tot, int mirror, vector<int> numr);
00598
00599
00600 static float ener(EMData* ave, vector<int> numr);
00601
00602 static float ener_tot(const vector<EMData*>& data, vector<int> numr, vector<float> tot);
00603
00605 static Dict min_dist_real(EMData* image, const vector<EMData*>& data);
00606
00608 static Dict min_dist_four(EMData* image, const vector<EMData*>& data);
00609
00616 static int k_means_cont_table_(int* group1, int* group2, int* stb, long int s1, long int s2, int flag);
00617
00618
00619
00629 static void bb_enumerate_(int* Parts, int* classDims, int nParts, int nClasses, int T, int n_guesses, int* levels);
00630
00634 static void initial_prune(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T);
00635
00639 static bool explore(vector <vector <int*> > & Parts, int* dimClasses, int nParts, int K, int T, int partref, int* curintx, int
00640 size_curintx, int* next, int size_next, int depth);
00641
00642 static int* branch(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* Levels, int nLevels, int curlevel, int n_guesses);
00643
00651 static int findTopLargest(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* matchlist, int max_num_matches, int*
00652 costlist, int n_guesses);
00653
00658 static int generatesubmax(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int n_guesses);
00659
00663 static void search2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int newT, int* curmax);
00664
00665 static int* explore2(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int newT, int* curintx, int size_curintx, int* next, int
00666 size_next, int depth);
00667
00672 static bool sanitycheck(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* output);
00673
00683 static vector<int> bb_enumerateMPI_(int* argParts, int* dimClasses, int nParts, int K, int T, int nTop, int n_guesses, bool doMPI, int* Levels);
00684
00690 static vector<int> branchMPIpy_(int* argParts, int* dimClasses, int nParts, int K, int T, int* Levels, int nLevels, int n_guesses, int nFirst, int* firstmatches);
00691
00697 static int* branchMPI(int* argParts, int* Indices, int* dimClasses, int nParts, int K, int T, int* Levels, int nLevels, int curlevel,int n_guesses, int nFirst, int* firstmatches);
00698
00699
00700
00701
00702 static vector<double> cml_weights(const vector<float>& cml);
00703
00705 static vector<int> cml_line_insino(vector<float> Rot, int i_prj, int n_prj);
00706
00708 static vector<int> cml_line_insino_all(vector<float> Rot, vector<int> seq, int n_prj, int n_lines);
00709
00711 static vector<double> cml_init_rot(vector<float> Ori);
00712
00714 static vector<float> cml_update_rot(vector<float> Rot, int iprj, float nph, float th, float nps);
00715
00717 static vector<double> cml_line_in3d(vector<float> Ori, vector<int> seq, int nprj, int nlines);
00718
00720 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);
00721
00723 static double cml_disc(const vector<EMData*>& data, vector<int> com, vector<int> seq, vector<float> weights, int n_lines);
00724
00730 static void set_line(EMData* img, int posline, EMData* line, int offset, int length);
00731
00739 static void cml_prepare_line(EMData* sino, EMData* line, int ilf, int ihf, int pos_line, int nblines);
00740
00741
00742
00743
00744
00745
00746
00747
00748 static EMData* decimate(EMData* img, int x_step,int y_step=1,int z_step=1);
00749
00750 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);
00751
00752 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");
00753
00754 static vector<float> histogram(EMData* image, EMData* mask, int nbins = 128, float hmin =0.0f, float hmax = 0.0f );
00755
00756 static Dict histc(EMData *ref,EMData *img,EMData *mask);
00757
00758 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);
00759
00760
00761
00762
00763
00764
00765
00766 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);
00767 static EMData *compress_image_mask(EMData* image, EMData* mask);
00768
00770 static EMData *reconstitute_image_mask(EMData *image,EMData *mask);
00771
00772 static vector<float> merge_peaks(vector<float> peak1, vector<float> peak2,float p_size);
00773 static vector<float> pw_extract(vector<float>pw, int n, int iswi,float ps);
00774 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);
00775 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);
00776 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
00777 int *iu, double *s);
00778 static float eval(char * images,EMData * img, vector<int> S,int N, int K,int size);
00779
00780
00781 static vector<double> vrdg(const vector<float>& ph, const vector<float>& th);
00782 static void hsortd(double *theta,double *phi,int *key,int len,int option);
00783 static void voronoidiag(double *theta,double *phi,double* weight,int n);
00784
00785
00786 static void voronoi(double *phi,double *theta,double *weight, int nt);
00787 static void disorder2(double *x,double *y,int *key,int len);
00788 static void ang_to_xyz(double *x,double *y,double *z,int len);
00789 static void flip23(double *x,double *y,double *z,int *key,int k,int len);
00790 struct tmpstruct{
00791 double theta1,phi1;
00792 int key1;
00793 };
00794 static bool cmp1(tmpstruct tmp1,tmpstruct tmp2);
00795 static bool cmp2(tmpstruct tmp1,tmpstruct tmp2);
00796
00797
00798 static int trmsh3_(int *n0, double *tol, double *x, double *y, double *z__, int *n, int *list, int *lptr,
00799 int *lend, int *lnew, int *indx, int *lcnt, int *near__, int *next, double *dist, int *ier);
00800 static double areav_(int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier);
00801
00802
00803
00804
00805
00806 static EMData* madn_scalar(EMData* img, EMData* img1, float scalar);
00807
00808 static EMData* mult_scalar(EMData* img, float scalar);
00809
00810 static EMData* addn_img(EMData* img, EMData* img1);
00811
00812 static EMData* subn_img(EMData* img, EMData* img1);
00813
00814 static EMData* muln_img(EMData* img, EMData* img1);
00815
00816 static EMData* divn_img(EMData* img, EMData* img1);
00817
00818 static EMData* divn_filter(EMData* img, EMData* img1);
00819
00820
00821 static void mad_scalar(EMData* img, EMData* img1, float scalar);
00822
00823 static void mul_scalar(EMData* img, float scalar);
00824
00825 static void add_img(EMData* img, EMData* img1);
00826
00827 static void add_img_abs(EMData* img, EMData* img1);
00828
00829 static void add_img2(EMData* img, EMData* img1);
00830
00831 static void sub_img(EMData* img, EMData* img1);
00832
00833 static void mul_img(EMData* img, EMData* img1);
00834
00835 static void div_img(EMData* img, EMData* img1);
00836
00837 static void div_filter(EMData* img, EMData* img1);
00838
00839 static EMData* pack_complex_to_real(EMData* img);
00840 private:
00841 static float ang_n(float peakp, string mode, int maxrin);
00842 public:
00843
00849 static vector<float> multiref_polar_ali_2d(EMData* image, const vector< EMData* >& crefim,
00850 float xrng, float yrng, float step, string mode,
00851 vector< int >numr, float cnx, float cny);
00852
00858 static vector<float> multiref_polar_ali_2d_delta(EMData* image, const vector< EMData* >& crefim,
00859 float xrng, float yrng, float step, string mode,
00860 vector< int >numr, float cnx, float cny, float delta_start, float delta);
00861
00867 static vector<float> multiref_polar_ali_2d_nom(EMData* image, const vector< EMData* >& crefim,
00868 float xrng, float yrng, float step, string mode,
00869 vector< int >numr, float cnx, float cny);
00870
00876 static vector<float> multiref_polar_ali_2d_local(EMData* image, const vector< EMData* >& crefim,
00877 float xrng, float yrng, float step, float ant, string mode,
00878 vector< int >numr, float cnx, float cny);
00879
00886 static vector<float> multiref_polar_ali_helical(EMData* image, const vector< EMData* >& crefim,
00887 float xrng, float yrng, float step, float psi_max, string mode,
00888 vector< int >numr, float cnx, float cny);
00889
00895 static vector<float> multiref_polar_ali_2d_local_psi(EMData* image, const vector< EMData* >& crefim,
00896 float xrng, float yrng, float step, float ant, float psi_max, string mode,
00897 vector< int >numr, float cnx, float cny);
00898
00906 static void multiref_peaks_ali2d(EMData* image, EMData* crefim,
00907 float xrng, float yrng, float step, string mode,
00908 vector< int >numr, float cnx, float cny, EMData* peaks, EMData* peakm);
00909
00917 static void multiref_peaks_compress_ali2d(EMData* image, EMData* crefim, float xrng, float yrng,
00918 float step, string mode, vector<int>numr, float cnx, float cny, EMData *peaks, EMData *peakm,
00919 EMData *peaks_compress, EMData *peakm_compress);
00920
00925 static vector<float> ali2d_ccf_list(EMData* image, EMData* crefim, float xrng, float yrng,
00926 float step, string mode, vector<int>numr, float cnx, float cny, double T);
00927
00928
00929
00930
00931
00932
00933 static vector<float> twoD_fine_ali(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00934
00935 static vector<float> twoD_fine_ali_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00936
00937 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);
00938
00939 static vector<float> twoD_fine_ali_SD(EMData* image, EMData *refim, EMData* mask, float ang, float sxs, float sys);
00940
00941 static float ccc_images(EMData *, EMData *, EMData *, float , float , float );
00942
00943 static vector<float> twoD_fine_ali_SD_G(EMData* image, EMData *refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sxs, float sys);
00944
00945 static float ccc_images_G(EMData* image, EMData* refim, EMData* mask, Util::KaiserBessel& kb, float ang, float sx, float sy);
00946
00947 static EMData* move_points(EMData* img, float qprob, int ri, int ro);
00948
00949 static EMData* get_biggest_cluster( EMData* mg );
00950
00951
00952 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);
00953
00954 static inline int mono(int k1, int k2) {
00955 #ifdef _WIN32
00956 int mk = _cpp_max(k1,k2);
00957 return _cpp_min(k1,k2) + mk*(mk-1)/2;
00958 #else
00959 int mk = std::max(k1,k2);
00960 return std::min(k1,k2) + mk*(mk-1)/2;
00961 #endif //_WIN32
00962 }
00963
00964 static inline int nint180(float arg) {
00965 int res = int(arg + 180.5) - 180;
00966 return res;
00967 }
00968
00969 static vector<float> cluster_pairwise(EMData* d, int K, float T, float F);
00970
00971 static vector<float> cluster_equalsize(EMData* d);
00972 static vector<float> vareas(EMData* d);
00973
00979 static EMData* get_slice(EMData *vol, int dim, int index);
00980
00981 static void image_mutation(EMData *img, float mutation_rate);
00982
00984 static void array_mutation(float* list, int len_list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
00985
00986 static vector<float> list_mutation(vector<float> list, float mutation_rate, float min_val, float max_val, int K, int is_mirror);
00987
00988
00989
00990 static inline float restrict1(float x, int nx) {
00991 while ( x < 0.0f ) x += nx;
00992 while ( x >= (float)(nx) ) x -= nx;
00993 return x;
00994 }
00995
00996 #endif //util__sparx_h__