00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00040 #ifndef emdata__core_h__
00041 #define emdata__core_h__
00042
00043 public:
00047 EMData *copy() const;
00048
00049
00053 EMData *copy_head() const;
00054
00055
00060 void add(float f,int keepzero=0);
00061
00062
00068 void add(const EMData & image);
00069
00075 void addsquare(const EMData & image);
00076
00077
00081 void sub(float f);
00082
00083
00088 void sub(const EMData & image);
00089
00095 void subsquare(const EMData & image);
00096
00097
00101 void mult(int n)
00102 {
00103 mult((float)n);
00104 }
00105
00106
00110 void mult(float f);
00111
00112
00120 void mult(const EMData & image, bool prevent_complex_multiplication=false);
00121
00122 void mult_complex_efficient(const EMData & em, const int radius);
00123
00127 void div(float f);
00128
00129
00135 void div(const EMData & image);
00136
00138 void to_zero();
00139
00140
00142 void to_one();
00143
00145 void to_value(const float& value);
00146
00157 float dot(EMData * with);
00158
00159
00166 EMData *get_row(int row_index) const;
00167
00168
00175 void set_row(const EMData * data, int row_index);
00176
00177
00184 EMData *get_col(int col_index) const;
00185
00186
00193 void set_col(const EMData * data, int col_index);
00194
00195
00204 inline float get_value_at(int x, int y, int z) const
00205 {
00206 return get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy];
00207 }
00208
00212 inline float get_value_at_index(int i)
00213 {
00214 return *(rdata + i);
00215 }
00216
00224 inline float get_value_at(int x, int y) const
00225 {
00226 return get_data()[x + y * nx];
00227 }
00228
00229
00237 inline float get_value_at(size_t i) const
00238 {
00239 return get_data()[i];
00240 }
00241
00253 std::complex<float> get_complex_at(const int &x,const int &y) const;
00254
00267 std::complex<float> get_complex_at(const int &x,const int &y,const int &z) const;
00268
00281 size_t get_complex_index(const int &x,const int &y,const int &z) const;
00282
00283 size_t get_complex_index(int x,int y,int z,const int &subx0,const int &suby0,const int &subz0,const int &fullnx,const int &fullny,const int &fullnz) const;
00284
00285 inline size_t get_complex_index_fast(const int &x,const int &y,const int &z) const {
00286
00287 if (x<0) {
00288 return (size_t)x*-2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
00289 }
00290 return x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
00291 }
00292
00305 void set_complex_at(const int &x,const int &y,const std::complex<float> &val);
00306
00320 void set_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val);
00321
00335 size_t add_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val);
00336
00337 inline size_t add_complex_at_fast(const int &x,const int &y,const int &z,const std::complex<float> &val) {
00338
00339
00340 size_t idx;
00341 if (x<0) {
00342 idx=-x*2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
00343 rdata[idx]+=(float)val.real();
00344 rdata[idx+1]-=(float)val.imag();
00345 return idx;
00346 }
00347
00348 idx=x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
00349 rdata[idx]+=(float)val.real();
00350 rdata[idx+1]+=(float)val.imag();
00351
00352 return idx;
00353 }
00354
00366 size_t add_complex_at(int x,int y,int z,const int &subx0,const int &suby0,const int &subz0,const int &fullnx,const int &fullny,const int &fullnz,const std::complex<float> &val);
00367
00377 float get_value_at_wrap(int x, int y, int z) const;
00378 float& get_value_at_wrap(int x, int y, int z);
00379
00388 float get_value_at_wrap(int x, int y) const;
00389 float& get_value_at_wrap(int x, int y);
00390
00398 float get_value_at_wrap(int x) const;
00399 float& get_value_at_wrap(int x);
00400
00410 float sget_value_at(int x, int y, int z) const;
00411
00412
00421 float sget_value_at(int x, int y) const;
00422
00423
00432 float sget_value_at(size_t i) const;
00433
00434
00442 float sget_value_at_interp(float x, float y) const;
00443
00444
00453 float sget_value_at_interp(float x, float y, float z) const;
00454
00455
00465 inline void set_value_at(int x, int y, int z, float v)
00466 {
00467 if( x>=nx || x<0 )
00468 {
00469 throw OutofRangeException(0, nx-1, x, "x dimension index");
00470 }
00471 else if( y>=ny || y<0 )
00472 {
00473 throw OutofRangeException(0, ny-1, y, "y dimension index");
00474 }
00475 else if( z>=nz || z<0 )
00476 {
00477 throw OutofRangeException(0, nz-1, z, "z dimension index");
00478 }
00479 else
00480 {
00481 get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy] = v;
00482 flags |= EMDATA_NEEDUPD;
00483 changecount++;
00484 }
00485 }
00486
00487
00497 inline void set_value_at_fast(int x, int y, int z, float v)
00498 {
00499 get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy] = v;
00500 flags |= EMDATA_NEEDUPD;
00501 changecount++;
00502 }
00503
00510 inline void set_value_at_index(int i, float v)
00511 {
00512 *(rdata + i) = v;
00513 }
00514
00523 inline void set_value_at(int x, int y, float v)
00524 {
00525 if( x>=nx || x<0 )
00526 {
00527 throw OutofRangeException(0, nx-1, x, "x dimension index");
00528 }
00529 else if( y>=ny || y<0 )
00530 {
00531 throw OutofRangeException(0, ny-1, y, "y dimension index");
00532 }
00533 else
00534 {
00535 get_data()[x + y * nx] = v;
00536 flags |= EMDATA_NEEDUPD;
00537 changecount++;
00538 }
00539 }
00540
00541
00549 inline void set_value_at_fast(int x, int y, float v)
00550 {
00551 get_data()[x + y * nx] = v;
00552 flags |= EMDATA_NEEDUPD;
00553 changecount++;
00554 }
00555
00556
00564 inline void set_value_at(int x, float v)
00565 {
00566 if( x>=nx || x<0 )
00567 {
00568 throw OutofRangeException(0, nx-1, x, "x dimension index");
00569 }
00570 else
00571 {
00572 get_data()[x] = v;
00573 flags |= EMDATA_NEEDUPD;
00574 changecount++;
00575 }
00576 }
00577
00584 inline void set_value_at_fast(int x, float v)
00585 {
00586 get_data()[x] = v;
00587 flags |= EMDATA_NEEDUPD;
00588 changecount++;
00589 }
00590
00591
00595 void free_memory();
00596
00600 void free_rdata();
00601
00602 EMData & operator+=(float n);
00603 EMData & operator-=(float n);
00604 EMData & operator*=(float n);
00605 EMData & operator/=(float n);
00606
00607 EMData & operator+=(const EMData & em);
00608 EMData & operator-=(const EMData & em);
00609 EMData & operator*=(const EMData & em);
00610 EMData & operator/=(const EMData & em);
00611
00612 bool operator==(const EMData& that) const;
00613
00615 inline float& operator()(const int ix, const int iy, const int iz) const {
00616 ptrdiff_t pos = (size_t)(ix-xoff) + ((iy-yoff) + (size_t)(iz-zoff)*ny)*nx;
00617 #ifdef BOUNDS_CHECKING
00618 if (pos < 0 || pos >= (size_t)nx*ny*nz) {
00619 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00620 }
00621 #endif // BOUNDS_CHECKING
00622 return *(get_data() + pos);
00623 }
00624
00625 inline float& operator()(const int ix, const int iy) const {
00626 ptrdiff_t pos = (ix - xoff) + (iy-yoff)*nx;
00627 #ifdef BOUNDS_CHECKING
00628 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00629 {
00630 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00631 }
00632 #endif // BOUNDS_CHECKING
00633 return *(get_data() + pos);
00634 }
00635
00636
00637 inline float& operator()(const size_t ix) const {
00638 ptrdiff_t pos = ix - xoff;
00639 #ifdef BOUNDS_CHECKING
00640 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00641 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00642 #endif // BOUNDS_CHECKING
00643 return *(get_data() + pos);
00644 }
00645
00646
00648 void set_array_offsets(const int xoff_=0, const int yoff_=0,
00649 const int zoff_=0) {
00650 xoff=xoff_; yoff=yoff_; zoff=zoff_;
00651 }
00652
00653
00654 void set_array_offsets(vector<int> offsets) {
00655 set_array_offsets(offsets[0],offsets[1],offsets[2]);
00656 }
00657
00658
00659 vector<int> get_array_offsets() {
00660 vector<int> offsets;
00661 offsets.push_back(xoff);
00662 offsets.push_back(yoff);
00663 offsets.push_back(zoff);
00664 return offsets;
00665 }
00666
00667
00669 std::complex<float>& cmplx(const int ix, const int iy, const int iz) {
00670 ptrdiff_t pos = 2*(ix-xoff)+((iy-yoff)+(iz-zoff)*ny)*(size_t)nx;
00671 #ifdef BOUNDS_CHECKING
00672 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00673 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00674 #endif // BOUNDS_CHECKING
00675 float* begin = get_data() + pos;
00676 return *(reinterpret_cast<std::complex<float>* >(begin));
00677 }
00678
00679
00680 std::complex<float>& cmplx(const int ix, const int iy) {
00681 ptrdiff_t pos = 2*(ix-xoff)+(iy-yoff)*nx;
00682 #ifdef BOUNDS_CHECKING
00683 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00684 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00685 #endif // BOUNDS_CHECKING
00686 float* begin = get_data() + pos;
00687 return *(reinterpret_cast<std::complex<float>* >(begin));
00688 }
00689
00690
00691 std::complex<float>& cmplx(const int ix) {
00692 ptrdiff_t pos = 2*(ix-xoff);
00693 #ifdef BOUNDS_CHECKING
00694 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00695 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00696 #endif // BOUNDS_CHECKING
00697 float* begin = get_data() + pos;
00698 return *(reinterpret_cast<std::complex<float>* >(begin));
00699 }
00700
00701
00707 EMData * power(int n) const;
00708
00713 EMData * sqrt() const;
00714
00715
00721 EMData * log() const;
00722
00723
00729 EMData * log10() const;
00730
00731
00736 EMData * real() const;
00737
00738
00744 EMData * imag() const;
00745
00751 EMData * absi() const;
00752
00753
00758 EMData * amplitude() const;
00759
00760
00765 EMData * phase() const;
00766
00767
00773 EMData * real2complex(float img = 0.0f) const;
00774
00775 #endif //emdata__core_h__