EMAN::Util Class Reference
[unit test in Python]

Util is a collection of utility functions. More...

#include <util.h>

List of all members.

Static Public Member Functions

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

For those util function developed by Pawel's group.

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)

Creates a Two D Test Pattern.

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

Given a tabulated function y of x (n unordered points), and Given the values of the m values xq to be interpolated This routine returns the interpolated array yq, PRB This function is called by splint.

static void spline (float *x, float *y, int n, float yp1, float ypn, float *y2)

Given a tabulated function y of x (unordered), and Given the values of the first derivatives at the end points This routine returns an array y2, that contains the second derivatives of the function at the tabulated points.

static void splint (float *xa, float *ya, float *y2a, int n, float *xq, float *yq, int m)

Given the arrays xa(ordered, ya of length n, which tabulate a function and given the array y2a which is the output of spline and an unordered array xq, this routine returns a cubic-spline interpolated array yq.

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

list the sorted lengths of the integer lattice sites of a square sided image of size Size.

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

Compute a vector containing quasi-evenly spaced Euler angles.

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

Quadratic interpolation (2D).

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

Quadratic interpolation (2D).

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 bilinear (float xold, float yold, int nsam, int nrow, float *xim)

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

Quadratic interpolation (3D).

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 alrl_ms (float *xim, int nsam, int nrow, float cns2, float cnr2, int *numr, float *circ, int lcirc, int nring, char mode)

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)

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)

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)

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)

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)

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)

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

static bool cmp1 (tmpstruct tmp1, tmpstruct tmp2)

static bool cmp2 (tmpstruct tmp1, tmpstruct tmp2)

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)

static EMData * madn_scalar (EMData *img, EMData *img1, float scalar)

static EMData * mult_scalar (EMData *img, float scalar)

static EMData * addn_img (EMData *img, EMData *img1)

static EMData * subn_img (EMData *img, EMData *img1)

static EMData * muln_img (EMData *img, EMData *img1)

static EMData * divn_img (EMData *img, EMData *img1)

static EMData * divn_filter (EMData *img, EMData *img1)

static void mad_scalar (EMData *img, EMData *img1, float scalar)

static void mul_scalar (EMData *img, float scalar)

static void add_img (EMData *img, EMData *img1)

static void add_img2 (EMData *img, EMData *img1)

static void sub_img (EMData *img, EMData *img1)

static void mul_img (EMData *img, EMData *img1)

static void div_img (EMData *img, EMData *img1)

static void div_filter (EMData *img, EMData *img1)

static EMData * pack_complex_to_real (EMData *img)

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 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, float cs, float wgh, float b_factor, float dza, float azz, float sign)

static int mono (int k1, int k2)

static int nint180 (float arg)

static vector< float > cluster_pairwise (EMData *d, int K, float T, float F)

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)

static float restrict1 (float x, int nx)

static void ap2ri (float *data, size_t n)

convert complex data array from Amplitude/Phase format into Real/Imaginary format.

static void flip_complex_phase (float *data, size_t n)

flip the phase of a complex data array.

static void rotate_phase_origin (float *data, size_t nx, size_t ny, size_t nz)

rotate data vertically by ny/2, to make the mrc phase origin compliance with EMAN1 and TNF reader

static int file_lock_wait (FILE *file)

lock a file.

static bool check_file_by_magic (const void *first_block, const char *magic)

check whether a file starts with certain magic string.

static bool is_file_exist (const string &filename)

check whether a file exists or not

static void flip_image (float *data, size_t nx, size_t ny)

Vertically flip the data of a 2D real image.

static vector< EMData * > svdcmp (const vector< EMData * > &data, int nvec)

Perform singular value decomposition on a set of images.

static string str_to_lower (const string &s)

Return a lower case version of the argument string.

static bool sstrncmp (const char *s1, const char *s2)

Safe string compare.

static string int2str (int n)

Get a string format of an integer, e.g.

static string get_line_from_string (char **str)

Extract a single line from a multi-line string.

static bool get_str_float (const char *str, const char *float_var, float *p_val)

Extract the float value from a variable=value string with format like "XYZ=1.1", where 'str' is "XYZ=1.1"; 'float_var' is "XYZ="; 'p_val' points to float number 1.1.

static bool get_str_float (const char *str, const char *float_var, float *p_v1, float *p_v2)

Extract the float values from a variable=value1,value2 string with format like "XYZ=1.1,1.2", where 'str' is "XYZ=1.1,1.2"; 'float_var' is "XYZ="; 'p_v1' points to 1.1; 'p_v2' points to 1.2.

static bool get_str_float (const char *str, const char *float_var, int *p_nvalues, float *p_v1, float *p_v2)

Extract number of values and the float values, if any, from a string whose format is either "variable=value1,value2 " or "variable".

static bool get_str_int (const char *str, const char *int_var, int *p_val)

Extract the int value from a variable=value string with format like "XYZ=1", where 'str' is "XYZ=1"; 'int_var' is "XYZ="; 'p_val' points to float number 1.

static bool get_str_int (const char *str, const char *int_var, int *p_v1, int *p_v2)

Extract the int value from a variable=value1,value2 string with format like "XYZ=1,2", where 'str' is "XYZ=1,2"; 'int_var' is "XYZ="; 'p_val' points to float number 1.

static bool get_str_int (const char *str, const char *int_var, int *p_nvalues, int *p_v1, int *p_v2)

Extract number of values and the int values, if any, from a string whose format is either "variable=value1,value2 " or "variable".

static string change_filename_ext (const string &old_filename, const string &new_ext)

Change a file's extension and return the new filename.

static string remove_filename_ext (const string &filename)

Remove a filename's extension and return the new filename.

static string get_filename_ext (const string &filename)

Get a filename's extension.

static string sbasename (const string &filename)

Get a filename's basename.

static void calc_least_square_fit (size_t nitems, const float *data_x, const float *data_y, float *p_slope, float *p_intercept, bool ignore_zero, float absmax=0)

calculate the least square fit value.

static Vec3f calc_bilinear_least_square (const vector< float > &points)

calculate bilinear least-square fit, z = a + b x + c y Takes a set of x,y,z vectors and produces an a,b,c vector does not accept error bars on z or return goodness of fit

static void save_data (const vector< float > &x_array, const vector< float > &y_array, const string &filename)

Save (x y) data array into a file.

static void save_data (float x0, float dx, const vector< float > &y_array, const string &filename)

Save x, y data into a file.

static void save_data (float x0, float dx, float *y_array, size_t array_size, const string &filename)

Save x, y data into a file.

static void sort_mat (float *left, float *right, int *leftPerm, int *rightPerm)

does a sort as in Matlab.

static unsigned long int get_randnum_seed ()

Get the seed for Randnum class.

static void set_randnum_seed (unsigned long int seed)

Set the seed for Randnum class.

static int get_irand (int low, int high)

Get an integer random number between low and high, [low, high].

static float get_frand (int low, int high)

Get a float random number between low and high, [low, high).

static float get_frand (float low, float high)

Get a float random number between low and high, [low, high).

static float get_frand (double low, double high)

Get a float random number between low and high, [low, high).

static float get_gauss_rand (float mean, float sigma)

Get a Gaussian random number.

static int round (float x)

Get ceiling round of a float number x.

static int round (double x)

Get ceiling round of a float number x.

static float linear_interpolate (float p1, float p2, float t)

Calculate linear interpolation.

static float bilinear_interpolate (float p1, float p2, float p3, float p4, float t, float u)

Calculate bilinear interpolation.

static float trilinear_interpolate (float p1, float p2, float p3, float p4, float p5, float p6, float p7, float p8, float t, float u, float v)

Calculate trilinear interpolation.

static void find_max (const float *data, size_t nitems, float *p_max_val, int *p_max_index=0)

Find the maximum value and (optional) its index in an array.

static void find_min_and_max (const float *data, size_t nitems, float *p_max_val, float *p_min_val, int *p_max_index=0, int *p_min_index=0)

Find the maximum value and (optional) its index, minimum value and (optional) its index in an array.

static Dict get_stats (const vector< float > &data)

Get the mean, standard deviation, skewness and kurtosis of the input data.

static Dict get_stats_cstyle (const vector< float > &data)

static int calc_best_fft_size (int low)

Search the best FFT size with good primes.

static EMData * calc_bessel (const int n, const float &x)

static int square (int n)

Calculate a number's square.

static float square (float x)

Calculate a number's square.

static float square (double x)

Calculate a number's square.

static float square_sum (float x, float y)

Calcuate (x*x + y*y).

static float hypot2 (float x, float y)

Euclidean distance function in 2D: f(x,y) = sqrt(x*x + y*y);.

static float hypot3 (int x, int y, int z)

Euclidean distance function in 3D: f(x,y,z) = sqrt(x*x + y*y + z*z);.

static float hypot3 (float x, float y, float z)

Euclidean distance function in 3D: f(x,y,z) = sqrt(x*x + y*y + z*z);.

static float hypot3 (double x, double y, double z)

Euclidean distance function in 3D: f(x,y,z) = sqrt(x*x + y*y + z*z);.

static float hypot_fast (int x, int y)

Euclidean distance in 2D for integers computed fast using a cached lookup table.

static short hypot_fast_int (int x, int y)

Euclidean distance in 2D for integers computed fast using a cached lookup table.

static int fast_floor (float x)

A fast way to calculate a floor, which is largest integral value not greater than argument.

static float agauss (float a, float dx, float dy, float dz, float d)

Calculate Gaussian value.

static int get_min (int f1, int f2)

Get the minimum of 2 numbers.

static int get_min (int f1, int f2, int f3)

Get the minimum of 3 numbers.

static float get_min (float f1, float f2)

Get the minimum of 2 numbers.

static float get_min (float f1, float f2, float f3)

Get the minimum of 3 numbers.

static float get_min (float f1, float f2, float f3, float f4)

Get the minimum of 4 numbers.

static float get_max (float f1, float f2)

Get the maximum of 2 numbers.

static float get_max (float f1, float f2, float f3)

Get the maximum of 3 numbers.

static float get_max (float f1, float f2, float f3, float f4)

Get the maximum of 4 numbers.

static float angle_sub_2pi (float x, float y)

Calculate the difference of 2 angles and makes the equivalent result to be less than Pi.

static float angle_sub_pi (float x, float y)

Calculate the difference of 2 angles and makes the equivalent result to be less than Pi/2.

static int goodf (const float *p_f)

Check whether a number is a good float.

static int goodf (const double *p_f)

static string get_time_label ()

Get the current time in a string with format "mm/dd/yyyy hh:mm".

static void set_log_level (int argc, char *argv[])

Set program logging level through command line option "-v N", where N is the level.

static float eman_copysign (float a, float b)

copy sign of a number.

static float eman_erfc (float x)

complementary error function.

static void equation_of_plane (const Vec3f &p1, const Vec3f &p2, const Vec3f &p3, float *plane)

Determine the equation of a plane that intersects 3 points in 3D space.

static bool point_is_in_convex_polygon_2d (const Vec2f &p1, const Vec2f &p2, const Vec2f &p3, const Vec2f &p4, const Vec2f &actual_point)

Determines if a point is in a 2D convex polygon described by 4 points using the Barycentric method, which is a fast way of performing the query.

static bool point_is_in_triangle_2d (const Vec2f &p1, const Vec2f &p2, const Vec2f &p3, const Vec2f &actual_point)

Determines if a point is in a 2D triangle using the Barycentric method, which is a fast way of performing the query Triangle points can be specified in any order.

static void printMatI3D (MIArray3D &mat, const string str=string(""), ostream &out=std::cout)

Print a 3D integer matrix to a file stream (std out by default).

template<class T>

static T sgn (T &val)

Sign function.

static float * getBaldwinGridWeights (const int &freq_cutoff, const float &P, const float &r, const float &dfreq=1, const float &alpha=0.5, const float &beta=0.2)

Get the isosurface value for 3D image.

static bool IsPower2 (int x)

Return true if an integer is positive and is power of 2.

static void apply_precision (float &value, const float &precision)

Static Private Member Functions

static float ang_n (float peakp, string mode, int maxrin)

Classes

class FakeKaiserBessel

class Gaussian

Gaussian function class. More...

class KaiserBessel

1-D Kaiser-Bessel window function class. More...

class sincBlackman

struct tmpstruct

Detailed Description

Util is a collection of utility functions.

Definition at line 80 of file util.h.

Member Function Documentation

int Util::coveig	(	int	n,
		float *	covmat,
		float *	eigval,
		float *	eigvec
	)			`[static]`

For those util function developed by Pawel's group.

$Id$ This file is a part of util.h, To use this file's functions, you should include "util.h" NEVER directly include this file

Definition at line 5915 of file util_sparx.cpp.

References ENTERFUNC, EXITFUNC, and ssyev_().

05916 {
05917         // n size of the covariance/correlation matrix
05918         // covmat --- covariance/correlation matrix (n by n)
05919         // eigval --- returns eigenvalues
05920         // eigvec --- returns eigenvectors
05921 
05922         ENTERFUNC;
05923 
05924         int i;
05925 
05926         // make a copy of covmat so that it will not be overwritten
05927         for ( i = 0 ; i < n * n ; i++ )   eigvec[i] = covmat[i];
05928 
05929         char NEEDV = 'V';
05930         char UPLO = 'U';
05931         int lwork = -1;
05932         int info = 0;
05933         float *work, wsize;
05934 
05935         //  query to get optimal workspace
05936         ssyev_(&NEEDV, &UPLO, &n, eigvec, &n, eigval, &wsize, &lwork, &info);
05937         lwork = (int)wsize;
05938 
05939         work = (float *)calloc(lwork, sizeof(float));
05940         //  calculate eigs
05941         ssyev_(&NEEDV, &UPLO, &n, eigvec, &n, eigval, work, &lwork, &info);
05942         free(work);
05943         EXITFUNC;
05944         return info;
05945 }

Dict Util::coveig_for_py	(	int	ncov,
		const vector< float > &	covmatpy
	)			`[static]`

Definition at line 5948 of file util_sparx.cpp.

References coveig(), covmat, eigval, eigvec, ENTERFUNC, EXITFUNC, and status.

05949 {
05950 
05951         ENTERFUNC;
05952         int len = covmatpy.size();
05953         float *eigvec;
05954         float *eigval;
05955         float *covmat;
05956         int status = 0;
05957         eigval = (float*)calloc(ncov,sizeof(float));
05958         eigvec = (float*)calloc(ncov*ncov,sizeof(float));
05959         covmat = (float*)calloc(ncov*ncov, sizeof(float));
05960 
05961         const float *covmat_ptr;
05962         covmat_ptr = &covmatpy[0];
05963         for(int i=0;i<len;i++){
05964             covmat[i] = covmat_ptr[i];
05965         }
05966 
05967         status = Util::coveig(ncov, covmat, eigval, eigvec);
05968 
05969         vector<float> eigval_py(ncov);
05970         const float *eigval_ptr;
05971         eigval_ptr = &eigval[0];
05972         for(int i=0;i<ncov;i++){
05973             eigval_py[i] = eigval_ptr[i];
05974         }
05975 
05976         vector<float> eigvec_py(ncov*ncov);
05977         const float *eigvec_ptr;
05978         eigvec_ptr = &eigvec[0];
05979         for(int i=0;i<ncov*ncov;i++){
05980             eigvec_py[i] = eigvec_ptr[i];
05981         }
05982 
05983         Dict res;
05984         res["eigval"] = eigval_py;
05985         res["eigvec"] = eigvec_py;
05986 
05987         EXITFUNC;
05988         return res;
05989 }

Dict Util::ExpMinus4YSqr	(	float	ymax,
		int	nsamples
	)			`[static]`

Definition at line 5746 of file util_sparx.cpp.

05747 {
05748         //exp(-16) is 1.0E-7 approximately)
05749         vector<float> expvect;
05750 
05751         double inc = double(ymax)/nsamples;
05752         double temp;
05753         for(int i =0;i<nsamples;i++) {
05754                 temp = exp((-4*(i*inc)*(i*inc)));
05755                 expvect.push_back(float(temp));
05756         }
05757         expvect.push_back(0.0);
05758         Dict lookupdict;
05759         lookupdict["table"]    = expvect;
05760         lookupdict["ymax"]     = ymax;
05761         lookupdict["nsamples"] = nsamples;
05762 
05763         return lookupdict;
05764 }

void Util::WTM	(	EMData *	PROJ,
		vector< float >	SS,
		int	DIAMETER,
		int	NUMP
	)			`[static]`

Definition at line 5620 of file util_sparx.cpp.

References AMAX1, AMIN1, CC, CP, EMAN::EMData::depad(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::norm_pad(), PROJ, RI, EMAN::EMData::set_size(), sqrt(), SS, EMAN::EMData::to_one(), EMAN::EMData::update(), VP, VV, and W.

05621 {
05622         float rad2deg =(180.0f/3.1415926f);
05623         float deg2rad = (3.1415926f/180.0f);
05624 
05625         int NSAM,NROW,NNNN,NR2,NANG,L,JY;
05626 
05627         NSAM = PROJ->get_xsize();
05628         NROW = PROJ->get_ysize();
05629         NNNN = NSAM+2-(NSAM%2);
05630         NR2  = NROW/2;
05631         NANG = int(SS.size())/6;
05632 
05633         float RI[9];
05634         RI(1,1)=SS(1,NUMP)*SS(3,NUMP)*SS(5,NUMP)-SS(2,NUMP)*SS(6,NUMP);
05635         RI(2,1)=-SS(1,NUMP)*SS(3,NUMP)*SS(6,NUMP)-SS(2,NUMP)*SS(5,NUMP);
05636         RI(3,1)=SS(1,NUMP)*SS(4,NUMP);
05637         RI(1,2)=SS(2,NUMP)*SS(3,NUMP)*SS(5,NUMP)+SS(1,NUMP)*SS(6,NUMP);
05638         RI(2,2)=-SS(2,NUMP)*SS(3,NUMP)*SS(6,NUMP)+SS(1,NUMP)*SS(5,NUMP);
05639         RI(3,2)=SS(2,NUMP)*SS(4,NUMP);
05640         RI(1,3)=-SS(4,NUMP)*SS(5,NUMP);
05641         RI(2,3)=SS(4,NUMP)*SS(6,NUMP);
05642         RI(3,3)=SS(3,NUMP);
05643 
05644         float THICK=static_cast<float>( NSAM)/DIAMETER/2.0f ;
05645 
05646         EMData* W = new EMData();
05647         int Wnx = NNNN/2;
05648         W->set_size(NNNN/2,NROW,1);
05649         W->to_one();
05650         float *Wptr = W->get_data();
05651 
05652         float ALPHA,TMP,FV,RT,FM,CCN,CC[3],CP[2],VP[2],VV[3];
05653 
05654         for (L=1; L<=NANG; L++) {
05655                 if (L != NUMP) {
05656                         CC(1)=SS(2,L)*SS(4,L)*SS(3,NUMP)-SS(3,L)*SS(2,NUMP)*SS(4,NUMP);
05657                         CC(2)=SS(3,L)*SS(1,NUMP)*SS(4,NUMP)-SS(1,L)*SS(4,L)*SS(3,NUMP);
05658                         CC(3)=SS(1,L)*SS(4,L)*SS(2,NUMP)*SS(4,NUMP)-SS(2,L)*SS(4,L)*SS(1,NUMP)*SS(4,NUMP);
05659 
05660                         TMP = sqrt(CC(1)*CC(1) +  CC(2)*CC(2) + CC(3)*CC(3));
05661                         CCN=static_cast<float>( AMAX1( AMIN1(TMP,1.0) ,-1.0) );
05662                         ALPHA=rad2deg*float(asin(CCN));
05663                         if (ALPHA>180.0f) ALPHA=ALPHA-180.0f;
05664                         if (ALPHA>90.0f) ALPHA=180.0f-ALPHA;
05665                         if(ALPHA<1.0E-6) {
05666                                 for(int J=1;J<=NROW;J++) for(int I=1;I<=NNNN/2;I++) W(I,J)+=1.0;
05667                         } else {
05668                                 FM=THICK/(fabs(sin(ALPHA*deg2rad)));
05669                                 CC(1)   = CC(1)/CCN;CC(2)   = CC(2)/CCN;CC(3)   = CC(3)/CCN;
05670                                 VV(1)= SS(2,L)*SS(4,L)*CC(3)-SS(3,L)*CC(2);
05671                                 VV(2)= SS(3,L)*CC(1)-SS(1,L)*SS(4,L)*CC(3);
05672                                 VV(3)= SS(1,L)*SS(4,L)*CC(2)-SS(2,L)*SS(4,L)*CC(1);
05673                                 CP(1)   = 0.0;CP(2) = 0.0;
05674                                 VP(1)   = 0.0;VP(2) = 0.0;
05675 
05676                                 CP(1) = CP(1) + RI(1,1)*CC(1) + RI(1,2)*CC(2) + RI(1,3)*CC(3);
05677                                 CP(2) = CP(2) + RI(2,1)*CC(1) + RI(2,2)*CC(2) + RI(2,3)*CC(3);
05678                                 VP(1) = VP(1) + RI(1,1)*VV(1) + RI(1,2)*VV(2) + RI(1,3)*VV(3);
05679                                 VP(2) = VP(2) + RI(2,1)*VV(1) + RI(2,2)*VV(2) + RI(2,3)*VV(3);
05680 
05681                                 TMP = CP(1)*VP(2)-CP(2)*VP(1);
05682 
05683                                 //     PREVENT TMP TO BE TOO SMALL, SIGN IS IRRELEVANT
05684                                 TMP = AMAX1(1.0E-4f,fabs(TMP));
05685                                 float tmpinv = 1.0f/TMP;
05686                                 for(int J=1;J<=NROW;J++) {
05687                                         JY = (J-1);
05688                                         if (JY>NR2)  JY=JY-NROW;
05689                                         for(int I=1;I<=NNNN/2;I++) {
05690                                                 FV     = fabs((JY*CP(1)-(I-1)*CP(2))*tmpinv);
05691                                                 RT     = 1.0f-FV/FM;
05692                                                 W(I,J) += ((RT>0.0f)*RT);
05693                                         }
05694                                 }
05695                         }
05696                 }
05697         }
05698 
05699         EMData* proj_in = PROJ;
05700 
05701         PROJ = PROJ->norm_pad( false, 2);
05702         PROJ->do_fft_inplace();
05703         PROJ->update();
05704         float *PROJptr = PROJ->get_data();
05705 
05706         int KX;
05707         float WW;
05708         for(int J=1; J<=NROW; J++)
05709                 for(int I=1; I<=NNNN; I+=2) {
05710                         KX          =  (I+1)/2;
05711                         WW          =  1.0f/W(KX,J);
05712                         PROJ(I,J)   = PROJ(I,J)*WW;
05713                         PROJ(I+1,J) = PROJ(I+1,J)*WW;
05714                 }
05715         delete W; W = 0;
05716         PROJ->do_ift_inplace();
05717         PROJ->depad();
05718 
05719         float* data_src = PROJ->get_data();
05720         float* data_dst = proj_in->get_data();
05721 
05722         int ntotal = NSAM*NROW;
05723         for( int i=0; i < ntotal; ++i )
05724         {
05725             data_dst[i] = data_src[i];
05726         }
05727 
05728         proj_in->update();
05729         delete PROJ;
05730 }

void Util::WTF	(	EMData *	PROJ,
		vector< float >	SS,
		float	SNR,
		int	K,
		vector< float >	exptable
	)			`[static]`

Definition at line 5528 of file util_sparx.cpp.

References EMAN::EMData::depad(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::norm_pad(), PROJ, SS, EMAN::EMData::update(), W, and Y.

05529 {
05530         int NSAM,NROW,NNNN,NR2,L,JY,KX,NANG;
05531         float WW,OX,OY,Y;
05532 
05533         NSAM = PROJ->get_xsize();
05534         NROW = PROJ->get_ysize();
05535         NNNN = NSAM+2-(NSAM%2);
05536         NR2  = NROW/2;
05537 
05538         NANG = int(SS.size())/6;
05539 
05540         EMData* W = new EMData();
05541         int Wnx = NNNN/2;
05542         W->set_size(Wnx,NROW,1);
05543         W->to_zero();
05544         float *Wptr = W->get_data();
05545         float *PROJptr = PROJ->get_data();
05546         float indcnst = 1000/2.0;
05547         // we create look-up table for 1001 uniformly distributed samples [0,2];
05548 
05549         for (L=1; L<=NANG; L++) {
05550                 OX = SS(6,K)*SS(4,L)*(-SS(1,L)*SS(2,K)+ SS(1,K)*SS(2,L)) + SS(5,K)*(-SS(3,L)*SS(4,K)+SS(3,K)*SS(4,L)*(SS(1,K)*SS(1,L) + SS(2,K)*SS(2,L)));
05551                 OY = SS(5,K)*SS(4,L)*(-SS(1,L)*SS(2,K)+ SS(1,K)*SS(2,L)) - SS(6,K)*(-SS(3,L)*SS(4,K)+SS(3,K)*SS(4,L)*(SS(1,K)*SS(1,L) + SS(2,K)*SS(2,L)));
05552 
05553                 if(OX != 0.0f || OY!=0.0f) {
05554                         //int count = 0;
05555                         for(int J=1;J<=NROW;J++) {
05556                                 JY = (J-1);
05557                                 if(JY > NR2) JY=JY-NROW;
05558                                 for(int I=1;I<=NNNN/2;I++) {
05559                                         Y =  fabs(OX * (I-1) + OY * JY);
05560                                         if(Y < 2.0f) W(I,J) += exptable[int(Y*indcnst)];//exp(-4*Y*Y);//
05561                                     //if(Y < 2.0f) Wptr[count++] += exp(-4*Y*Y);//exptable[int(Y*indcnst)];//
05562                                 }
05563                         }
05564                 } else {
05565                         for(int J=1;J<=NROW;J++) for(int I=1;I<=NNNN/2;I++)  W(I,J) += 1.0f;
05566                 }
05567         }
05568 
05569         EMData* proj_in = PROJ;
05570 
05571         PROJ = PROJ->norm_pad( false, 2);
05572         PROJ->do_fft_inplace();
05573         PROJ->update();
05574         PROJptr = PROJ->get_data();
05575 
05576         float WNRMinv,temp;
05577         float osnr = 1.0f/SNR;
05578         WNRMinv = 1/W(1,1);
05579         for(int J=1;J<=NROW;J++)
05580                 for(int I=1;I<=NNNN;I+=2) {
05581                         KX          = (I+1)/2;
05582                         temp    = W(KX,J)*WNRMinv;
05583                         WW      = temp/(temp*temp + osnr);
05584                         PROJ(I,J)       *= WW;
05585                         PROJ(I+1,J) *= WW;
05586                 }
05587         delete W; W = 0;
05588         PROJ->do_ift_inplace();
05589         PROJ->depad();
05590 
05591         float* data_src = PROJ->get_data();
05592         float* data_dst = proj_in->get_data();
05593 
05594         int ntotal = NSAM*NROW;
05595         for( int i=0; i < ntotal; ++i )
05596         {
05597             data_dst[i] = data_src[i];
05598         }
05599 
05600         proj_in->update();
05601 
05602         delete PROJ;
05603 }

Dict Util::CANG	(	float	PHI,
		float	THETA,
		float	PSI
	)			`[static]`

Definition at line 5418 of file util_sparx.cpp.

References DGR_TO_RAD, DM, and SS.

05419 {
05420         double CPHI,SPHI,CTHE,STHE,CPSI,SPSI;
05421         vector<float>   DM,SS;
05422 
05423         for(int i =0;i<9;i++) DM.push_back(0);
05424 
05425         for(int i =0;i<6;i++) SS.push_back(0);
05426 
05427         CPHI = cos(double(PHI)*DGR_TO_RAD);
05428         SPHI = sin(double(PHI)*DGR_TO_RAD);
05429         CTHE = cos(double(THETA)*DGR_TO_RAD);
05430         STHE = sin(double(THETA)*DGR_TO_RAD);
05431         CPSI = cos(double(PSI)*DGR_TO_RAD);
05432         SPSI = sin(double(PSI)*DGR_TO_RAD);
05433 
05434         SS(1) = float(CPHI);
05435         SS(2) = float(SPHI);
05436         SS(3) = float(CTHE);
05437         SS(4) = float(STHE);
05438         SS(5) = float(CPSI);
05439         SS(6) = float(SPSI);
05440 
05441         DM(1) = float(CPHI*CTHE*CPSI-SPHI*SPSI);
05442         DM(2) = float(SPHI*CTHE*CPSI+CPHI*SPSI);
05443         DM(3) = float(-STHE*CPSI);
05444         DM(4) = float(-CPHI*CTHE*SPSI-SPHI*CPSI);
05445         DM(5) = float(-SPHI*CTHE*SPSI+CPHI*CPSI);
05446         DM(6) = float(STHE*SPSI);
05447         DM(7) = float(STHE*CPHI);
05448         DM(8) = float(STHE*SPHI);
05449         DM(9) = float(CTHE);
05450 
05451         Dict DMnSS;
05452         DMnSS["DM"] = DM;
05453         DMnSS["SS"] = SS;
05454 
05455         return(DMnSS);
05456 }

void Util::BPCQ	(	EMData *	B,
		EMData *	CUBE,
		vector< float >	DM
	)			`[static]`

Definition at line 5466 of file util_sparx.cpp.

References B, CUBE, DM, EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::Transform::get_params(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), and t.

05467 {
05468 
05469         float  *Bptr = B->get_data();
05470         float  *CUBEptr = CUBE->get_data();
05471 
05472         int NSAM,NROW,NX3D,NY3D,NZC,KZ,IQX,IQY,LDPX,LDPY,LDPZ,LDPNMX,LDPNMY,NZ1;
05473         float DIPX,DIPY,XB,YB,XBB,YBB;
05474 
05475         Transform * t = B->get_attr("xform.projection");
05476         Dict d = t->get_params("spider");
05477         if(t) {delete t; t=0;}
05478         //  Unsure about sign of shifts, check later PAP 06/28/09
05479         float x_shift = d[ "tx" ];
05480         float y_shift = d[ "ty" ];
05481         x_shift = -x_shift;
05482         y_shift = -y_shift;
05483 
05484         NSAM = B->get_xsize();
05485         NROW = B->get_ysize();
05486         NX3D = CUBE->get_xsize();
05487         NY3D = CUBE->get_ysize();
05488         NZC  = CUBE->get_zsize();
05489 
05490 
05491         LDPX   = NX3D/2 +1;
05492         LDPY   = NY3D/2 +1;
05493         LDPZ   = NZC/2 +1;
05494         LDPNMX = NSAM/2 +1;
05495         LDPNMY = NROW/2 +1;
05496         NZ1    = 1;
05497 
05498         for(int K=1;K<=NZC;K++) {
05499                 KZ=K-1+NZ1;
05500                 for(int J=1;J<=NY3D;J++) {
05501                         XBB = (1-LDPX)*DM(1)+(J-LDPY)*DM(2)+(KZ-LDPZ)*DM(3);
05502                         YBB = (1-LDPX)*DM(4)+(J-LDPY)*DM(5)+(KZ-LDPZ)*DM(6);
05503                         for(int I=1;I<=NX3D;I++) {
05504                                 XB  = (I-1)*DM(1)+XBB-x_shift;
05505                                 IQX = int(XB+float(LDPNMX));
05506                                 if (IQX <1 || IQX >= NSAM) continue;
05507                                 YB  = (I-1)*DM(4)+YBB-y_shift;
05508                                 IQY = int(YB+float(LDPNMY));
05509                                 if (IQY<1 || IQY>=NROW)  continue;
05510                                 DIPX = XB+LDPNMX-IQX;
05511                                 DIPY = YB+LDPNMY-IQY;
05512 
05513                                 CUBE(I,J,K) = CUBE(I,J,K)+B(IQX,IQY)+DIPY*(B(IQX,IQY+1)-B(IQX,IQY))+DIPX*(B(IQX+1,IQY)-B(IQX,IQY)+DIPY*(B(IQX+1,IQY+1)-B(IQX+1,IQY)-B(IQX,IQY+1)+B(IQX,IQY)));
05514                         }
05515                 }
05516         }
05517 }

vector< float > Util::infomask	(	EMData *	Vol,
		EMData *	mask,
		bool	flip = `false`
	)			`[static]`

Definition at line 62 of file util_sparx.cpp.

References ENTERFUNC, EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, ImageFormatException, LOGERR, nx, ny, and sqrt().

00063              :  find statistics under the mask (mask >0.5)
00064 //  flip false: find statistics ourside the mask (mask <0.5)
00065 {
00066         ENTERFUNC;
00067         vector<float> stats;
00068         float *Volptr, *maskptr,MAX,MIN;
00069         long double Sum1,Sum2;
00070         long count;
00071 
00072         MAX = -FLT_MAX;
00073         MIN =  FLT_MAX;
00074         count = 0L;
00075         Sum1 = 0.0L;
00076         Sum2 = 0.0L;
00077 
00078         if (mask == NULL) {
00079            //Vol->update_stat();
00080            stats.push_back(Vol->get_attr("mean"));
00081            stats.push_back(Vol->get_attr("sigma"));
00082            stats.push_back(Vol->get_attr("minimum"));
00083            stats.push_back(Vol->get_attr("maximum"));
00084            return stats;
00085         }
00086 
00087         /* Check if the sizes of the mask and image are same */
00088 
00089         size_t nx = Vol->get_xsize();
00090         size_t ny = Vol->get_ysize();
00091         size_t nz = Vol->get_zsize();
00092 
00093         size_t mask_nx = mask->get_xsize();
00094         size_t mask_ny = mask->get_ysize();
00095         size_t mask_nz = mask->get_zsize();
00096 
00097         if  (nx != mask_nx || ny != mask_ny || nz != mask_nz )
00098                 throw ImageDimensionException("The dimension of the image does not match the dimension of the mask!");
00099 
00100  /*       if (nx != mask_nx ||
00101             ny != mask_ny ||
00102             nz != mask_nz  ) {
00103            // should throw an exception here!!! (will clean it up later CY)
00104            fprintf(stderr, "The dimension of the image does not match the dimension of the mask!\n");
00105            fprintf(stderr, " nx = %d, mask_nx = %d\n", nx, mask_nx);
00106            fprintf(stderr, " ny = %d, mask_ny = %d\n", ny, mask_ny);
00107            fprintf(stderr, " nz = %d, mask_nz = %d\n", nz, mask_nz);
00108            exit(1);
00109         }
00110  */
00111         Volptr = Vol->get_data();
00112         maskptr = mask->get_data();
00113 
00114         for (size_t i = 0; i < nx*ny*nz; i++) {
00115               if (maskptr[i]>0.5f == flip) {
00116                 Sum1 += Volptr[i];
00117                 Sum2 += Volptr[i]*Volptr[i];
00118                 MAX = (MAX < Volptr[i])?Volptr[i]:MAX;
00119                 MIN = (MIN > Volptr[i])?Volptr[i]:MIN;
00120                 count++;
00121               }
00122         }
00123 
00124         if (count == 0) {
00125                 LOGERR("Invalid mask");
00126                 throw ImageFormatException( "Invalid mask");
00127         }
00128 
00129        float avg = static_cast<float>(Sum1/count);
00130        float sig2 = static_cast<float>(Sum2 - count*avg*avg)/(count-1);
00131        float sig = sqrt(sig2);
00132 
00133        stats.push_back(avg);
00134        stats.push_back(sig);
00135        stats.push_back(MIN);
00136        stats.push_back(MAX);
00137 
00138        return stats;
00139 }

void Util::colreverse	(	float *	beg,
		float *	end,
		int	nx
	)			`[static]`

Definition at line 5132 of file util_sparx.cpp.

Referenced by cyclicshift(), and slicereverse().

05132                                                     {
05133         float* tmp = new float[nx];
05134         int n = (end - beg)/nx;
05135         int nhalf = n/2;
05136         for (int i = 0; i < nhalf; i++) {
05137                 // swap col i and col n-1-i
05138                 memcpy(tmp, beg+i*nx, nx*sizeof(float));
05139                 memcpy(beg+i*nx, beg+(n-1-i)*nx, nx*sizeof(float));
05140                 memcpy(beg+(n-1-i)*nx, tmp, nx*sizeof(float));
05141         }
05142         delete[] tmp;
05143 }

void Util::slicereverse	(	float *	beg,
		float *	end,
		int	nx,
		int	ny
	)			`[static]`

Definition at line 5145 of file util_sparx.cpp.

References colreverse().

Referenced by cyclicshift().

05146 {
05147         int nxy = nx*ny;
05148         colreverse(beg, end, nxy);
05149 }

void Util::cyclicshift	(	EMData *	image,
		Dict	params
	)			`[static]`

Definition at line 5152 of file util_sparx.cpp.

References colreverse(), data, EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageFormatException, EMAN::EMData::is_complex(), nx, ny, slicereverse(), and EMAN::EMData::update().

05152                                                  {
05153 /*
05154  Performs inplace integer cyclic shift as specified by the "dx","dy","dz" parameters on a 3d volume.
05155  Implements the inplace swapping using reversals as descibed in  also:
05156     http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Intro/Eg01/
05157 
05158 
05159 * @author  Phani Ivatury
05160 * @date 18-2006
05161 * @see http://www.csse.monash.edu.au/~lloyd/tildeAlgDS/Intro/Eg01/
05162 *
05163 *
05164 * A[0] A[1] A[2] A[3] A[4] A[5] A[6] A[7] A[8] A[9]
05165 *
05166 * 10   20   30   40   50   60   70   80   90   100
05167 * ------------
05168 *    m = 3 (shift left three places)
05169 *
05170 * Reverse the items from 0..m-1 and m..N-1:
05171 *
05172 * 30   20   10   100  90   80   70   60   50   40
05173 *
05174 * Now reverse the entire sequence:
05175 *
05176 * 40   50   60   70   80   90   100  10   20   30
05177 
05178 
05179     cycl_shift() in libpy/fundementals.py calls this function
05180 
05181     Usage:
05182     EMData *im1 = new EMData();
05183     im1->set_size(70,80,85);
05184     im1->to_one();
05185     Dict params; params["dx"] = 10;params["dy"] = 10000;params["dz"] = -10;
05186     Utils::cyclicshift(im1,params);
05187     im1.peak_search(1,1)
05188 */
05189 
05190         if (image->is_complex()) throw ImageFormatException("Real image required for IntegerCyclicShift2DProcessor");
05191 
05192         int dx = params["dx"];
05193         int dy = params["dy"];
05194         int dz = params["dz"];
05195 
05196         // The reverse trick we're using shifts to the left (a negative shift)
05197         int nx = image->get_xsize();
05198         dx %= nx;
05199         if (dx < 0) dx += nx;
05200         int ny = image->get_ysize();
05201         dy %= ny;
05202         if (dy < 0) dy += ny;
05203         int nz = image->get_zsize();
05204         dz %= nz;
05205         if (dz < 0) dz += nz;
05206 
05207         int mx = -(dx - nx);
05208         int my = -(dy - ny);
05209         int mz = -(dz - nz);
05210 
05211         float* data = image->get_data();
05212         // x-reverses
05213         if (mx != 0) {
05214                 for (int iz = 0; iz < nz; iz++)
05215                        for (int iy = 0; iy < ny; iy++) {
05216                                 // reverses for column iy
05217                                 int offset = nx*iy + nx*ny*iz; // starting location for column iy in slice iz
05218                                 reverse(&data[offset],&data[offset+mx]);
05219                                 reverse(&data[offset+mx],&data[offset+nx]);
05220                                 reverse(&data[offset],&data[offset+nx]);
05221                         }
05222         }
05223         // y-reverses
05224         if (my != 0) {
05225                 for (int iz = 0; iz < nz; iz++) {
05226                         int offset = nx*ny*iz;
05227                         colreverse(&data[offset], &data[offset + my*nx], nx);
05228                         colreverse(&data[offset + my*nx], &data[offset + ny*nx], nx);
05229                         colreverse(&data[offset], &data[offset + ny*nx], nx);
05230                 }
05231         }
05232         if (mz != 0) {
05233                 slicereverse(&data[0], &data[mz*ny*nx], nx, ny);
05234                 slicereverse(&data[mz*ny*nx], &data[nz*ny*nx], nx, ny);
05235                 slicereverse(&data[0], &data[nz*ny*nx], nx ,ny);
05236         }
05237         image->update();
05238 }

Dict Util::im_diff	(	EMData *	V1,
		EMData *	V2,
		EMData *	mask = `0`
	)			`[static]`

Definition at line 144 of file util_sparx.cpp.

References B, ENTERFUNC, EXITFUNC, EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageFormatException, EMAN::EMUtil::is_same_size(), LOGERR, nx, ny, EMAN::EMData::set_size(), EMAN::EMData::to_one(), and EMAN::EMData::update().

00145 {
00146         ENTERFUNC;
00147 
00148         if (!EMUtil::is_same_size(V1, V2)) {
00149                 LOGERR("images not same size");
00150                 throw ImageFormatException( "images not same size");
00151         }
00152 
00153         size_t nx = V1->get_xsize();
00154         size_t ny = V1->get_ysize();
00155         size_t nz = V1->get_zsize();
00156         size_t size = nx*ny*nz;
00157 
00158         EMData *BD = new EMData();
00159         BD->set_size(nx, ny, nz);
00160 
00161         float *params = new float[2];
00162 
00163         float *V1ptr, *V2ptr, *MASKptr, *BDptr, A, B;
00164         long double S1=0.L,S2=0.L,S3=0.L,S4=0.L;
00165         int nvox = 0L;
00166 
00167         V1ptr = V1->get_data();
00168         V2ptr = V2->get_data();
00169         BDptr = BD->get_data();
00170 
00171 
00172         if(!mask){
00173                 EMData * Mask = new EMData();
00174                 Mask->set_size(nx,ny,nz);
00175                 Mask->to_one();
00176                 MASKptr = Mask->get_data();
00177         } else {
00178                 if (!EMUtil::is_same_size(V1, mask)) {
00179                         LOGERR("mask not same size");
00180                         throw ImageFormatException( "mask not same size");
00181                 }
00182 
00183                 MASKptr = mask->get_data();
00184         }
00185 
00186 
00187 
00188 //       calculation of S1,S2,S3,S3,nvox
00189 
00190         for (size_t i = 0L;i < size; i++) {
00191               if (MASKptr[i]>0.5f) {
00192                S1 += V1ptr[i]*V2ptr[i];
00193                S2 += V1ptr[i]*V1ptr[i];
00194                S3 += V2ptr[i];
00195                S4 += V1ptr[i];
00196                nvox ++;
00197               }
00198         }
00199 
00200         if ((nvox*S1 - S3*S4) == 0. || (nvox*S2 - S4*S4) == 0) {
00201                 A =1.0f ;
00202         } else {
00203                 A = static_cast<float>( (nvox*S1 - S3*S4)/(nvox*S2 - S4*S4) );
00204         }
00205         B = static_cast<float> (A*S4  -  S3)/nvox;
00206 
00207         // calculation of the difference image
00208 
00209         for (size_t i = 0L;i < size; i++) {
00210              if (MASKptr[i]>0.5f) {
00211                BDptr[i] = A*V1ptr[i] -  B  - V2ptr[i];
00212              }  else  {
00213                BDptr[i] = 0.f;
00214              }
00215         }
00216 
00217         BD->update();
00218 
00219         params[0] = A;
00220         params[1] = B;
00221 
00222         Dict BDnParams;
00223         BDnParams["imdiff"] = BD;
00224         BDnParams["A"] = params[0];
00225         BDnParams["B"] = params[1];
00226 
00227         EXITFUNC;
00228         return BDnParams;
00229  }

EMData * Util::TwoDTestFunc	(	int	Size,
		float	p,
		float	q,
		float	a,
		float	b,
		int	flag = `0`,
		float	alphaDeg = `0`
	)			`[static]`

Creates a Two D Test Pattern.

Parameters:

`[in]`	Size	must be odd
`[in]`	p	the x frequency
`[in]`	q	the y frequency
`[in]`	a	the x falloff
	b	the y falloff
	flag
`[in]`	alphaDeg	the projection angle

Returns:: The 2D test pattern in real space, fourier space, or the projection in real or fourier space or the FH of the pattern

Definition at line 235 of file util_sparx.cpp.

References C, ENTERFUNC, EXITFUNC, pi, real2FH(), set_complex(), EMAN::EMData::set_complex(), EMAN::EMData::set_fftodd(), set_ri(), EMAN::EMData::set_ri(), EMAN::EMData::set_shuffled(), EMAN::EMData::set_size(), sqrt(), EMAN::EMData::to_zero(), update(), EMAN::EMData::update(), x, and y.

00236 {
00237         ENTERFUNC;
00238         int Mid= (Size+1)/2;
00239 
00240         if (flag==0) { // This is the real function
00241                 EMData* ImBW = new EMData();
00242                 ImBW->set_size(Size,Size,1);
00243                 ImBW->to_zero();
00244 
00245                 float tempIm;
00246                 float x,y;
00247 
00248                 for (int ix=(1-Mid);  ix<Mid; ix++){
00249                         for (int iy=(1-Mid);  iy<Mid; iy++){
00250                                 x = (float)ix;
00251                                 y = (float)iy;
00252                         tempIm= static_cast<float>( (1/(2*M_PI)) * cos(p*x)* cos(q*y) * exp(-.5*x*x/(a*a))* exp(-.5*y*y/(b*b)) );
00253                                 (*ImBW)(ix+Mid-1,iy+Mid-1) = tempIm * exp(.5f*p*p*a*a)* exp(.5f*q*q*b*b);
00254                         }
00255                 }
00256                 ImBW->update();
00257                 ImBW->set_complex(false);
00258                 ImBW->set_ri(true);
00259 
00260                 return ImBW;
00261         }
00262         else if (flag==1) {  // This is the Fourier Transform
00263                 EMData* ImBWFFT = new EMData();
00264                 ImBWFFT ->set_size(2*Size,Size,1);
00265                 ImBWFFT ->to_zero();
00266 
00267                 float r,s;
00268 
00269                 for (int ir=(1-Mid);  ir<Mid; ir++){
00270                         for (int is=(1-Mid);  is<Mid; is++){
00271                                 r = (float)ir;
00272                                 s = (float)is;
00273                         (*ImBWFFT)(2*(ir+Mid-1),is+Mid-1)= cosh(p*r*a*a) * cosh(q*s*b*b) *
00274                                 exp(-.5f*r*r*a*a)* exp(-.5f*s*s*b*b);
00275                         }
00276                 }
00277                 ImBWFFT->update();
00278                 ImBWFFT->set_complex(true);
00279                 ImBWFFT->set_ri(true);
00280                 ImBWFFT->set_shuffled(true);
00281                 ImBWFFT->set_fftodd(true);
00282 
00283                 return ImBWFFT;
00284         }
00285         else if (flag==2 || flag==3) { //   This is the projection in Real Space
00286                 float alpha = static_cast<float>( alphaDeg*M_PI/180.0 );
00287                 float C=cos(alpha);
00288                 float S=sin(alpha);
00289                 float D= sqrt(S*S*b*b + C*C*a*a);
00290                 //float D2 = D*D;   PAP - to get rid of warning
00291 
00292                 float P = p * C *a*a/D ;
00293                 float Q = q * S *b*b/D ;
00294 
00295                 if (flag==2) {
00296                         EMData* pofalpha = new EMData();
00297                         pofalpha ->set_size(Size,1,1);
00298                         pofalpha ->to_zero();
00299 
00300                         float Norm0 =  D*(float)sqrt(2*pi);
00301                         float Norm1 =  exp( .5f*(P+Q)*(P+Q)) / Norm0 ;
00302                         float Norm2 =  exp( .5f*(P-Q)*(P-Q)) / Norm0 ;
00303                         float sD;
00304 
00305                         for (int is=(1-Mid);  is<Mid; is++){
00306                                 sD = is/D ;
00307                                 (*pofalpha)(is+Mid-1) =  Norm1 * exp(-.5f*sD*sD)*cos(sD*(P+Q))
00308                          + Norm2 * exp(-.5f*sD*sD)*cos(sD*(P-Q));
00309                         }
00310                         pofalpha-> update();
00311                         pofalpha-> set_complex(false);
00312                         pofalpha-> set_ri(true);
00313 
00314                         return pofalpha;
00315                 }
00316                 if (flag==3) { // This is the projection in Fourier Space
00317                         float vD;
00318 
00319                         EMData* pofalphak = new EMData();
00320                         pofalphak ->set_size(2*Size,1,1);
00321                         pofalphak ->to_zero();
00322 
00323                         for (int iv=(1-Mid);  iv<Mid; iv++){
00324                                 vD = iv*D ;
00325                                 (*pofalphak)(2*(iv+Mid-1)) =  exp(-.5f*vD*vD)*(cosh(vD*(P+Q)) + cosh(vD*(P-Q)) );
00326                         }
00327                         pofalphak-> update();
00328                         pofalphak-> set_complex(false);
00329                         pofalphak-> set_ri(true);
00330 
00331                         return pofalphak;
00332                 }
00333         }
00334         else if (flag==4) {
00335                 cout <<" FH under construction";
00336                 EMData* OutFT= TwoDTestFunc(Size, p, q, a, b, 1);
00337                 EMData* TryFH= OutFT -> real2FH(4.0);
00338                 return TryFH;
00339         } else {
00340                 cout <<" flag must be 0,1,2,3, or 4";
00341         }
00342 
00343         EXITFUNC;
00344         return 0;
00345 }

void Util::spline_mat	(	float *	x,
		float *	y,
		int	n,
		float *	xq,
		float *	yq,
		int	m
	)			`[static]`

Given a tabulated function y of x (n unordered points), and Given the values of the m values xq to be interpolated This routine returns the interpolated array yq, PRB This function is called by splint.

Parameters:

	x
`[in]`	y	of x is the tabulated function of length n
	n
`[in]`	xq	is the x values to be splined: has m points.
`[out]`	yq	are the splined values
	m

Definition at line 348 of file util_sparx.cpp.

References spline(), and splint().

00349 {
00350 
00351         float x0= x[0];
00352         float x1= x[1];
00353         float x2= x[2];
00354         float y0= y[0];
00355         float y1= y[1];
00356         float y2= y[2];
00357         float yp1 =  (y1-y0)/(x1-x0) +  (y2-y0)/(x2-x0) - (y2-y1)/(x2-x1)  ;
00358         float xn  = x[n];
00359         float xnm1= x[n-1];
00360         float xnm2= x[n-2];
00361         float yn  = y[n];
00362         float ynm1= y[n-1];
00363         float ynm2= y[n-2];
00364         float ypn=  (yn-ynm1)/(xn-xnm1) +  (yn-ynm2)/(xn-xnm2) - (ynm1-ynm2)/(xnm1-xnm2) ;
00365         float *y2d = new float[n];
00366         Util::spline(x,y,n,yp1,ypn,y2d);
00367         Util::splint(x,y,y2d,n,xq,yq,m); //PRB
00368         delete [] y2d;
00369         return;
00370 }

void Util::spline	(	float *	x,
		float *	y,
		int	n,
		float	yp1,
		float	ypn,
		float *	y2
	)			`[static]`

Given a tabulated function y of x (unordered), and Given the values of the first derivatives at the end points This routine returns an array y2, that contains the second derivatives of the function at the tabulated points.

PRB This function is called by splint

Parameters:

	x
`[in]`	y	of x is the tabulated function of length n
	n
`[in]`	yp1	: the derivatives of the first point.
`[in]`	ypn	: the derivatives of the last point.
`[out]`	y2	is the value of the second derivatives

Definition at line 373 of file util_sparx.cpp.

00374 {
00375         int i,k;
00376         float p, qn, sig, un, *u;
00377         u = new float[n-1];
00378 
00379         if (yp1 > .99e30){
00380                 y2[0]=u[0]=0.0;
00381         } else {
00382                 y2[0]=-.5f;
00383                 u[0] =(3.0f/ (x[1] -x[0]))*( (y[1]-y[0])/(x[1]-x[0]) -yp1);
00384         }
00385 
00386         for (i=1; i < n-1; i++) {
00387                 sig= (x[i] - x[i-1])/(x[i+1] - x[i-1]);
00388                 p = sig*y2[i-1] + 2.0f;
00389                 y2[i]  = (sig-1.0f)/p;
00390                 u[i] = (y[i+1] - y[i] )/(x[i+1]-x[i] ) -  (y[i] - y[i-1] )/(x[i] -x[i-1]);
00391                 u[i] = (6.0f*u[i]/ (x[i+1]-x[i-1]) - sig*u[i-1])/p;
00392         }
00393 
00394         if (ypn>.99e30){
00395                 qn=0; un=0;
00396         } else {
00397                 qn= .5f;
00398                 un= (3.0f/(x[n-1] -x[n-2])) * (ypn -  (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
00399         }
00400         y2[n-1]= (un - qn*u[n-2])/(qn*y2[n-2]+1.0f);
00401         for (k=n-2; k>=0; k--){
00402                 y2[k]=y2[k]*y2[k+1]+u[k];
00403         }
00404         delete [] u;
00405 }

void Util::splint	(	float *	xa,
		float *	ya,
		float *	y2a,
		int	n,
		float *	xq,
		float *	yq,
		int	m
	)			`[static]`

Given the arrays xa(ordered, ya of length n, which tabulate a function and given the array y2a which is the output of spline and an unordered array xq, this routine returns a cubic-spline interpolated array yq.

Parameters:

	xa
`[in]`	ya	of x is the tabulated function of length n
`[in]`	y2a	is returned from spline: second derivs
	n
`[in]`	xq	is the x values to be splined: has m points.
`[out]`	yq	are the splined values
	m

Definition at line 408 of file util_sparx.cpp.

References b.

00409 {
00410         int klo, khi, k;
00411         float h, b, a;
00412 
00413 //      klo=0; // can try to put here
00414         for (int j=0; j<m;j++){
00415                 klo=0;
00416                 khi=n-1;
00417                 while (khi-klo >1) {
00418                         k=(khi+klo) >>1;
00419                         if  (xa[k]>xq[j]){ khi=k;}
00420                         else { klo=k;}
00421                 }
00422                 h=xa[khi]- xa[klo];
00423                 if (h==0.0) printf("Bad XA input to routine SPLINT \n");
00424                 a =(xa[khi]-xq[j])/h;
00425                 b=(xq[j]-xa[klo])/h;
00426                 yq[j]=a*ya[klo] + b*ya[khi]
00427                         + ((a*a*a-a)*y2a[klo]
00428                              +(b*b*b-b)*y2a[khi]) *(h*h)/6.0f;
00429         }
00430 //      printf("h=%f, a = %f, b=%f, ya[klo]=%f, ya[khi]=%f , yq=%f\n",h, a, b, ya[klo], ya[khi],yq[0]);
00431 }

void Util::Radialize	(	int *	PermMatTr,
		float *	kValsSorted,
		float *	weightofkvalsSorted,
		int	Size,
		int *	SizeReturned
	)			`[static]`

list the sorted lengths of the integer lattice sites of a square sided image of size Size.

PRB

Parameters:

`[out]`	PermMatTr	The matrix telling the ordering of the lattice sites wrt the array
`[out]`	kValsSorted
`[out]`	weightofkvalsSorted	the number of sites at that distance
`[in]`	Size	The length of the image
`[out]`	SizeReturned

Definition at line 434 of file util_sparx.cpp.

References sort_mat().

00436 {
00437         int iMax = (int) floor( (Size-1.0)/2 +.01);
00438         int CountMax = (iMax+2)*(iMax+1)/2;
00439         int Count=-1;
00440         float *kVals     = new float[CountMax];
00441         float *weightMat = new float[CountMax];
00442         int *PermMat     = new   int[CountMax];
00443         SizeReturned[0] = CountMax;
00444 
00445 //      printf("Aa \n");        fflush(stdout);
00446         for (int jkx=0; jkx< iMax+1; jkx++) {
00447                 for (int jky=0; jky< jkx+1; jky++) {
00448                         Count++;
00449                         kVals[Count] = sqrtf((float) (jkx*jkx +jky*jky));
00450                         weightMat[Count]=  1.0;
00451                         if (jkx!=0)  { weightMat[Count] *=2;}
00452                         if (jky!=0)  { weightMat[Count] *=2;}
00453                         if (jkx!=jky){ weightMat[Count] *=2;}
00454                         PermMat[Count]=Count+1;
00455                 }
00456         }
00457 
00458         int lkVals = Count+1;
00459 //      printf("Cc \n");fflush(stdout);
00460 
00461         sort_mat(&kVals[0],&kVals[Count],
00462              &PermMat[0],  &PermMat[Count]);  //PermMat is
00463                                 //also returned as well as kValsSorted
00464         fflush(stdout);
00465 
00466         int newInd;
00467 
00468         for (int iP=0; iP < lkVals ; iP++ ) {
00469                 newInd =  PermMat[iP];
00470                 PermMatTr[newInd-1] = iP+1;
00471         }
00472 
00473 //      printf("Ee \n"); fflush(stdout);
00474 
00475         int CountA=-1;
00476         int CountB=-1;
00477 
00478         while (CountB< (CountMax-1)) {
00479                 CountA++;
00480                 CountB++;
00481 //              printf("CountA=%d ; CountB=%d \n", CountA,CountB);fflush(stdout);
00482                 kValsSorted[CountA] = kVals[CountB] ;
00483                 if (CountB<(CountMax-1) ) {
00484                         while (fabs(kVals[CountB] -kVals[CountB+1])<.0000001  ) {
00485                                 SizeReturned[0]--;
00486                                 for (int iP=0; iP < lkVals; iP++){
00487 //                                      printf("iP=%d \n", iP);fflush(stdout);
00488                                         if  (PermMatTr[iP]>CountA+1) {
00489                                                 PermMatTr[iP]--;
00490                                         }
00491                                 }
00492                                 CountB++;
00493                         }
00494                 }
00495         }
00496 
00497 
00498         for (int CountD=0; CountD < CountMax; CountD++) {
00499             newInd = PermMatTr[CountD];
00500             weightofkValsSorted[newInd-1] += weightMat[CountD];
00501         }
00502 
00503 }

vector< float > Util::even_angles	(	float	delta,
		float	t1 = `0`,
		float	t2 = `90`,
		float	p1 = `0`,
		float	p2 = `359.999`
	)			`[static]`

Compute a vector containing quasi-evenly spaced Euler angles.

The order of angles in the vector is phi, theta, psi.

Parameters:

`[in]`	delta	Delta theta (spacing in theta).
`[in]`	t1	Starting (min) value of theta in degrees, default = 0.
`[in]`	t2	Ending (max) value of theta in degrees, default = 90.
`[in]`	p1	Starting (min) value of phi in degrees, default = 0.
`[in]`	p2	Ending (max) value of phi in degrees, default = 359.9.

Returns:: Vector of angles as a flat list of phi_0, theta_0, psi_0, ..., phi_N, theta_N, psi_N.

Definition at line 507 of file util_sparx.cpp.

References angles, dgr_to_rad, phi, and theta.

00508 {
00509         vector<float> angles;
00510         float psi = 0.0;
00511         if ((0.0 == t1)&&(0.0 == t2)||(t1 >= t2)) {
00512                 t1 = 0.0f;
00513                 t2 = 90.0f;
00514         }
00515         if ((0.0 == p1)&&(0.0 == p2)||(p1 >= p2)) {
00516                 p1 = 0.0f;
00517                 p2 = 359.9f;
00518         }
00519         bool skip = ((t1 < 90.0)&&(90.0 == t2)&&(0.0 == p1)&&(p2 > 180.0));
00520         for (float theta = t1; theta <= t2; theta += delta) {
00521                 float detphi;
00522                 int lt;
00523                 if ((0.0 == theta)||(180.0 == theta)) {
00524                         detphi = 360.0f;
00525                         lt = 1;
00526                 } else {
00527                         detphi = delta/sin(theta*static_cast<float>(dgr_to_rad));
00528                         lt = int((p2 - p1)/detphi)-1;
00529                         if (lt < 1) lt = 1;
00530                         detphi = (p2 - p1)/lt;
00531                 }
00532                 for (int i = 0; i < lt; i++) {
00533                         float phi = p1 + i*detphi;
00534                         if (skip&&(90.0 == theta)&&(phi > 180.0)) continue;
00535                         angles.push_back(phi);
00536                         angles.push_back(theta);
00537                         angles.push_back(psi);
00538                 }
00539         }
00540         return angles;
00541 }

float Util::quadri	(	float	x,
		float	y,
		int	nx,
		int	ny,
		float *	image
	)			`[static]`

Quadratic interpolation (2D).

Note: This routine starts counting from 1, not 0!

This routine uses six image points for interpolation:

See also:: M. Abramowitz & I.E. Stegun, Handbook of Mathematical Functions (Dover, New York, 1964), Sec. 25.2.67. http://www.math.sfu.ca/~cbm/aands/page_882.htm
http://www.cl.cam.ac.uk/users/nad/pubs/quad.pdf

                f3    fc
                |
                | x
         f2-----f0----f1
                |
                |
                f4
		 *

f0 - f4 are image values near the interpolated point X. f0 is the interior mesh point nearest x.

Coords:

f0 = (x0, y0)
f1 = (xb, y0)
f2 = (xa, y0)
f3 = (x0, yb)
f4 = (x0, ya)
fc = (xc, yc)

Mesh spacings:

hxa -- x- mesh spacing to the left of f0
hxb -- x- mesh spacing to the right of f0
hyb -- y- mesh spacing above f0
hya -- y- mesh spacing below f0

Interpolant: f = f0 + c1*(x-x0) + c2*(x-x0)*(x-x1) + c3*(y-y0) + c4*(y-y0)*(y-y1) + c5*(x-x0)*(y-y0)

Parameters:

`[in]`	x	x-coord value
`[in]`	y	y-coord value
	nx
	ny
`[in]`	image	Image object (pointer)

Returns:: Interpolated value

Definition at line 646 of file util_sparx.cpp.

References fdata, x, and y.

Referenced by alrl_ms(), Polar2D(), Polar2Dm(), quadri_background(), and EMAN::EMData::rot_scale_trans2D().

00647 {
00648 //  purpose: quadratic interpolation
00649 //  Optimized for speed, circular closer removed, checking of ranges removed
00650         float  x, y, dx0, dy0, f0, c1, c2, c3, c4, c5, dxb, dyb;
00651         float  quadri;
00652         int    i, j, ip1, im1, jp1, jm1, ic, jc, hxc, hyc;
00653 
00654         x = xx;
00655         y = yy;
00656 
00657         //     any xx and yy
00658         while ( x < 1.0 )                 x += nxdata;
00659         while ( x >= (float)(nxdata+1) )  x -= nxdata;
00660         while ( y < 1.0 )                 y += nydata;
00661         while ( y >= (float)(nydata+1) )  y -= nydata;
00662 
00663         i   = (int) x;
00664         j   = (int) y;
00665 
00666         dx0 = x - i;
00667         dy0 = y - j;
00668 
00669         ip1 = i + 1;
00670         im1 = i - 1;
00671         jp1 = j + 1;
00672         jm1 = j - 1;
00673 
00674         if (ip1 > nxdata) ip1 -= nxdata;
00675         if (im1 < 1)      im1 += nxdata;
00676         if (jp1 > nydata) jp1 -= nydata;
00677         if (jm1 < 1)      jm1 += nydata;
00678 
00679         f0  = fdata(i,j);
00680         c1  = fdata(ip1,j) - f0;
00681         c2  = (c1 - f0 + fdata(im1,j)) * 0.5f;
00682         c3  = fdata(i,jp1) - f0;
00683         c4  = (c3 - f0 + fdata(i,jm1)) * 0.5f;
00684 
00685         dxb = dx0 - 1;
00686         dyb = dy0 - 1;
00687 
00688         // hxc & hyc are either 1 or -1
00689         if (dx0 >= 0) hxc = 1; else hxc = -1;
00690         if (dy0 >= 0) hyc = 1; else hyc = -1;
00691 
00692         ic  = i + hxc;
00693         jc  = j + hyc;
00694 
00695         if (ic > nxdata) ic -= nxdata;  else if (ic < 1) ic += nxdata;
00696         if (jc > nydata) jc -= nydata;  else if (jc < 1) jc += nydata;
00697 
00698         c5  =  ( (fdata(ic,jc) - f0 - hxc * c1 - (hxc * (hxc - 1.0f)) * c2
00699                 - hyc * c3 - (hyc * (hyc - 1.0f)) * c4) * (hxc * hyc));
00700 
00701 
00702         quadri = f0 + dx0 * (c1 + dxb * c2 + dy0 * c5) + dy0 * (c3 + dyb * c4);
00703 
00704         return quadri;
00705 }

float Util::quadri_background	(	float	x,
		float	y,
		int	nx,
		int	ny,
		float *	image,
		int	xnew,
		int	ynew
	)			`[static]`

Quadratic interpolation (2D).

This is identical to quadri except the wrap around is not done circulantly.

Note: This routine starts counting from 1, not 0!

This routine uses six image points for interpolation:

See also:: M. Abramowitz & I.E. Stegun, Handbook of Mathematical Functions (Dover, New York, 1964), Sec. 25.2.67. http://www.math.sfu.ca/~cbm/aands/page_882.htm
http://www.cl.cam.ac.uk/users/nad/pubs/quad.pdf

                f3    fc
                |
                | x
         f2-----f0----f1
                |
                |
                f4
		 *

f0 - f4 are image values near the interpolated point X. f0 is the interior mesh point nearest x.

Coords:

f0 = (x0, y0)
f1 = (xb, y0)
f2 = (xa, y0)
f3 = (x0, yb)
f4 = (x0, ya)
fc = (xc, yc)

Mesh spacings:

hxa -- x- mesh spacing to the left of f0
hxb -- x- mesh spacing to the right of f0
hyb -- y- mesh spacing above f0
hya -- y- mesh spacing below f0

Interpolant: f = f0 + c1*(x-x0) + c2*(x-x0)*(x-x1) + c3*(y-y0) + c4*(y-y0)*(y-y1) + c5*(x-x0)*(y-y0)

Parameters:

`[in]`	x	x-coord value
`[in]`	y	y-coord value
	nx
	ny
`[in]`	image	Image object (pointer)

Returns:: Interpolated value

Definition at line 710 of file util_sparx.cpp.

References fdata, quadri(), x, and y.

Referenced by EMAN::EMData::rot_scale_trans2D_background().

00711 {
00712 //  purpose: quadratic interpolation
00713 //  Optimized for speed, circular closer removed, checking of ranges removed
00714         float  x, y, dx0, dy0, f0, c1, c2, c3, c4, c5, dxb, dyb;
00715         float  quadri;
00716         int    i, j, ip1, im1, jp1, jm1, ic, jc, hxc, hyc;
00717 
00718         x = xx;
00719         y = yy;
00720 
00721         // wrap around is not done circulantly; if (x,y) is not in the image, then x = xnew and y = ynew
00722         if ( (x < 1.0) || ( x >= (float)(nxdata+1) ) || ( y < 1.0 ) || ( y >= (float)(nydata+1) )){
00723               x = xnew;
00724                   y = ynew;
00725      }
00726             
00727 
00728         i   = (int) x;
00729         j   = (int) y;
00730 
00731         dx0 = x - i;
00732         dy0 = y - j;
00733 
00734         ip1 = i + 1;
00735         im1 = i - 1;
00736         jp1 = j + 1;
00737         jm1 = j - 1;
00738 
00739         if (ip1 > nxdata) ip1 -= nxdata;
00740         if (im1 < 1)      im1 += nxdata;
00741         if (jp1 > nydata) jp1 -= nydata;
00742         if (jm1 < 1)      jm1 += nydata;
00743 
00744         f0  = fdata(i,j);
00745         c1  = fdata(ip1,j) - f0;
00746         c2  = (c1 - f0 + fdata(im1,j)) * 0.5f;
00747         c3  = fdata(i,jp1) - f0;
00748         c4  = (c3 - f0 + fdata(i,jm1)) * 0.5f;
00749 
00750         dxb = dx0 - 1;
00751         dyb = dy0 - 1;
00752 
00753         // hxc & hyc are either 1 or -1
00754         if (dx0 >= 0) hxc = 1; else hxc = -1;
00755         if (dy0 >= 0) hyc = 1; else hyc = -1;
00756 
00757         ic  = i + hxc;
00758         jc  = j + hyc;
00759 
00760         if (ic > nxdata) ic -= nxdata;  else if (ic < 1) ic += nxdata;
00761         if (jc > nydata) jc -= nydata;  else if (jc < 1) jc += nydata;
00762 
00763         c5  =  ( (fdata(ic,jc) - f0 - hxc * c1 - (hxc * (hxc - 1.0f)) * c2
00764                 - hyc * c3 - (hyc * (hyc - 1.0f)) * c4) * (hxc * hyc));
00765 
00766 
00767         quadri = f0 + dx0 * (c1 + dxb * c2 + dy0 * c5) + dy0 * (c3 + dyb * c4);
00768 
00769         return quadri;
00770 }

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

Definition at line 775 of file util_sparx.cpp.

References EMAN::Util::KaiserBessel::get_window_size(), EMAN::Util::KaiserBessel::i0win_tab(), restrict1(), and round().

Referenced by Polar2Dmi(), and EMAN::EMData::rot_scale_conv_new().

00775                                                                                                                              {
00776 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
00777 
00778 /* Commented by Zhengfan Yang on 04/20/07
00779 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
00780 I made the following changes to get_pixel_conv():
00781 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
00782 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
00783 3. Unfold the 'for' loop
00784 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)
00785 
00786 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
00787 size, say N=5, you can easily modify it by referring my code.
00788 */
00789         int K = kb.get_window_size();
00790         int kbmin = -K/2;
00791         int kbmax = -kbmin;
00792         int kbc = kbmax+1;
00793 
00794         float pixel =0.0f;
00795         float w=0.0f;
00796 
00797         delx = restrict1(delx, nx);
00798         int inxold = int(round(delx));
00799         if ( ny < 2 ) {  //1D
00800                 float tablex1 = kb.i0win_tab(delx-inxold+3);
00801                 float tablex2 = kb.i0win_tab(delx-inxold+2);
00802                 float tablex3 = kb.i0win_tab(delx-inxold+1);
00803                 float tablex4 = kb.i0win_tab(delx-inxold);
00804                 float tablex5 = kb.i0win_tab(delx-inxold-1);
00805                 float tablex6 = kb.i0win_tab(delx-inxold-2);
00806                 float tablex7 = kb.i0win_tab(delx-inxold-3);
00807 
00808                 int x1, x2, x3, x4, x5, x6, x7;
00809 
00810                 if ( inxold <= kbc || inxold >=nx-kbc-2 )  {
00811                         x1 = (inxold-3+nx)%nx;
00812                         x2 = (inxold-2+nx)%nx;
00813                         x3 = (inxold-1+nx)%nx;
00814                         x4 = (inxold  +nx)%nx;
00815                         x5 = (inxold+1+nx)%nx;
00816                         x6 = (inxold+2+nx)%nx;
00817                         x7 = (inxold+3+nx)%nx;
00818                 } else {
00819                         x1 = inxold-3;
00820                         x2 = inxold-2;
00821                         x3 = inxold-1;
00822                         x4 = inxold;
00823                         x5 = inxold+1;
00824                         x6 = inxold+2;
00825                         x7 = inxold+3;
00826                 }
00827 
00828                 pixel = data[x1]*tablex1 + data[x2]*tablex2 + data[x3]*tablex3 +
00829                         data[x4]*tablex4 + data[x5]*tablex5 + data[x6]*tablex6 +
00830                         data[x7]*tablex7 ;
00831 
00832                 w = tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7;
00833         } else if ( nz < 2 ) {  // 2D
00834                 dely = restrict1(dely, ny);
00835                 int inyold = int(round(dely));
00836                 float tablex1 = kb.i0win_tab(delx-inxold+3);
00837                 float tablex2 = kb.i0win_tab(delx-inxold+2);
00838                 float tablex3 = kb.i0win_tab(delx-inxold+1);
00839                 float tablex4 = kb.i0win_tab(delx-inxold);
00840                 float tablex5 = kb.i0win_tab(delx-inxold-1);
00841                 float tablex6 = kb.i0win_tab(delx-inxold-2);
00842                 float tablex7 = kb.i0win_tab(delx-inxold-3);
00843 
00844                 float tabley1 = kb.i0win_tab(dely-inyold+3);
00845                 float tabley2 = kb.i0win_tab(dely-inyold+2);
00846                 float tabley3 = kb.i0win_tab(dely-inyold+1);
00847                 float tabley4 = kb.i0win_tab(dely-inyold);
00848                 float tabley5 = kb.i0win_tab(dely-inyold-1);
00849                 float tabley6 = kb.i0win_tab(dely-inyold-2);
00850                 float tabley7 = kb.i0win_tab(dely-inyold-3);
00851 
00852                 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;
00853 
00854                 if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 )  {
00855                         x1 = (inxold-3+nx)%nx;
00856                         x2 = (inxold-2+nx)%nx;
00857                         x3 = (inxold-1+nx)%nx;
00858                         x4 = (inxold  +nx)%nx;
00859                         x5 = (inxold+1+nx)%nx;
00860                         x6 = (inxold+2+nx)%nx;
00861                         x7 = (inxold+3+nx)%nx;
00862 
00863                         y1 = ((inyold-3+ny)%ny)*nx;
00864                         y2 = ((inyold-2+ny)%ny)*nx;
00865                         y3 = ((inyold-1+ny)%ny)*nx;
00866                         y4 = ((inyold  +ny)%ny)*nx;
00867                         y5 = ((inyold+1+ny)%ny)*nx;
00868                         y6 = ((inyold+2+ny)%ny)*nx;
00869                         y7 = ((inyold+3+ny)%ny)*nx;
00870                 } else {
00871                         x1 = inxold-3;
00872                         x2 = inxold-2;
00873                         x3 = inxold-1;
00874                         x4 = inxold;
00875                         x5 = inxold+1;
00876                         x6 = inxold+2;
00877                         x7 = inxold+3;
00878 
00879                         y1 = (inyold-3)*nx;
00880                         y2 = (inyold-2)*nx;
00881                         y3 = (inyold-1)*nx;
00882                         y4 = inyold*nx;
00883                         y5 = (inyold+1)*nx;
00884                         y6 = (inyold+2)*nx;
00885                         y7 = (inyold+3)*nx;
00886                 }
00887 
00888                 pixel    = ( data[x1+y1]*tablex1 + data[x2+y1]*tablex2 + data[x3+y1]*tablex3 +
00889                              data[x4+y1]*tablex4 + data[x5+y1]*tablex5 + data[x6+y1]*tablex6 +
00890                              data[x7+y1]*tablex7 ) * tabley1 +
00891                            ( data[x1+y2]*tablex1 + data[x2+y2]*tablex2 + data[x3+y2]*tablex3 +
00892                              data[x4+y2]*tablex4 + data[x5+y2]*tablex5 + data[x6+y2]*tablex6 +
00893                              data[x7+y2]*tablex7 ) * tabley2 +
00894                            ( data[x1+y3]*tablex1 + data[x2+y3]*tablex2 + data[x3+y3]*tablex3 +
00895                              data[x4+y3]*tablex4 + data[x5+y3]*tablex5 + data[x6+y3]*tablex6 +
00896                              data[x7+y3]*tablex7 ) * tabley3 +
00897                            ( data[x1+y4]*tablex1 + data[x2+y4]*tablex2 + data[x3+y4]*tablex3 +
00898                              data[x4+y4]*tablex4 + data[x5+y4]*tablex5 + data[x6+y4]*tablex6 +
00899                              data[x7+y4]*tablex7 ) * tabley4 +
00900                            ( data[x1+y5]*tablex1 + data[x2+y5]*tablex2 + data[x3+y5]*tablex3 +
00901                              data[x4+y5]*tablex4 + data[x5+y5]*tablex5 + data[x6+y5]*tablex6 +
00902                              data[x7+y5]*tablex7 ) * tabley5 +
00903                            ( data[x1+y6]*tablex1 + data[x2+y6]*tablex2 + data[x3+y6]*tablex3 +
00904                              data[x4+y6]*tablex4 + data[x5+y6]*tablex5 + data[x6+y6]*tablex6 +
00905                              data[x7+y6]*tablex7 ) * tabley6 +
00906                            ( data[x1+y7]*tablex1 + data[x2+y7]*tablex2 + data[x3+y7]*tablex3 +
00907                              data[x4+y7]*tablex4 + data[x5+y7]*tablex5 + data[x6+y7]*tablex6 +
00908                              data[x7+y7]*tablex7 ) * tabley7;
00909 
00910                 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
00911                     (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
00912         } else {  //  3D
00913                 dely = restrict1(dely, ny);
00914                 int inyold = int(Util::round(dely));
00915                 delz = restrict1(delz, nz);
00916                 int inzold = int(Util::round(delz));
00917 
00918                 float tablex1 = kb.i0win_tab(delx-inxold+3);
00919                 float tablex2 = kb.i0win_tab(delx-inxold+2);
00920                 float tablex3 = kb.i0win_tab(delx-inxold+1);
00921                 float tablex4 = kb.i0win_tab(delx-inxold);
00922                 float tablex5 = kb.i0win_tab(delx-inxold-1);
00923                 float tablex6 = kb.i0win_tab(delx-inxold-2);
00924                 float tablex7 = kb.i0win_tab(delx-inxold-3);
00925 
00926                 float tabley1 = kb.i0win_tab(dely-inyold+3);
00927                 float tabley2 = kb.i0win_tab(dely-inyold+2);
00928                 float tabley3 = kb.i0win_tab(dely-inyold+1);
00929                 float tabley4 = kb.i0win_tab(dely-inyold);
00930                 float tabley5 = kb.i0win_tab(dely-inyold-1);
00931                 float tabley6 = kb.i0win_tab(dely-inyold-2);
00932                 float tabley7 = kb.i0win_tab(dely-inyold-3);
00933 
00934                 float tablez1 = kb.i0win_tab(delz-inzold+3);
00935                 float tablez2 = kb.i0win_tab(delz-inzold+2);
00936                 float tablez3 = kb.i0win_tab(delz-inzold+1);
00937                 float tablez4 = kb.i0win_tab(delz-inzold);
00938                 float tablez5 = kb.i0win_tab(delz-inzold-1);
00939                 float tablez6 = kb.i0win_tab(delz-inzold-2);
00940                 float tablez7 = kb.i0win_tab(delz-inzold-3);
00941 
00942                 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7, z1, z2, z3, z4, z5, z6, z7;
00943 
00944                 if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >= nz-kbc-2 )  {
00945                         x1 = (inxold-3+nx)%nx;
00946                         x2 = (inxold-2+nx)%nx;
00947                         x3 = (inxold-1+nx)%nx;
00948                         x4 = (inxold  +nx)%nx;
00949                         x5 = (inxold+1+nx)%nx;
00950                         x6 = (inxold+2+nx)%nx;
00951                         x7 = (inxold+3+nx)%nx;
00952 
00953                         y1 = ((inyold-3+ny)%ny)*nx;
00954                         y2 = ((inyold-2+ny)%ny)*nx;
00955                         y3 = ((inyold-1+ny)%ny)*nx;
00956                         y4 = ((inyold  +ny)%ny)*nx;
00957                         y5 = ((inyold+1+ny)%ny)*nx;
00958                         y6 = ((inyold+2+ny)%ny)*nx;
00959                         y7 = ((inyold+3+ny)%ny)*nx;
00960 
00961                         z1 = ((inzold-3+nz)%nz)*nx*ny;
00962                         z2 = ((inzold-2+nz)%nz)*nx*ny;
00963                         z3 = ((inzold-1+nz)%nz)*nx*ny;
00964                         z4 = ((inzold  +nz)%nz)*nx*ny;
00965                         z5 = ((inzold+1+nz)%nz)*nx*ny;
00966                         z6 = ((inzold+2+nz)%nz)*nx*ny;
00967                         z7 = ((inzold+3+nz)%nz)*nx*ny;
00968                 } else {
00969                         x1 = inxold-3;
00970                         x2 = inxold-2;
00971                         x3 = inxold-1;
00972                         x4 = inxold;
00973                         x5 = inxold+1;
00974                         x6 = inxold+2;
00975                         x7 = inxold+3;
00976 
00977                         y1 = (inyold-3)*nx;
00978                         y2 = (inyold-2)*nx;
00979                         y3 = (inyold-1)*nx;
00980                         y4 = inyold*nx;
00981                         y5 = (inyold+1)*nx;
00982                         y6 = (inyold+2)*nx;
00983                         y7 = (inyold+3)*nx;
00984 
00985                         z1 = (inzold-3)*nx*ny;
00986                         z2 = (inzold-2)*nx*ny;
00987                         z3 = (inzold-1)*nx*ny;
00988                         z4 = inzold*nx*ny;
00989                         z5 = (inzold+1)*nx*ny;
00990                         z6 = (inzold+2)*nx*ny;
00991                         z7 = (inzold+3)*nx*ny;
00992                 }
00993 
00994                 pixel  = ( ( data[x1+y1+z1]*tablex1 + data[x2+y1+z1]*tablex2 + data[x3+y1+z1]*tablex3 +
00995                              data[x4+y1+z1]*tablex4 + data[x5+y1+z1]*tablex5 + data[x6+y1+z1]*tablex6 +
00996                              data[x7+y1+z1]*tablex7 ) * tabley1 +
00997                            ( data[x1+y2+z1]*tablex1 + data[x2+y2+z1]*tablex2 + data[x3+y2+z1]*tablex3 +
00998                              data[x4+y2+z1]*tablex4 + data[x5+y2+z1]*tablex5 + data[x6+y2+z1]*tablex6 +
00999                              data[x7+y2+z1]*tablex7 ) * tabley2 +
01000                            ( data[x1+y3+z1]*tablex1 + data[x2+y3+z1]*tablex2 + data[x3+y3+z1]*tablex3 +
01001                              data[x4+y3+z1]*tablex4 + data[x5+y3+z1]*tablex5 + data[x6+y3+z1]*tablex6 +
01002                              data[x7+y3+z1]*tablex7 ) * tabley3 +
01003                            ( data[x1+y4+z1]*tablex1 + data[x2+y4+z1]*tablex2 + data[x3+y4+z1]*tablex3 +
01004                              data[x4+y4+z1]*tablex4 + data[x5+y4+z1]*tablex5 + data[x6+y4+z1]*tablex6 +
01005                              data[x7+y4+z1]*tablex7 ) * tabley4 +
01006                            ( data[x1+y5+z1]*tablex1 + data[x2+y5+z1]*tablex2 + data[x3+y5+z1]*tablex3 +
01007                              data[x4+y5+z1]*tablex4 + data[x5+y5+z1]*tablex5 + data[x6+y5+z1]*tablex6 +
01008                              data[x7+y5+z1]*tablex7 ) * tabley5 +
01009                            ( data[x1+y6+z1]*tablex1 + data[x2+y6+z1]*tablex2 + data[x3+y6+z1]*tablex3 +
01010                              data[x4+y6+z1]*tablex4 + data[x5+y6+z1]*tablex5 + data[x6+y6+z1]*tablex6 +
01011                              data[x7+y6+z1]*tablex7 ) * tabley6 +
01012                            ( data[x1+y7+z1]*tablex1 + data[x2+y7+z1]*tablex2 + data[x3+y7+z1]*tablex3 +
01013                              data[x4+y7+z1]*tablex4 + data[x5+y7+z1]*tablex5 + data[x6+y7+z1]*tablex6 +
01014                              data[x7+y7+z1]*tablex7 ) * tabley7 ) *tablez1 +
01015                          ( ( data[x1+y1+z2]*tablex1 + data[x2+y1+z2]*tablex2 + data[x3+y1+z2]*tablex3 +
01016                              data[x4+y1+z2]*tablex4 + data[x5+y1+z2]*tablex5 + data[x6+y1+z2]*tablex6 +
01017                              data[x7+y1+z2]*tablex7 ) * tabley1 +
01018                            ( data[x1+y2+z2]*tablex1 + data[x2+y2+z2]*tablex2 + data[x3+y2+z2]*tablex3 +
01019                              data[x4+y2+z2]*tablex4 + data[x5+y2+z2]*tablex5 + data[x6+y2+z2]*tablex6 +
01020                              data[x7+y2+z2]*tablex7 ) * tabley2 +
01021                            ( data[x1+y3+z2]*tablex1 + data[x2+y3+z2]*tablex2 + data[x3+y3+z2]*tablex3 +
01022                              data[x4+y3+z2]*tablex4 + data[x5+y3+z2]*tablex5 + data[x6+y3+z2]*tablex6 +
01023                              data[x7+y3+z2]*tablex7 ) * tabley3 +
01024                            ( data[x1+y4+z2]*tablex1 + data[x2+y4+z2]*tablex2 + data[x3+y4+z2]*tablex3 +
01025                              data[x4+y4+z2]*tablex4 + data[x5+y4+z2]*tablex5 + data[x6+y4+z2]*tablex6 +
01026                              data[x7+y4+z2]*tablex7 ) * tabley4 +
01027                            ( data[x1+y5+z2]*tablex1 + data[x2+y5+z2]*tablex2 + data[x3+y5+z2]*tablex3 +
01028                              data[x4+y5+z2]*tablex4 + data[x5+y5+z2]*tablex5 + data[x6+y5+z2]*tablex6 +
01029                              data[x7+y5+z2]*tablex7 ) * tabley5 +
01030                            ( data[x1+y6+z2]*tablex1 + data[x2+y6+z2]*tablex2 + data[x3+y6+z2]*tablex3 +
01031                              data[x4+y6+z2]*tablex4 + data[x5+y6+z2]*tablex5 + data[x6+y6+z2]*tablex6 +
01032                              data[x7+y6+z2]*tablex7 ) * tabley6 +
01033                            ( data[x1+y7+z2]*tablex1 + data[x2+y7+z2]*tablex2 + data[x3+y7+z2]*tablex3 +
01034                              data[x4+y7+z2]*tablex4 + data[x5+y7+z2]*tablex5 + data[x6+y7+z2]*tablex6 +
01035                              data[x7+y7+z2]*tablex7 ) * tabley7 ) *tablez2 +
01036                          ( ( data[x1+y1+z3]*tablex1 + data[x2+y1+z3]*tablex2 + data[x3+y1+z3]*tablex3 +
01037                              data[x4+y1+z3]*tablex4 + data[x5+y1+z3]*tablex5 + data[x6+y1+z3]*tablex6 +
01038                              data[x7+y1+z3]*tablex7 ) * tabley1 +
01039                            ( data[x1+y2+z3]*tablex1 + data[x2+y2+z3]*tablex2 + data[x3+y2+z3]*tablex3 +
01040                              data[x4+y2+z3]*tablex4 + data[x5+y2+z3]*tablex5 + data[x6+y2+z3]*tablex6 +
01041                              data[x7+y2+z3]*tablex7 ) * tabley2 +
01042                            ( data[x1+y3+z3]*tablex1 + data[x2+y3+z3]*tablex2 + data[x3+y3+z3]*tablex3 +
01043                              data[x4+y3+z3]*tablex4 + data[x5+y3+z3]*tablex5 + data[x6+y3+z3]*tablex6 +
01044                              data[x7+y3+z3]*tablex7 ) * tabley3 +
01045                            ( data[x1+y4+z3]*tablex1 + data[x2+y4+z3]*tablex2 + data[x3+y4+z3]*tablex3 +
01046                              data[x4+y4+z3]*tablex4 + data[x5+y4+z3]*tablex5 + data[x6+y4+z3]*tablex6 +
01047                              data[x7+y4+z3]*tablex7 ) * tabley4 +
01048                            ( data[x1+y5+z3]*tablex1 + data[x2+y5+z3]*tablex2 + data[x3+y5+z3]*tablex3 +
01049                              data[x4+y5+z3]*tablex4 + data[x5+y5+z3]*tablex5 + data[x6+y5+z3]*tablex6 +
01050                              data[x7+y5+z3]*tablex7 ) * tabley5 +
01051                            ( data[x1+y6+z3]*tablex1 + data[x2+y6+z3]*tablex2 + data[x3+y6+z3]*tablex3 +
01052                              data[x4+y6+z3]*tablex4 + data[x5+y6+z3]*tablex5 + data[x6+y6+z3]*tablex6 +
01053                              data[x7+y6+z3]*tablex7 ) * tabley6 +
01054                            ( data[x1+y7+z3]*tablex1 + data[x2+y7+z3]*tablex2 + data[x3+y7+z3]*tablex3 +
01055                              data[x4+y7+z3]*tablex4 + data[x5+y7+z3]*tablex5 + data[x6+y7+z3]*tablex6 +
01056                              data[x7+y7+z3]*tablex7 ) * tabley7 ) *tablez3 +
01057                          ( ( data[x1+y1+z4]*tablex1 + data[x2+y1+z4]*tablex2 + data[x3+y1+z4]*tablex3 +
01058                              data[x4+y1+z4]*tablex4 + data[x5+y1+z4]*tablex5 + data[x6+y1+z4]*tablex6 +
01059                              data[x7+y1+z4]*tablex7 ) * tabley1 +
01060                            ( data[x1+y2+z4]*tablex1 + data[x2+y2+z4]*tablex2 + data[x3+y2+z4]*tablex3 +
01061                              data[x4+y2+z4]*tablex4 + data[x5+y2+z4]*tablex5 + data[x6+y2+z4]*tablex6 +
01062                              data[x7+y2+z4]*tablex7 ) * tabley2 +
01063                            ( data[x1+y3+z4]*tablex1 + data[x2+y3+z4]*tablex2 + data[x3+y3+z4]*tablex3 +
01064                              data[x4+y3+z4]*tablex4 + data[x5+y3+z4]*tablex5 + data[x6+y3+z4]*tablex6 +
01065                              data[x7+y3+z4]*tablex7 ) * tabley3 +
01066                            ( data[x1+y4+z4]*tablex1 + data[x2+y4+z4]*tablex2 + data[x3+y4+z4]*tablex3 +
01067                              data[x4+y4+z4]*tablex4 + data[x5+y4+z4]*tablex5 + data[x6+y4+z4]*tablex6 +
01068                              data[x7+y4+z4]*tablex7 ) * tabley4 +
01069                            ( data[x1+y5+z4]*tablex1 + data[x2+y5+z4]*tablex2 + data[x3+y5+z4]*tablex3 +
01070                              data[x4+y5+z4]*tablex4 + data[x5+y5+z4]*tablex5 + data[x6+y5+z4]*tablex6 +
01071                              data[x7+y5+z4]*tablex7 ) * tabley5 +
01072                            ( data[x1+y6+z4]*tablex1 + data[x2+y6+z4]*tablex2 + data[x3+y6+z4]*tablex3 +
01073                              data[x4+y6+z4]*tablex4 + data[x5+y6+z4]*tablex5 + data[x6+y6+z4]*tablex6 +
01074                              data[x7+y6+z4]*tablex7 ) * tabley6 +
01075                            ( data[x1+y7+z4]*tablex1 + data[x2+y7+z4]*tablex2 + data[x3+y7+z4]*tablex3 +
01076                              data[x4+y7+z4]*tablex4 + data[x5+y7+z4]*tablex5 + data[x6+y7+z4]*tablex6 +
01077                              data[x7+y7+z4]*tablex7 ) * tabley7 ) *tablez4 +
01078                          ( ( data[x1+y1+z5]*tablex1 + data[x2+y1+z5]*tablex2 + data[x3+y1+z5]*tablex3 +
01079                              data[x4+y1+z5]*tablex4 + data[x5+y1+z5]*tablex5 + data[x6+y1+z5]*tablex6 +
01080                              data[x7+y1+z5]*tablex7 ) * tabley1 +
01081                            ( data[x1+y2+z5]*tablex1 + data[x2+y2+z5]*tablex2 + data[x3+y2+z5]*tablex3 +
01082                              data[x4+y2+z5]*tablex4 + data[x5+y2+z5]*tablex5 + data[x6+y2+z5]*tablex6 +
01083                              data[x7+y2+z5]*tablex7 ) * tabley2 +
01084                            ( data[x1+y3+z5]*tablex1 + data[x2+y3+z5]*tablex2 + data[x3+y3+z5]*tablex3 +
01085                              data[x4+y3+z5]*tablex4 + data[x5+y3+z5]*tablex5 + data[x6+y3+z5]*tablex6 +
01086                              data[x7+y3+z5]*tablex7 ) * tabley3 +
01087                            ( data[x1+y4+z5]*tablex1 + data[x2+y4+z5]*tablex2 + data[x3+y4+z5]*tablex3 +
01088                              data[x4+y4+z5]*tablex4 + data[x5+y4+z5]*tablex5 + data[x6+y4+z5]*tablex6 +
01089                              data[x7+y4+z5]*tablex7 ) * tabley4 +
01090                            ( data[x1+y5+z5]*tablex1 + data[x2+y5+z5]*tablex2 + data[x3+y5+z5]*tablex3 +
01091                              data[x4+y5+z5]*tablex4 + data[x5+y5+z5]*tablex5 + data[x6+y5+z5]*tablex6 +
01092                              data[x7+y5+z5]*tablex7 ) * tabley5 +
01093                            ( data[x1+y6+z5]*tablex1 + data[x2+y6+z5]*tablex2 + data[x3+y6+z5]*tablex3 +
01094                              data[x4+y6+z5]*tablex4 + data[x5+y6+z5]*tablex5 + data[x6+y6+z5]*tablex6 +
01095                              data[x7+y6+z5]*tablex7 ) * tabley6 +
01096                            ( data[x1+y7+z5]*tablex1 + data[x2+y7+z5]*tablex2 + data[x3+y7+z5]*tablex3 +
01097                              data[x4+y7+z5]*tablex4 + data[x5+y7+z5]*tablex5 + data[x6+y7+z5]*tablex6 +
01098                              data[x7+y7+z5]*tablex7 ) * tabley7 ) *tablez5 +
01099                          ( ( data[x1+y1+z6]*tablex1 + data[x2+y1+z6]*tablex2 + data[x3+y1+z6]*tablex3 +
01100                              data[x4+y1+z6]*tablex4 + data[x5+y1+z6]*tablex5 + data[x6+y1+z6]*tablex6 +
01101                              data[x7+y1+z6]*tablex7 ) * tabley1 +
01102                            ( data[x1+y2+z6]*tablex1 + data[x2+y2+z6]*tablex2 + data[x3+y2+z6]*tablex3 +
01103                              data[x4+y2+z6]*tablex4 + data[x5+y2+z6]*tablex5 + data[x6+y2+z6]*tablex6 +
01104                              data[x7+y2+z6]*tablex7 ) * tabley2 +
01105                            ( data[x1+y3+z6]*tablex1 + data[x2+y3+z6]*tablex2 + data[x3+y3+z6]*tablex3 +
01106                              data[x4+y3+z6]*tablex4 + data[x5+y3+z6]*tablex5 + data[x6+y3+z6]*tablex6 +
01107                              data[x7+y3+z6]*tablex7 ) * tabley3 +
01108                            ( data[x1+y4+z6]*tablex1 + data[x2+y4+z6]*tablex2 + data[x3+y4+z6]*tablex3 +
01109                              data[x4+y4+z6]*tablex4 + data[x5+y4+z6]*tablex5 + data[x6+y4+z6]*tablex6 +
01110                              data[x7+y4+z6]*tablex7 ) * tabley4 +
01111                            ( data[x1+y5+z6]*tablex1 + data[x2+y5+z6]*tablex2 + data[x3+y5+z6]*tablex3 +
01112                              data[x4+y5+z6]*tablex4 + data[x5+y5+z6]*tablex5 + data[x6+y5+z6]*tablex6 +
01113                              data[x7+y5+z6]*tablex7 ) * tabley5 +
01114                            ( data[x1+y6+z6]*tablex1 + data[x2+y6+z6]*tablex2 + data[x3+y6+z6]*tablex3 +
01115                              data[x4+y6+z6]*tablex4 + data[x5+y6+z6]*tablex5 + data[x6+y6+z6]*tablex6 +
01116                              data[x7+y6+z6]*tablex7 ) * tabley6 +
01117                            ( data[x1+y7+z6]*tablex1 + data[x2+y7+z6]*tablex2 + data[x3+y7+z6]*tablex3 +
01118                              data[x4+y7+z6]*tablex4 + data[x5+y7+z6]*tablex5 + data[x6+y7+z6]*tablex6 +
01119                              data[x7+y7+z6]*tablex7 ) * tabley7 ) *tablez6 +
01120                          ( ( data[x1+y1+z7]*tablex1 + data[x2+y1+z7]*tablex2 + data[x3+y1+z7]*tablex3 +
01121                              data[x4+y1+z7]*tablex4 + data[x5+y1+z7]*tablex5 + data[x6+y1+z7]*tablex6 +
01122                              data[x7+y1+z7]*tablex7 ) * tabley1 +
01123                            ( data[x1+y2+z7]*tablex1 + data[x2+y2+z7]*tablex2 + data[x3+y2+z7]*tablex3 +
01124                              data[x4+y2+z7]*tablex4 + data[x5+y2+z7]*tablex5 + data[x6+y2+z7]*tablex6 +
01125                              data[x7+y2+z7]*tablex7 ) * tabley2 +
01126                            ( data[x1+y3+z7]*tablex1 + data[x2+y3+z7]*tablex2 + data[x3+y3+z7]*tablex3 +
01127                              data[x4+y3+z7]*tablex4 + data[x5+y3+z7]*tablex5 + data[x6+y3+z7]*tablex6 +
01128                              data[x7+y3+z7]*tablex7 ) * tabley3 +
01129                            ( data[x1+y4+z7]*tablex1 + data[x2+y4+z7]*tablex2 + data[x3+y4+z7]*tablex3 +
01130                              data[x4+y4+z7]*tablex4 + data[x5+y4+z7]*tablex5 + data[x6+y4+z7]*tablex6 +
01131                              data[x7+y4+z7]*tablex7 ) * tabley4 +
01132                            ( data[x1+y5+z7]*tablex1 + data[x2+y5+z7]*tablex2 + data[x3+y5+z7]*tablex3 +
01133                              data[x4+y5+z7]*tablex4 + data[x5+y5+z7]*tablex5 + data[x6+y5+z7]*tablex6 +
01134                              data[x7+y5+z7]*tablex7 ) * tabley5 +
01135                            ( data[x1+y6+z7]*tablex1 + data[x2+y6+z7]*tablex2 + data[x3+y6+z7]*tablex3 +
01136                              data[x4+y6+z7]*tablex4 + data[x5+y6+z7]*tablex5 + data[x6+y6+z7]*tablex6 +
01137                              data[x7+y6+z7]*tablex7 ) * tabley6 +
01138                            ( data[x1+y7+z7]*tablex1 + data[x2+y7+z7]*tablex2 + data[x3+y7+z7]*tablex3 +
01139                              data[x4+y7+z7]*tablex4 + data[x5+y7+z7]*tablex5 + data[x6+y7+z7]*tablex6 +
01140                              data[x7+y7+z7]*tablex7 ) * tabley7 ) *tablez7;
01141 
01142                 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
01143                     (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7) *
01144                     (tablez1+tablez2+tablez3+tablez4+tablez5+tablez6+tablez7);
01145         }
01146         return pixel/w;
01147 }

float Util::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]`

Definition at line 1149 of file util_sparx.cpp.

References EMAN::Util::KaiserBessel::get_window_size(), EMAN::Util::KaiserBessel::i0win_tab(), restrict1(), and round().

Referenced by EMAN::EMData::rot_scale_conv_new_background().

01149                                                                                                                                                             {
01150 // Here counting is in C style, so coordinates of the pixel delx should be [0-nx-1]
01151 
01152 /* Commented by Zhengfan Yang on 04/20/07
01153 This function is written to replace get_pixel_conv(), which is too slow to use in practice.
01154 I made the following changes to get_pixel_conv():
01155 1. Use the same data passing scheme as quadri() and move the function from emdata_sparx.cpp to util_sparx.cpp
01156 2. Reduce usage of i0win_tab (from 98 calls to 14 calls in 2D case, from 1029 calls to 21 calls in 3D case!)
01157 3. Unfold the 'for' loop
01158 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)
01159 
01160 The shortcoming of this routine is that it only works for window size N=7. In case you want to use other window
01161 size, say N=5, you can easily modify it by referring my code.
01162 */
01163         int K = kb.get_window_size();
01164         int kbmin = -K/2;
01165         int kbmax = -kbmin;
01166         int kbc = kbmax+1;
01167 
01168         float pixel =0.0f;
01169         float w=0.0f;
01170 
01171     float argdelx = delx; // adding this for 2D case where the wrap around is not done circulantly using restrict1.
01172         delx = restrict1(delx, nx);
01173         int inxold = int(round(delx));
01174         if ( ny < 2 ) {  //1D
01175                 float tablex1 = kb.i0win_tab(delx-inxold+3);
01176                 float tablex2 = kb.i0win_tab(delx-inxold+2);
01177                 float tablex3 = kb.i0win_tab(delx-inxold+1);
01178                 float tablex4 = kb.i0win_tab(delx-inxold);
01179                 float tablex5 = kb.i0win_tab(delx-inxold-1);
01180                 float tablex6 = kb.i0win_tab(delx-inxold-2);
01181                 float tablex7 = kb.i0win_tab(delx-inxold-3);
01182 
01183                 int x1, x2, x3, x4, x5, x6, x7;
01184 
01185                 if ( inxold <= kbc || inxold >=nx-kbc-2 )  {
01186                         x1 = (inxold-3+nx)%nx;
01187                         x2 = (inxold-2+nx)%nx;
01188                         x3 = (inxold-1+nx)%nx;
01189                         x4 = (inxold  +nx)%nx;
01190                         x5 = (inxold+1+nx)%nx;
01191                         x6 = (inxold+2+nx)%nx;
01192                         x7 = (inxold+3+nx)%nx;
01193                 } else {
01194                         x1 = inxold-3;
01195                         x2 = inxold-2;
01196                         x3 = inxold-1;
01197                         x4 = inxold;
01198                         x5 = inxold+1;
01199                         x6 = inxold+2;
01200                         x7 = inxold+3;
01201                 }
01202 
01203                 pixel = data[x1]*tablex1 + data[x2]*tablex2 + data[x3]*tablex3 +
01204                         data[x4]*tablex4 + data[x5]*tablex5 + data[x6]*tablex6 +
01205                         data[x7]*tablex7 ;
01206 
01207                 w = tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7;
01208         } else if ( nz < 2 ) {  // 2D
01209                 
01210                 delx = argdelx;
01211                 // the wrap around is not done circulantly for 2D case; if (argdelx, argdely) is not in the image, then make them (xnew, ynew) which is definitely in the image
01212                 if ((delx < 0.0f) || (delx >= (float) (nx)) || (dely < 0.0f) || (dely >= (float) (ny)) ){
01213                 delx = (float)xnew*2.0;
01214                 dely = (float)ynew*2.0;
01215                 }
01216                 
01217                 int inxold = int(round(delx));
01218                 int inyold = int(round(dely));
01219                 
01220                 float tablex1 = kb.i0win_tab(delx-inxold+3);
01221                 float tablex2 = kb.i0win_tab(delx-inxold+2);
01222                 float tablex3 = kb.i0win_tab(delx-inxold+1);
01223                 float tablex4 = kb.i0win_tab(delx-inxold);
01224                 float tablex5 = kb.i0win_tab(delx-inxold-1);
01225                 float tablex6 = kb.i0win_tab(delx-inxold-2);
01226                 float tablex7 = kb.i0win_tab(delx-inxold-3);
01227 
01228                 float tabley1 = kb.i0win_tab(dely-inyold+3);
01229                 float tabley2 = kb.i0win_tab(dely-inyold+2);
01230                 float tabley3 = kb.i0win_tab(dely-inyold+1);
01231                 float tabley4 = kb.i0win_tab(dely-inyold);
01232                 float tabley5 = kb.i0win_tab(dely-inyold-1);
01233                 float tabley6 = kb.i0win_tab(dely-inyold-2);
01234                 float tabley7 = kb.i0win_tab(dely-inyold-3);
01235 
01236                 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7;
01237 
01238                 if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 )  {
01239                         x1 = (inxold-3+nx)%nx;
01240                         x2 = (inxold-2+nx)%nx;
01241                         x3 = (inxold-1+nx)%nx;
01242                         x4 = (inxold  +nx)%nx;
01243                         x5 = (inxold+1+nx)%nx;
01244                         x6 = (inxold+2+nx)%nx;
01245                         x7 = (inxold+3+nx)%nx;
01246 
01247                         y1 = ((inyold-3+ny)%ny)*nx;
01248                         y2 = ((inyold-2+ny)%ny)*nx;
01249                         y3 = ((inyold-1+ny)%ny)*nx;
01250                         y4 = ((inyold  +ny)%ny)*nx;
01251                         y5 = ((inyold+1+ny)%ny)*nx;
01252                         y6 = ((inyold+2+ny)%ny)*nx;
01253                         y7 = ((inyold+3+ny)%ny)*nx;
01254                 } else {
01255                         x1 = inxold-3;
01256                         x2 = inxold-2;
01257                         x3 = inxold-1;
01258                         x4 = inxold;
01259                         x5 = inxold+1;
01260                         x6 = inxold+2;
01261                         x7 = inxold+3;
01262 
01263                         y1 = (inyold-3)*nx;
01264                         y2 = (inyold-2)*nx;
01265                         y3 = (inyold-1)*nx;
01266                         y4 = inyold*nx;
01267                         y5 = (inyold+1)*nx;
01268                         y6 = (inyold+2)*nx;
01269                         y7 = (inyold+3)*nx;
01270                 }
01271 
01272                 pixel    = ( data[x1+y1]*tablex1 + data[x2+y1]*tablex2 + data[x3+y1]*tablex3 +
01273                              data[x4+y1]*tablex4 + data[x5+y1]*tablex5 + data[x6+y1]*tablex6 +
01274                              data[x7+y1]*tablex7 ) * tabley1 +
01275                            ( data[x1+y2]*tablex1 + data[x2+y2]*tablex2 + data[x3+y2]*tablex3 +
01276                              data[x4+y2]*tablex4 + data[x5+y2]*tablex5 + data[x6+y2]*tablex6 +
01277                              data[x7+y2]*tablex7 ) * tabley2 +
01278                            ( data[x1+y3]*tablex1 + data[x2+y3]*tablex2 + data[x3+y3]*tablex3 +
01279                              data[x4+y3]*tablex4 + data[x5+y3]*tablex5 + data[x6+y3]*tablex6 +
01280                              data[x7+y3]*tablex7 ) * tabley3 +
01281                            ( data[x1+y4]*tablex1 + data[x2+y4]*tablex2 + data[x3+y4]*tablex3 +
01282                              data[x4+y4]*tablex4 + data[x5+y4]*tablex5 + data[x6+y4]*tablex6 +
01283                              data[x7+y4]*tablex7 ) * tabley4 +
01284                            ( data[x1+y5]*tablex1 + data[x2+y5]*tablex2 + data[x3+y5]*tablex3 +
01285                              data[x4+y5]*tablex4 + data[x5+y5]*tablex5 + data[x6+y5]*tablex6 +
01286                              data[x7+y5]*tablex7 ) * tabley5 +
01287                            ( data[x1+y6]*tablex1 + data[x2+y6]*tablex2 + data[x3+y6]*tablex3 +
01288                              data[x4+y6]*tablex4 + data[x5+y6]*tablex5 + data[x6+y6]*tablex6 +
01289                              data[x7+y6]*tablex7 ) * tabley6 +
01290                            ( data[x1+y7]*tablex1 + data[x2+y7]*tablex2 + data[x3+y7]*tablex3 +
01291                              data[x4+y7]*tablex4 + data[x5+y7]*tablex5 + data[x6+y7]*tablex6 +
01292                              data[x7+y7]*tablex7 ) * tabley7;
01293 
01294                 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
01295                     (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7);
01296         } else {  //  3D
01297                 dely = restrict1(dely, ny);
01298                 int inyold = int(Util::round(dely));
01299                 delz = restrict1(delz, nz);
01300                 int inzold = int(Util::round(delz));
01301 
01302                 float tablex1 = kb.i0win_tab(delx-inxold+3);
01303                 float tablex2 = kb.i0win_tab(delx-inxold+2);
01304                 float tablex3 = kb.i0win_tab(delx-inxold+1);
01305                 float tablex4 = kb.i0win_tab(delx-inxold);
01306                 float tablex5 = kb.i0win_tab(delx-inxold-1);
01307                 float tablex6 = kb.i0win_tab(delx-inxold-2);
01308                 float tablex7 = kb.i0win_tab(delx-inxold-3);
01309 
01310                 float tabley1 = kb.i0win_tab(dely-inyold+3);
01311                 float tabley2 = kb.i0win_tab(dely-inyold+2);
01312                 float tabley3 = kb.i0win_tab(dely-inyold+1);
01313                 float tabley4 = kb.i0win_tab(dely-inyold);
01314                 float tabley5 = kb.i0win_tab(dely-inyold-1);
01315                 float tabley6 = kb.i0win_tab(dely-inyold-2);
01316                 float tabley7 = kb.i0win_tab(dely-inyold-3);
01317 
01318                 float tablez1 = kb.i0win_tab(delz-inzold+3);
01319                 float tablez2 = kb.i0win_tab(delz-inzold+2);
01320                 float tablez3 = kb.i0win_tab(delz-inzold+1);
01321                 float tablez4 = kb.i0win_tab(delz-inzold);
01322                 float tablez5 = kb.i0win_tab(delz-inzold-1);
01323                 float tablez6 = kb.i0win_tab(delz-inzold-2);
01324                 float tablez7 = kb.i0win_tab(delz-inzold-3);
01325 
01326                 int x1, x2, x3, x4, x5, x6, x7, y1, y2, y3, y4, y5, y6, y7, z1, z2, z3, z4, z5, z6, z7;
01327 
01328                 if ( inxold <= kbc || inxold >=nx-kbc-2 || inyold <= kbc || inyold >=ny-kbc-2 || inzold <= kbc || inzold >= nz-kbc-2 )  {
01329                         x1 = (inxold-3+nx)%nx;
01330                         x2 = (inxold-2+nx)%nx;
01331                         x3 = (inxold-1+nx)%nx;
01332                         x4 = (inxold  +nx)%nx;
01333                         x5 = (inxold+1+nx)%nx;
01334                         x6 = (inxold+2+nx)%nx;
01335                         x7 = (inxold+3+nx)%nx;
01336 
01337                         y1 = ((inyold-3+ny)%ny)*nx;
01338                         y2 = ((inyold-2+ny)%ny)*nx;
01339                         y3 = ((inyold-1+ny)%ny)*nx;
01340                         y4 = ((inyold  +ny)%ny)*nx;
01341                         y5 = ((inyold+1+ny)%ny)*nx;
01342                         y6 = ((inyold+2+ny)%ny)*nx;
01343                         y7 = ((inyold+3+ny)%ny)*nx;
01344 
01345                         z1 = ((inzold-3+nz)%nz)*nx*ny;
01346                         z2 = ((inzold-2+nz)%nz)*nx*ny;
01347                         z3 = ((inzold-1+nz)%nz)*nx*ny;
01348                         z4 = ((inzold  +nz)%nz)*nx*ny;
01349                         z5 = ((inzold+1+nz)%nz)*nx*ny;
01350                         z6 = ((inzold+2+nz)%nz)*nx*ny;
01351                         z7 = ((inzold+3+nz)%nz)*nx*ny;
01352                 } else {
01353                         x1 = inxold-3;
01354                         x2 = inxold-2;
01355                         x3 = inxold-1;
01356                         x4 = inxold;
01357                         x5 = inxold+1;
01358                         x6 = inxold+2;
01359                         x7 = inxold+3;
01360 
01361                         y1 = (inyold-3)*nx;
01362                         y2 = (inyold-2)*nx;
01363                         y3 = (inyold-1)*nx;
01364                         y4 = inyold*nx;
01365                         y5 = (inyold+1)*nx;
01366                         y6 = (inyold+2)*nx;
01367                         y7 = (inyold+3)*nx;
01368 
01369                         z1 = (inzold-3)*nx*ny;
01370                         z2 = (inzold-2)*nx*ny;
01371                         z3 = (inzold-1)*nx*ny;
01372                         z4 = inzold*nx*ny;
01373                         z5 = (inzold+1)*nx*ny;
01374                         z6 = (inzold+2)*nx*ny;
01375                         z7 = (inzold+3)*nx*ny;
01376                 }
01377 
01378                 pixel  = ( ( data[x1+y1+z1]*tablex1 + data[x2+y1+z1]*tablex2 + data[x3+y1+z1]*tablex3 +
01379                              data[x4+y1+z1]*tablex4 + data[x5+y1+z1]*tablex5 + data[x6+y1+z1]*tablex6 +
01380                              data[x7+y1+z1]*tablex7 ) * tabley1 +
01381                            ( data[x1+y2+z1]*tablex1 + data[x2+y2+z1]*tablex2 + data[x3+y2+z1]*tablex3 +
01382                              data[x4+y2+z1]*tablex4 + data[x5+y2+z1]*tablex5 + data[x6+y2+z1]*tablex6 +
01383                              data[x7+y2+z1]*tablex7 ) * tabley2 +
01384                            ( data[x1+y3+z1]*tablex1 + data[x2+y3+z1]*tablex2 + data[x3+y3+z1]*tablex3 +
01385                              data[x4+y3+z1]*tablex4 + data[x5+y3+z1]*tablex5 + data[x6+y3+z1]*tablex6 +
01386                              data[x7+y3+z1]*tablex7 ) * tabley3 +
01387                            ( data[x1+y4+z1]*tablex1 + data[x2+y4+z1]*tablex2 + data[x3+y4+z1]*tablex3 +
01388                              data[x4+y4+z1]*tablex4 + data[x5+y4+z1]*tablex5 + data[x6+y4+z1]*tablex6 +
01389                              data[x7+y4+z1]*tablex7 ) * tabley4 +
01390                            ( data[x1+y5+z1]*tablex1 + data[x2+y5+z1]*tablex2 + data[x3+y5+z1]*tablex3 +
01391                              data[x4+y5+z1]*tablex4 + data[x5+y5+z1]*tablex5 + data[x6+y5+z1]*tablex6 +
01392                              data[x7+y5+z1]*tablex7 ) * tabley5 +
01393                            ( data[x1+y6+z1]*tablex1 + data[x2+y6+z1]*tablex2 + data[x3+y6+z1]*tablex3 +
01394                              data[x4+y6+z1]*tablex4 + data[x5+y6+z1]*tablex5 + data[x6+y6+z1]*tablex6 +
01395                              data[x7+y6+z1]*tablex7 ) * tabley6 +
01396                            ( data[x1+y7+z1]*tablex1 + data[x2+y7+z1]*tablex2 + data[x3+y7+z1]*tablex3 +
01397                              data[x4+y7+z1]*tablex4 + data[x5+y7+z1]*tablex5 + data[x6+y7+z1]*tablex6 +
01398                              data[x7+y7+z1]*tablex7 ) * tabley7 ) *tablez1 +
01399                          ( ( data[x1+y1+z2]*tablex1 + data[x2+y1+z2]*tablex2 + data[x3+y1+z2]*tablex3 +
01400                              data[x4+y1+z2]*tablex4 + data[x5+y1+z2]*tablex5 + data[x6+y1+z2]*tablex6 +
01401                              data[x7+y1+z2]*tablex7 ) * tabley1 +
01402                            ( data[x1+y2+z2]*tablex1 + data[x2+y2+z2]*tablex2 + data[x3+y2+z2]*tablex3 +
01403                              data[x4+y2+z2]*tablex4 + data[x5+y2+z2]*tablex5 + data[x6+y2+z2]*tablex6 +
01404                              data[x7+y2+z2]*tablex7 ) * tabley2 +
01405                            ( data[x1+y3+z2]*tablex1 + data[x2+y3+z2]*tablex2 + data[x3+y3+z2]*tablex3 +
01406                              data[x4+y3+z2]*tablex4 + data[x5+y3+z2]*tablex5 + data[x6+y3+z2]*tablex6 +
01407                              data[x7+y3+z2]*tablex7 ) * tabley3 +
01408                            ( data[x1+y4+z2]*tablex1 + data[x2+y4+z2]*tablex2 + data[x3+y4+z2]*tablex3 +
01409                              data[x4+y4+z2]*tablex4 + data[x5+y4+z2]*tablex5 + data[x6+y4+z2]*tablex6 +
01410                              data[x7+y4+z2]*tablex7 ) * tabley4 +
01411                            ( data[x1+y5+z2]*tablex1 + data[x2+y5+z2]*tablex2 + data[x3+y5+z2]*tablex3 +
01412                              data[x4+y5+z2]*tablex4 + data[x5+y5+z2]*tablex5 + data[x6+y5+z2]*tablex6 +
01413                              data[x7+y5+z2]*tablex7 ) * tabley5 +
01414                            ( data[x1+y6+z2]*tablex1 + data[x2+y6+z2]*tablex2 + data[x3+y6+z2]*tablex3 +
01415                              data[x4+y6+z2]*tablex4 + data[x5+y6+z2]*tablex5 + data[x6+y6+z2]*tablex6 +
01416                              data[x7+y6+z2]*tablex7 ) * tabley6 +
01417                            ( data[x1+y7+z2]*tablex1 + data[x2+y7+z2]*tablex2 + data[x3+y7+z2]*tablex3 +
01418                              data[x4+y7+z2]*tablex4 + data[x5+y7+z2]*tablex5 + data[x6+y7+z2]*tablex6 +
01419                              data[x7+y7+z2]*tablex7 ) * tabley7 ) *tablez2 +
01420                          ( ( data[x1+y1+z3]*tablex1 + data[x2+y1+z3]*tablex2 + data[x3+y1+z3]*tablex3 +
01421                              data[x4+y1+z3]*tablex4 + data[x5+y1+z3]*tablex5 + data[x6+y1+z3]*tablex6 +
01422                              data[x7+y1+z3]*tablex7 ) * tabley1 +
01423                            ( data[x1+y2+z3]*tablex1 + data[x2+y2+z3]*tablex2 + data[x3+y2+z3]*tablex3 +
01424                              data[x4+y2+z3]*tablex4 + data[x5+y2+z3]*tablex5 + data[x6+y2+z3]*tablex6 +
01425                              data[x7+y2+z3]*tablex7 ) * tabley2 +
01426                            ( data[x1+y3+z3]*tablex1 + data[x2+y3+z3]*tablex2 + data[x3+y3+z3]*tablex3 +
01427                              data[x4+y3+z3]*tablex4 + data[x5+y3+z3]*tablex5 + data[x6+y3+z3]*tablex6 +
01428                              data[x7+y3+z3]*tablex7 ) * tabley3 +
01429                            ( data[x1+y4+z3]*tablex1 + data[x2+y4+z3]*tablex2 + data[x3+y4+z3]*tablex3 +
01430                              data[x4+y4+z3]*tablex4 + data[x5+y4+z3]*tablex5 + data[x6+y4+z3]*tablex6 +
01431                              data[x7+y4+z3]*tablex7 ) * tabley4 +
01432                            ( data[x1+y5+z3]*tablex1 + data[x2+y5+z3]*tablex2 + data[x3+y5+z3]*tablex3 +
01433                              data[x4+y5+z3]*tablex4 + data[x5+y5+z3]*tablex5 + data[x6+y5+z3]*tablex6 +
01434                              data[x7+y5+z3]*tablex7 ) * tabley5 +
01435                            ( data[x1+y6+z3]*tablex1 + data[x2+y6+z3]*tablex2 + data[x3+y6+z3]*tablex3 +
01436                              data[x4+y6+z3]*tablex4 + data[x5+y6+z3]*tablex5 + data[x6+y6+z3]*tablex6 +
01437                              data[x7+y6+z3]*tablex7 ) * tabley6 +
01438                            ( data[x1+y7+z3]*tablex1 + data[x2+y7+z3]*tablex2 + data[x3+y7+z3]*tablex3 +
01439                              data[x4+y7+z3]*tablex4 + data[x5+y7+z3]*tablex5 + data[x6+y7+z3]*tablex6 +
01440                              data[x7+y7+z3]*tablex7 ) * tabley7 ) *tablez3 +
01441                          ( ( data[x1+y1+z4]*tablex1 + data[x2+y1+z4]*tablex2 + data[x3+y1+z4]*tablex3 +
01442                              data[x4+y1+z4]*tablex4 + data[x5+y1+z4]*tablex5 + data[x6+y1+z4]*tablex6 +
01443                              data[x7+y1+z4]*tablex7 ) * tabley1 +
01444                            ( data[x1+y2+z4]*tablex1 + data[x2+y2+z4]*tablex2 + data[x3+y2+z4]*tablex3 +
01445                              data[x4+y2+z4]*tablex4 + data[x5+y2+z4]*tablex5 + data[x6+y2+z4]*tablex6 +
01446                              data[x7+y2+z4]*tablex7 ) * tabley2 +
01447                            ( data[x1+y3+z4]*tablex1 + data[x2+y3+z4]*tablex2 + data[x3+y3+z4]*tablex3 +
01448                              data[x4+y3+z4]*tablex4 + data[x5+y3+z4]*tablex5 + data[x6+y3+z4]*tablex6 +
01449                              data[x7+y3+z4]*tablex7 ) * tabley3 +
01450                            ( data[x1+y4+z4]*tablex1 + data[x2+y4+z4]*tablex2 + data[x3+y4+z4]*tablex3 +
01451                              data[x4+y4+z4]*tablex4 + data[x5+y4+z4]*tablex5 + data[x6+y4+z4]*tablex6 +
01452                              data[x7+y4+z4]*tablex7 ) * tabley4 +
01453                            ( data[x1+y5+z4]*tablex1 + data[x2+y5+z4]*tablex2 + data[x3+y5+z4]*tablex3 +
01454                              data[x4+y5+z4]*tablex4 + data[x5+y5+z4]*tablex5 + data[x6+y5+z4]*tablex6 +
01455                              data[x7+y5+z4]*tablex7 ) * tabley5 +
01456                            ( data[x1+y6+z4]*tablex1 + data[x2+y6+z4]*tablex2 + data[x3+y6+z4]*tablex3 +
01457                              data[x4+y6+z4]*tablex4 + data[x5+y6+z4]*tablex5 + data[x6+y6+z4]*tablex6 +
01458                              data[x7+y6+z4]*tablex7 ) * tabley6 +
01459                            ( data[x1+y7+z4]*tablex1 + data[x2+y7+z4]*tablex2 + data[x3+y7+z4]*tablex3 +
01460                              data[x4+y7+z4]*tablex4 + data[x5+y7+z4]*tablex5 + data[x6+y7+z4]*tablex6 +
01461                              data[x7+y7+z4]*tablex7 ) * tabley7 ) *tablez4 +
01462                          ( ( data[x1+y1+z5]*tablex1 + data[x2+y1+z5]*tablex2 + data[x3+y1+z5]*tablex3 +
01463                              data[x4+y1+z5]*tablex4 + data[x5+y1+z5]*tablex5 + data[x6+y1+z5]*tablex6 +
01464                              data[x7+y1+z5]*tablex7 ) * tabley1 +
01465                            ( data[x1+y2+z5]*tablex1 + data[x2+y2+z5]*tablex2 + data[x3+y2+z5]*tablex3 +
01466                              data[x4+y2+z5]*tablex4 + data[x5+y2+z5]*tablex5 + data[x6+y2+z5]*tablex6 +
01467                              data[x7+y2+z5]*tablex7 ) * tabley2 +
01468                            ( data[x1+y3+z5]*tablex1 + data[x2+y3+z5]*tablex2 + data[x3+y3+z5]*tablex3 +
01469                              data[x4+y3+z5]*tablex4 + data[x5+y3+z5]*tablex5 + data[x6+y3+z5]*tablex6 +
01470                              data[x7+y3+z5]*tablex7 ) * tabley3 +
01471                            ( data[x1+y4+z5]*tablex1 + data[x2+y4+z5]*tablex2 + data[x3+y4+z5]*tablex3 +
01472                              data[x4+y4+z5]*tablex4 + data[x5+y4+z5]*tablex5 + data[x6+y4+z5]*tablex6 +
01473                              data[x7+y4+z5]*tablex7 ) * tabley4 +
01474                            ( data[x1+y5+z5]*tablex1 + data[x2+y5+z5]*tablex2 + data[x3+y5+z5]*tablex3 +
01475                              data[x4+y5+z5]*tablex4 + data[x5+y5+z5]*tablex5 + data[x6+y5+z5]*tablex6 +
01476                              data[x7+y5+z5]*tablex7 ) * tabley5 +
01477                            ( data[x1+y6+z5]*tablex1 + data[x2+y6+z5]*tablex2 + data[x3+y6+z5]*tablex3 +
01478                              data[x4+y6+z5]*tablex4 + data[x5+y6+z5]*tablex5 + data[x6+y6+z5]*tablex6 +
01479                              data[x7+y6+z5]*tablex7 ) * tabley6 +
01480                            ( data[x1+y7+z5]*tablex1 + data[x2+y7+z5]*tablex2 + data[x3+y7+z5]*tablex3 +
01481                              data[x4+y7+z5]*tablex4 + data[x5+y7+z5]*tablex5 + data[x6+y7+z5]*tablex6 +
01482                              data[x7+y7+z5]*tablex7 ) * tabley7 ) *tablez5 +
01483                          ( ( data[x1+y1+z6]*tablex1 + data[x2+y1+z6]*tablex2 + data[x3+y1+z6]*tablex3 +
01484                              data[x4+y1+z6]*tablex4 + data[x5+y1+z6]*tablex5 + data[x6+y1+z6]*tablex6 +
01485                              data[x7+y1+z6]*tablex7 ) * tabley1 +
01486                            ( data[x1+y2+z6]*tablex1 + data[x2+y2+z6]*tablex2 + data[x3+y2+z6]*tablex3 +
01487                              data[x4+y2+z6]*tablex4 + data[x5+y2+z6]*tablex5 + data[x6+y2+z6]*tablex6 +
01488                              data[x7+y2+z6]*tablex7 ) * tabley2 +
01489                            ( data[x1+y3+z6]*tablex1 + data[x2+y3+z6]*tablex2 + data[x3+y3+z6]*tablex3 +
01490                              data[x4+y3+z6]*tablex4 + data[x5+y3+z6]*tablex5 + data[x6+y3+z6]*tablex6 +
01491                              data[x7+y3+z6]*tablex7 ) * tabley3 +
01492                            ( data[x1+y4+z6]*tablex1 + data[x2+y4+z6]*tablex2 + data[x3+y4+z6]*tablex3 +
01493                              data[x4+y4+z6]*tablex4 + data[x5+y4+z6]*tablex5 + data[x6+y4+z6]*tablex6 +
01494                              data[x7+y4+z6]*tablex7 ) * tabley4 +
01495                            ( data[x1+y5+z6]*tablex1 + data[x2+y5+z6]*tablex2 + data[x3+y5+z6]*tablex3 +
01496                              data[x4+y5+z6]*tablex4 + data[x5+y5+z6]*tablex5 + data[x6+y5+z6]*tablex6 +
01497                              data[x7+y5+z6]*tablex7 ) * tabley5 +
01498                            ( data[x1+y6+z6]*tablex1 + data[x2+y6+z6]*tablex2 + data[x3+y6+z6]*tablex3 +
01499                              data[x4+y6+z6]*tablex4 + data[x5+y6+z6]*tablex5 + data[x6+y6+z6]*tablex6 +
01500                              data[x7+y6+z6]*tablex7 ) * tabley6 +
01501                            ( data[x1+y7+z6]*tablex1 + data[x2+y7+z6]*tablex2 + data[x3+y7+z6]*tablex3 +
01502                              data[x4+y7+z6]*tablex4 + data[x5+y7+z6]*tablex5 + data[x6+y7+z6]*tablex6 +
01503                              data[x7+y7+z6]*tablex7 ) * tabley7 ) *tablez6 +
01504                          ( ( data[x1+y1+z7]*tablex1 + data[x2+y1+z7]*tablex2 + data[x3+y1+z7]*tablex3 +
01505                              data[x4+y1+z7]*tablex4 + data[x5+y1+z7]*tablex5 + data[x6+y1+z7]*tablex6 +
01506                              data[x7+y1+z7]*tablex7 ) * tabley1 +
01507                            ( data[x1+y2+z7]*tablex1 + data[x2+y2+z7]*tablex2 + data[x3+y2+z7]*tablex3 +
01508                              data[x4+y2+z7]*tablex4 + data[x5+y2+z7]*tablex5 + data[x6+y2+z7]*tablex6 +
01509                              data[x7+y2+z7]*tablex7 ) * tabley2 +
01510                            ( data[x1+y3+z7]*tablex1 + data[x2+y3+z7]*tablex2 + data[x3+y3+z7]*tablex3 +
01511                              data[x4+y3+z7]*tablex4 + data[x5+y3+z7]*tablex5 + data[x6+y3+z7]*tablex6 +
01512                              data[x7+y3+z7]*tablex7 ) * tabley3 +
01513                            ( data[x1+y4+z7]*tablex1 + data[x2+y4+z7]*tablex2 + data[x3+y4+z7]*tablex3 +
01514                              data[x4+y4+z7]*tablex4 + data[x5+y4+z7]*tablex5 + data[x6+y4+z7]*tablex6 +
01515                              data[x7+y4+z7]*tablex7 ) * tabley4 +
01516                            ( data[x1+y5+z7]*tablex1 + data[x2+y5+z7]*tablex2 + data[x3+y5+z7]*tablex3 +
01517                              data[x4+y5+z7]*tablex4 + data[x5+y5+z7]*tablex5 + data[x6+y5+z7]*tablex6 +
01518                              data[x7+y5+z7]*tablex7 ) * tabley5 +
01519                            ( data[x1+y6+z7]*tablex1 + data[x2+y6+z7]*tablex2 + data[x3+y6+z7]*tablex3 +
01520                              data[x4+y6+z7]*tablex4 + data[x5+y6+z7]*tablex5 + data[x6+y6+z7]*tablex6 +
01521                              data[x7+y6+z7]*tablex7 ) * tabley6 +
01522                            ( data[x1+y7+z7]*tablex1 + data[x2+y7+z7]*tablex2 + data[x3+y7+z7]*tablex3 +
01523                              data[x4+y7+z7]*tablex4 + data[x5+y7+z7]*tablex5 + data[x6+y7+z7]*tablex6 +
01524                              data[x7+y7+z7]*tablex7 ) * tabley7 ) *tablez7;
01525 
01526                 w = (tablex1+tablex2+tablex3+tablex4+tablex5+tablex6+tablex7) *
01527                     (tabley1+tabley2+tabley3+tabley4+tabley5+tabley6+tabley7) *
01528                     (tablez1+tablez2+tablez3+tablez4+tablez5+tablez6+tablez7);
01529         }
01530         return pixel/w;
01531 }

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

Definition at line 1711 of file util_sparx.cpp.

References EMAN::EMData::cmplx(), EMAN::Util::KaiserBessel::i0win_tab(), ImageDimensionException, and round().

01711                                                                                                                    {
01712 
01713         int nxreal = nx - 2;
01714         if (nxreal != ny)
01715                 throw ImageDimensionException("extractpoint requires ny == nx");
01716         int nhalf = nxreal/2;
01717         bool flip = (nuxnew < 0.f);
01718         if (flip) {
01719                 nuxnew *= -1;
01720                 nuynew *= -1;
01721         }
01722         if (nuynew >= nhalf-0.5)  {
01723                 nuynew -= nxreal;
01724         } else if (nuynew < -nhalf-0.5) {
01725                 nuynew += nxreal;
01726         }
01727 
01728         // put (xnew,ynew) on a grid.  The indices will be wrong for
01729         // the Fourier elements in the image, but the grid sizing will
01730         // be correct.
01731         int ixn = int(Util::round(nuxnew));
01732         int iyn = int(Util::round(nuynew));
01733 
01734         // set up some temporary weighting arrays
01735         static float wy[7];
01736         static float wx[7];
01737 
01738         float iynn = nuynew - iyn;
01739         wy[0] = kb.i0win_tab(iynn+3);
01740         wy[1] = kb.i0win_tab(iynn+2);
01741         wy[2] = kb.i0win_tab(iynn+1);
01742         wy[3] = kb.i0win_tab(iynn);
01743         wy[4] = kb.i0win_tab(iynn-1);
01744         wy[5] = kb.i0win_tab(iynn-2);
01745         wy[6] = kb.i0win_tab(iynn-3);
01746 
01747         float ixnn = nuxnew - ixn;
01748         wx[0] = kb.i0win_tab(ixnn+3);
01749         wx[1] = kb.i0win_tab(ixnn+2);
01750         wx[2] = kb.i0win_tab(ixnn+1);
01751         wx[3] = kb.i0win_tab(ixnn);
01752         wx[4] = kb.i0win_tab(ixnn-1);
01753         wx[5] = kb.i0win_tab(ixnn-2);
01754         wx[6] = kb.i0win_tab(ixnn-3);
01755 
01756         float wsum = (wx[0]+wx[1]+wx[2]+wx[3]+wx[4]+wx[5]+wx[6])*(wy[0]+wy[1]+wy[2]+wy[3]+wy[4]+wy[5]+wy[6]);
01757 
01758         complex<float> result(0.f,0.f);
01759         if ((ixn >= 3) && (ixn <= nhalf-3) && (iyn >= -nhalf+3) && (iyn <= nhalf-4)) {
01760                 // (xin,yin) not within window border from the edge
01761                 for (int iy = 0; iy < 7; iy++) {
01762                         int iyp = iyn + iy - 3 ;
01763                         for (int ix = 0; ix < 7; ix++) {
01764                                 int ixp = ixn + ix - 3;
01765                                 float w = wx[ix]*wy[iy];
01766                                 complex<float> val = fimage->cmplx(ixp,iyp);
01767                                 result += val*w;
01768                         }
01769                 }
01770         } else {
01771                 // points that "stick out"
01772                 for (int iy = 0; iy < 7; iy++) {
01773                         int iyp = iyn + iy - 3;
01774                         for (int ix = 0; ix < 7; ix++) {
01775                                 int ixp = ixn + ix - 3;
01776                                 bool mirror = false;
01777                                 int ixt = ixp, iyt = iyp;
01778                                 if (ixt < 0) {
01779                                         ixt = -ixt;
01780                                         iyt = -iyt;
01781                                         mirror = !mirror;
01782                                 }
01783                                 if (ixt > nhalf) {
01784                                         ixt = nxreal - ixt;
01785                                         iyt = -iyt;
01786                                         mirror = !mirror;
01787                                 }
01788                                 if (iyt > nhalf-1)  iyt -= nxreal;
01789                                 if (iyt < -nhalf)   iyt += nxreal;
01790                                 float w = wx[ix]*wy[iy];
01791                                 complex<float> val = fimage->cmplx(ixt,iyt);
01792                                 if (mirror)  result += conj(val)*w;
01793                                 else         result += val*w;
01794                         }
01795                 }
01796         }
01797         if (flip)  result = conj(result)/wsum;
01798         else result /= wsum;
01799         return result;
01800 }

float Util::bilinear	(	float	xold,
		float	yold,
		int	nsam,
		int	nrow,
		float *	xim
	)			`[static]`

Definition at line 2337 of file util_sparx.cpp.

References xim.

02338 {
02339 /*
02340 c  purpose: linear interpolation
02341   Optimized for speed, circular closer removed, checking of ranges removed
02342 */
02343     float bilinear;
02344     int   ixold, iyold;
02345 
02346 /*
02347         float xdif, ydif, xrem, yrem;
02348         ixold   = (int) floor(xold);
02349         iyold   = (int) floor(yold);
02350         ydif = yold - iyold;
02351         yrem = 1.0f - ydif;
02352 
02353         //  May want to insert if?
02354 //              IF ((IYOLD .GE. 1 .AND. IYOLD .LE. NROW-1) .AND.
02355 //     &            (IXOLD .GE. 1 .AND. IXOLD .LE. NSAM-1)) THEN
02356 //c                INSIDE BOUNDARIES OF OUTPUT IMAGE
02357         xdif = xold - ixold;
02358         xrem = 1.0f- xdif;
02359 //                 RBUF(K) = YDIF*(BUF(NADDR+NSAM)*XREM
02360 //     &                    +BUF(NADDR+NSAM+1)*XDIF)
02361 //     &                    +YREM*(BUF(NADDR)*XREM + BUF(NADDR+1)*XDIF)
02362         bilinear = ydif*(xim(ixold,iyold+1)*xrem + xim(ixold+1,iyold+1)*xdif) +
02363                                         yrem*(xim(ixold,iyold)*xrem+xim(ixold+1,iyold)*xdif);
02364 
02365     return bilinear;
02366 }
02367 */
02368         float xdif, ydif;
02369 
02370         ixold   = (int) xold;
02371         iyold   = (int) yold;
02372         ydif = yold - iyold;
02373 
02374         //  May want to insert it?
02375 //              IF ((IYOLD .GE. 1 .AND. IYOLD .LE. NROW-1) .AND.
02376 //     &            (IXOLD .GE. 1 .AND. IXOLD .LE. NSAM-1)) THEN
02377 //c                INSIDE BOUNDARIES OF OUTPUT IMAGE
02378         xdif = xold - ixold;
02379         bilinear = xim(ixold, iyold) + ydif* (xim(ixold, iyold+1) - xim(ixold, iyold)) +
02380                    xdif* (xim(ixold+1, iyold) - xim(ixold, iyold) +
02381                            ydif* (xim(ixold+1, iyold+1) - xim(ixold+1, iyold) - xim(ixold, iyold+1) + xim(ixold, iyold)) );
02382 
02383         return bilinear;
02384 }

float Util::triquad	(	float	r,
		float	s,
		float	t,
		float *	fdata
	)			`[static]`

Quadratic interpolation (3D).

Parameters:

	r
	s
	t
	fdata

Returns:: Interpolated value

Definition at line 1930 of file util_sparx.cpp.

01931 {
01932 
01933     const float C2 = 0.5f;    //1.0 / 2.0;
01934     const float C4 = 0.25f;   //1.0 / 4.0;
01935     const float C8 = 0.125f;  //1.0 / 8.0;
01936 
01937     float  RS   = R * S;
01938     float  ST   = S * T;
01939     float  RT   = R * T;
01940     float  RST  = R * ST;
01941 
01942     float  RSQ  = 1-R*R;
01943     float  SSQ  = 1-S*S;
01944     float  TSQ  = 1-T*T;
01945 
01946     float  RM1  = (1-R);
01947     float  SM1  = (1-S);
01948     float  TM1  = (1-T);
01949 
01950     float  RP1  = (1+R);
01951     float  SP1  = (1+S);
01952     float  TP1  = (1+T);
01953 
01954     float triquad =
01955     (-C8) * RST * RM1  * SM1  * TM1 * fdata[0] +
01956         ( C4) * ST  * RSQ  * SM1  * TM1 * fdata[1] +
01957         ( C8) * RST * RP1  * SM1  * TM1 * fdata[2] +
01958         ( C4) * RT  * RM1  * SSQ  * TM1 * fdata[3] +
01959         (-C2) * T   * RSQ  * SSQ  * TM1 * fdata[4] +
01960         (-C4) * RT  * RP1  * SSQ  * TM1 * fdata[5] +
01961         ( C8) * RST * RM1  * SP1  * TM1 * fdata[6] +
01962         (-C4) * ST  * RSQ  * SP1  * TM1 * fdata[7] +
01963         (-C8) * RST * RP1  * SP1  * TM1 * fdata[8] +
01964 //
01965         ( C4) * RS  * RM1  * SM1  * TSQ * fdata[9]  +
01966         (-C2) * S   * RSQ  * SM1  * TSQ * fdata[10] +
01967         (-C4) * RS  * RP1  * SM1  * TSQ * fdata[11] +
01968         (-C2) * R   * RM1  * SSQ  * TSQ * fdata[12] +
01969                       RSQ  * SSQ  * TSQ * fdata[13] +
01970         ( C2) * R   * RP1  * SSQ  * TSQ * fdata[14] +
01971         (-C4) * RS  * RM1  * SP1  * TSQ * fdata[15] +
01972         ( C2) * S   * RSQ  * SP1  * TSQ * fdata[16] +
01973         ( C4) * RS  * RP1  * SP1  * TSQ * fdata[17] +
01974  //
01975         ( C8) * RST * RM1  * SM1  * TP1 * fdata[18] +
01976         (-C4) * ST  * RSQ  * SM1  * TP1 * fdata[19] +
01977         (-C8) * RST * RP1  * SM1  * TP1 * fdata[20] +
01978         (-C4) * RT  * RM1  * SSQ  * TP1 * fdata[21] +
01979         ( C2) * T   * RSQ  * SSQ  * TP1 * fdata[22] +
01980         ( C4) * RT  * RP1  * SSQ  * TP1 * fdata[23] +
01981         (-C8) * RST * RM1  * SP1  * TP1 * fdata[24] +
01982         ( C4) * ST  * RSQ  * SP1  * TP1 * fdata[25] +
01983         ( C8) * RST * RP1  * SP1  * TP1 * fdata[26]   ;
01984      return triquad;
01985 }

EMData * Util::Polar2D	(	EMData *	image,
		vector< int >	numr,
		string	mode
	)			`[static]`

Definition at line 2163 of file util_sparx.cpp.

References circ, EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), numr, quadri(), EMAN::EMData::set_size(), xim, and y.

02163                                                                   {
02164         int nsam = image->get_xsize();
02165         int nrow = image->get_ysize();
02166         int nring = numr.size()/3;
02167         int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
02168         EMData* out = new EMData();
02169         out->set_size(lcirc,1,1);
02170         char mode = (cmode == "F" || cmode == "f") ? 'f' : 'h';
02171         float *xim  = image->get_data();
02172         float *circ = out->get_data();
02173 /*   alrq(image->get_data(), nsam, nrow, &numr[0], out->get_data(), lcirc, nring, cmode);
02174    return out;
02175 }
02176 void Util::alrq(float *xim,  int nsam , int nrow , int *numr,
02177           float *circ, int lcirc, int nring, char mode)
02178 {*/
02179 /*
02180 c
02181 c  purpose:
02182 c
02183 c  resmaple to polar coordinates
02184 c
02185 */
02186         //  dimension         xim(nsam,nrow),circ(lcirc)
02187         //  integer           numr(3,nring)
02188 
02189         double dfi, dpi;
02190         int    ns2, nr2, i, inr, l, nsim, kcirc, lt, j;
02191         float  yq, xold, yold, fi, x, y;
02192 
02193         ns2 = nsam/2+1;
02194         nr2 = nrow/2+1;
02195         dpi = 2.0*atan(1.0);
02196 
02197         for (i=1;i<=nring;i++) {
02198                 // radius of the ring
02199                 inr = numr(1,i);
02200                 yq  = static_cast<float>(inr);
02201                 l   = numr(3,i);
02202                 if (mode == 'h' || mode == 'H')  lt = l/2;
02203                 else                             lt = l/4;
02204 
02205                 nsim           = lt-1;
02206                 dfi            = dpi/(nsim+1);
02207                 kcirc          = numr(2,i);
02208                 xold           = 0.0f;
02209                 yold           = static_cast<float>(inr);
02210                 circ(kcirc)    = quadri(xold+(float)ns2,yold+(float)nr2,nsam,nrow,xim);
02211                 xold           = static_cast<float>(inr);
02212                 yold           = 0.0f;
02213                 circ(lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02214 
02215                 if (mode == 'f' || mode == 'F') {
02216                         xold              = 0.0f;
02217                         yold              = static_cast<float>(-inr);
02218                         circ(lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02219                         xold              = static_cast<float>(-inr);
02220                         yold              = 0.0f;
02221                         circ(lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02222                 }
02223 
02224                 for (j=1;j<=nsim;j++) {
02225                         fi               = static_cast<float>(dfi*j);
02226                         x                = sin(fi)*yq;
02227                         y                = cos(fi)*yq;
02228                         xold             = x;
02229                         yold             = y;
02230                         circ(j+kcirc)    = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02231                         xold             =  y;
02232                         yold             = -x;
02233                         circ(j+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02234 
02235                         if (mode == 'f' || mode == 'F')  {
02236                                 xold                = -x;
02237                                 yold                = -y;
02238                                 circ(j+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02239                                 xold                = -y;
02240                                 yold                =  x;
02241                                 circ(j+lt+lt+lt+kcirc) = quadri(xold+ns2,yold+nr2,nsam,nrow,xim);
02242                         }
02243                 }
02244         }
02245         return  out;
02246 }

EMData * Util::Polar2Dm	(	EMData *	image,
		float	cns2,
		float	cnr2,
		vector< int >	numr,
		string	cmode
	)			`[static]`

Definition at line 2248 of file util_sparx.cpp.

References Assert, circ, EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), numr, quadri(), EMAN::EMData::set_size(), xim, and y.

Referenced by ali2d_ccf_list(), multiref_peaks_ali2d(), multiref_peaks_compress_ali2d(), multiref_polar_ali_2d(), multiref_polar_ali_2d_local(), multiref_polar_ali_2d_local_psi(), and multiref_polar_ali_2d_nom().

02248                                                                                            {
02249         int nsam = image->get_xsize();
02250         int nrow = image->get_ysize();
02251         int nring = numr.size()/3;
02252         int lcirc = numr[3*nring-2]+numr[3*nring-1]-1;
02253         EMData* out = new EMData();
02254         out->set_size(lcirc,1,1);
02255         char mode = (cmode == "F" || cmode == "f") ? 'f' : 'h';
02256         float *xim  = image->get_data();
02257         float *circ = out->get_data();
02258         double dpi, dfi;
02259         int    it, jt, inr, l, nsim, kcirc, lt;
02260         float  xold, yold, fi, x, y;
02261 
02262         //     cns2 and cnr2 are predefined centers
02263         //     no need to set to zero, all elements are defined
02264         dpi = 2*atan(1.0);
02265         for (it=1; it<=nring; it++) {
02266                 // radius of the ring
02267                 inr = numr(1,it);
02268 
02269                 // "F" means a full circle interpolation
02270                 // "H" means a half circle interpolation
02271 
02272                 l = numr(3,it);
02273                 if ( mode == 'h' || mode == 'H' ) lt = l / 2;
02274                 else                              lt = l / 4;
02275 
02276                 nsim  = lt - 1;
02277                 dfi   = dpi / (nsim+1);
02278                 kcirc = numr(2,it);
02279                 xold  = 0.0f+cns2;
02280                 yold  = inr+cnr2;
02281 
02282                 Assert( kcirc <= lcirc );
02283                 circ(kcirc) = quadri(xold,yold,nsam,nrow,xim);    // Sampling on 90 degree
02284 
02285                 xold  = inr+cns2;
02286                 yold  = 0.0f+cnr2;
02287                 Assert( lt+kcirc <= lcirc );
02288                 circ(lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on 0 degree
02289 
02290                 if ( mode == 'f' || mode == 'F' ) {
02291                         xold = 0.0f+cns2;
02292                         yold = -inr+cnr2;
02293                         Assert( lt+lt+kcirc <= lcirc );
02294                         circ(lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on 270 degree
02295 
02296                         xold = -inr+cns2;
02297                         yold = 0.0f+cnr2;
02298                         Assert(lt+lt+lt+kcirc <= lcirc );
02299                         circ(lt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim); // Sampling on 180 degree
02300                 }
02301 
02302                 for (jt=1; jt<=nsim; jt++) {
02303                         fi   = static_cast<float>(dfi * jt);
02304                         x    = sin(fi) * inr;
02305                         y    = cos(fi) * inr;
02306 
02307                         xold = x+cns2;
02308                         yold = y+cnr2;
02309 
02310                         Assert( jt+kcirc <= lcirc );
02311                         circ(jt+kcirc) = quadri(xold,yold,nsam,nrow,xim);      // Sampling on the first quadrant
02312 
02313                         xold = y+cns2;
02314                         yold = -x+cnr2;
02315 
02316                         Assert( jt+lt+kcirc <= lcirc );
02317                         circ(jt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);    // Sampling on the fourth quadrant
02318 
02319                         if ( mode == 'f' || mode == 'F' ) {
02320                                 xold = -x+cns2;
02321                                 yold = -y+cnr2;
02322 
02323                                 Assert( jt+lt+lt+kcirc <= lcirc );
02324                                 circ(jt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim); // Sampling on the third quadrant
02325 
02326                                 xold = -y+cns2;
02327                                 yold = x+cnr2;
02328 
02329                                 Assert( jt+lt+lt+lt+kcirc <= lcirc );
02330                                 circ(jt+lt+lt+lt+kcirc) = quadri(xold,yold,nsam,nrow,xim);  // Sampling on the second quadrant
02331                         }
02332                 } // end for jt
02333         } //end for it
02334         return out;
02335 }

void Util::alrl_ms	(	float *	xim,
		int	nsam,
		int	nrow,
		float	cns2,
		float	cnr2,
		int *	numr,
		float *	circ,
		int	lcirc,
		int	nring,
		char	mode
	)			`[static]`

Definition at line 2386 of file util_sparx.cpp.

References circ, numr, quadri(), and y.

02387                                                                     {
02388         double dpi, dfi;
02389         int    it, jt, inr, l, nsim, kcirc, lt;
02390         float   xold, yold, fi, x, y;
02391 
02392         //     cns2 and cnr2 are predefined centers
02393         //     no need to set to zero, all elements are defined
02394 
02395         dpi = 2*atan(1.0);
02396         for (it=1; it<=nring; it++) {
02397                 // radius of the ring
02398                 inr = numr(1,it);
02399 
02400                 l = numr(3,it);
02401

EMAN::Util Class Reference [unit test in Python]

Static Public Member Functions

Static Private Member Functions

Classes

Detailed Description

Member Function Documentation

EMAN::Util Class Reference
[unit test in Python]