EMAN2
Classes | Public Member Functions | Static Public Attributes | Private Types | Private Member Functions | Private Attributes | Friends | List of all members
EMAN::EMData Class Reference
tested » unit test in Python

EMData stores an image's data and defines core image processing routines. More...

#include <emdata.h>

Collaboration diagram for EMAN::EMData:
Collaboration graph
[legend]

Classes

class  ClipInplaceVariables
 

Public Member Functions

 EMData ()
 For all image I/O. More...
 
 ~EMData ()
 
 EMData (const string &filename, int image_index=0)
 Construct from an image file. More...
 
 EMData (int nx, int ny, int nz=1, bool is_real=true)
 
 EMData (float *data, const int nx, const int ny, const int nz, const Dict &attr_dict=Dict())
 Construction from a data pointer, dimensions must be supplied. More...
 
 EMData (float *data, float *cudadata, const int nx, const int ny, const int nz, const Dict &attr_dict=Dict())
 Construction from a data pointer for usage in cuda, dimensions must be supplied. More...
 
 EMData (const EMData &that)
 Construct from an EMData (copy constructor). More...
 
EMData & operator= (const EMData &that)
 EMData assignment operator Performs a deep copy. More...
 
EMData * get_clip (const Region &area, const float fill=0) const
 Get an inclusive clip. More...
 
void clip_inplace (const Region &area, const float &fill_value=0)
 Clip the image inplace - clipping region must be smaller than the current region internally memory is reallocated. More...
 
EMData * get_top_half () const
 Get the top half of this 3D image. More...
 
EMData * get_rotated_clip (const Transform &xform, const IntSize &size, float scale=1.0)
 This will extract an arbitrarily oriented and sized region from the image. More...
 
EMData * window_center (int l)
 Window the center of an image. More...
 
float * setup4slice (bool redo=true)
 Set up for fftslice operations. More...
 
void scale (float scale_factor)
 scale the image by a factor. More...
 
void translate (float dx, float dy, float dz)
 Translate this image. More...
 
void translate (const Vec3f &translation)
 Translate this image. More...
 
void translate (int dx, int dy, int dz)
 Translate this image. More...
 
void translate (const Vec3i &translation)
 Translate this image. More...
 
void rotate (const Transform &t)
 Rotate this image. More...
 
float max_3D_pixel_error (const Transform &t1, const Transform &t2, float r)
 
void rotate (float az, float alt, float phi)
 Rotate this image. More...
 
void transform (const Transform &t)
 Rotate then translate the image. More...
 
void rotate_translate (const Transform &t)
 Apply a transformation to the image. More...
 
void rotate_translate (float az, float alt, float phi, float dx, float dy, float dz)
 Rotate then translate the image. More...
 
void rotate_translate (float az, float alt, float phi, float dx, float dy, float dz, float pdx, float pdy, float pdz)
 Rotate then translate the image. More...
 
void rotate_x (int dx)
 This performs a translation of each line along x with wraparound. More...
 
void rotate_180 ()
 Fast rotation by 180 degrees. More...
 
double dot_rotate_translate (EMData *with, float dx, float dy, float da, const bool mirror=false)
 dot product of 2 images. More...
 
EMData * little_big_dot (EMData *little_img, bool do_sigma=false)
 This does a normalized dot product of a little image with a big image using real-space methods. More...
 
EMData * do_radon ()
 Radon Transform: an algorithm that transforms an original image into a series of equiangular projections. More...
 
EMData * calc_ccf (EMData *with=0, fp_flag fpflag=CIRCULANT, bool center=false)
 Calculate Cross-Correlation Function (CCF). More...
 
EMData * calc_ccf_masked (EMData *with, EMData *withsq=0, EMData *mask=0)
 Calculate cross correlation between this and with where this is assumed to have had a mask applied to it (zero value in masked-out regions). More...
 
void zero_corner_circulant (const int radius=0)
 Zero the pixels in the bottom left corner of the image If radius is greater than 1, than circulant zeroing occurs assuming that the center of operation starts in the bottom left corner and proceed outwards to the NE and backwards in a circulant fashion towards the SW. More...
 
EMData * calc_ccfx (EMData *const with, int y0=0, int y1=-1, bool nosum=false, bool flip=false, bool usez=false)
 Calculate Cross-Correlation Function (CCF) in the x-direction and adds them up, result in 1D. More...
 
EMData * make_rotational_footprint (bool unwrap=true)
 Makes a 'rotational footprint', which is an 'unwound' autocorrelation function. More...
 
EMData * make_rotational_footprint_e1 (bool unwrap=true)
 
EMData * make_rotational_footprint_cmc (bool unwrap=true)
 
EMData * make_footprint (int type=0)
 Makes a 'footprint' for the current image. More...
 
EMData * calc_mutual_correlation (EMData *with, bool tocorner=false, EMData *filter=0)
 Calculates mutual correlation function (MCF) between 2 images. More...
 
EMData * unwrap (int r1=-1, int r2=-1, int xs=-1, int dx=0, int dy=0, bool do360=false, bool weight_radial=true) const
 Maps to polar coordinates from Cartesian coordinates. More...
 
EMData * unwrap_largerR (int r1, int r2, int xs, float rmax_f)
 
EMData * oneDfftPolar (int size, float rmax, float MAXR)
 
void apply_radial_func (float x0, float dx, vector< float >array, bool interp=true)
 multiplies by a radial function in fourier space. More...
 
vector< float > calc_radial_dist (int n, float x0, float dx, int inten)
 calculates radial distribution. More...
 
vector< float > calc_radial_dist (int n, float x0, float dx, int nwedge, float offset, bool inten)
 calculates radial distribution subdivided by angle. More...
 
void cconj ()
 Replace the image its complex conjugate. More...
 
void add_incoherent (EMData *obj)
 Adds 'obj' to 'this' incoherently. More...
 
vector< float > calc_hist (int hist_size=128, float hist_min=0, float hist_max=0, const float &brt=0.0f, const float &cont=1.0f)
 Calculates the histogram of 'this' image. More...
 
vector< float > calc_az_dist (int n, float a0, float da, float rmin, float rmax)
 Caculates the azimuthal distributions. More...
 
float calc_dist (EMData *second_img, int y_index=0) const
 Calculates the distance between 2 vectors. More...
 
EMData * calc_flcf (EMData *with)
 Calculates the cross correlation with local normalization between 2 images. More...
 
EMData * calc_fast_sigma_image (EMData *mask)
 Calculates the local standard deviation (sigma) image using the given mask image. More...
 
EMData * convolute (EMData *with)
 Convolutes 2 data sets. More...
 
void common_lines (EMData *image1, EMData *image2, int mode=0, int steps=180, bool horizontal=false)
 Finds common lines between 2 complex images. More...
 
void common_lines_real (EMData *image1, EMData *image2, int steps=180, bool horizontal=false)
 Finds common lines between 2 real images. More...
 
void cut_slice (const EMData *const map, const Transform &tr, bool interpolate=true)
 cut a 2D slice out of a real 3D map. More...
 
void uncut_slice (EMData *const map, const Transform &tr) const
 Opposite of the cut_slice(). More...
 
EMData * extract_box (const Transform &cs, const Region &r)
 Extract a box from EMData in an abritrary orrientation. More...
 
int getResolution () const
 function for MarchingCubes, for 3D image display More...
 
void debug_print_parms ()
 Printing EMData params for debugging purpose. More...
 
void set_xyz_origin (float origin_x, float origin_y, float origin_z)
 Set the x,y,z origin of the image. More...
 
EMData * compute_missingwedge (float wedgeangle, float start=0.05, float stop=0.5)
 Find the mean and variance of voxels in the missing wedge. More...
 

Static Public Attributes

static int totalalloc =0
 

Private Types

enum  EMDataFlags {
  EMDATA_BUSY = 1 << 3 , EMDATA_HASCTFF = 1 << 4 , EMDATA_NEEDUPD = 1 << 5 , EMDATA_FLIP = 1 << 7 ,
  EMDATA_PAD = 1 << 8 , EMDATA_FFTODD = 1 << 9 , EMDATA_SHUFFLE = 1 << 10 , EMDATA_FH = 1 << 11 ,
  EMDATA_CPU_NEEDS_UPDATE = 1 << 12 , EMDATA_GPU_NEEDS_UPDATE = 1 << 13 , EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14
}
 This EMDataFlags is deprecated. More...
 

Private Member Functions

void update_stat () const
 
void save_byteorder_to_dict (ImageIO *imageio)
 

Private Attributes

Dict attr_dict
 to store all image header info More...
 
float * rdata
 image real data More...
 
float * supp
 supplementary data array More...
 
int flags
 CTF data All CTF data become attribute ctf(vector<float>) in attr_dict –Grant Tang. More...
 
int changecount
 
int nx
 image size More...
 
int ny
 
int nz
 
int nxy
 
size_t nxyz
 
int xoff
 array index offsets More...
 
int yoff
 
int zoff
 
Vec3f all_translation
 translation from the original location More...
 
string path
 
int pathnum
 
EMData * rot_fp
 This is a cached rotational footprint, can save much time. More...
 

Friends

class GLUtil
 

Detailed Description

EMData stores an image's data and defines core image processing routines.

The image is 1D, 2D or 3D, in real space or fourier space (complex image). 

Data are ordered with x increasing fastest, then y, then z.

Definition at line 81 of file emdata.h.

Member Enumeration Documentation

◆ EMDataFlags

enum EMAN::EMData::EMDataFlags
private

This EMDataFlags is deprecated.

For anything which is currently handled by setting a bit in 'flags', instead, set or read an appropriately named attribute in the attributes dictionary. While there is a small overhead in the string lookup, none of these things should be called in the innermost loop anywhere, so it should be fine. –Grant

Enumerator
EMDATA_BUSY 
EMDATA_HASCTFF 
EMDATA_NEEDUPD 
EMDATA_FLIP 
EMDATA_PAD 
EMDATA_FFTODD 
EMDATA_SHUFFLE 
EMDATA_FH 
EMDATA_CPU_NEEDS_UPDATE 
EMDATA_GPU_NEEDS_UPDATE 
EMDATA_GPU_RO_NEEDS_UPDATE 

Definition at line 811 of file emdata.h.

811 {
812// EMDATA_COMPLEX = 1 << 1,
813// EMDATA_RI = 1 << 2, // real/imaginary or amp/phase
814 EMDATA_BUSY = 1 << 3, // someone is modifying data, NO LONGER USED
815 EMDATA_HASCTFF = 1 << 4, // has CTF info in the image file
816 EMDATA_NEEDUPD = 1 << 5, // needs a real update
817// EMDATA_COMPLEXX = 1 << 6, // 1D fft's in X
818 EMDATA_FLIP = 1 << 7, // is the image flipped
819 EMDATA_PAD = 1 << 8, // is the image fft padded
820 EMDATA_FFTODD = 1 << 9, // is the (real-space) nx odd
821 EMDATA_SHUFFLE = 1 << 10, // fft been shuffled? (so O is centered) PRB
822 EMDATA_FH = 1 << 11, // is the complex image a FH image
823 EMDATA_CPU_NEEDS_UPDATE = 1 << 12, // CUDA related: is the CPU version of the image out out data
824 EMDATA_GPU_NEEDS_UPDATE = 1 << 13, // CUDA related: is the GPU version of the image out out data
825 EMDATA_GPU_RO_NEEDS_UPDATE = 1 << 14 // // CUDA related: is the GPU RO version of the image out out data
826 };
EMAN::EMData::EMDATA_GPU_RO_NEEDS_UPDATE
@ EMDATA_GPU_RO_NEEDS_UPDATE
Definition: emdata.h:825
EMAN::EMData::EMDATA_NEEDUPD
@ EMDATA_NEEDUPD
Definition: emdata.h:816
EMAN::EMData::EMDATA_SHUFFLE
@ EMDATA_SHUFFLE
Definition: emdata.h:821
EMAN::EMData::EMDATA_HASCTFF
@ EMDATA_HASCTFF
Definition: emdata.h:815
EMAN::EMData::EMDATA_PAD
@ EMDATA_PAD
Definition: emdata.h:819
EMAN::EMData::EMDATA_GPU_NEEDS_UPDATE
@ EMDATA_GPU_NEEDS_UPDATE
Definition: emdata.h:824
EMAN::EMData::EMDATA_CPU_NEEDS_UPDATE
@ EMDATA_CPU_NEEDS_UPDATE
Definition: emdata.h:823
EMAN::EMData::EMDATA_BUSY
@ EMDATA_BUSY
Definition: emdata.h:814
EMAN::EMData::EMDATA_FLIP
@ EMDATA_FLIP
Definition: emdata.h:818
EMAN::EMData::EMDATA_FFTODD
@ EMDATA_FFTODD
Definition: emdata.h:820

Constructor & Destructor Documentation

◆ EMData() [1/3]

EMData::EMData ( )

For all image I/O.

For anything read/set image's information For modular class functions, process, align, etc. For fft, wavelet, insert data For get/set values, basic math operations, operators This is the header of EMData stay in sparx directory This is the header of EMData stay in sphire directory Construct an empty EMData instance. It has no image data.

Definition at line 70 of file emdata.cpp.

70 :
71#ifdef EMAN2_USING_CUDA
72 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
73#endif //EMAN2_USING_CUDA
74#ifdef FFT_CACHING
75 fftcache(0),
76#endif //FFT_CACHING
77 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0),
78 zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
79
80{
81 ENTERFUNC;
82
83 attr_dict["apix_x"] = 1.0f;
84 attr_dict["apix_y"] = 1.0f;
85 attr_dict["apix_z"] = 1.0f;
86
87 attr_dict["is_complex"] = int(0);
88 attr_dict["is_complex_x"] = int(0);
89 attr_dict["is_complex_ri"] = int(1);
90
91 attr_dict["datatype"] = (int)EMUtil::EM_FLOAT;
92
93 EMData::totalalloc++;
94#ifdef MEMDEBUG2
95 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
96#endif
97
98 EXITFUNC;
99}
EMAN::EMData::pathnum
int pathnum
Definition: emdata.h:858
EMAN::EMData::rot_fp
EMData * rot_fp
This is a cached rotational footprint, can save much time.
Definition: emdata.h:861
EMAN::EMData::yoff
int yoff
Definition: emdata.h:851
EMAN::EMData::totalalloc
static int totalalloc
Definition: emdata.h:803
EMAN::EMData::supp
float * supp
supplementary data array
Definition: emdata.h:837
EMAN::EMData::changecount
int changecount
Definition: emdata.h:846
EMAN::EMData::rdata
float * rdata
image real data
Definition: emdata.h:835
EMAN::EMData::attr_dict
Dict attr_dict
to store all image header info
Definition: emdata.h:833
EMAN::EMData::path
string path
Definition: emdata.h:857
EMAN::EMData::nx
int nx
image size
Definition: emdata.h:848
EMAN::EMData::flags
int flags
CTF data All CTF data become attribute ctf(vector<float>) in attr_dict –Grant Tang.
Definition: emdata.h:844
EMAN::EMData::all_translation
Vec3f all_translation
translation from the original location
Definition: emdata.h:854
EMAN::EMData::nxyz
size_t nxyz
Definition: emdata.h:849
EMAN::EMData::xoff
int xoff
array index offsets
Definition: emdata.h:851
EMAN::EMData::zoff
int zoff
Definition: emdata.h:851
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49

References attr_dict, EMAN::EMUtil::EM_FLOAT, ENTERFUNC, EXITFUNC, and totalalloc.

Referenced by absi(), bispecRotTransInvDirect(), bispecRotTransInvN(), calc_ccfx(), calc_fast_sigma_image(), calc_flcf(), compute_missingwedge(), copy_head(), do_radon(), EMData(), extract_box(), get_circle_mean(), get_clip(), get_col(), get_rotated_clip(), get_row(), get_top_half(), make_footprint(), make_rotational_footprint(), make_rotational_footprint_cmc(), make_rotational_footprint_e1(), oneDfftPolar(), operator=(), phase(), read_binedimage(), real2complex(), unwrap(), and unwrap_largerR().

◆ ~EMData()

EMData::~EMData ( )

Definition at line 331 of file emdata.cpp.

332{
333 ENTERFUNC;
334#ifdef FFT_CACHING
335 if (fftcache!=0) { delete fftcache; fftcache=0;}
336#endif //FFT_CACHING
337 free_memory();
338
339#ifdef EMAN2_USING_CUDA
340 if(cudarwdata){rw_free();}
341 if(cudarodata){ro_free();}
342#endif // EMAN2_USING_CUDA
343 EMData::totalalloc--;
344#ifdef MEMDEBUG2
345 printf("EMDATA- %4d %p\n",EMData::totalalloc,this);
346#endif
347 EXITFUNC;
348}
free_memory
void free_memory()
Free memory associated with this EMData Called in destructor and in assignment operator.

References ENTERFUNC, EXITFUNC, free_memory(), and totalalloc.

◆ EMData() [2/3]

EMAN::EMData::EMData ( const string &  filename,
int  image_index = 0 
)
explicit

Construct from an image file.

Parameters
filenamethe image file name
image_indexthe image index for stack image file, default 0

◆ EMData() [3/3]

EMData::EMData ( float *  data,
const int  nx,
const int  ny,
const int  nz,
const Dict &  attr_dict = Dict() 
)

Construction from a data pointer, dimensions must be supplied.

Takes possession of the pointer. data pointer must be allocated using malloc!

Parameters
dataa pointer to the pixel data which is stored in memory. Takes possession
nxthe number of pixels in the x direction
nythe number of pixels in the y direction
nzthe number of pixels in the z direction
attr_dictattribute dictionary for this image

Definition at line 276 of file emdata.cpp.

276 :
277#ifdef EMAN2_USING_CUDA
278 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
279#endif //EMAN2_USING_CUDA
280#ifdef FFT_CACHING
281 fftcache(0),
282#endif //FFT_CACHING
283 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
284 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
285{
286 ENTERFUNC;
287 // used to replace cube 'pixel'
288 attr_dict["apix_x"] = 1.0f;
289 attr_dict["apix_y"] = 1.0f;
290 attr_dict["apix_z"] = 1.0f;
291
292 EMData::totalalloc++;
293#ifdef MEMDEBUG2
294 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
295#endif
296
297 update();
298 EXITFUNC;
299}
update
void update()
Mark EMData as changed, statistics, etc will be updated at need.
Definition: emdata_metadata.h:138
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517

References attr_dict, ENTERFUNC, EXITFUNC, totalalloc, and update().

Member Function Documentation

◆ add_incoherent()

void EMData::add_incoherent ( EMData *  obj)

Adds 'obj' to 'this' incoherently.

'obj' and 'this' should be same size. Both images should be complex images.

Parameters
objThe image added to 'this' image.
Exceptions
ImageFormatExceptionIf the 2 images are not same size; or if the 2 images are not complex images.

Definition at line 3364 of file emdata.cpp.

3365{
3366 ENTERFUNC;
3367
3368 if (!obj) {
3369 LOGERR("NULL image");
3370 throw NullPointerException("NULL image");
3371 }
3372
3373 if (!obj->is_complex() || !is_complex()) {
3374 throw ImageFormatException("complex images only");
3375 }
3376
3377 if (!EMUtil::is_same_size(this, obj)) {
3378 throw ImageFormatException("images not same size");
3379 }
3380
3381 ri2ap();
3382 obj->ri2ap();
3383
3384 float *dest = get_data();
3385 float *src = obj->get_data();
3386 size_t size = (size_t)nx * ny * nz;
3387 for (size_t j = 0; j < size; j += 2) {
3388 dest[j] = (float) hypot(src[j], dest[j]);
3389 dest[j + 1] = 0;
3390 }
3391
3392 obj->update();
3393 update();
3394 EXITFUNC;
3395}
EMAN::EMUtil::is_same_size
static bool is_same_size(const EMData *image1, const EMData *image2)
Check whether two EMData images are of the same size.
Definition: emutil.cpp:1224
is_complex
bool is_complex() const
Is this a complex image?
Definition: emdata_metadata.h:648
get_data
float * get_data() const
Get the image pixel density data in a 1D float array.
Definition: emdata_metadata.h:78
ri2ap
void ri2ap()
convert the complex image from real/imaginary to amplitude/phase
Definition: emdata_transform.cpp:999
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
NullPointerException
#define NullPointerException(desc)
Definition: exception.h:241
LOGERR
#define LOGERR
Definition: log.h:51

References ENTERFUNC, EXITFUNC, get_data(), ImageFormatException, is_complex(), EMAN::EMUtil::is_same_size(), LOGERR, NullPointerException, nx, ny, nz, ri2ap(), and update().

◆ apply_radial_func()

void EMData::apply_radial_func ( float  x0,
float  dx,
vector< float >  array,
bool  interp = true 
)

multiplies by a radial function in fourier space.

Parameters
x0starting point x coordinate.
dxstep of x.
arrayradial function data array.
interpDo the interpolation or not.

Definition at line 2677 of file emdata.cpp.

2678{
2679 ENTERFUNC;
2680
2681 if (!is_complex()) throw ImageFormatException("apply_radial_func requires a complex image");
2682
2683 int n = (int)(array.size());
2684
2685 if (n*step>2.0) printf("Warning, apply_radial_func takes x0 and step with respect to Nyquist (0.5)\n");
2686
2687// printf("%f %f %f\n",array[0],array[25],array[50]);
2688
2689 ap2ri();
2690
2691 size_t ndims = get_ndim();
2692 float * data = get_data();
2693 if (ndims == 2) {
2694 int k = 0;
2695 for (int j = 0; j < ny; j++) {
2696 for (int i = 0; i < nx; i += 2, k += 2) {
2697 float r;
2698 if (j<ny/2) r = (float)hypot(i/(float)(nx*2), j/(float)ny);
2699 else r = (float)hypot(i/(float)(nx*2), (ny-j)/(float)ny);
2700 r = (r - x0) / step;
2701
2702 int l = 0;
2703 if (interp) {
2704 l = (int) floor(r);
2705 }
2706 else {
2707 l = (int) floor(r + 1);
2708 }
2709
2710
2711 float f = 0;
2712 if (l >= n - 2) {
2713 f = array[n - 1];
2714 }
2715 else {
2716 if (interp) {
2717 r -= l;
2718 f = (array[l] * (1.0f - r) + array[l + 1] * r);
2719 }
2720 else {
2721 f = array[l];
2722 }
2723 }
2724
2725 data[k] *= f;
2726 data[k + 1] *= f;
2727 }
2728 }
2729 }
2730 else if (ndims == 3) {
2731 int k = 0;
2732 for (int m = 0; m < nz; m++) {
2733 float mnz;
2734 if (m<nz/2) mnz=m*m/(float)(nz*nz);
2735 else { mnz=(nz-m)/(float)nz; mnz*=mnz; }
2736
2737 for (int j = 0; j < ny; j++) {
2738 float jny;
2739 if (j<ny/2) jny= j*j/(float)(ny*ny);
2740 else { jny=(ny-j)/(float)ny; jny*=jny; }
2741
2742 for (int i = 0; i < nx; i += 2, k += 2) {
2743 float r = std::sqrt((i * i / (nx*nx*4.0f)) + jny + mnz);
2744 r = (r - x0) / step;
2745
2746 int l = 0;
2747 if (interp) {
2748 l = (int) floor(r);
2749 }
2750 else {
2751 l = (int) floor(r + 1);
2752 }
2753
2754 float f = 0;
2755 if (l >= n - 2) {
2756 f = array[n - 1];
2757 }
2758 else {
2759 if (interp) {
2760 r -= l;
2761 f = (array[l] * (1.0f - r) + array[l + 1] * r);
2762 }
2763 else {
2764 f = array[l];
2765 }
2766 }
2767
2768//if (k%5000==0) printf("%d %d %d %f\n",i,j,m,f);
2769 data[k] *= f;
2770 data[k + 1] *= f;
2771 }
2772 }
2773 }
2774
2775 }
2776
2777 update();
2778 EXITFUNC;
2779}
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
get_ndim
int get_ndim() const
Get image dimension.
Definition: emdata_metadata.h:596
ap2ri
void ap2ri()
convert the complex image from amplitude/phase to real/imaginary

References ap2ri(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageFormatException, is_complex(), nx, ny, nz, sqrt(), and update().

Referenced by EMAN::OptVarianceCmp::cmp(), EMAN::CtfSimProcessor::process(), EMAN::FSCFourierProcessor::process(), EMAN::FourierProcessor::process_inplace(), EMAN::FourierAnlProcessor::process_inplace(), EMAN::SNRProcessor::process_inplace(), and EMAN::FileFourierProcessor::process_inplace().

◆ calc_az_dist()

vector< float > EMData::calc_az_dist ( int  n,
float  a0,
float  da,
float  rmin,
float  rmax 
)

Caculates the azimuthal distributions.

works for real or complex images, 2D only.

Parameters
nNumber of elements.
a0Starting angle.
daAngle step.
rminMinimum radius.
rmaxMaximum radius.
Exceptions
ImageDimensionExceptionIf image is 3D.
Returns
Float array to store the data.

Definition at line 2472 of file emdata.cpp.

2473{
2474 ENTERFUNC;
2475
2476 a0=a0*M_PI/180.0f;
2477 da=da*M_PI/180.0f;
2478
2479 if (get_ndim() > 2) {
2480 throw ImageDimensionException("no 3D image");
2481 }
2482
2483 float *yc = new float[n];
2484
2485 vector<float> vd(n);
2486 for (int i = 0; i < n; i++) {
2487 vd[i]=0;
2488 yc[i] = 0.00001f;
2489 }
2490
2491 int isri=is_ri();
2492
2493 float * data = get_data();
2494 if (is_complex()) {
2495 int c = 0;
2496 for (int y = 0; y < ny; y++) {
2497 for (int x = 0; x < nx; x += 2, c += 2) {
2498 int x1 = x / 2;
2499 int y1 = y<ny/2?y:y-ny;
2500 float r = (float)Util::hypot_fast(x1,y1);
2501
2502 if (r >= rmin && r <= rmax) {
2503 float a = 0;
2504
2505 if (y != ny / 2 || x != 0) {
2506 a = (atan2((float)y1, (float)x1) - a0) / da;
2507 }
2508
2509 int i = (int)(floor(a));
2510 a -= i;
2511
2512 float h=0;
2513 if (isri) {
2514 h = (float)hypot(data[c], data[c + 1]);
2515 }
2516 else h = data[c];
2517 if (i == 0) {
2518 vd[0] += h * (1.0f - a);
2519 yc[0] += (1.0f - a);
2520 }
2521 else if (i == n - 1) {
2522 vd[n - 1] += h * a;
2523 yc[n - 1] += a;
2524 }
2525 else if (i > 0 && i < (n - 1)) {
2526 vd[i] += h * (1.0f - a);
2527 yc[i] += (1.0f - a);
2528 vd[i + 1] += h * a;
2529 yc[i + 1] += a;
2530 }
2531 }
2532 }
2533 }
2534 }
2535 else {
2536 int c = 0;
2537 float half_nx = (nx - 1) / 2.0f;
2538 float half_ny = (ny - 1) / 2.0f;
2539
2540 for (int y = 0; y < ny; y++) {
2541 for (int x = 0; x < nx; x++, c++) {
2542 float y1 = y - half_ny;
2543 float x1 = x - half_nx;
2544 float r = (float)hypot(x1, y1);
2545
2546 if (r >= rmin && r <= rmax) {
2547 float a = 0;
2548 if (x1 != 0 || y1 != 0) {
2549 a = atan2(y1, x1);
2550 if (a < 0) {
2551 a += M_PI * 2;
2552 }
2553 }
2554
2555 a = (a - a0) / da;
2556 int i = (int)(floor(a));
2557 a -= i;
2558
2559 if (i == 0) {
2560 vd[0] += data[c] * (1.0f - a);
2561 yc[0] += (1.0f - a);
2562 }
2563 else if (i == n - 1) {
2564 vd[n - 1] += data[c] * a;
2565 yc[n - 1] += (a);
2566 }
2567 else if (i > 0 && i < (n - 1)) {
2568 vd[i] += data[c] * (1.0f - a);
2569 yc[i] += (1.0f - a);
2570 vd[i + 1] += data[c] * a;
2571 yc[i + 1] += a;
2572 }
2573 }
2574 }
2575 }
2576 }
2577
2578
2579 for (int i = 0; i < n; i++) {
2580 if (vd[i]<0||yc[i]<0) printf("%d vd %f yc %f\n",i,vd[i],yc[i]);
2581 vd[i] /= yc[i];
2582 }
2583
2584 if( yc )
2585 {
2586 delete[]yc;
2587 yc = 0;
2588 }
2589
2590 return vd;
2591
2592 EXITFUNC;
2593}
EMAN::Util::hypot_fast
static float hypot_fast(int x, int y)
Euclidean distance in 2D for integers computed fast using a cached lookup table.
Definition: util.cpp:742
is_ri
bool is_ri() const
Is this image a real/imaginary format complex image?
Definition: emdata_metadata.h:797
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166

References ENTERFUNC, EXITFUNC, get_data(), get_ndim(), EMAN::Util::hypot_fast(), ImageDimensionException, is_complex(), is_ri(), nx, ny, x, and y.

◆ calc_ccf_masked()

EMData * EMData::calc_ccf_masked ( EMData *  with,
EMData *  withsq = 0,
EMData *  mask = 0 
)

Calculate cross correlation between this and with where this is assumed to have had a mask applied to it (zero value in masked-out regions).

The mask will be used to compute a local normalization to prevent regions of high brightness from being preferentially aligned to the reference. Note that this local normalization is applied only to one of the two images, so the correlation values remain on arbitrary scale. If FFT of the square of with and the fft of the mask have been precalculated, passing them in will save time. Arguments should be either all-real or all comples. If mask is not provided, it will be computed from this based on zero values

Definition at line 1600 of file emdata.cpp.

1601{
1602 if ((withsq==0 && mask!=0)||(withsq!=0 && mask==0))
1603 throw NullPointerException("calc_ccf_masked error, both or neither of withsq and mask must be specified");
1604
1605 bool needfree=0;
1606 if (withsq==0) {
1607 withsq=with->process("math.squared");
1608 mask=this->process("threshold.notzero");
1609 needfree=1;
1610 }
1611
1612 EMData *c1=this->calc_ccf(with);
1613 EMData *c2=mask->calc_ccf(withsq);
1614
1615
1616 c2->process_inplace("math.reciprocal",Dict("zero_to",0.0f));
1617 c2->process_inplace("math.sqrt");
1618 c1->mult(*c2);
1619 delete c2;
1620
1621 if (needfree) {
1622 delete withsq;
1623 delete mask;
1624 }
1625 return c1;
1626}
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
process
EMData * process(const string &processorname, const Dict &params=Dict()) const
Apply a processor with its parameters on a copy of this image, return result as a a new image.
EMAN::EMData::calc_ccf
EMData * calc_ccf(EMData *with=0, fp_flag fpflag=CIRCULANT, bool center=false)
Calculate Cross-Correlation Function (CCF).
Definition: emdata.cpp:1499

References calc_ccf(), NullPointerException, and process().

Referenced by EMAN::RT3DLocalTreeAligner::testort(), and EMAN::RT3DLocalTreeAligner::xform_align_nbest().

◆ calc_dist()

float EMData::calc_dist ( EMData *  second_img,
int  y_index = 0 
) const

Calculates the distance between 2 vectors.

'this' image is 1D, which contains a vector; 'second_img' may be nD. One of its row is used as the second vector. 'second_img' and 'this' must have the same x size.

Parameters
second_imgThe image used to caculate the distance.
y_indexSpecifies which row in 'second_img' is used to do the caculation.
Exceptions
ImageDimensionExceptionIf 'this' image is not 1D.
ImageFormatExceptionIf the 2 images don't have same xsize.
Returns
The distance between 2 vectors.

Definition at line 3398 of file emdata.cpp.

3399{
3400 ENTERFUNC;
3401
3402 if (get_ndim() != 1) {
3403 throw ImageDimensionException("'this' image is 1D only");
3404 }
3405
3406 if (second_img->get_xsize() != nx || ny != 1) {
3407 throw ImageFormatException("image xsize not same");
3408 }
3409
3410 if (y_index > second_img->get_ysize() || y_index < 0) {
3411 return -1;
3412 }
3413
3414 float ret = 0;
3415 float *d1 = get_data();
3416 float *d2 = second_img->get_data() + second_img->get_xsize() * y_index;
3417
3418 for (int i = 0; i < nx; i++) {
3419 ret += Util::square(d1[i] - d2[i]);
3420 }
3421 EXITFUNC;
3422 return std::sqrt(ret);
3423}
EMAN::Util::square
static int square(int n)
Calculate a number's square.
Definition: util.h:736

References ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, ImageFormatException, nx, ny, sqrt(), and EMAN::Util::square().

◆ calc_fast_sigma_image()

EMData * EMData::calc_fast_sigma_image ( EMData *  mask)

Calculates the local standard deviation (sigma) image using the given mask image.

The mask image is typically much smaller than this image, and consists of ones, or is a small circle consisting of ones. The extent of the non zero neighborhood explicitly defines the range over which the local standard deviation is determined. Fourier convolution is used to do the math, ala Roseman (2003, Ultramicroscopy) However, Roseman was just working on methods Van Heel had presented earlier. The normalize flag causes the mask image to be processed so that it has a unit sum. Works in 1,2 and 3D

Parameters
maskthe image that will be used to define the neighborhood for determine the local standard deviation
Returns
the sigma image, the phase origin is at the corner (not the center)
Exceptions
ImageDimensionExceptionif the dimensions of with do not match those of this
ImageDimensionExceptionif any of the dimensions sizes of with exceed of this image's.
Author
David Woolford
Date
April 2008

Definition at line 3426 of file emdata.cpp.

3427{
3428 ENTERFUNC;
3429
3430 bool maskflag = false;
3431 if (mask == 0) {
3432 mask = new EMData(nx,ny,nz);
3433 mask->process_inplace("testimage.circlesphere");
3434 maskflag = true;
3435 }
3436
3437 if (get_ndim() != mask->get_ndim() ) throw ImageDimensionException("The dimensions do not match");
3438
3439 int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
3440
3441 if ( mnx > nx || mny > ny || mnz > nz)
3442 throw ImageDimensionException("Can not calculate variance map using an image that is larger than this image");
3443
3444 size_t P = 0;
3445 for(size_t i = 0; i < mask->get_size(); ++i){
3446 if (mask->get_value_at(i) != 0){
3447 ++P;
3448 }
3449 }
3450 float normfac = 1.0f/(float)P;
3451
3452// bool undoclip = false;
3453
3454 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
3455// if ( mnx < nx || mny < ny || mnz < nz) {
3456 Region r;
3457 if (ny == 1) r = Region((mnx-nxc)/2,nxc);
3458 else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
3459 else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3460 mask->clip_inplace(r,0.0);
3461 //Region r((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3462 //mask->clip_inplace(r);
3463 //undoclip = true;
3464 //}
3465
3466 // Here we generate the local average of the squares
3467 Region r2;
3468 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
3469 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
3470 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
3471 EMData* squared = get_clip(r2,get_edge_mean());
3472
3473 EMData* tmp = squared->copy();
3474 Dict pow;
3475 pow["pow"] = 2.0f;
3476 squared->process_inplace("math.pow",pow);
3477 EMData* s = mask->convolute(squared);//ming, mask squared exchange
3478 squared->mult(normfac);
3479
3480 EMData* m = mask->convolute(tmp);//ming, tmp mask exchange
3481 m->mult(normfac);
3482 m->process_inplace("math.pow",pow);
3483 delete tmp; tmp = 0;
3484 s->sub(*m);
3485 // Here we finally generate the standard deviation image
3486 s->process_inplace("math.sqrt");
3487
3488// if ( undoclip ) {
3489// Region r((nx-mnx)/2, (ny-mny)/2, (nz-mnz)/2,mnx,mny,mnz);
3490// mask->clip_inplace(r);
3491// }
3492
3493 if (maskflag) {
3494 delete mask;
3495 mask = 0;
3496 } else {
3497 Region r;
3498 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
3499 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
3500 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
3501 mask->clip_inplace(r);
3502 }
3503
3504 delete squared;
3505 delete m;
3506
3507 s->process_inplace("xform.phaseorigin.tocenter");
3508 Region r3;
3509 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
3510 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
3511 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
3512 s->clip_inplace(r3);
3513 EXITFUNC;
3514 return s;
3515}
EMAN::EMData::EMData
EMData()
For all image I/O.
Definition: emdata.cpp:70
EMAN::EMData::get_clip
EMData * get_clip(const Region &area, const float fill=0) const
Get an inclusive clip.
Definition: emdata.cpp:592
EMAN::EMData::convolute
EMData * convolute(EMData *with)
Convolutes 2 data sets.
Definition: emdata.cpp:3567
EMAN::EMData::clip_inplace
void clip_inplace(const Region &area, const float &fill_value=0)
Clip the image inplace - clipping region must be smaller than the current region internally memory is...
Definition: emdata.cpp:350
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
get_edge_mean
float get_edge_mean() const
Calculates the mean pixel values around the (1 pixel) edge of the image.
Definition: emdata_metadata.cpp:713

References clip_inplace(), convolute(), EMData(), ENTERFUNC, EXITFUNC, get_clip(), get_edge_mean(), get_ndim(), ImageDimensionException, nx, ny, and nz.

Referenced by calc_flcf().

◆ calc_flcf()

EMData * EMData::calc_flcf ( EMData *  with)

Calculates the cross correlation with local normalization between 2 images.

This is a faster version of local correlation that make use of Fourier convolution and correlation. With is the template - the thing that you wish to find in the this image. It should not be unecessarily padded. The function uses the size of with to determine the extent of the local neighborhood used in the local normalization (for technical details, see calc_fast_sigma_image). Note that this function circularly masks the template at its radius so the calling function need not do this beforehand. Works in 1,2 and 3D.

Parameters
withThe image used to calculate cross correlation (the template)
Returns
the local normalized cross correlation image - the phase origin is at the corner of the image
Author
David Woolford
Date
April 2008

Definition at line 3519 of file emdata.cpp.

3520{
3521 ENTERFUNC;
3522 EMData *this_copy=this;
3523 this_copy=copy();
3524
3525 int mnx = with->get_xsize(); int mny = with->get_ysize(); int mnz = with->get_zsize();
3526 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
3527
3528 // Ones is a circular/spherical mask, consisting of 1s.
3529 EMData* ones = new EMData(mnx,mny,mnz);
3530 ones->process_inplace("testimage.circlesphere");
3531
3532 // Get a copy of with, we will eventually resize it
3533 EMData* with_resized = with->copy();
3534 with_resized->process_inplace("normalize");
3535 with_resized->mult(*ones);
3536
3537 EMData* s = calc_fast_sigma_image(ones);// Get the local sigma image
3538
3539 Region r1;
3540 if (ny == 1) r1 = Region((mnx-nxc)/2,nxc);
3541 else if (nz == 1) r1 = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
3542 else r1 = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3543 with_resized->clip_inplace(r1,0.0);
3544
3545 Region r2;
3546 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
3547 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
3548 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
3549 this_copy->clip_inplace(r2,0.0);
3550
3551 EMData* corr = this_copy->calc_ccf(with_resized); // the ccf results should have same size as sigma
3552
3553 corr->process_inplace("xform.phaseorigin.tocenter");
3554 Region r3;
3555 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
3556 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
3557 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
3558 corr->clip_inplace(r3);
3559
3560 corr->div(*s);
3561
3562 delete with_resized; delete ones; delete this_copy; delete s;
3563 EXITFUNC;
3564 return corr;
3565}
EMAN::EMData::calc_fast_sigma_image
EMData * calc_fast_sigma_image(EMData *mask)
Calculates the local standard deviation (sigma) image using the given mask image.
Definition: emdata.cpp:3426
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....

References calc_ccf(), calc_fast_sigma_image(), clip_inplace(), copy(), EMData(), ENTERFUNC, EXITFUNC, nx, ny, and nz.

◆ calc_hist()

vector< float > EMData::calc_hist ( int  hist_size = 128,
float  hist_min = 0,
float  hist_max = 0,
const float &  brt = 0.0f,
const float &  cont = 1.0f 
)

Calculates the histogram of 'this' image.

The result is stored in float array 'hist'. If hist_min = hist_max, use image data min as hist_min; use image data max as hist_max.

Parameters
hist_sizeHistogram array's size.
hist_minMinimum histogram value.
hist_maxMaximum histogram value.
brt
cont
Returns
histogram array of this image.

Definition at line 2365 of file emdata.cpp.

2366{
2367 ENTERFUNC;
2368
2369// static size_t prime[] = { 1, 3, 7, 11, 17, 23, 37, 59, 127, 253, 511 };
2370
2371 size_t size = (size_t)nx * ny * nz;
2372 if (histmin == histmax) {
2373 histmin = get_attr("minimum");
2374 histmax = get_attr("maximum");
2375 }
2376 int n=hist_size;
2377 float w = (float)(n-1) / (histmax - histmin);
2378 float * data = get_data();
2379
2380 vector <float> hist(hist_size, 0.0);
2381
2382 if (cont != 1.0f || brt != 0) {
2383 for(size_t i=0; i<size; i++) {
2384 float val = cont*(data[i]+brt);
2385 int j = Util::round((val - histmin) * w);
2386
2387 // Outliers now go in the edge bins
2388 if (j>=n) j=n-1;
2389 if (j<0) j=0;
2390 hist[j]+=1;
2391 }
2392 }
2393 else {
2394 for(size_t i=0; i<size; i++) {
2395 float val = data[i];
2396 int j = Util::round((val - histmin) * w);
2397 if (j>=n) j=n-1;
2398 if (j<0) j=0;
2399 hist[j]+=1;
2400 }
2401 }
2402
2403// 20 years ago, it was expensive to compute histograms of large images
2404// so there was a complicated strategy to approximate them
2405// int p0 = 0;
2406// int p1 = 0;
2407// size_t size = (size_t)nx * ny * nz;
2408// if (size < 300000) {
2409// p0 = 0;
2410// p1 = 0;
2411// }
2412// else if (size < 2000000) {
2413// p0 = 2;
2414// p1 = 3;
2415// }
2416// else if (size < 8000000) {
2417// p0 = 4;
2418// p1 = 6;
2419// }
2420// else {
2421// p0 = 7;
2422// p1 = 9;
2423// }
2424//
2425// if (is_complex() && p0 > 0) {
2426// p0++;
2427// p1++;
2428// }
2429
2430// size_t di = 0;
2431// // float norm = 0;
2432// size_t n = hist.size();
2433//
2434// float * data = get_data();
2435// for (int k = p0; k <= p1; ++k) {
2436// if (is_complex()) {
2437// di = prime[k] * 2;
2438// }
2439// else {
2440// di = prime[k];
2441// }
2442//
2443// // norm += (float)size / (float) di;
2444// float w = (float)n / (histmax - histmin);
2445//
2446// for(size_t i=0; i<=size-di; i += di) {
2447// float val;
2448// if (cont != 1.0f || brt != 0)val = cont*(data[i]+brt);
2449// else val = data[i];
2450// int j = Util::round((val - histmin) * w);
2451// if (j >= 0 && j < (int) n) {
2452// hist[j] += 1;
2453// }
2454// }
2455// }
2456// /*
2457// for (size_t i = 0; i < hist.size(); ++i) {
2458// if (norm != 0) {
2459// hist[i] = hist[i] / norm;
2460// }
2461// }
2462// */
2463 return hist;
2464
2465 EXITFUNC;
2466}
EMAN::Util::round
static int round(float x)
Get ceiling round of a float number x.
Definition: util.h:465
get_attr
EMObject get_attr(const string &attr_name) const
The generic way to get any image header information given a header attribute name.
Definition: emdata_metadata.cpp:1006

References ENTERFUNC, EXITFUNC, get_attr(), get_data(), nx, ny, nz, and EMAN::Util::round().

Referenced by EMAN::ToCenterProcessor::process_inplace().

◆ calc_mutual_correlation()

EMData * EMData::calc_mutual_correlation ( EMData *  with,
bool  tocorner = false,
EMData *  filter = 0 
)

Calculates mutual correlation function (MCF) between 2 images.

If 'with' is NULL, this does mirror ACF.

Parameters
withThe image used to calculate MCF.
tocornerSet whether to translate the result image to the corner.
filterThe filter image used in calculating MCF.
Exceptions
ImageFormatExceptionIf 'with' is not NULL and it doesn't have the same size to 'this' image.
NullPointerExceptionIf FFT returns NULL image.
Returns
Mutual correlation function image.

Definition at line 2209 of file emdata.cpp.

2210{
2211 ENTERFUNC;
2212
2213 if (with && !EMUtil::is_same_size(this, with)) {
2214 LOGERR("images not same size");
2215 throw ImageFormatException( "images not same size");
2216 }
2217
2218#ifdef EMAN2_USING_CUDA
2219 if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata)
2220 {
2221
2222 EMData* this_fft = do_fft_cuda();
2223
2224 EMData *cf = 0;
2225 if (with && with != this) {
2226 cf = with->do_fft_cuda();
2227 }else{
2228 cf = this_fft->copy();
2229 }
2230
2231 if (filter) {
2232 if (!EMUtil::is_same_size(filter, cf)) {
2233 LOGERR("improperly sized filter");
2234 throw ImageFormatException("improperly sized filter");
2235 }
2236 mult_complex_efficient_cuda(cf->cudarwdata, filter->cudarwdata, cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), 1);
2237 mult_complex_efficient_cuda(this_fft->cudarwdata, filter->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize(), 1);
2238 }
2239
2240 mcf_cuda(this_fft->cudarwdata, cf->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize());
2241
2242 EMData *f2 = cf->do_ift_cuda();
2243
2244 if (tocenter) {
2245 f2->process_inplace("xform.phaseorigin.tocenter");
2246 }
2247
2248 if( cf )
2249 {
2250 delete cf;
2251 cf = 0;
2252 }
2253
2254 if( this_fft )
2255 {
2256 delete this_fft;
2257 this_fft = 0;
2258 }
2259
2260 f2->set_attr("label", "MCF");
2261 f2->set_path("/tmp/eman.mcf");
2262 f2->update();
2263
2264 EXITFUNC;
2265 return f2;
2266 }
2267#endif
2268
2269 EMData *this_fft = 0;
2270 this_fft = do_fft();
2271
2272 if (!this_fft) {
2273
2274 LOGERR("FFT returns NULL image");
2275 throw NullPointerException("FFT returns NULL image");
2276 }
2277
2278 this_fft->ap2ri(); //this is not needed!
2279 EMData *cf = 0;
2280
2281 if (with && with != this) {
2282 cf = with->do_fft();
2283 if (!cf) {
2284 LOGERR("FFT returns NULL image");
2285 throw NullPointerException("FFT returns NULL image");
2286 }
2287 cf->ap2ri(); //nor is this!
2288 }
2289 else {
2290 cf = this_fft->copy();
2291 }
2292
2293 if (filter) {
2294 if (!EMUtil::is_same_size(filter, cf)) {
2295 LOGERR("improperly sized filter");
2296 throw ImageFormatException("improperly sized filter");
2297 }
2298
2299 cf->mult_complex_efficient(*filter,true); //insanely this is required....
2300 this_fft->mult(*filter,true);
2301 //cf->mult_complex_efficient(*filter,7); // takes advantage of the fact that the filter is 1 everywhere but near the origin
2302 //this_fft->mult_complex_efficient(*filter,7);
2303 /*cf->mult_complex_efficient(*filter,5);
2304 this_fft->mult_complex_efficient(*filter,5);*/
2305 }
2306
2307 float *rdata1 = this_fft->get_data();
2308 float *rdata2 = cf->get_data();
2309 size_t this_fft_size = (size_t)this_fft->get_xsize() * this_fft->get_ysize() * this_fft->get_zsize();
2310
2311 if (with == this) {
2312 for (size_t i = 0; i < this_fft_size; i += 2) {
2313 rdata2[i] = std::sqrt(rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
2314 rdata2[i + 1] = 0;
2315 }
2316
2317 this_fft->update();
2318 cf->update();
2319 }
2320 else {
2321 for (size_t i = 0; i < this_fft_size; i += 2) {
2322 rdata2[i] = (rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
2323 rdata2[i + 1] = (rdata1[i + 1] * rdata2[i] - rdata1[i] * rdata2[i + 1]);
2324 }
2325
2326 //This seems like a bug, but it probably is never used....
2327 for (size_t i = 0; i < this_fft_size; i += 2) {
2328 float t = Util::square(rdata2[i]) + Util::square(rdata2[i + 1]);
2329 if (t != 0) {
2330 t = pow(t, 0.25f);
2331 rdata2[i] /= t;
2332 rdata2[i + 1] /= t;
2333 }
2334 }
2335 this_fft->update();
2336 cf->update();
2337 }
2338
2339 EMData *f2 = cf->do_ift();
2340
2341 if (tocenter) {
2342 f2->process_inplace("xform.phaseorigin.tocenter");
2343 }
2344
2345 if( cf )
2346 {
2347 delete cf;
2348 cf = 0;
2349 }
2350
2351 if( this_fft )
2352 {
2353 delete this_fft;
2354 this_fft = 0;
2355 }
2356
2357 f2->set_attr("label", "MCF");
2358 f2->set_path("/tmp/eman.mcf");
2359
2360 EXITFUNC;
2361 return f2;
2362}
mcf_cuda
void mcf_cuda(const float *data1, float *data2, const int nx, const int ny, const int nz)
mult_complex_efficient_cuda
void mult_complex_efficient_cuda(float *data, const float *src_data, const int nx, const int ny, const int nz, const int radius)
do_fft
EMData * do_fft() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
Definition: emdata_transform.cpp:55

References do_fft(), ENTERFUNC, EXITFUNC, ImageFormatException, EMAN::EMUtil::is_same_size(), LOGERR, mcf_cuda(), mult_complex_efficient_cuda(), NullPointerException, sqrt(), and EMAN::Util::square().

Referenced by make_rotational_footprint_cmc(), and make_rotational_footprint_e1().

◆ calc_radial_dist() [1/2]

vector< float > EMData::calc_radial_dist ( int  n,
float  x0,
float  dx,
int  inten 
)

calculates radial distribution.

works for real and imaginary images. Returns mean radial amplitude, or intensity if inten is set. Note that the complex origin is at (0,0), with periodic boundaries. Note that the inten option is NOT equivalent to returning amplitude and squaring the result.

Parameters
nnumber of points.
x0starting point x coordinate.
dxstep of x.
inten0 -> amp, 1-> inten (amp^2), 2->min, 3-> max, 4-> sigma or rmsd intensity from zero if complex, 5-> for 1/2 images along x in real space, origin at x=0, mean amplitude
Returns
The radial distribution in an array.

Definition at line 2781 of file emdata.cpp.

2782{
2783 ENTERFUNC;
2784
2785 vector<double>ret(n);
2786 vector<double>norm(n);
2787 vector<double>count(n);
2788
2789 int x,y,z,i;
2790 int step=is_complex()?2:1;
2791 int isinten=get_attr_default("is_intensity",0);
2792 int isri=is_ri();
2793
2794 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
2795
2796 switch (inten){
2797 case 0:
2798 case 1:
2799 case 4:
2800 case 5:
2801 for (i=0; i<n; i++) ret[i]=norm[i]=count[i]=0.0;
2802 break;
2803 case 2:
2804 for (i=0; i<n; i++) ret[i]=1.0e27;
2805 break;
2806 case 3:
2807 for (i=0; i<n; i++) ret[i]=-1.0e27;
2808 break;
2809 }
2810
2811 float * data = get_data();
2812
2813 // We do 2D separately to avoid the hypot3 call
2814 if (nz==1) {
2815 for (y=i=0; y<ny; y++) {
2816 for (x=0; x<nx; x+=step,i+=step) {
2817 float r,v;
2818 int f;
2819 if (step==2) { //complex
2820 if (x==0 && y>ny/2) continue;
2821 r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y)); // origin at 0,0; periodic
2822 r=(r-x0)/dx;
2823 f=int(r); // safe truncation, so floor isn't needed
2824 if (f<0 || f>=n) continue;
2825 switch (inten) {
2826 case 0:
2827 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2828 else v=data[i]; // amp/phase, just get amp
2829 break;
2830 case 1:
2831 if (isinten) v=data[i];
2832 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2833 else v=data[i]*data[i];
2834 break;
2835 case 2:
2836 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2837 else v=data[i];
2838 if (v<ret[f]) ret[f]=v;
2839 break;
2840 case 3:
2841 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2842 else v=data[i];
2843 if (v>ret[f]) ret[f]=v;
2844 break;
2845 case 4:
2846 if (isinten) v=data[i];
2847 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2848 else v=data[i]*data[i];
2849 ret[f]+=std::sqrt(v);
2850 norm[f]+=v;
2851 count[f]+=1.0;
2852 break;
2853 }
2854 }
2855 else {
2856 if (inten==5) r=(float)(Util::hypot_fast(x,y-ny/2));
2857 else r=(float)(Util::hypot_fast(x-nx/2,y-ny/2));
2858 r=(r-x0)/dx;
2859 f=int(r); // safe truncation, so floor isn't needed
2860 if (f<0 || f>=n) continue;
2861 switch (inten) {
2862 case 0:
2863 case 5:
2864 v=data[i];
2865 break;
2866 case 1:
2867 v=data[i]*data[i];
2868 break;
2869 case 2:
2870 if (data[i]<ret[f]) ret[f]=data[i];
2871 break;
2872 case 3:
2873 if (data[i]>ret[f]) ret[f]=data[i];
2874 break;
2875 case 4:
2876 ret[f]+=data[i];
2877 norm[f]+=data[i]*data[i];
2878 count[f]+=1.0;
2879 break;
2880 }
2881 }
2882
2883 if (inten<2||inten==5) {
2884 r-=float(f); // r is now the fractional spacing between bins
2885 // printf("%d\t%d\t%d\t%1.3f\t%d\t%1.3f\t%1.4g\n",x,y,f,r,step,Util::hypot_fast(x/2,y<ny/2?y:ny-y),v);
2886 ret[f]+=v*(1.0f-r);
2887 norm[f]+=(1.0f-r);
2888 if (f<n-1) {
2889 ret[f+1]+=v*r;
2890 norm[f+1]+=r;
2891 }
2892 }
2893 }
2894 }
2895 }
2896 else {
2897// FILE *out = fopen("x.txt","w");
2898 size_t i; //3D file may have >2G size
2899 for (z=i=0; z<nz; ++z) {
2900 for (y=0; y<ny; ++y) {
2901 for (x=0; x<nx; x+=step,i+=step) {
2902 float r,v;
2903 int f;
2904 if (step==2) { //complex
2905 if (x==0 && z>nz/2) continue;
2906 if (x==0 && z==nz/2 && y>ny/2) continue;
2907 r=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z); // origin at 0,0; periodic
2908 r=(r-x0)/dx;
2909 f=int(r); // safe truncation, so floor isn't needed
2910 if (f<0 || f>=n) continue;
2911 switch(inten) {
2912 case 0:
2913 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2914 else v=data[i]; // amp/phase, just get amp
2915 break;
2916 case 1:
2917 if (isinten) v=data[i];
2918 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2919 else v=data[i]*data[i];
2920 break;
2921 case 2:
2922 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2923 else v=data[i]; // amp/phase, just get amp
2924 if (v<ret[f]) ret[f]=v;
2925 break;
2926 case 3:
2927 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
2928 else v=data[i]; // amp/phase, just get amp
2929 if (v>ret[f]) ret[f]=v;
2930 break;
2931 case 4:
2932 if (isinten) v=data[i];
2933 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
2934 else v=data[i]*data[i];
2935 ret[f]+=std::sqrt(v);
2936 norm[f]+=v;
2937 count[f]+=1.0;
2938// if (f==100) fprintf(out,"%1.4g\n",std::sqrt(v));
2939 break;
2940 }
2941 }
2942 else {
2943 if (inten==5) r=Util::hypot3(x,y-ny/2,z-nz/2);
2944 else r=Util::hypot3(x-nx/2,y-ny/2,z-nz/2);
2945 r=(r-x0)/dx;
2946 f=int(r); // safe truncation, so floor isn't needed
2947 if (f<0 || f>=n) continue;
2948 switch(inten) {
2949 case 0:
2950 case 5:
2951 v=data[i];
2952 break;
2953 case 1:
2954 v=data[i]*data[i];
2955 break;
2956 case 2:
2957 if (data[i]<ret[f]) ret[f]=data[i];
2958 break;
2959 case 3:
2960 if (data[i]>ret[f]) ret[f]=data[i];
2961 break;
2962 case 4:
2963 ret[f]+=data[i];
2964 norm[f]+=data[i]*data[i];
2965 count[f]+=1.0;
2966 break;
2967 }
2968 }
2969
2970 if (inten<2||inten==5) {
2971// ret[f]+=v;
2972// norm[f]+=1.0;
2973
2974 r-=float(f); // r is now the fractional spacing between bins
2975 ret[f]+=v*(1.0f-r);
2976 norm[f]+=(1.0f-r);
2977 if (f<n-1) {
2978 ret[f+1]+=v*r;
2979 norm[f+1]+=r;
2980 }
2981 }
2982 }
2983 }
2984 }
2985// fclose(out);
2986 }
2987
2988 if (inten<2||inten==5) {
2989 for (i=0; i<n; i++) ret[i]/=(norm[i]==0?1.0f:norm[i]); // Normalize
2990 }
2991 else if (inten==4) {
2992 for (i=0; i<n; i++) {
2993 ret[i]/=count[i]; // becomes mean
2994 norm[i]/=count[i]; // avg amp^2
2995 ret[i]*=ret[i];
2996 if (norm[i]<=ret[i]) ret[i]=0.0;
2997 else ret[i]=std::sqrt(norm[i]-ret[i]); // sigma
2998 }
2999 }
3000
3001 EXITFUNC;
3002
3003 return vector<float>(ret.begin(),ret.end());
3004}
EMAN::Util::hypot3
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);.
Definition: util.h:827
get_attr_default
EMObject get_attr_default(const string &attr_name, const EMObject &em_obj=EMObject()) const
The generic way to get any image header information given a header attribute name.
InvalidParameterException
#define InvalidParameterException(desc)
Definition: exception.h:361

References ENTERFUNC, EXITFUNC, get_attr_default(), get_data(), EMAN::Util::hypot3(), EMAN::Util::hypot_fast(), InvalidParameterException, is_complex(), is_ri(), nx, ny, nz, sqrt(), x, and y.

Referenced by EMAN::TomoAverager::add_image(), EMAN::TomoWedgeFscCmp::cmp(), EMAN::OptVarianceCmp::cmp(), EMAN::FRCCmp::cmp(), EMAN::MatchSFProcessor::create_radial_func(), EMAN::CtfCAutoAverager::finish(), EMAN::CtfCWautoAverager::finish(), EMAN::FourierAnlProcessor::process_inplace(), EMAN::NonConvexProcessor::process_inplace(), EMAN::FFTPeakProcessor::process_inplace(), EMAN::WedgeFillProcessor::process_inplace(), EMAN::RotationalAverageProcessor::process_inplace(), EMAN::RotationalSubstractProcessor::process_inplace(), EMAN::CTFSNRWeightProcessor::process_inplace(), EMAN::RT3DTreeAligner::xform_align_nbest(), and EMAN::RT3DLocalTreeAligner::xform_align_nbest().

◆ calc_radial_dist() [2/2]

vector< float > EMData::calc_radial_dist ( int  n,
float  x0,
float  dx,
int  nwedge,
float  offset,
bool  inten 
)

calculates radial distribution subdivided by angle.

works for real and imaginary images. 2-D only. The first returns a single vector of n*nwedge points, with radius varying first. That is, the first n points represent the radial profile in the first wedge.

Parameters
nnumber of points.
x0starting x coordinate.
dxstep of x.
nwedgeint number of wedges to divide the circle into
nwedgefloat angular offset in radians for start of first bin
intenreturns intensity (amp^2) rather than amplitude if set
Exceptions
ImageDimensionExceptionIf 'this' image is not 2D.
Returns
nwedge radial distributions packed into a single vector<float>

Definition at line 3006 of file emdata.cpp.

3007{
3008 ENTERFUNC;
3009
3010 if (nz > 1) {
3011 LOGERR("2D images only.");
3012 throw ImageDimensionException("2D images only");
3013 }
3014 int isinten=get_attr_default("is_intensity",0);
3015
3016 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
3017
3018
3019 vector<double>ret(n*nwedge);
3020 vector<double>norm(n*nwedge);
3021
3022 int x,y,i;
3023 int isri = is_ri(); // this has become expensive!
3024 int step=is_complex()?2:1;
3025 float astep=(float)(M_PI*2.0/nwedge);
3026 if (is_complex()) astep/=2; // Since we only have the right 1/2 of Fourier space
3027 float* data = get_data();
3028 for (i=0; i<n*nwedge; i++) ret[i]=norm[i]=0.0;
3029
3030 // We do 2D separately to avoid the hypot3 call
3031 for (y=i=0; y<ny; y++) {
3032 for (x=0; x<nx; x+=step,i+=step) {
3033 float r,v,a;
3034 int bin;
3035 if (is_complex()) {
3036 r=(float)(hypot(x/2.0,y<ny/2?y:ny-y)); // origin at 0,0; periodic
3037 a=atan2(float(y<ny/2?y:y-ny),x/2.0f);
3038 if (!inten) {
3039 if (isri) v=(float)(hypot(data[i],data[i+1])); // real/imag, compute amplitude
3040 else v=data[i]; // amp/phase, just get amp
3041 } else {
3042 if (isinten) v=data[i];
3043 else if (isri) v=data[i]*data[i]+data[i+1]*data[i+1];
3044 else v=data[i]*data[i];
3045 }
3046 bin=n*int(floor((a+M_PI/2.0f+offset)/astep));
3047 }
3048 else {
3049 r=(float)(hypot(x-nx/2,y-ny/2));
3050 a=atan2(float(y-ny/2),float(x-nx/2));
3051 if (inten) v=data[i]*data[i];
3052 else v=data[i];
3053 bin=n*int(floor((a+M_PI+offset)/astep));
3054 }
3055 if (bin>=nwedge*n) bin-=nwedge*n;
3056 if (bin<0) bin+=nwedge*n;
3057 r=(r-x0)/dx;
3058 int f=int(r); // safe truncation, so floor isn't needed
3059 r-=float(f); // r is now the fractional spacing between bins
3060// printf("%d %d %d %d %1.3f %1.3f\n",x,y,bin,f,r,a);
3061 if (f>=0 && f<n) {
3062 ret[f+bin]+=v*(1.0f-r);
3063 norm[f+bin]+=(1.0f-r);
3064 if (f<n-1) {
3065 ret[f+1+bin]+=v*r;
3066 norm[f+1+bin]+=r;
3067 }
3068 }
3069 }
3070 }
3071
3072 for (i=0; i<n*nwedge; i++) ret[i]/=norm[i]?norm[i]:1.0f; // Normalize
3073 EXITFUNC;
3074
3075 return vector<float>(ret.begin(),ret.end());
3076}

References ENTERFUNC, EXITFUNC, get_attr_default(), get_data(), ImageDimensionException, InvalidParameterException, is_complex(), is_ri(), LOGERR, nx, ny, nz, x, and y.

◆ cconj()

void EMData::cconj ( )

Replace the image its complex conjugate.

Exceptions
ImageFormatExceptionImage must be complex (and RI)

Definition at line 3078 of file emdata.cpp.

3078 {
3079 ENTERFUNC;
3080 if (!is_complex() || !is_ri())
3081 throw ImageFormatException("EMData::conj requires a complex, ri image");
3082 int nxreal = nx -2 + int(is_fftodd());
3083 int nxhalf = nxreal/2;
3084 for (int iz = 0; iz < nz; iz++)
3085 for (int iy = 0; iy < ny; iy++)
3086 for (int ix = 0; ix <= nxhalf; ix++)
3087 cmplx(ix,iy,iz) = conj(cmplx(ix,iy,iz));
3088 EXITFUNC;
3089}
cmplx
std::complex< float > & cmplx(const int ix, const int iy, const int iz)
Return reference to complex elements.
Definition: emdata_core.h:781
is_fftodd
bool is_fftodd() const
Does this image correspond to a (real-space) odd nx?
Definition: emdata_metadata.h:864

References cmplx(), ENTERFUNC, EXITFUNC, ImageFormatException, is_complex(), is_fftodd(), is_ri(), nx, ny, and nz.

◆ clip_inplace()

void EMData::clip_inplace ( const Region &  area,
const float &  fill_value = 0 
)

Clip the image inplace - clipping region must be smaller than the current region internally memory is reallocated.

Exceptions
ImageDimensionExceptionif any of the dimensions of the argument region are negative
Parameters
areaThe clip area, can be 2D/3D.
fill_valuethe value that new region

Definition at line 350 of file emdata.cpp.

351{
352 // Added by d.woolford
353 ENTERFUNC;
354
355// printf("cip %d %d %d %d %d %d %f %d %d %d\n",area.origin[0],area.origin[1],area.origin[2],area.size[0],area.size[1],area.size[2],fill_value,nx,ny,nz);
356 // Store the current dimension values
357 int prev_nx = nx, prev_ny = ny, prev_nz = nz;
358 size_t prev_size = (size_t)nx*ny*nz;
359
360 // Get the zsize, ysize and xsize of the final area, these are the new dimension sizes of the pixel data
361 int new_nz = ( area.size[2]==0 ? 1 : (int)area.size[2]);
362 int new_ny = ( area.size[1]==0 ? 1 : (int)area.size[1]);
363 int new_nx = (int)area.size[0];
364
365 if ( new_nz < 0 || new_ny < 0 || new_nx < 0 )
366 {
367 // Negative image dimensions were never tested nor considered when creating this implementation
368 throw ImageDimensionException("New image dimensions are negative - this is not supported in the clip_inplace operation");
369 }
370
371 size_t new_size = (size_t)new_nz*new_ny*new_nx;
372
373 // Get the translation values, they are used to construct the ClipInplaceVariables object
374 int x0 = (int) area.origin[0];
375 int y0 = (int) area.origin[1];
376 int z0 = (int) area.origin[2];
377
378 // Get a object that calculates all the interesting variables associated with the clip inplace operation
379 ClipInplaceVariables civ(prev_nx, prev_ny, prev_nz, new_nx, new_ny, new_nz, x0, y0, z0);
380
381 get_data(); // Do this here to make sure rdata is up to date, applicable if GPU stuff is occurring
382 // Double check to see if any memory shifting even has to occur
383 if ( x0 > prev_nx || y0 > prev_ny || z0 > prev_nz || civ.x_iter == 0 || civ.y_iter == 0 || civ.z_iter == 0)
384 {
385 // In this case the volume has been shifted beyond the location of the pixel rdata and
386 // the client should expect to see a volume with nothing in it.
387
388 // Set size calls realloc,
389 set_size(new_nx, new_ny, new_nz);
390
391 // Set pixel memory to zero - the client should expect to see nothing
392 EMUtil::em_memset(rdata, 0, (size_t)new_nx*new_ny*new_nz);
393
394 return;
395 }
396
397 // Resize the volume before memory shifting occurs if the new volume is larger than the previous volume
398 // All of the pixel rdata is guaranteed to be at the start of the new volume because realloc (called in set size)
399 // guarantees this.
400 if ( new_size > prev_size )
401 set_size(new_nx, new_ny, new_nz);
402
403 // Store the clipped row size.
404 size_t clipped_row_size = (civ.x_iter) * sizeof(float);
405
406 // Get the new sector sizes to save multiplication later.
407 size_t new_sec_size = new_nx * new_ny;
408 size_t prev_sec_size = prev_nx * prev_ny;
409
410 // Determine the memory locations of the source and destination pixels - at the point nearest
411 // to the beginning of the volume (rdata)
412 size_t src_it_begin = civ.prv_z_bottom*prev_sec_size + civ.prv_y_front*prev_nx + civ.prv_x_left;
413 size_t dst_it_begin = civ.new_z_bottom*new_sec_size + civ.new_y_front*new_nx + civ.new_x_left;
414
415 // This loop is in the forward direction (starting at points nearest to the beginning of the volume)
416 // it copies memory only when the destination pointer is less the source pointer - therefore
417 // ensuring that no memory "copied to" is yet to be "copied from"
418 for (int i = 0; i < civ.z_iter; ++i) {
419 for (int j = 0; j < civ.y_iter; ++j) {
420
421 // Determine the memory increments as dependent on i and j
422 // This could be optimized so that not so many multiplications are occurring...
423 size_t dst_inc = dst_it_begin + j*new_nx + i*new_sec_size;
424 size_t src_inc = src_it_begin + j*prev_nx + i*prev_sec_size;
425 float* local_dst = rdata + dst_inc;
426 float* local_src = rdata + src_inc;
427
428 if ( dst_inc >= src_inc )
429 {
430 // this is fine, it will happen now and then and it will be necessary to continue.
431 // the tempatation is to break, but you can't do that (because the point where memory intersects
432 // could be in this slice - and yes, this aspect could be optimized).
433 continue;
434 }
435
436 // Asserts are compiled only in debug mode
437 // This situation not encountered in testing thus far
438 Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
439
440 // Finally copy the memory
441 memmove(local_dst, local_src, clipped_row_size);
442 }
443 }
444
445 // Determine the memory locations of the source and destination pixels - at the point nearest
446 // to the end of the volume (rdata+new_size)
447 size_t src_it_end = prev_size - civ.prv_z_top*prev_sec_size - civ.prv_y_back*prev_nx - prev_nx + civ.prv_x_left;
448 size_t dst_it_end = new_size - civ.new_z_top*new_sec_size - civ.new_y_back*new_nx - new_nx + civ.new_x_left;
449
450 // This loop is in the reverse direction (starting at points nearest to the end of the volume).
451 // It copies memory only when the destination pointer is greater than the source pointer therefore
452 // ensuring that no memory "copied to" is yet to be "copied from"
453 for (int i = 0; i < civ.z_iter; ++i) {
454 for (int j = 0; j < civ.y_iter; ++j) {
455
456 // Determine the memory increments as dependent on i and j
457 size_t dst_inc = dst_it_end - j*new_nx - i*new_sec_size;
458 size_t src_inc = src_it_end - j*prev_nx - i*prev_sec_size;
459 float* local_dst = rdata + dst_inc;
460 float* local_src = rdata + src_inc;
461
462 if (dst_inc <= (src_inc + civ.x_iter ))
463 {
464 // Overlap
465 if ( dst_inc > src_inc )
466 {
467 // Because the memcpy operation is the forward direction, and this "reverse
468 // direction" loop is proceeding in a backwards direction, it is possible
469 // that memory copied to is yet to be copied from (because memcpy goes forward).
470 // In this scenario pixel memory "copied to" is yet to be "copied from"
471 // i.e. there is overlap
472
473 // memmove handles overlapping cases.
474 // memmove could use a temporary buffer, or could go just go backwards
475 // the specification doesn't say how the function behaves...
476 // If memmove creates a temporary buffer is clip_inplace no longer inplace?
477 memmove(local_dst, local_src, clipped_row_size);
478 }
479 continue;
480 }
481
482 // This situation not encountered in testing thus far
483 Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
484
485 // Perform the memory copy
486 EMUtil::em_memcpy(local_dst, local_src, clipped_row_size);
487 }
488 }
489
490 // Resize the volume after memory shifting occurs if the new volume is smaller than the previous volume
491 // set_size calls realloc, guaranteeing that the pixel rdata is in the right location.
492 if ( new_size < prev_size )
493 set_size(new_nx, new_ny, new_nz);
494
495 // Now set all the edges to zero
496
497 // Set the extra bottom z slices to the fill_value
498 if ( z0 < 0 )
499 {
500 //EMUtil::em_memset(rdata, 0, (-z0)*new_sec_size*sizeof(float));
501 size_t inc = (-z0)*new_sec_size;
502 std::fill(rdata,rdata+inc,fill_value);
503 }
504
505 // Set the extra top z slices to the fill_value
506 if ( civ.new_z_top > 0 )
507 {
508 float* begin_pointer = rdata + (new_nz-civ.new_z_top)*new_sec_size;
509 //EMUtil::em_memset(begin_pointer, 0, (civ.new_z_top)*new_sec_size*sizeof(float));
510 float* end_pointer = begin_pointer+(civ.new_z_top)*new_sec_size;
511 std::fill(begin_pointer,end_pointer,fill_value);
512 }
513
514 // Next deal with x and y edges by iterating through each slice
515 for ( int i = civ.new_z_bottom; i < civ.new_z_bottom + civ.z_iter; ++i )
516 {
517 // Set the extra front y components to the fill_value
518 if ( y0 < 0 )
519 {
520 float* begin_pointer = rdata + i*new_sec_size;
521 //EMUtil::em_memset(begin_pointer, 0, (-y0)*new_nx*sizeof(float));
522 float* end_pointer = begin_pointer+(-y0)*new_nx;
523 std::fill(begin_pointer,end_pointer,fill_value);
524 }
525
526 // Set the extra back y components to the fill_value
527 if ( civ.new_y_back > 0 )
528 {
529 float* begin_pointer = rdata + i*new_sec_size + (new_ny-civ.new_y_back)*new_nx;
530 //EMUtil::em_memset(begin_pointer, 0, (civ.new_y_back)*new_nx*sizeof(float));
531 float* end_pointer = begin_pointer+(civ.new_y_back)*new_nx;
532 std::fill(begin_pointer,end_pointer,fill_value);
533 }
534
535 // Iterate through the y to set each correct x component to the fill_value
536 for (int j = civ.new_y_front; j <civ.new_y_front + civ.y_iter; ++j)
537 {
538 // Set the extra left x components to the fill_value
539 if ( x0 < 0 )
540 {
541 float* begin_pointer = rdata + i*new_sec_size + j*new_nx;
542 //EMUtil::em_memset(begin_pointer, 0, (-x0)*sizeof(float));
543 float* end_pointer = begin_pointer+(-x0);
544 std::fill(begin_pointer,end_pointer,fill_value);
545 }
546
547 // Set the extra right x components to the fill_value
548 if ( civ.new_x_right > 0 )
549 {
550 float* begin_pointer = rdata + i*new_sec_size + j*new_nx + (new_nx - civ.new_x_right);
551 //EMUtil::em_memset(begin_pointer, 0, (civ.new_x_right)*sizeof(float));
552 float* end_pointer = begin_pointer+(civ.new_x_right);
553 std::fill(begin_pointer,end_pointer,fill_value);
554 }
555
556 }
557 }
558
559// These couts may be useful
560// cout << "start starts " << civ.prv_x_left << " " << civ.prv_y_front << " " << civ.prv_z_bottom << endl;
561// cout << "start ends " << civ.prv_x_right << " " << civ.prv_y_back << " " << civ.prv_z_top << endl;
562// cout << "dst starts " << civ.new_x_left << " " << civ.new_y_front << " " << civ.new_z_bottom << endl;
563// cout << "dst ends " << civ.new_x_right << " " << civ.new_y_back << " " << civ.new_z_top << endl;
564// cout << "total iter z - " << civ.z_iter << " y - " << civ.y_iter << " x - " << civ.x_iter << endl;
565// cout << "=====" << endl;
566// cout << "dst_end is " << dst_it_end << " src end is " << src_it_end << endl;
567// cout << "dst_begin is " << dst_it_begin << " src begin is " << src_it_begin << endl;
568
569 // Update appropriate attributes (Copied and pasted from get_clip)
570 if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
571 attr_dict.has_key("origin_z") )
572 {
573 float xorigin = attr_dict["origin_x"];
574 float yorigin = attr_dict["origin_y"];
575 float zorigin = attr_dict["origin_z"];
576
577 float apix_x = attr_dict["apix_x"];
578 float apix_y = attr_dict["apix_y"];
579 float apix_z = attr_dict["apix_z"];
580
581 set_xyz_origin(xorigin + apix_x * area.origin[0],
582 yorigin + apix_y * area.origin[1],
583 zorigin + apix_z * area.origin[2]);
584 }
585
586 // Set the update flag because the size of the image has changed and stats should probably be recalculated if requested.
587 update();
588
589 EXITFUNC;
590}
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMData::set_xyz_origin
void set_xyz_origin(float origin_x, float origin_y, float origin_z)
Set the x,y,z origin of the image.
Definition: emdata.cpp:3286
EMAN::EMUtil::em_memset
static void em_memset(void *data, const int value, const size_t size)
Definition: emutil.h:377
EMAN::EMUtil::em_memcpy
static void em_memcpy(void *dst, const void *const src, const size_t size)
Definition: emutil.h:384
EMAN::Region::size
FloatSize size
Definition: geometry.h:655
EMAN::Region::origin
FloatPoint origin
Definition: geometry.h:654
Assert
#define Assert(s)
Define Assert() function that is effective only when -DDEBUG is used.
Definition: emassert.h:42
set_size
void set_size(int nx, int ny=1, int nz=1, bool noalloc=false)
Resize this EMData's main board memory pointer.

References Assert, attr_dict, EMAN::EMUtil::em_memcpy(), EMAN::EMUtil::em_memset(), ENTERFUNC, EXITFUNC, get_data(), EMAN::Dict::has_key(), ImageDimensionException, nx, ny, nz, EMAN::Region::origin, rdata, set_size(), set_xyz_origin(), EMAN::Region::size, and update().

Referenced by calc_ccf(), calc_fast_sigma_image(), calc_flcf(), EMAN::FFTResampleProcessor::fft_resample(), EMAN::ScaleTransformProcessor::process(), EMAN::XGradientProcessor::process_inplace(), EMAN::YGradientProcessor::process_inplace(), EMAN::ZGradientProcessor::process_inplace(), EMAN::IntTranslateProcessor::process_inplace(), EMAN::ScaleTransformProcessor::process_inplace(), and EMAN::FlattenBackgroundProcessor::process_inplace().

◆ compute_missingwedge()

EMData * EMData::compute_missingwedge ( float  wedgeangle,
float  start = 0.05,
float  stop = 0.5 
)

Find the mean and variance of voxels in the missing wedge.

Parameters
wedgeanglethe angle of the missing wedge
influncethe region of influnce in fourier space. This is a fudge factor between 0 and 0.5
wedgedirectionthe direction of the wedge, so far only a wedge along Z is supported (set wedgedirection to 0)

Definition at line 4369 of file emdata.cpp.

4370{
4371 EMData* test = new EMData();
4372 test->set_size(nx,ny,nz);
4373
4374 float ratio = tan((90.0f-wedgeangle)*M_PI/180.0f);
4375
4376 int offset_i = 2*int(start*nz/2);
4377 int offset_f = int(stop*nz/2);
4378
4379 int step = 0;
4380 float sum = 0.0;
4381 double square_sum = 0.0;
4382 for (int j = 0; j < offset_f; j++){
4383 for (int k = offset_i; k < offset_f; k++) {
4384 for (int i = 0; i < nx; i+=2) {
4385 if (i < int(k*ratio)) {
4386 test->set_value_at(i, j, k, 1.0);
4387 float v = std::sqrt(pow(get_value_at_wrap(i, j, k),2) + pow(get_value_at_wrap(i+1, j, k),2));
4388 sum += v;
4389 square_sum += v * (double)(v);
4390 step++;
4391 }
4392 }
4393 }
4394 }
4395
4396 float mean = sum / step;
4397 double var = (square_sum - sum*mean) / (step-1);
4398 float sigma = var >= 0.0 ? (float) std::sqrt(var) : 0.0;
4399
4400 cout << "Mean sqr wedge amp " << mean << " Sigma Squ wedge Amp " << sigma << endl;
4401 set_attr("spt_wedge_mean", mean);
4402 set_attr("spt_wedge_sigma", sigma);
4403
4404 return test;
4405}
get_value_at_wrap
float get_value_at_wrap(int x, int y, int z) const
Get the pixel density value at coordinates (x,y,z).
Definition: emdata_core.cpp:839
set_attr
void set_attr(const string &key, EMObject val)
Set a header attribute's value.

References EMData(), get_value_at_wrap(), nx, ny, nz, set_attr(), and sqrt().

◆ convolute()

EMData * EMData::convolute ( EMData *  with)

Convolutes 2 data sets.

The 2 images must be of the same size.

Parameters
withOne data set. 'this' image is the other data set.
Exceptions
NullPointerExceptionIf FFT resturns NULL image.
Returns
The result image.

Definition at line 3567 of file emdata.cpp.

3568{
3569 ENTERFUNC;
3570
3571 EMData *f1 = do_fft();
3572 if (!f1) {
3573 LOGERR("FFT returns NULL image");
3574 throw NullPointerException("FFT returns NULL image");
3575 }
3576
3577 f1->ap2ri();
3578
3579 EMData *cf = 0;
3580 if (with) {
3581 if (with-> is_real())
3582 cf = with->do_fft();
3583 else
3584 cf = with->copy();
3585 if (!cf) {
3586 LOGERR("FFT returns NULL image");
3587 throw NullPointerException("FFT returns NULL image");
3588 }
3589 cf->ap2ri();
3590 }
3591 else {
3592 cf = f1->copy();
3593 }
3594 //printf("cf_x=%d, f1y=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f1->get_ysize(),this->get_xsize(),with->get_xsize());
3595 if (with && !EMUtil::is_same_size(f1, cf)) {
3596 LOGERR("images not same size");
3597 throw ImageFormatException("images not same size");
3598 }
3599
3600 float *rdata1 = f1->get_data();
3601 float *rdata2 = cf->get_data();
3602 size_t cf_size = (size_t)cf->get_xsize() * cf->get_ysize() * cf->get_zsize();
3603
3604 float re,im;
3605
3606 for (size_t i = 0; i < cf_size; i += 2) {
3607 re = rdata1[i] * rdata2[i] - rdata1[i + 1] * rdata2[i + 1];
3608 im = rdata1[i + 1] * rdata2[i] + rdata1[i] * rdata2[i + 1];
3609 rdata2[i]=re;
3610 rdata2[i+1]=im;
3611 }
3612 cf->update();
3613 EMData *f2 = cf->do_ift();//ming change cf to cf_temp
3614 //printf("cf_x=%d, f2x=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f2->get_xsize(),this->get_xsize(),with->get_xsize());
3615 if( cf )
3616 {
3617 delete cf;
3618 cf = 0;
3619 }
3620
3621 if( f1 )
3622 {
3623 delete f1;
3624 f1=0;
3625 }
3626
3627 EXITFUNC;
3628 return f2;
3629}
is_real
bool is_real() const
Is this a real image?
Definition: emdata_metadata.h:667

References do_fft(), ENTERFUNC, EXITFUNC, ImageFormatException, is_real(), EMAN::EMUtil::is_same_size(), LOGERR, and NullPointerException.

Referenced by calc_fast_sigma_image(), and EMAN::FlattenBackgroundProcessor::process_inplace().

◆ cut_slice()

void EMData::cut_slice ( const EMData *const  map,
const Transform &  tr,
bool  interpolate = true 
)

cut a 2D slice out of a real 3D map.

Put slice into 'this' image.

Parameters
mapThe real 3D map.
trorientation of the slice as encapsulated in a Transform object.
interpolateDo interpolation or not.
Exceptions
NullPointerExceptionIf map is NULL.
ImageDimensionExceptionIf this image is not 2D.
ImageDimensionExceptionIf map image is not 3D.
ImageFormatExceptionIf this image is complex
ImageFormatExceptionIf map is complex
Author
David Woolford (adapted from an original version by Steve Ludtke)
Date
Feb 2008

Definition at line 4018 of file emdata.cpp.

4019{
4020 ENTERFUNC;
4021
4022 if (!map) throw NullPointerException("NULL image");
4023 // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
4024 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
4025 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
4026 // Now check for complex images - this is really just being thorough
4027 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
4028 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
4029
4030
4031 float *sdata = map->get_data();
4032 float *ddata = get_data();
4033
4034 int map_nx = map->get_xsize();
4035 int map_ny = map->get_ysize();
4036 int map_nz = map->get_zsize();
4037 int map_nxy = map_nx * map_ny;
4038
4039 int ymax = ny/2;
4040 if ( ny % 2 == 1 ) ymax += 1;
4041 int xmax = nx/2;
4042 if ( nx % 2 == 1 ) xmax += 1;
4043 for (int y = -ny/2; y < ymax; y++) {
4044 for (int x = -nx/2; x < xmax; x++) {
4045 Vec3f coord(x,y,0);
4046 Vec3f soln = transform*coord;
4047
4048// float xx = (x+pretrans[0]) * (*ort)[0][0] + (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
4049// float yy = (x+pretrans[0]) * (*ort)[1][0] + (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
4050// float zz = (x+pretrans[0]) * (*ort)[2][0] + (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
4051
4052
4053// xx += map_nx/2;
4054// yy += map_ny/2;
4055// zz += map_nz/2;
4056
4057 float xx = soln[0]+map_nx/2;
4058 float yy = soln[1]+map_ny/2;
4059 float zz = soln[2]+map_nz/2;
4060
4061 int l = (x+nx/2) + (y+ny/2) * nx;
4062
4063 float t = xx - floor(xx);
4064 float u = yy - floor(yy);
4065 float v = zz - floor(zz);
4066
4067 if (xx < 0 || yy < 0 || zz < 0 ) {
4068 ddata[l] = 0;
4069 continue;
4070 }
4071 if (interpolate) {
4072 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
4073 ddata[l] = 0;
4074 continue;
4075 }
4076 int k = (int) (Util::fast_floor(xx) + Util::fast_floor(yy) * map_nx + Util::fast_floor(zz) * map_nxy);
4077
4078
4079 if (xx < (map_nx - 1) && yy < (map_ny - 1) && zz < (map_nz - 1)) {
4080 ddata[l] = Util::trilinear_interpolate(sdata[k],
4081 sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
4082 sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
4083 sdata[k + map_nx + map_nxy + 1],t, u, v);
4084 }
4085 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) && zz == (map_nz - 1) ) {
4086 ddata[l] += sdata[k];
4087 }
4088 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) ) {
4089 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nxy],v);
4090 }
4091 else if ( xx == (map_nx - 1) && zz == (map_nz - 1) ) {
4092 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nx],u);
4093 }
4094 else if ( yy == (map_ny - 1) && zz == (map_nz - 1) ) {
4095 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + 1],t);
4096 }
4097 else if ( xx == (map_nx - 1) ) {
4098 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + map_nx], sdata[k + map_nxy], sdata[k + map_nxy + map_nx],u,v);
4099 }
4100 else if ( yy == (map_ny - 1) ) {
4101 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nxy], sdata[k + map_nxy + 1],t,v);
4102 }
4103 else if ( zz == (map_nz - 1) ) {
4104 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nx], sdata[k + map_nx + 1],t,u);
4105 }
4106
4107// if (k >= map->get_size()) {
4108// cout << xx << " " << yy << " " << zz << " " << endl;
4109// cout << k << " " << get_size() << endl;
4110// cout << get_xsize() << " " << get_ysize() << " " << get_zsize() << endl;
4111// throw;
4112// }
4113//
4114// ddata[l] = Util::trilinear_interpolate(sdata[k],
4115// sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
4116// sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
4117// sdata[k + map_nx + map_nxy + 1],t, u, v);
4118 }
4119 else {
4120 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
4121 ddata[l] = 0;
4122 continue;
4123 }
4124 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
4125 ddata[l] = sdata[k];
4126 }
4127
4128 }
4129 }
4130
4131 update();
4132
4133 EXITFUNC;
4134}
EMAN::EMData::transform
void transform(const Transform &t)
Rotate then translate the image.
Definition: emdata.h:295
EMAN::Util::linear_interpolate
static float linear_interpolate(float p1, float p2, float t)
Calculate linear interpolation.
Definition: util.h:518
EMAN::Util::fast_floor
static int fast_floor(float x)
A fast way to calculate a floor, which is largest integral value not greater than argument.
Definition: util.h:874
EMAN::Util::trilinear_interpolate
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.
Definition: util.h:619
EMAN::Util::bilinear_interpolate
static float bilinear_interpolate(float p1, float p2, float p3, float p4, float t, float u)
Calculate bilinear interpolation.
Definition: util.h:543

References EMAN::Util::bilinear_interpolate(), ENTERFUNC, EXITFUNC, EMAN::Util::fast_floor(), get_data(), get_ndim(), ImageDimensionException, ImageFormatException, is_complex(), EMAN::Util::linear_interpolate(), NullPointerException, nx, ny, EMAN::Util::round(), transform(), EMAN::Util::trilinear_interpolate(), update(), x, and y.

◆ debug_print_parms()

void EMAN::EMData::debug_print_parms ( )
inline

Printing EMData params for debugging purpose.

Definition at line 779 of file emdata.h.

780 {
781 std::cout << "Printing EMData params" << std::endl;
782 for ( Dict::const_iterator it = attr_dict.begin(); it != attr_dict.end(); ++it )
783 {
784 std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
785 }
786 std::cout << "Done printing EMData params" << std::endl;
787 }
EMAN::Dict::end
iterator end()
Definition: emobject.cpp:1061
EMAN::Dict::begin
iterator begin()
Definition: emobject.cpp:1045

References attr_dict, EMAN::Dict::begin(), and EMAN::Dict::end().

◆ do_radon()

EMData * EMData::do_radon ( )

Radon Transform: an algorithm that transforms an original image into a series of equiangular projections.

When applied to a 2D object, the output of the Radon transform is a series of 1D lines.

Do radon transformation on this image. This image must be 2D square.

Exceptions
ImageFormatExceptionIf the image is not square.
ImageDimensionExceptionIf the image is not 2D.
Returns
Radon transform image in square.

Definition at line 1413 of file emdata.cpp.

1414{
1415 ENTERFUNC;
1416
1417 if (get_ndim() != 2) {
1418 throw ImageDimensionException("2D only");
1419 }
1420
1421 if (nx != ny) {
1422 throw ImageFormatException("square image only");
1423 }
1424
1425 EMData *result = new EMData();
1426 result->set_size(nx, ny, 1);
1427 result->to_zero();
1428 float *result_data = result->get_data();
1429
1430 EMData *this_copy = this;
1431 this_copy = copy();
1432
1433 for (int i = 0; i < nx; i++) {
1434 Transform t(Dict("type","2d","alpha",(float) M_PI * 2.0f * i / nx));
1435 this_copy->transform(t);
1436
1437 float *copy_data = this_copy->get_data();
1438
1439 for (int y = 0; y < nx; y++) {
1440 for (int x = 0; x < nx; x++) {
1441 if (Util::square(x - nx / 2) + Util::square(y - nx / 2) <= nx * nx / 4) {
1442 result_data[i + y * nx] += copy_data[x + y * nx];
1443 }
1444 }
1445 }
1446
1447 this_copy->update();
1448 }
1449
1450 result->update();
1451
1452 if( this_copy )
1453 {
1454 delete this_copy;
1455 this_copy = 0;
1456 }
1457
1458 EXITFUNC;
1459 return result;
1460}
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75

References copy(), EMData(), ENTERFUNC, EXITFUNC, get_ndim(), ImageDimensionException, ImageFormatException, nx, ny, EMAN::Util::square(), transform(), x, and y.

◆ dot_rotate_translate()

double EMData::dot_rotate_translate ( EMData *  with,
float  dx,
float  dy,
float  da,
const bool  mirror = false 
)

dot product of 2 images.

Then 'this' image is rotated/translated. It is much faster than Rotate/Translate then dot product. 2D images only.

Parameters
withThe image used to do the dot product.
dxTranslation distance in x direction.
dyTranslation distance in y direction.
daRotation euler angle in degrees
mirror
Exceptions
ImageFormatExceptionIf the 2 images are not the same size.
ImageDimensionExceptionIf the image is 3D.
Returns

Definition at line 1239 of file emdata.cpp.

1240{
1241 ENTERFUNC;
1242
1243 if (!EMUtil::is_same_size(this, with)) {
1244 LOGERR("images not same size");
1245 throw ImageFormatException("images not same size");
1246 }
1247
1248 if (get_ndim() == 3) {
1249 LOGERR("1D/2D Images only");
1250 throw ImageDimensionException("1D/2D only");
1251 }
1252
1253 float *this_data = 0;
1254
1255 this_data = get_data();
1256
1257 float da_rad = da*(float)M_PI/180.0f;
1258
1259 float *with_data = with->get_data();
1260 float mx0 = cos(da_rad);
1261 float mx1 = sin(da_rad);
1262 float y = -ny / 2.0f;
1263 float my0 = mx0 * (-nx / 2.0f - 1.0f) + nx / 2.0f - dx;
1264 float my1 = -mx1 * (-nx / 2.0f - 1.0f) + ny / 2.0f - dy;
1265 double result = 0;
1266
1267 for (int j = 0; j < ny; j++) {
1268 float x2 = my0 + mx1 * y;
1269 float y2 = my1 + mx0 * y;
1270
1271 int ii = Util::fast_floor(x2);
1272 int jj = Util::fast_floor(y2);
1273 float t = x2 - ii;
1274 float u = y2 - jj;
1275
1276 for (int i = 0; i < nx; i++) {
1277 t += mx0;
1278 u -= mx1;
1279
1280 if (t >= 1.0f) {
1281 ii++;
1282 t -= 1.0f;
1283 }
1284
1285 if (u >= 1.0f) {
1286 jj++;
1287 u -= 1.0f;
1288 }
1289
1290 if (t < 0) {
1291 ii--;
1292 t += 1.0f;
1293 }
1294
1295 if (u < 0) {
1296 jj--;
1297 u += 1.0f;
1298 }
1299
1300 if (ii >= 0 && ii <= nx - 2 && jj >= 0 && jj <= ny - 2) {
1301 int k0 = ii + jj * nx;
1302 int k1 = k0 + 1;
1303 int k2 = k0 + nx + 1;
1304 int k3 = k0 + nx;
1305
1306 float tt = 1 - t;
1307 float uu = 1 - u;
1308 int idx = i + j * nx;
1309 if (mirror) idx = nx-1-i+j*nx; // mirroring of Transforms is always about the y axis
1310 result += (this_data[k0] * tt * uu + this_data[k1] * t * uu +
1311 this_data[k2] * t * u + this_data[k3] * tt * u) * with_data[idx];
1312 }
1313 }
1314 y += 1.0f;
1315 }
1316
1317 EXITFUNC;
1318 return result;
1319}

References ENTERFUNC, EXITFUNC, EMAN::Util::fast_floor(), get_data(), get_ndim(), ImageDimensionException, ImageFormatException, EMAN::EMUtil::is_same_size(), LOGERR, nx, ny, and y.

Referenced by refalifnfast().

◆ extract_box()

EMData * EMData::extract_box ( const Transform &  cs,
const Region &  r 
)

Extract a box from EMData in an abritrary orrientation.

Used for extracting helix boxes from tomograms

EMData extract_box(const Transform& cs, const int ix, const int fx, const int iy, const int yf, const int zi, const int zf);

Parameters
csTransform describing the coordinate system of the box realative to the standard coord system
rRegion describe the volume to extract, in the local coordinate system
Author
John Flanagan
Date
Aug 2011

Definition at line 4301 of file emdata.cpp.

4302{
4303 vector<float> cs_matrix = cs.get_matrix();
4304
4305 EMData* box = new EMData();
4306 box->set_size((r.get_width()-r.x_origin()), (r.get_height()-r.y_origin()), (r.get_depth()-r.z_origin()));
4307 int box_nx = box->get_xsize();
4308 int box_ny = box->get_ysize();
4309 int box_nxy = box_nx*box_ny;
4310 float* bdata = box->get_data();
4311 float* ddata = get_data();
4312
4313 for (int x = r.x_origin(); x < r.get_width(); x++) {
4314 for (int y = r.y_origin(); y < r.get_height(); y++) {
4315 for (int z = r.z_origin(); z < r.get_depth(); z++) {
4316 //float xb = cs_matrix[0]*x + cs_matrix[1]*y + cs_matrix[2]*z + cs_matrix[3];
4317 //float yb = cs_matrix[4]*x + cs_matrix[5]*y + cs_matrix[6]*z + cs_matrix[7];
4318 //float zb = cs_matrix[8]*x + cs_matrix[9]*y + cs_matrix[10]*z + cs_matrix[11];
4319 float xb = cs_matrix[0]*x + y*cs_matrix[4] + z*cs_matrix[8] + cs_matrix[3];
4320 float yb = cs_matrix[1]*x + y*cs_matrix[5] + z*cs_matrix[9] + cs_matrix[7];
4321 float zb = cs_matrix[2]*x + y*cs_matrix[6] + z*cs_matrix[10] + cs_matrix[11];
4322 float t = xb - Util::fast_floor(xb);
4323 float u = yb - Util::fast_floor(yb);
4324 float v = zb - Util::fast_floor(zb);
4325
4326 //cout << x << " " << y << " " << z << " Box " << xb << " " << yb << " " << zb << endl;
4327 int l = (x - r.x_origin()) + (y - r.y_origin())*box_nx + (z - r.z_origin())*box_nxy;
4328 int k = (int) (Util::fast_floor(xb) + Util::fast_floor(yb) * nx + Util::fast_floor(zb) * nxy);
4329 //cout << k << " " << l << endl;
4330 if ( xb > nx - 1 || yb > ny - 1 || zb > nz - 1) {
4331 bdata[l] = 0;
4332 continue;
4333 }
4334 if (xb < 0 || yb < 0 || zb < 0){
4335 bdata[l] = 0;
4336 continue;
4337 }
4338
4339 if (xb < (nx - 1) && yb < (ny - 1) && zb < (nz - 1)) {
4340 bdata[l] = Util::trilinear_interpolate(ddata[k], ddata[k + 1], ddata[k + nx],ddata[k + nx + 1], ddata[k + nxy], ddata[k + nxy + 1], ddata[k + nx + nxy], ddata[k + nx + nxy + 1],t, u, v);
4341 }
4342 }
4343 }
4344 }
4345
4346 return box;
4347}
EMAN::Region::get_width
float get_width() const
get the width
Definition: geometry.h:606
EMAN::Region::y_origin
float y_origin() const
get the y element of the origin
Definition: geometry.h:622
EMAN::Region::x_origin
float x_origin() const
get the x element of the origin
Definition: geometry.h:620
EMAN::Region::get_depth
float get_depth() const
get the depth
Definition: geometry.h:610
EMAN::Region::z_origin
float z_origin() const
get the z element of the origin
Definition: geometry.h:624
EMAN::Region::get_height
float get_height() const
get the height
Definition: geometry.h:608
EMAN::Transform::get_matrix
vector< float > get_matrix() const
Get the transformation matrix using a vector.
Definition: transform.cpp:181

References EMData(), EMAN::Util::fast_floor(), get_data(), EMAN::Region::get_depth(), EMAN::Region::get_height(), EMAN::Transform::get_matrix(), EMAN::Region::get_width(), nx, nxy, ny, nz, EMAN::Util::trilinear_interpolate(), x, EMAN::Region::x_origin(), y, EMAN::Region::y_origin(), and EMAN::Region::z_origin().

◆ get_clip()

EMData * EMData::get_clip ( const Region &  area,
const float  fill = 0 
) const

Get an inclusive clip.

Pads to fill if larger than this image.

Parameters
areaThe clip area, can be 2D/3D
fillthe value to assign new pixels outside the area of the original image
Exceptions
ImageDimensionExceptionif any of the dimensions of the argument region are negative
Returns
The clip image.

Definition at line 592 of file emdata.cpp.

593{
594 ENTERFUNC;
595 // Steve: removed check 10/1/17. Prevented some sane operations, like a 3d clip which was actually 2d from a 2d volume
596// if (get_ndim() != area.get_ndim()) {
597// LOGERR("cannot get %dD clip out of %dD image", area.get_ndim(),get_ndim());
598// return 0;
599// }
600
601 EMData *result = new EMData();
602
603 // Ensure that all of the metadata of this is stored in the new object
604 // Originally added to ensure that euler angles were retained when preprocessing (zero padding) images
605 // prior to insertion into the 3D for volume in the reconstruction phase (see reconstructor.cpp/h).
606 result->attr_dict = this->attr_dict;
607 int zsize = (int)area.size[2];
608 if (zsize == 0 && nz <= 1) {
609 zsize = 1;
610 }
611 int ysize = (ny<=1 && (int)area.size[1]==0 ? 1 : (int)area.size[1]);
612
613 if ( (int)area.size[0] < 0 || ysize < 0 || zsize < 0 )
614 {
615 // Negative image dimensions not supported - added retrospectively by d.woolford (who didn't write get_clip but wrote clip_inplace)
616 throw ImageDimensionException("New image dimensions are negative - this is not supported in the the get_clip operation");
617 }
618
619//#ifdef EMAN2_USING_CUDA
620 // Strategy is always to prefer using the GPU if possible
621// bool use_gpu = false;
622// if ( gpu_operation_preferred() ) {
623// result->set_size_cuda((int)area.size[0], ysize, zsize);
624 //CudaDataLock lock(this); // Just so we never have to recopy this data to and from the GPU
625// result->get_cuda_data(); // Force the allocation - set_size_cuda is lazy
626 // Setting the value is necessary seeing as cuda data is not automatically zeroed
627// result->to_value(fill); // This will automatically use the GPU.
628// use_gpu = true;
629// } else { // cpu == True
630// result->set_size((int)area.size[0], ysize, zsize);
631// if (fill != 0.0) { result->to_value(fill); };
632// }
633//#else
634 result->set_size((int)area.size[0], ysize, zsize);
635 if (fill != 0.0) { result->to_value(fill); };
636//#endif //EMAN2_USING_CUDA
637
638 int x0 = (int) area.origin[0];
639 x0 = x0 < 0 ? 0 : x0;
640
641 int y0 = (int) area.origin[1];
642 y0 = y0 < 0 ? 0 : y0;
643
644 int z0 = (int) area.origin[2];
645 z0 = z0 < 0 ? 0 : z0;
646
647 int x1 = (int) (area.origin[0] + area.size[0]);
648 x1 = x1 > nx ? nx : x1;
649
650 int y1 = (int) (area.origin[1] + area.size[1]);
651 y1 = y1 > ny ? ny : y1;
652
653 int z1 = (int) (area.origin[2] + area.size[2]);
654 z1 = z1 > nz ? nz : z1;
655 if (z1 <= 0) {
656 z1 = 1;
657 }
658
659 result->insert_clip(this,-((IntPoint)area.origin));
660
661 if( attr_dict.has_key("apix_x") && attr_dict.has_key("apix_y") &&
662 attr_dict.has_key("apix_z") )
663 {
664 float apix_x = attr_dict["apix_x"];
665 float apix_y = attr_dict["apix_y"];
666 float apix_z = attr_dict["apix_z"];
667 result->set_attr("apix_x",apix_x);
668 result->set_attr("apix_y",apix_y);
669 result->set_attr("apix_z",apix_z);
670
671 if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
672 attr_dict.has_key("origin_z") )
673 {
674 float xorigin = attr_dict["origin_x"];
675 float yorigin = attr_dict["origin_y"];
676 float zorigin = attr_dict["origin_z"];
677
678
679 result->set_xyz_origin(xorigin + apix_x * area.origin[0],
680 yorigin + apix_y * area.origin[1],
681 zorigin + apix_z * area.origin[2]);
682 }
683 else {
684 result->set_xyz_origin(apix_x * area.origin[0],
685 apix_y * area.origin[1],
686 apix_z * area.origin[2]);
687 }
688
689 }
690
691//#ifdef EMAN2_USING_CUDA
692// if (use_gpu) result->gpu_update();
693// else result->update();
694//#else
695 result->update();
696//#endif // EMAN2_USING_CUDA
697
698
699 result->set_path(path);
700 result->set_pathnum(pathnum);
701
702 EXITFUNC;
703 return result;
704}
EMAN::IntPoint
IntPoint defines an integer-coordinate point in a 1D/2D/3D space.
Definition: geometry.h:192

References attr_dict, EMData(), ENTERFUNC, EXITFUNC, EMAN::Dict::has_key(), ImageDimensionException, nx, ny, nz, EMAN::Region::origin, path, pathnum, set_xyz_origin(), and EMAN::Region::size.

Referenced by calc_fast_sigma_image(), make_footprint(), EMAN::HarmonicProcessor::process(), EMAN::CtfSimProcessor::process(), EMAN::IntTranslateProcessor::process(), EMAN::ScaleTransformProcessor::process(), EMAN::CutoffBlockProcessor::process_inplace(), EMAN::SNRProcessor::process_inplace(), window_center(), EMAN::RT2DTreeAligner::xform_align_nbest(), EMAN::RT2Dto3DTreeAligner::xform_align_nbest(), EMAN::RT3DTreeAligner::xform_align_nbest(), and EMAN::RT3DLocalTreeAligner::xform_align_nbest().

◆ get_rotated_clip()

EMData * EMData::get_rotated_clip ( const Transform &  xform,
const IntSize &  size,
float  scale = 1.0 
)

This will extract an arbitrarily oriented and sized region from the image.

Parameters
xformThe transformation of the region.
sizeSize of the clip.
scaleScaling put on the returned image.
Returns
The clip image.

Definition at line 733 of file emdata.cpp.

735{
736 EMData *result = new EMData();
737 result->set_size(size[0],size[1],size[2]);
738
739 if (nz==1) {
740 for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
741 for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
742 Vec3f xv=xform.transform(Vec3f((float)x,(float)y,0.0f));
743 float v = 0;
744
745 if (xv[0]<0||xv[1]<0||xv[0]>nx-2||xv[1]>ny-2) v=0.;
746 else v=sget_value_at_interp(xv[0],xv[1]);
747 result->set_value_at(x+size[0]/2,y+size[1]/2,v);
748 }
749 }
750 }
751 else {
752 for (int z=-size[2]/2; z<(size[2]+1)/2; z++) {
753 for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
754 for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
755 Vec3f xv=xform.transform(Vec3f((float)x,(float)y,(float)z));
756 float v = 0;
757
758 if (xv[0]<0||xv[1]<0||xv[2]<0||xv[0]>nx-2||xv[1]>ny-2||xv[2]>nz-2) v=0.;
759 else v=sget_value_at_interp(xv[0],xv[1],xv[2]);
760 result->set_value_at(x+size[0]/2,y+size[1]/2,z+size[2]/2,v);
761 }
762 }
763 }
764 }
765
766 result->attr_dict["apix_x"] = attr_dict["apix_x"];
767 result->attr_dict["apix_y"] = attr_dict["apix_y"];
768 result->attr_dict["apix_z"] = attr_dict["apix_z"];
769
770 result->update();
771
772 return result;
773}
EMAN::Transform::transform
Vec2f transform(const float &x, const float &y) const
Transform 2D coordinates using the internal transformation matrix.
Definition: transform.h:417
sget_value_at_interp
float sget_value_at_interp(float x, float y) const
Get pixel density value at interpolation of (x,y).
Definition: emdata_core.cpp:914
EMAN::Vec3f
Vec3< float > Vec3f
Definition: vec3.h:693

References attr_dict, EMData(), nx, ny, nz, sget_value_at_interp(), EMAN::Transform::transform(), x, and y.

◆ get_top_half()

EMData * EMData::get_top_half ( ) const

Get the top half of this 3D image.

Exceptions
ImageDimensionExceptionIf this image is not 3D.
Returns
The top half of this image.

Definition at line 707 of file emdata.cpp.

708{
709 ENTERFUNC;
710
711 if (get_ndim() != 3) {
712 throw ImageDimensionException("3D only");
713 }
714
715 EMData *half = new EMData();
716 half->attr_dict = attr_dict;
717 half->set_size(nx, ny, nz / 2);
718
719 float *half_data = half->get_data();
720 EMUtil::em_memcpy(half_data, &(get_data()[(size_t)nz / 2 * (size_t)nx * (size_t)ny]), sizeof(float) * (size_t)nx * (size_t)ny * (size_t)nz / 2lu);
721
722 float apix_z = attr_dict["apix_z"];
723 float origin_z = attr_dict["origin_z"];
724 origin_z += apix_z * nz / 2;
725 half->attr_dict["origin_z"] = origin_z;
726 half->update();
727
728 EXITFUNC;
729 return half;
730}

References attr_dict, EMAN::EMUtil::em_memcpy(), EMData(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, nx, ny, and nz.

◆ getResolution()

int EMAN::EMData::getResolution ( ) const
inline

function for MarchingCubes, for 3D image display

Returns
the resolution

Definition at line 766 of file emdata.h.

766 {
767 int resolution = 0;
768 int num = 1;
769 while(num < get_xsize()) {
770 resolution++;
771 num = 1 << resolution;
772 }
773
774 return resolution;
775 }
get_xsize
int get_xsize() const
Get the image x-dimensional size.
Definition: emdata_metadata.h:544

References get_xsize().

◆ little_big_dot()

EMData * EMData::little_big_dot ( EMData *  little_img,
bool  do_sigma = false 
)

This does a normalized dot product of a little image with a big image using real-space methods.

The result is the same size as 'this', but a border 1/2 the size of 'little_img' will be zero. This routine is only efficient when 'little_img' is fairly small.

Parameters
little_imgA small image.
do_sigmaCalculate sigma or not.
Exceptions
ImageDimensionExceptionIf the image is not 1D/2D.
Returns
normalized dot product image.

Definition at line 1322 of file emdata.cpp.

1323{
1324 ENTERFUNC;
1325
1326 if (get_ndim() > 2) {
1327 throw ImageDimensionException("1D/2D only");
1328 }
1329
1330 EMData *ret = copy_head();
1331 ret->set_size(nx,ny,nz);
1332 ret->to_zero();
1333
1334 int nx2 = with->get_xsize();
1335 int ny2 = with->get_ysize();
1336 float em = with->get_edge_mean();
1337
1338 float *data = get_data();
1339 float *with_data = with->get_data();
1340 float *ret_data = ret->get_data();
1341
1342 float sum2 = (Util::square((float)with->get_attr("sigma")) +
1343 Util::square((float)with->get_attr("mean")));
1344 if (do_sigma) {
1345 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
1346 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
1347 float sum = 0;
1348 float sum1 = 0;
1349 float summ = 0;
1350 int k = 0;
1351
1352 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1353 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1354 int l = ii + jj * nx;
1355 sum1 += Util::square(data[l]);
1356 summ += data[l];
1357 sum += data[l] * with_data[k];
1358 k++;
1359 }
1360 }
1361 float tmp_f1 = (sum1 / 2.0f - sum) / (nx2 * ny2);
1362 float tmp_f2 = Util::square((float)with->get_attr("mean") -
1363 summ / (nx2 * ny2));
1364 ret_data[i + j * nx] = sum2 + tmp_f1 - tmp_f2;
1365 }
1366 }
1367 }
1368 else {
1369 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
1370 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
1371 float eml = 0;
1372 float dot = 0;
1373 float dot2 = 0;
1374
1375 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1376 eml += data[ii + (j - ny2 / 2) * nx] + data[ii + (j + ny2 / 2 - 1) * nx];
1377 }
1378
1379 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1380 eml += data[i - nx2 / 2 + jj * nx] + data[i + nx2 / 2 - 1 + jj * nx];
1381 }
1382
1383 eml /= (nx2 + ny2) * 2.0f;
1384 int k = 0;
1385
1386 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
1387 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
1388 dot += (data[ii + jj * nx] - eml) * (with_data[k] - em);
1389 dot2 += Util::square(data[ii + jj * nx] - eml);
1390 k++;
1391 }
1392 }
1393
1394 dot2 = std::sqrt(dot2);
1395
1396 if (dot2 == 0) {
1397 ret_data[i + j * nx] = 0;
1398 }
1399 else {
1400 ret_data[i + j * nx] = dot / (nx2 * ny2 * dot2 * (float)with->get_attr("sigma"));
1401 }
1402 }
1403 }
1404 }
1405
1406 ret->update();
1407
1408 EXITFUNC;
1409 return ret;
1410}
copy_head
EMData * copy_head() const
Make an image with a copy of the current image's header.
Definition: emdata_core.cpp:100
EMAN::dot
double dot(const Vector3 &w, const Vector3 &v)
Definition: vecmath.h:305

References copy_head(), EMAN::dot(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, nx, ny, nz, sqrt(), and EMAN::Util::square().

◆ make_footprint()

EMData * EMData::make_footprint ( int  type = 0)

Makes a 'footprint' for the current image.

This is image containing a rotational & translational invariant of the parent image. The size of the resulting image depends on the selected type.

type 0- The original, default footprint derived from the rotational footprint types 1-6 - bispectrum-based types 1,3,5 - returns Fouier-like images types 2,4,6 - returns real-space-like images type 1,2 - simple r1,r2, 2-D footprints type 3,4 - r1,r2,anle 3D footprints type 5,6 - same as 1,2 but with the cube root of the final products used

Parameters
typeSelect one of several possible algorithms for producing the invariants
Exceptions
ImageFormatExceptionIf image size is not even.
Returns
The footprint image.

Definition at line 2035 of file emdata.cpp.

2036{
2037// printf("Make fp %d\n",type);
2038 if (type==0) {
2039 EMData *un=make_rotational_footprint_e1(); // Use EMAN1's footprint strategy
2040 if (un->get_ysize() <= 6) {
2041 throw UnexpectedBehaviorException("In EMData::make_footprint. The rotational footprint is too small");
2042 }
2043 EMData *tmp=un->get_clip(Region(0,4,un->get_xsize(),un->get_ysize()-6)); // 4 and 6 are empirical
2044 EMData *cx=tmp->calc_ccfx(tmp,0,-1,1);
2045 EMData *fp=cx->get_clip(Region(0,0,cx->get_xsize()/2,cx->get_ysize()));
2046 delete un;
2047 delete tmp;
2048 delete cx;
2049 return fp;
2050 }
2051 else if (type==1 || type==2 ||type==5 || type==6) {
2052 int i,j,kx,ky,lx,ly;
2053
2054 EMData *fft=do_fft();
2055
2056 // map for x,y -> radius for speed
2057 int rmax=(get_xsize()+1)/2;
2058 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
2059 for (i=0; i<rmax; i++) {
2060 for (j=0; j<rmax; j++) {
2061 rmap[i+j*rmax]=hypot((float)i,(float)j);
2062// printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
2063 }
2064 }
2065
2066 EMData *fp=new EMData(rmax*2+2,rmax*2,1);
2067 fp->set_complex(1);
2068 fp->to_zero();
2069
2070 // Two vectors in to complex space (kx,ky) and (lx,ly)
2071 // We are computing the bispectrum, f(k).f(l).f*(k+l)
2072 // but integrating out two dimensions, leaving |k|,|l|
2073 for (kx=-rmax+1; kx<rmax; kx++) {
2074 for (ky=-rmax+1; ky<rmax; ky++) {
2075 for (lx=-rmax+1; lx<rmax; lx++) {
2076 for (ly=-rmax+1; ly<rmax; ly++) {
2077 int ax=kx+lx;
2078 int ay=ky+ly;
2079 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
2080 int r1=(int)floor(.5+rmap[abs(kx)+rmax*abs(ky)]);
2081 int r2=(int)floor(.5+rmap[abs(lx)+rmax*abs(ly)]);
2082// if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
2083// float r3=rmap[ax+rmax*ay];
2084 if (r1+r2>=rmax) continue;
2085
2086 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
2087 fp->set_value_at(r1*2,r2,p.real()+fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2088// fp->set_value_at(r1*2,rmax*2-r2-1, fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2089// fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2)); // We keep only the real component in anticipation of zero phase sum
2090 fp->set_value_at(r1*2+1,r2,fp->get_value_at(r1*2+1,r2)+1); // a normalization counter
2091 }
2092 }
2093 }
2094 }
2095
2096 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
2097 if (type==5 || type==6) {
2098 for (i=0; i<rmax*2; i+=2) {
2099 for (j=0; j<rmax; j++) {
2100 float norm=fp->get_value_at(i+1,j);
2101 fp->set_value_at(i,rmax*2-j-1,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
2102 fp->set_value_at(i,j,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
2103 fp->set_value_at(i+1,j,0.0);
2104 }
2105 }
2106 }
2107 else {
2108 for (i=0; i<rmax*2; i+=2) {
2109 for (j=0; j<rmax; j++) {
2110 float norm=fp->get_value_at(i+1,j);
2111 fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
2112 fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
2113 fp->set_value_at(i+1,j,0.0);
2114 }
2115 }
2116 }
2117
2118 free(rmap);
2119 if (type==2||type==6) {
2120 EMData *f2=fp->do_ift();
2121 if (f2->get_value_at(0,0)<0) f2->mult(-1.0f);
2122 f2->process_inplace("xform.phaseorigin.tocorner");
2123 delete fp;
2124 return f2;
2125 }
2126 return fp;
2127 }
2128 else if (type==3 || type==4) {
2129 int h,i,j,kx,ky,lx,ly;
2130
2131 EMData *fft=do_fft();
2132
2133 // map for x,y -> radius for speed
2134 int rmax=(get_xsize()+1)/2;
2135 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
2136 for (i=0; i<rmax; i++) {
2137 for (j=0; j<rmax; j++) {
2138 rmap[i+j*rmax]=hypot((float)i,(float)j);
2139// printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
2140 }
2141 }
2142
2143 EMData *fp=new EMData(rmax*2+2,rmax*2,16);
2144
2145 fp->set_complex(1);
2146 fp->to_zero();
2147
2148 // Two vectors in to complex space (kx,ky) and (lx,ly)
2149 // We are computing the bispectrum, f(k).f(l).f*(k+l)
2150 // but integrating out two dimensions, leaving |k|,|l|
2151 for (kx=-rmax+1; kx<rmax; kx++) {
2152 for (ky=-rmax+1; ky<rmax; ky++) {
2153 for (lx=-rmax+1; lx<rmax; lx++) {
2154 for (ly=-rmax+1; ly<rmax; ly++) {
2155 int ax=kx+lx;
2156 int ay=ky+ly;
2157 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
2158 float rr1=rmap[abs(kx)+rmax*abs(ky)];
2159 float rr2=rmap[abs(lx)+rmax*abs(ly)];
2160 int r1=(int)floor(.5+rr1);
2161 int r2=(int)floor(.5+rr2);
2162// if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
2163// float r3=rmap[ax+rmax*ay];
2164 if (r1+r2>=rmax || rr1==0 ||rr2==0) continue;
2165
2166 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
2167 int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5); // projection of k on l 0-31
2168 if (dot<0) dot=16+dot;
2169// int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5+8.0); // projection of k on l 0-15
2170 fp->set_value_at(r1*2,r2,dot,p.real()+fp->get_value_at(r1*2,r2,dot)); // We keep only the real component in anticipation of zero phase sum
2171// fp->set_value_at(r1*2,rmax*2-r2-1, fp->get_value_at(r1*2,r2)); // We keep only the real component in anticipation of zero phase sum
2172// fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2)); // We keep only the real component in anticipation of zero phase sum
2173 fp->set_value_at(r1*2+1,r2,dot,fp->get_value_at(r1*2+1,r2,dot)+1); // a normalization counter
2174 }
2175 }
2176 }
2177 }
2178
2179 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
2180 for (i=0; i<rmax*2; i+=2) {
2181 for (j=0; j<rmax; j++) {
2182 for (h=0; h<16; h++) {
2183 float norm=fp->get_value_at(i+1,j,h);
2184// fp->set_value_at(i,rmax*2-j-1,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
2185// fp->set_value_at(i,j,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
2186 fp->set_value_at(i,rmax*2-j-1,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
2187 fp->set_value_at(i,j,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
2188 // fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0?1.0:norm));
2189 // fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0?1.0:norm));
2190 fp->set_value_at(i+1,j,h,0.0);
2191 }
2192 }
2193 }
2194
2195 free(rmap);
2196 if (type==4) {
2197 EMData *f2=fp->do_ift();
2198 if (f2->get_value_at(0,0,0)<0) f2->mult(-1.0f);
2199 f2->process_inplace("xform.phaseorigin.tocorner");
2200 delete fp;
2201 return f2;
2202 }
2203 return fp;
2204 }
2205 throw UnexpectedBehaviorException("There is not implementation for the parameters you specified");
2206}
EMAN::EMData::make_rotational_footprint_e1
EMData * make_rotational_footprint_e1(bool unwrap=true)
Definition: emdata.cpp:1898
UnexpectedBehaviorException
#define UnexpectedBehaviorException(desc)
Definition: exception.h:400
EMAN::EMData::calc_ccfx
EMData * calc_ccfx(EMData *const with, int y0=0, int y1=-1, bool nosum=false, bool flip=false, bool usez=false)
Calculate Cross-Correlation Function (CCF) in the x-direction and adds them up, result in 1D.
Definition: emdata.cpp:1628

References calc_ccfx(), do_fft(), EMAN::dot(), EMData(), get_clip(), get_xsize(), make_rotational_footprint_e1(), and UnexpectedBehaviorException.

◆ make_rotational_footprint_cmc()

EMData * EMData::make_rotational_footprint_cmc ( bool  unwrap = true)

Definition at line 1795 of file emdata.cpp.

1795 {
1796 ENTERFUNC;
1797 update_stat();
1798 // Note that rotational_footprint caching saves a large amount of time
1799 // but this is at the expense of memory. Note that a policy is hardcoded here,
1800 // that is that caching is only employed when premasked is false and unwrap
1801 // is true - this is probably going to be what is used in most scenarios
1802 // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
1803 // generated by e2speedtest.
1804 if ( rot_fp != 0 && unwrap == true) {
1805 return new EMData(*rot_fp);
1806 }
1807
1808 static EMData obj_filt;
1809 static int updsize = 0; // for threadsafety
1810 EMData* filt = &obj_filt;
1811 filt->set_complex(true);
1812 int delfilt = 0;
1813
1814 // The filter object is nothing more than a cached high pass filter
1815 // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
1816 // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
1817 // set to true, which is used for speed reasons.
1818 if (filt->get_xsize() != nx+2-(nx%2) || filt->get_ysize() != ny ||
1819 filt->get_zsize() != nz ) {
1820 // In case of a threading conflict, this will ignore the cache for this specific computation
1821 if (updsize) {
1822 delfilt=1;
1823 printf("FAILSAFE!!!!\n");
1824 filt=new EMData(nx+2-(nx%2), ny, nz);
1825 filt->to_one();
1826 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1827 }
1828 else {
1829 updsize=1;
1830 filt->set_size(nx+2-(nx%2), ny, nz);
1831 filt->to_one();
1832
1833 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1834 }
1835 }
1836
1837 EMData *ccf = this->calc_mutual_correlation(this, true,filt);
1838
1839 ccf->sub(ccf->get_edge_mean());
1840 EMData *result = ccf->unwrap();
1841 delete ccf; ccf = 0;
1842 updsize=0;
1843 if (delfilt) delete filt;
1844
1845 EXITFUNC;
1846 if ( unwrap == true)
1847 {
1848 // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
1849
1850// Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
1851// to throw any exception
1852// if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
1853
1854// Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
1855// The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
1856 rot_fp = result;
1857 return new EMData(*rot_fp);
1858 }
1859 else return result;
1860}
EMAN::EMData::update_stat
void update_stat() const
Definition: emdata.cpp:3091
EMAN::EMData::calc_mutual_correlation
EMData * calc_mutual_correlation(EMData *with, bool tocorner=false, EMData *filter=0)
Calculates mutual correlation function (MCF) between 2 images.
Definition: emdata.cpp:2209
EMAN::EMData::unwrap
EMData * unwrap(int r1=-1, int r2=-1, int xs=-1, int dx=0, int dy=0, bool do360=false, bool weight_radial=true) const
Maps to polar coordinates from Cartesian coordinates.
Definition: emdata.cpp:2596

References calc_mutual_correlation(), EMData(), ENTERFUNC, EXITFUNC, nx, ny, nz, rot_fp, unwrap(), and update_stat().

Referenced by EMAN::RotationalAligner::align_180_ambiguous().

◆ make_rotational_footprint_e1()

EMData * EMData::make_rotational_footprint_e1 ( bool  unwrap = true)

Definition at line 1898 of file emdata.cpp.

1899{
1900 ENTERFUNC;
1901
1902 update_stat();
1903 // Note that rotational_footprint caching saves a large amount of time
1904 // but this is at the expense of memory. Note that a policy is hardcoded here,
1905 // that is that caching is only employed when premasked is false and unwrap
1906 // is true - this is probably going to be what is used in most scenarios
1907 // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
1908 // generated by e2speedtest.
1909 if ( rot_fp != 0 && unwrap == true) {
1910 return new EMData(*rot_fp);
1911 }
1912
1913 static EMData obj_filt;
1914 EMData* filt = &obj_filt;
1915 filt->set_complex(true);
1916// Region filt_region;
1917
1918// if (nx & 1) {
1919// LOGERR("even image xsize only"); throw ImageFormatException("even image xsize only");
1920// }
1921
1922 int cs = (((nx * 7 / 4) & 0xfffff8) - nx) / 2; // this pads the image to 1 3/4 * size with result divis. by 8
1923
1924 static EMData big_clip;
1925 int big_x = nx+2*cs;
1926 int big_y = ny+2*cs;
1927 int big_z = 1;
1928 if ( nz != 1 ) {
1929 big_z = nz+2*cs;
1930 }
1931
1932
1933 if ( big_clip.get_xsize() != big_x || big_clip.get_ysize() != big_y || big_clip.get_zsize() != big_z ) {
1934 big_clip.set_size(big_x,big_y,big_z);
1935 }
1936 // It is important to set all newly established pixels around the boundaries to the mean
1937 // If this is not done then the associated rotational alignment routine breaks, in fact
1938 // everythin just goes foo.
1939
1940 big_clip.to_value(get_edge_mean());
1941
1942 if (nz != 1) {
1943 big_clip.insert_clip(this,IntPoint(cs,cs,cs));
1944 } else {
1945 big_clip.insert_clip(this,IntPoint(cs,cs,0));
1946 }
1947
1948 // The filter object is nothing more than a cached high pass filter
1949 // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
1950 // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
1951 // set to true, which is used for speed reasons.
1952 if (filt->get_xsize() != big_clip.get_xsize() +2-(big_clip.get_xsize()%2) || filt->get_ysize() != big_clip.get_ysize() ||
1953 filt->get_zsize() != big_clip.get_zsize()) {
1954 filt->set_size(big_clip.get_xsize() + 2-(big_clip.get_xsize()%2), big_clip.get_ysize(), big_clip.get_zsize());
1955 filt->to_one();
1956 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
1957#ifdef EMAN2_USING_CUDA
1958 /*
1959 if(EMData::usecuda == 1 && big_clip.cudarwdata)
1960 {
1961 filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
1962 }
1963 */
1964#endif
1965 }
1966#ifdef EMAN2_USING_CUDA
1967 /*
1968 if(EMData::usecuda == 1 && big_clip.cudarwdata && !filt->cudarwdata)
1969 {
1970 filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
1971 }
1972 */
1973#endif
1974
1975 EMData *mc = big_clip.calc_mutual_correlation(&big_clip, true,filt);
1976 mc->sub(mc->get_edge_mean());
1977
1978 static EMData sml_clip;
1979 int sml_x = nx * 3 / 2;
1980 int sml_y = ny * 3 / 2;
1981 int sml_z = 1;
1982 if ( nz != 1 ) {
1983 sml_z = nz * 3 / 2;
1984 }
1985
1986 if ( sml_clip.get_xsize() != sml_x || sml_clip.get_ysize() != sml_y || sml_clip.get_zsize() != sml_z ) {
1987 sml_clip.set_size(sml_x,sml_y,sml_z); }
1988 if (nz != 1) {
1989 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,-cs+nz/4));
1990 } else {
1991 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,0)); // same as padding change above
1992 }
1993
1994 delete mc; mc = 0;
1995 EMData * result = NULL;
1996
1997 if (nz == 1) {
1998 if (!unwrap) {
1999#ifdef EMAN2_USING_CUDA
2000 //if(EMData::usecuda == 1 && sml_clip.cudarwdata) throw UnexpectedBehaviorException("shap masking is not yet supported by CUDA");
2001#endif
2002 result = sml_clip.process("mask.sharp", Dict("outer_radius", -1, "value", 0));
2003
2004 }
2005 else {
2006 result = sml_clip.unwrap();
2007 }
2008 }
2009 else {
2010 // I am not sure why there is any consideration of non 2D images, but it was here
2011 // in the first port so I kept when I cleaned this function up (d.woolford)
2012// result = clipped_mc;
2013 result = new EMData(sml_clip);
2014 }
2015
2016#ifdef EMAN2_USING_CUDA
2017 //if (EMData::usecuda == 1) sml_clip.roneedsanupdate(); //If we didn't do this then unwrap would use data from the previous call of this function, happens b/c sml_clip is static
2018#endif
2019 EXITFUNC;
2020 if ( unwrap == true)
2021 { // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
2022
2023 // Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
2024 // to throw any exception
2025 if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
2026
2027 // Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
2028 // The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
2029 rot_fp = result;
2030 return new EMData(*rot_fp);
2031 }
2032 else return result;
2033}

References calc_mutual_correlation(), EMData(), ENTERFUNC, EXITFUNC, get_edge_mean(), nx, ny, nz, rot_fp, UnexpectedBehaviorException, unwrap(), and update_stat().

Referenced by EMAN::RotationalAligner::align_180_ambiguous(), and make_footprint().

◆ max_3D_pixel_error()

float EMData::max_3D_pixel_error ( const Transform &  t1,
const Transform &  t2,
float  r 
)

Definition at line 977 of file emdata.cpp.

977 {
978
979 Transform t;
980 int r0 = (int)r;
981 float ddmax = 0.0f;
982
983 t = t2*t1.inverse();
984
985 for (int i=0; i<int(2*M_PI*r0+0.5); i++) {
986 float ang = (float)i/r;
987 Vec3f v = Vec3f(r0*cos(ang), r0*sin(ang), 0.0f);
988 Vec3f d = t*v-v;
989 ddmax = Util::get_max(ddmax, d[0]*d[0]+d[1]*d[1]+d[2]*d[2]);
990 }
991
992 return std::sqrt(ddmax);
993}
EMAN::Transform::inverse
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
EMAN::Util::get_max
static float get_max(float f1, float f2)
Get the maximum of 2 numbers.
Definition: util.h:998

References EMAN::Util::get_max(), EMAN::Transform::inverse(), and sqrt().

◆ oneDfftPolar()

EMData * EMData::oneDfftPolar ( int  size,
float  rmax,
float  MAXR 
)

Definition at line 4213 of file emdata.cpp.

4213 { // sent MAXR value here later!!
4214 float *pcs=get_data();
4215 EMData *imagepcsfft = new EMData;
4216 imagepcsfft->set_size((size+2), (int)MAXR+1, 1);
4217 float *d=imagepcsfft->get_data();
4218
4219 EMData *data_in=new EMData;
4220 data_in->set_size(size,1,1);
4221 float *in=data_in->get_data();
4222
4223 for(int row=0; row<=(int)MAXR; ++row){
4224 if(row<=(int)rmax) {
4225 for(int i=0; i<size;++i) in[i] = pcs[i+row*size]; // ming
4226 data_in->set_complex(false);
4227 data_in->do_fft_inplace();
4228 for(int j=0;j<size+2;j++) d[j+row*(size+2)]=in[j];
4229 }
4230 else for(int j=0;j<size+2;j++) d[j+row*(size+2)]=0.0;
4231 }
4232 imagepcsfft->update();
4233 delete data_in;
4234 return imagepcsfft;
4235}

References EMData(), and get_data().

◆ rotate() [1/2]

void EMData::rotate ( const Transform &  t)

Rotate this image.

DEPRECATED USE EMData::transform()

Parameters
tTransformation rotation.

Definition at line 971 of file emdata.cpp.

972{
973 cout << "Deprecation warning in EMData::rotate. Please consider using EMData::transform() instead " << endl;
974 transform(t);
975}

References transform().

◆ rotate() [2/2]

void EMData::rotate ( float  az,
float  alt,
float  phi 
)

Rotate this image.

DEPRECATED USE EMData::Transform

Parameters
azRotation euler angle az in EMAN convention.
altRotation euler angle alt in EMAN convention.
phiRotation euler angle phi in EMAN convention.

Definition at line 959 of file emdata.cpp.

960{
961 Dict d("type","eman");
962 d["az"] = az;
963 d["alt"] = alt;
964 d["phi"] = phi;
965 Transform t(d);
966 transform(t);
967}

References transform().

◆ rotate_180()

void EMAN::EMData::rotate_180 ( )
inline

Fast rotation by 180 degrees.

Square 2D image only.

Exceptions
ImageFormatExceptionIf the image is not square.
ImageDimensionExceptionIf the image is not 2D.

Definition at line 350 of file emdata.h.

350 {
351 ENTERFUNC;
352 process_inplace("math.rotate.180",Dict());
353 EXITFUNC;
354 }
process_inplace
void process_inplace(const string &processorname, const Dict &params=Dict())
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....

References ENTERFUNC, EXITFUNC, and process_inplace().

◆ rotate_translate() [1/3]

void EMAN::EMData::rotate_translate ( const Transform &  t)
inline

Apply a transformation to the image.

DEPRECATED USE EMData::Transform

Parameters
ttransform object that describes the transformation to be applied to the image.

Definition at line 306 of file emdata.h.

306 {
307 cout << "Deprecation warning. Please consider using EMData::transform() instead " << endl;
308 transform(t); }

References transform().

Referenced by rotate_translate().

◆ rotate_translate() [2/3]

void EMData::rotate_translate ( float  az,
float  alt,
float  phi,
float  dx,
float  dy,
float  dz 
)

Rotate then translate the image.

DEPRECATED USE EMData::Transform

Parameters
azRotation euler angle az in EMAN convention.
altRotation euler angle alt in EMAN convention.
phiRotation euler angle phi in EMAN convention.
dxTranslation distance in x direction.
dyTranslation distance in y direction.
dzTranslation distance in z direction.

Definition at line 995 of file emdata.cpp.

996{
997 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
998// Transform3D t( az, alt, phi,Vec3f(dx, dy, dz));
999 Transform t;
1000 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
1001 t.set_trans(dx, dy, dz);
1002 rotate_translate(t);
1003}
EMAN::EMData::rotate_translate
void rotate_translate(const Transform &t)
Apply a transformation to the image.
Definition: emdata.h:306
EMAN::Transform::set_rotation
void set_rotation(const Dict &rotation)
Set a rotation using a specific Euler type and the dictionary interface Works for all Euler types.
Definition: transform.cpp:519
EMAN::Transform::set_trans
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036

References rotate_translate(), EMAN::Transform::set_rotation(), and EMAN::Transform::set_trans().

◆ rotate_translate() [3/3]

void EMData::rotate_translate ( float  az,
float  alt,
float  phi,
float  dx,
float  dy,
float  dz,
float  pdx,
float  pdy,
float  pdz 
)

Rotate then translate the image.

DEPRECATED USE EMData::Transform

Parameters
azRotation euler angle az in EMAN convention.
altRotation euler angle alt in EMAN convention.
phiRotation euler angle phi in EMAN convention.
dxTranslation distance in x direction.
dyTranslation distance in y direction.
dzTranslation distance in z direction.
pdxPretranslation distance in x direction.
pdyPretranslation distance in y direction.
pdzPretranslation distance in z direction.

Definition at line 1006 of file emdata.cpp.

1008{
1009 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
1010// Transform3D t(Vec3f(dx, dy, dz), az, alt, phi, Vec3f(pdx,pdy,pdz));
1011// rotate_translate(t);
1012
1013 Transform t;
1014 t.set_pre_trans(Vec3f(dx, dy, dz));
1015 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
1016 t.set_trans(pdx, pdy, pdz);
1017 rotate_translate(t);
1018}
EMAN::Transform::set_pre_trans
void set_pre_trans(const type &v)
Set the translational component of the matrix as though it was MSRT_ not MTSR, where T_ is the pre tr...
Definition: transform.h:556

References rotate_translate(), EMAN::Transform::set_pre_trans(), EMAN::Transform::set_rotation(), and EMAN::Transform::set_trans().

◆ rotate_x()

void EMData::rotate_x ( int  dx)

This performs a translation of each line along x with wraparound.

This is equivalent to a rotation when performed on 'unwrapped' maps.

Parameters
dxTranslation distance align x direction.
Exceptions
ImageDimensionExceptionIf the image is 3D.

Definition at line 1209 of file emdata.cpp.

1210{
1211 ENTERFUNC;
1212
1213 if (get_ndim() > 2) {
1214 throw ImageDimensionException("no 3D image");
1215 }
1216
1217
1218 size_t row_size = nx * sizeof(float);
1219 float *tmp = (float*)EMUtil::em_malloc(row_size);
1220 float * data = get_data();
1221
1222 for (int y = 0; y < ny; y++) {
1223 int y_nx = y * nx;
1224 for (int x = 0; x < nx; x++) {
1225 tmp[x] = data[y_nx + (x + dx) % nx];
1226 }
1227 EMUtil::em_memcpy(&data[y_nx], tmp, row_size);
1228 }
1229
1230 update();
1231 if( tmp )
1232 {
1233 delete[]tmp;
1234 tmp = 0;
1235 }
1236 EXITFUNC;
1237}
EMAN::EMUtil::em_malloc
static void * em_malloc(const size_t size)
Definition: emutil.h:366

References EMAN::EMUtil::em_malloc(), EMAN::EMUtil::em_memcpy(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, nx, ny, update(), x, and y.

◆ save_byteorder_to_dict()

void EMData::save_byteorder_to_dict ( ImageIO *  imageio)
private

Definition at line 4349 of file emdata.cpp.

4350{
4351 string image_endian = "ImageEndian";
4352 string host_endian = "HostEndian";
4353
4354 if (imageio->is_image_big_endian()) {
4355 attr_dict[image_endian] = "big";
4356 }
4357 else {
4358 attr_dict[image_endian] = "little";
4359 }
4360
4361 if (ByteOrder::is_host_big_endian()) {
4362 attr_dict[host_endian] = "big";
4363 }
4364 else {
4365 attr_dict[host_endian] = "little";
4366 }
4367}
EMAN::ByteOrder::is_host_big_endian
static bool is_host_big_endian()
Definition: byteorder.cpp:40
EMAN::ImageIO::is_image_big_endian
virtual bool is_image_big_endian()=0
Is this image in big endian or not.

References attr_dict, EMAN::ByteOrder::is_host_big_endian(), and EMAN::ImageIO::is_image_big_endian().

Referenced by read_binedimage().

◆ scale()

void EMData::scale ( float  scale_factor)

scale the image by a factor.

Parameters
scale_factorscale factor.

Definition at line 886 of file emdata.cpp.

887{
888 ENTERFUNC;
889 Transform t;
890 t.set_scale(s);
891 transform(t);
892 EXITFUNC;
893}
EMAN::Transform::set_scale
void set_scale(const float &scale)
Set the scale.
Definition: transform.cpp:1123

References ENTERFUNC, EXITFUNC, EMAN::Transform::set_scale(), and transform().

Referenced by do_ift(), do_ift_inplace(), render_amp24(), render_ap24(), and scale_pixel().

◆ set_xyz_origin()

void EMData::set_xyz_origin ( float  origin_x,
float  origin_y,
float  origin_z 
)

Set the x,y,z origin of the image.

Parameters
origin_xthe x origin
origin_ythe y origin
origin_zthe z origin

Definition at line 3286 of file emdata.cpp.

3287{
3288 attr_dict["origin_x"] = origin_x;
3289 attr_dict["origin_y"] = origin_y;
3290 attr_dict["origin_z"] = origin_z;
3291}

References attr_dict.

Referenced by clip_inplace(), and get_clip().

◆ setup4slice()

float * EMData::setup4slice ( bool  redo = true)

Set up for fftslice operations.

When interpolating in fourier space there is a little problem when we get close to x=0, since f(-x,-y,-z) = f(x,y,z)* . So this makes a supplementary array that allows for up to +-2 point interpolation all the way to the origin in x.

3D only; complex image only

Parameters
redoIf true, recalculate the supplementary array.
Exceptions
ImageFormatExceptionIf the image is not a complex image.
ImageDimensionExceptionIf the image is not 3D.
Returns
The supplementary array.

Definition at line 822 of file emdata.cpp.

823{
824 ENTERFUNC;
825
826 if (!is_complex()) {
827 throw ImageFormatException("complex image only");
828 }
829
830 if (get_ndim() != 3) {
831 throw ImageDimensionException("3D only");
832 }
833
834 if (supp) {
835 if (redo) {
836 EMUtil::em_free(supp);
837 supp = 0;
838 }
839 else {
840 EXITFUNC;
841 return supp;
842 }
843 }
844
845 const int SUPP_ROW_SIZE = 8;
846 const int SUPP_ROW_OFFSET = 4;
847 const int supp_size = SUPP_ROW_SIZE + SUPP_ROW_OFFSET;
848
849 supp = (float *) EMUtil::em_calloc(supp_size * ny * nz, sizeof(float));
850 int nxy = nx * ny;
851 int supp_xy = supp_size * ny;
852 float * data = get_data();
853
854 for (int z = 0; z < nz; z++) {
855 size_t cur_z1 = z * nxy;
856 size_t cur_z2 = z * supp_xy;
857
858 for (int y = 0; y < ny; y++) {
859 size_t cur_y1 = y * nx;
860 size_t cur_y2 = y * supp_size;
861
862 for (int x = 0; x < SUPP_ROW_SIZE; x++) {
863 size_t k = (x + SUPP_ROW_OFFSET) + cur_y2 + cur_z2;
864 supp[k] = data[x + cur_y1 + cur_z1];
865 }
866 }
867 }
868
869 for (int z = 1, zz = nz - 1; z < nz; z++, zz--) {
870 size_t cur_z1 = zz * nxy;
871 size_t cur_z2 = z * supp_xy;
872
873 for (int y = 1, yy = ny - 1; y < ny; y++, yy--) {
874 supp[y * 12 + cur_z2] = data[4 + yy * nx + cur_z1];
875 supp[1 + y * 12 + cur_z2] = -data[5 + yy * nx + cur_z1];
876 supp[2 + y * 12 + cur_z2] = data[2 + yy * nx + cur_z1];
877 supp[3 + y * 12 + cur_z2] = -data[3 + yy * nx + cur_z1];
878 }
879 }
880
881 EXITFUNC;
882 return supp;
883}
EMAN::EMUtil::em_calloc
static void * em_calloc(const size_t nmemb, const size_t size)
Definition: emutil.h:370
EMAN::EMUtil::em_free
static void em_free(void *data)
Definition: emutil.h:380

References EMAN::EMUtil::em_calloc(), EMAN::EMUtil::em_free(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, ImageFormatException, is_complex(), nx, nxy, ny, nz, supp, x, and y.

Referenced by EMAN::GaussFFTProjector::interp_ft_3d().

◆ transform()

void EMAN::EMData::transform ( const Transform &  t)
inline

Rotate then translate the image.

DEPRECATED USE EMData::Transform

Parameters
tThe rotation and translation transformation to be done. Transform the image
tthe transform object that describes the transformation to be applied to the image.

Definition at line 295 of file emdata.h.

295 {
296 ENTERFUNC;
297 process_inplace("xform",Dict("transform",(Transform*)(&t)));
298 //update(); no need, process_inplace did it
299 EXITFUNC;
300 }

References ENTERFUNC, EXITFUNC, and process_inplace().

Referenced by cut_slice(), do_radon(), EMAN::BackProjectionReconstructor::insert_slice(), EMAN::TestImageSinewave::process_inplace(), rotate(), rotate_translate(), scale(), and uncut_slice().

◆ translate() [1/4]

void EMData::translate ( const Vec3f &  translation)

Translate this image.

Parameters
translationThe translation distance vector.

Definition at line 940 of file emdata.cpp.

941{
942 ENTERFUNC;
943
944 if( translation[0] == 0.0f && translation[1] == 0.0f && translation[2] == 0.0f ) {
945 EXITFUNC;
946 return;
947 }
948
949 Transform* t = new Transform();
950 t->set_trans(translation);
951 process_inplace("xform",Dict("transform",t));
952 delete t;
953
954 all_translation += translation;
955 EXITFUNC;
956}

References all_translation, ENTERFUNC, EXITFUNC, process_inplace(), and EMAN::Transform::set_trans().

◆ translate() [2/4]

void EMData::translate ( const Vec3i &  translation)

Translate this image.

integer only translation could be done faster, without interpolation.

Parameters
translationThe translation distance vector.

Definition at line 920 of file emdata.cpp.

921{
922 ENTERFUNC;
923
924 //if traslation is 0, do nothing
925 if( translation[0] == 0 && translation[1] == 0 && translation[2] == 0) {
926 EXITFUNC;
927 return;
928 }
929
930 Dict params("trans",static_cast< vector<int> >(translation));
931 process_inplace("xform.translate.int",params);
932
933 // update() - clip_inplace does the update
934 all_translation += translation;
935
936 EXITFUNC;
937}

References all_translation, ENTERFUNC, EXITFUNC, and process_inplace().

◆ translate() [3/4]

void EMData::translate ( float  dx,
float  dy,
float  dz 
)

Translate this image.

Parameters
dxTranslation distance in x direction.
dyTranslation distance in y direction.
dzTranslation distance in z direction.

Definition at line 904 of file emdata.cpp.

905{
906 ENTERFUNC;
907 int dx_ = Util::round(dx);
908 int dy_ = Util::round(dy);
909 int dz_ = Util::round(dz);
910 if( ( (dx-dx_) == 0 ) && ( (dy-dy_) == 0 ) && ( (dz-dz_) == 0 )) {
911 translate(dx_, dy_, dz_);
912 }
913 else {
914 translate(Vec3f(dx, dy, dz));
915 }
916 EXITFUNC;
917}
EMAN::EMData::translate
void translate(float dx, float dy, float dz)
Translate this image.
Definition: emdata.cpp:904

References ENTERFUNC, EXITFUNC, EMAN::Util::round(), and translate().

Referenced by EMAN::PhaseToMassCenterProcessor::process_inplace(), EMAN::ToCenterProcessor::process_inplace(), EMAN::ToMassCenterProcessor::process_inplace(), EMAN::ACFCenterProcessor::process_inplace(), EMAN::GaussFFTProjector::project3d(), EMAN::RT2Dto3DTreeAligner::testort(), and translate().

◆ translate() [4/4]

void EMData::translate ( int  dx,
int  dy,
int  dz 
)

Translate this image.

integer only translation could be done faster, without interpolation.

Parameters
dxTranslation distance in x direction.
dyTranslation distance in y direction.
dzTranslation distance in z direction.

Definition at line 896 of file emdata.cpp.

897{
898 ENTERFUNC;
899 translate(Vec3i(dx, dy, dz));
900 EXITFUNC;
901}
EMAN::Vec3i
Vec3< int > Vec3i
Definition: vec3.h:694

References ENTERFUNC, EXITFUNC, and translate().

◆ uncut_slice()

void EMData::uncut_slice ( EMData *const  map,
const Transform &  tr 
) const

Opposite of the cut_slice().

It will take a slice and insert the data into a real 3D map. It does not interpolate, it uses the nearest neighbor.

Parameters
mapThe real 3D map.
trOrientation of the slice.
Exceptions
NullPointerExceptionIf map is NULL.
ImageDimensionExceptionIf this image is not 2D.
ImageDimensionExceptionIf map image is not 3D.
ImageFormatExceptionIf this image is complex
ImageFormatExceptionIf map is complex
Author
David Woolford (adapted from an original version by Steve Ludtke)
Date
Feb 2008

Definition at line 4237 of file emdata.cpp.

4238{
4239 ENTERFUNC;
4240
4241 if (!map) throw NullPointerException("NULL image");
4242 // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
4243 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
4244 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
4245 // Now check for complex images - this is really just being thorough
4246 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
4247 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
4248
4249// Transform3D r( 0, 0, 0); // EMAN by default
4250// if (!ort) {
4251// ort = &r;
4252// }
4253
4254 float *ddata = map->get_data();
4255 float *sdata = get_data();
4256
4257 int map_nx = map->get_xsize();
4258 int map_ny = map->get_ysize();
4259 int map_nz = map->get_zsize();
4260 int map_nxy = map_nx * map_ny;
4261 float map_nz_round_limit = (float) map_nz-0.5f;
4262 float map_ny_round_limit = (float) map_ny-0.5f;
4263 float map_nx_round_limit = (float) map_nx-0.5f;
4264/*
4265 Vec3f posttrans = ort->get_posttrans();
4266 Vec3f pretrans = ort->get_pretrans();*/
4267
4268 int ymax = ny/2;
4269 if ( ny % 2 == 1 ) ymax += 1;
4270 int xmax = nx/2;
4271 if ( nx % 2 == 1 ) xmax += 1;
4272 for (int y = -ny/2; y < ymax; y++) {
4273 for (int x = -nx/2; x < xmax; x++) {
4274 Vec3f coord(x,y,0);
4275 Vec3f soln = transform*coord;
4276// float xx = (x+pretrans[0]) * (*ort)[0][0] + (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
4277// float yy = (x+pretrans[0]) * (*ort)[1][0] + (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
4278// float zz = (x+pretrans[0]) * (*ort)[2][0] + (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
4279//
4280// xx += map_nx/2;
4281// yy += map_ny/2;
4282// zz += map_nz/2;
4283//
4284 float xx = soln[0]+map_nx/2;
4285 float yy = soln[1]+map_ny/2;
4286 float zz = soln[2]+map_nz/2;
4287
4288 // These 0.5 offsets are here because the round function rounds to the nearest whole number.
4289 if (xx < -0.5 || yy < -0.5 || zz < -0.5 || xx >= map_nx_round_limit || yy >= map_ny_round_limit || zz >= map_nz_round_limit) continue;
4290
4291 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
4292 int l = (x+nx/2) + (y+ny/2) * nx;
4293 ddata[k] = sdata[l];
4294 }
4295 }
4296
4297 map->update();
4298 EXITFUNC;
4299}

References ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, ImageFormatException, is_complex(), NullPointerException, nx, ny, EMAN::Util::round(), transform(), x, and y.

◆ unwrap_largerR()

EMData * EMData::unwrap_largerR ( int  r1,
int  r2,
int  xs,
float  rmax_f 
)

Definition at line 4136 of file emdata.cpp.

4136 {
4137 float *d,*dd;
4138 int do360=2;
4139 int rmax = (int)(rmax_f+0.5f);
4140 unsigned long i;
4141 unsigned int nvox=get_xsize()*get_ysize();//ming
4142 float maxmap=0.0f, minmap=0.0f;
4143 float temp=0.0f, diff_den=0.0f, norm=0.0f;
4144 float cut_off_va =0.0f;
4145
4146 d=get_data();
4147 maxmap=-1000000.0f;
4148 minmap=1000000.0f;
4149 for (i=0;i<nvox;i++){
4150 if(d[i]>maxmap) maxmap=d[i];
4151 if(d[i]<minmap) minmap=d[i];
4152 }
4153 diff_den = maxmap-minmap;
4154 for(i=0;i<nvox;i++) {
4155 temp = (d[i]-minmap)/diff_den;
4156 if(cut_off_va != 0.0) { // cut off the lowerset ?% noisy information
4157 if(temp < cut_off_va)
4158 d[i] = 0.0f; // set the empty part density=0.0
4159 else
4160 d[i] = temp-cut_off_va;
4161 }
4162 else d[i] = temp;
4163 }
4164
4165 for(i=0;i<nvox;i++) {
4166 temp=d[i];
4167 norm += temp*temp;
4168 }
4169 for(i=0;i<nvox;i++) d[i] /= norm; // y' = y/norm(y)
4170
4171 if (xs<1) {
4172 xs = (int) floor(do360*M_PI*get_ysize()/4); // ming
4173 xs=Util::calc_best_fft_size(xs); // ming
4174 }
4175 if (r1<0) r1=0;
4176 float maxext=ceil(0.6f*std::sqrt((float)(get_xsize()*get_xsize()+get_ysize()*get_ysize())));// ming add std::
4177
4178 if (r2<r1) r2=(int)maxext;
4179 EMData *ret = new EMData;
4180
4181 ret->set_size(xs,r2+1,1);
4182
4183 dd=ret->get_data();
4184
4185 for (int i=0; i<xs; i++) {
4186 float si=sin(i*M_PI*2/xs);
4187 float co=cos(i*M_PI*2/xs);
4188 for (int r=0; r<=maxext; r++) {
4189 float x=(r+r1)*co+get_xsize()/2; // ming
4190 float y=(r+r1)*si+get_ysize()/2; // ming
4191 if(x<0.0 || x>=get_xsize()-1.0 || y<0.0 || y>=get_ysize()-1.0 || r>rmax){ //Ming , ~~~~ rmax need pass here
4192 for(;r<=r2;r++) // here r2=MAXR
4193 dd[i+r*xs]=0.0;
4194 break;
4195 }
4196 int x1=(int)floor(x);
4197 int y1=(int)floor(y);
4198 float t=x-x1;
4199 float u=y-y1;
4200 float f11= d[x1+y1*get_xsize()]; // ming
4201 float f21= d[(x1+1)+y1*get_xsize()]; // ming
4202 float f12= d[x1+(y1+1)*get_xsize()]; // ming
4203 float f22= d[(x1+1)+(y1+1)*get_xsize()]; // ming
4204 dd[i+r*xs] = (1-t)*(1-u)*f11+t*(1-u)*f21+t*u*f22+(1-t)*u*f12;
4205 }
4206 }
4207 update();
4208 ret->update();
4209 return ret;
4210}
EMAN::Util::calc_best_fft_size
static int calc_best_fft_size(int low)
Search the best FFT size with good primes.
Definition: util.cpp:1021
get_ysize
int get_ysize() const
Get the image y-dimensional size.
Definition: emdata_metadata.h:553

References EMAN::Util::calc_best_fft_size(), EMData(), get_data(), get_xsize(), get_ysize(), sqrt(), update(), x, and y.

◆ update_stat()

void EMData::update_stat ( ) const
private

Definition at line 3091 of file emdata.cpp.

3092{
3093 ENTERFUNC;
3094// printf("update stat %f %d\n",(float)attr_dict["mean"],flags);
3095 if (!(flags & EMDATA_NEEDUPD))
3096 {
3097 EXITFUNC;
3098 return;
3099 }
3100 if (rdata==0) return;
3101
3102 float* data = get_data();
3103 float max = -FLT_MAX;
3104 float min = -max;
3105
3106 double sum = 0;
3107 double square_sum = 0;
3108
3109 int step = 1;
3110 if (is_complex() && !is_ri()) {
3111 step = 2;
3112 }
3113
3114 int n_nonzero = 0;
3115
3116 size_t size = (size_t)nx*ny*nz;
3117
3118 for (size_t i = 0; i < size; i += step) {
3119 float v = data[i];
3120 max = Util::get_max(max, v);
3121 min = Util::get_min(min, v);
3122 sum += v;
3123 square_sum += v * (double)(v);
3124 if (v != 0) n_nonzero++;
3125 }
3126
3127 size_t n = size / step;
3128 double mean = sum / n;
3129 double var = (square_sum - sum*sum / n) / (n-1);
3130 double sigma = var >= 0.0 ? std::sqrt(var) : 0.0;
3131 n_nonzero = Util::get_max(1, n_nonzero);
3132 double varn = (square_sum - sum*sum / n_nonzero) / (n_nonzero-1);
3133 double sigma_nonzero = varn >= 0.0 ? std::sqrt(varn) : 0.0;
3134 double mean_nonzero = sum / n_nonzero; // previous version overcounted! G2
3135
3136 attr_dict["minimum"] = min;
3137 attr_dict["maximum"] = max;
3138 attr_dict["mean"] = (float)(mean);
3139 attr_dict["sigma"] = (float)(sigma);
3140 attr_dict["square_sum"] = (float)(square_sum);
3141 attr_dict["mean_nonzero"] = (float)(mean_nonzero);
3142 attr_dict["sigma_nonzero"] = (float)(sigma_nonzero);
3143 attr_dict["is_complex"] = (int) is_complex();
3144 attr_dict["is_complex_ri"] = (int) is_ri();
3145
3146 flags &= ~EMDATA_NEEDUPD;
3147
3148 if (rot_fp != 0)
3149 {
3150 delete rot_fp; rot_fp = 0;
3151 }
3152
3153 EXITFUNC;
3154// printf("done stat %f %f %f\n",(float)mean,(float)max,(float)sigma);
3155}
EMAN::Util::get_min
static int get_min(int f1, int f2)
Get the minimum of 2 numbers.
Definition: util.h:922

References attr_dict, EMDATA_NEEDUPD, ENTERFUNC, EXITFUNC, flags, get_data(), EMAN::Util::get_max(), EMAN::Util::get_min(), is_complex(), is_ri(), nx, ny, nz, rdata, rot_fp, and sqrt().

Referenced by get_attr(), get_attr_dict(), make_rotational_footprint(), make_rotational_footprint_cmc(), and make_rotational_footprint_e1().

◆ window_center()

EMData * EMData::window_center ( int  l)

Window the center of an image.

Often an image is padded with zeros for fourier interpolation. In that case the desired lxlxl volume (or lxl area) lies in the center of a larger volume (or area). This routine creates a new object that contains only the desired window. (This routine is a thin wrapper around get_clip.)

Parameters
lLength of the window.
Returns
An image object that has been windowed.

Definition at line 776 of file emdata.cpp.

776 {
777 ENTERFUNC;
778 // sanity checks
779 int n = nx;
780 if (is_complex()) {
781 LOGERR("Need real-space data for window_center()");
782 throw ImageFormatException(
783 "Complex input image; real-space expected.");
784 }
785 if (is_fftpadded()) {
786 // image has been fft-padded, compute the real-space size
787 n -= (2 - int(is_fftodd()));
788 }
789 int corner = n/2 - l/2;
790 int ndim = get_ndim();
791 EMData* ret;
792 switch (ndim) {
793 case 3:
794 if ((n != ny) || (n != nz)) {
795 LOGERR("Need the real-space image to be cubic.");
796 throw ImageFormatException(
797 "Need cubic real-space image.");
798 }
799 ret = get_clip(Region(corner, corner, corner, l, l, l));
800 break;
801 case 2:
802 if (n != ny) {
803 LOGERR("Need the real-space image to be square.");
804 throw ImageFormatException(
805 "Need square real-space image.");
806 }
807 //cout << "Using corner " << corner << endl;
808 ret = get_clip(Region(corner, corner, l, l));
809 break;
810 case 1:
811 ret = get_clip(Region(corner, l));
812 break;
813 default:
814 throw ImageDimensionException(
815 "window_center only supports 1-d, 2-d, and 3-d images");
816 }
817 return ret;
818 EXITFUNC;
819}
is_fftpadded
bool is_fftpadded() const
Is this image already extended along x for ffts?
Definition: emdata_metadata.h:831

References ENTERFUNC, EXITFUNC, get_clip(), get_ndim(), ImageDimensionException, ImageFormatException, is_complex(), is_fftodd(), is_fftpadded(), LOGERR, nx, ny, and nz.

Referenced by EMAN::FourierGriddingProjector::project3d().

◆ zero_corner_circulant()

void EMData::zero_corner_circulant ( const int  radius = 0)

Zero the pixels in the bottom left corner of the image If radius is greater than 1, than circulant zeroing occurs assuming that the center of operation starts in the bottom left corner and proceed outwards to the NE and backwards in a circulant fashion towards the SW.

Intended to zero the area corresponding to the middle of the image, as generated by calc_ccf

Parameters
radiusthe radius to the zeroing operation
Exceptions
ImageDimensionExceptionif nx > 1 and nx < 2*radius + 1
ImageDimensionExceptionif ny > 1 and ny < 2*radius + 1
ImageDimensionExceptionif nz > 1 and nz < 2*radius + 1

Definition at line 1462 of file emdata.cpp.

1463{
1464 if ( nz > 1 && nz < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nz is too small");
1465 if ( ny > 1 && ny < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - ny is too small");
1466 if ( nx > 1 && nx < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nx is too small");
1467
1468 int it_z = radius;
1469 int it_y = radius;
1470 int it_x = radius;
1471
1472 if ( nz == 1 ) it_z = 0;
1473 if ( ny == 1 ) it_y = 0;
1474 if ( nx == 1 ) it_z = 0;
1475
1476 if ( nz == 1 && ny == 1 )
1477 {
1478 for ( int x = -it_x; x <= it_x; ++x )
1479 get_value_at_wrap(x) = 0;
1480
1481 }
1482 else if ( nz == 1 )
1483 {
1484 for ( int y = -it_y; y <= it_y; ++y)
1485 for ( int x = -it_x; x <= it_x; ++x )
1486 get_value_at_wrap(x,y) = 0;
1487 }
1488 else
1489 {
1490 for( int z = -it_z; z <= it_z; ++z )
1491 for ( int y = -it_y; y <= it_y; ++y)
1492 for ( int x = -it_x; x < it_x; ++x )
1493 get_value_at_wrap(x,y,z) = 0;
1494
1495 }
1496
1497}

References get_value_at_wrap(), ImageDimensionException, nx, ny, nz, x, and y.

Friends And Related Function Documentation

◆ GLUtil

friend class GLUtil
friend

Definition at line 83 of file emdata.h.

Member Data Documentation

◆ all_translation

Vec3f EMAN::EMData::all_translation
private

translation from the original location

Definition at line 854 of file emdata.h.

Referenced by copy_head(), operator=(), and translate().

◆ attr_dict

Dict EMAN::EMData::attr_dict
mutableprivate

to store all image header info

Definition at line 833 of file emdata.h.

Referenced by clip_inplace(), copy_head(), debug_print_parms(), del_attr(), EMData(), get_attr(), get_attr_dict(), get_clip(), get_ctf(), get_rotated_clip(), get_top_half(), operator=(), read_binedimage(), save_byteorder_to_dict(), scale_pixel(), set_xyz_origin(), update_stat(), and write_lst().

◆ changecount

int EMAN::EMData::changecount
private

Definition at line 846 of file emdata.h.

Referenced by get_attr(), get_attr_dict(), and operator=().

◆ flags

int EMAN::EMData::flags
mutableprivate

CTF data All CTF data become attribute ctf(vector<float>) in attr_dict –Grant Tang.

flags

Definition at line 844 of file emdata.h.

Referenced by copy_head(), get_attr(), operator=(), render_ap24(), and update_stat().

◆ nx

int EMAN::EMData::nx
private

image size

Definition at line 848 of file emdata.h.

Referenced by absi(), add_complex_at(), add_incoherent(), apply_radial_func(), calc_az_dist(), calc_ccf(), calc_ccfx(), calc_center_density(), calc_center_of_mass(), calc_dist(), calc_fast_sigma_image(), calc_flcf(), calc_highest_locations(), calc_hist(), calc_max_index(), calc_max_location(), calc_max_location_wrap(), calc_min_index(), calc_min_location(), calc_n_highest_locations(), calc_radial_dist(), calc_sigma_diff(), cconj(), clip_inplace(), common_lines(), compute_missingwedge(), copy_head(), cut_slice(), do_fft(), do_fft_inplace(), do_ift(), do_ift_inplace(), do_radon(), dot_rotate_translate(), EMData(), extract_box(), find_pixels_with_value(), get_2dcview(), get_2dview(), get_3dcview(), get_3dview(), get_attr(), get_attr_dict(), get_circle_mean(), get_clip(), get_col(), get_complex_at(), get_complex_index(), get_data_pickle(), get_edge_mean(), get_fft_amplitude(), get_fft_amplitude2D(), get_fft_phase(), get_rotated_clip(), get_row(), get_top_half(), get_value_at_wrap(), little_big_dot(), make_rotational_footprint_cmc(), make_rotational_footprint_e1(), operator=(), phase(), read_binedimage(), real2complex(), render_amp24(), render_ap24(), ri2ap(), ri2inten(), rotate_x(), set_complex_at(), set_data_pickle(), setup4slice(), sget_value_at(), to_value(), uncut_slice(), unwrap(), update_stat(), window_center(), and zero_corner_circulant().

◆ nxy

int EMAN::EMData::nxy
private

Definition at line 848 of file emdata.h.

Referenced by add_complex_at(), calc_center_of_mass(), calc_highest_locations(), calc_max_location(), calc_min_location(), calc_n_highest_locations(), extract_box(), get_complex_at(), get_complex_index(), get_edge_mean(), get_value_at_wrap(), set_complex_at(), setup4slice(), and sget_value_at().

◆ nxyz

size_t EMAN::EMData::nxyz
private

Definition at line 849 of file emdata.h.

Referenced by add_complex_at(), do_fft(), get_complex_index(), log(), log10(), and sqrt().

◆ ny

int EMAN::EMData::ny
private

Definition at line 848 of file emdata.h.

Referenced by absi(), add_complex_at(), add_incoherent(), apply_radial_func(), calc_az_dist(), calc_ccf(), calc_ccfx(), calc_center_density(), calc_center_of_mass(), calc_dist(), calc_fast_sigma_image(), calc_flcf(), calc_highest_locations(), calc_hist(), calc_max_index(), calc_max_location(), calc_max_location_wrap(), calc_min_index(), calc_min_location(), calc_n_highest_locations(), calc_radial_dist(), calc_sigma_diff(), cconj(), clip_inplace(), common_lines(), compute_missingwedge(), copy_head(), cut_slice(), do_fft(), do_fft_inplace(), do_ift(), do_ift_inplace(), do_radon(), dot_rotate_translate(), EMData(), extract_box(), find_pixels_with_value(), get_2dcview(), get_2dview(), get_3dcview(), get_3dview(), get_attr(), get_attr_dict(), get_circle_mean(), get_clip(), get_col(), get_complex_at(), get_complex_index(), get_data_pickle(), get_edge_mean(), get_fft_amplitude(), get_fft_amplitude2D(), get_fft_phase(), get_rotated_clip(), get_top_half(), get_value_at_wrap(), little_big_dot(), make_rotational_footprint_cmc(), make_rotational_footprint_e1(), operator=(), phase(), read_binedimage(), real2complex(), render_amp24(), render_ap24(), ri2ap(), ri2inten(), rotate_x(), set_complex_at(), set_data_pickle(), setup4slice(), sget_value_at(), to_value(), uncut_slice(), unwrap(), update_stat(), window_center(), and zero_corner_circulant().

◆ nz

int EMAN::EMData::nz
private

Definition at line 848 of file emdata.h.

Referenced by absi(), add_complex_at(), add_incoherent(), apply_radial_func(), calc_ccf(), calc_ccfx(), calc_center_density(), calc_center_of_mass(), calc_fast_sigma_image(), calc_flcf(), calc_highest_locations(), calc_hist(), calc_max_location(), calc_max_location_wrap(), calc_min_location(), calc_n_highest_locations(), calc_radial_dist(), calc_sigma_diff(), cconj(), clip_inplace(), compute_missingwedge(), copy_head(), do_fft(), do_fft_inplace(), do_ift(), do_ift_inplace(), EMData(), extract_box(), find_pixels_with_value(), get_3dcview(), get_3dview(), get_attr(), get_attr_dict(), get_circle_mean(), get_clip(), get_complex_at(), get_complex_index(), get_data_pickle(), get_edge_mean(), get_fft_amplitude(), get_fft_amplitude2D(), get_fft_phase(), get_rotated_clip(), get_top_half(), get_value_at_wrap(), little_big_dot(), make_rotational_footprint_cmc(), make_rotational_footprint_e1(), operator=(), phase(), read_binedimage(), real2complex(), ri2ap(), ri2inten(), set_complex_at(), set_data_pickle(), setup4slice(), sget_value_at(), to_value(), update_stat(), window_center(), and zero_corner_circulant().

◆ path

string EMAN::EMData::path
private

Definition at line 857 of file emdata.h.

Referenced by copy_head(), get_clip(), and operator=().

◆ pathnum

int EMAN::EMData::pathnum
private

Definition at line 858 of file emdata.h.

Referenced by copy_head(), get_clip(), and operator=().

◆ rdata

float* EMAN::EMData::rdata
mutableprivate

image real data

Definition at line 835 of file emdata.h.

Referenced by add_complex_at(), clip_inplace(), do_fft_inplace(), do_ift(), EMData(), free_rdata(), get_attr_dict(), get_complex_at(), operator=(), set_complex_at(), and update_stat().

◆ rot_fp

EMData* EMAN::EMData::rot_fp
mutableprivate

This is a cached rotational footprint, can save much time.

Definition at line 861 of file emdata.h.

Referenced by EMData(), make_rotational_footprint(), make_rotational_footprint_cmc(), make_rotational_footprint_e1(), operator=(), and update_stat().

◆ supp

float* EMAN::EMData::supp
private

supplementary data array

Definition at line 837 of file emdata.h.

Referenced by set_supp_pickle(), and setup4slice().

◆ totalalloc

int EMData::totalalloc =0
static

Definition at line 803 of file emdata.h.

Referenced by EMData(), and ~EMData().

◆ xoff

int EMAN::EMData::xoff
private

array index offsets

Definition at line 851 of file emdata.h.

Referenced by operator=().

◆ yoff

int EMAN::EMData::yoff
private

Definition at line 851 of file emdata.h.

Referenced by operator=().

◆ zoff

int EMAN::EMData::zoff
private

Definition at line 851 of file emdata.h.

Referenced by operator=().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2