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;
00614 bool equal(const EMData& that) const;
00615
00617 inline float& operator()(const int ix, const int iy, const int iz) const {
00618 ptrdiff_t pos = (size_t)(ix-xoff) + ((iy-yoff) + (size_t)(iz-zoff)*ny)*nx;
00619 #ifdef BOUNDS_CHECKING
00620 if (pos < 0 || pos >= (size_t)nx*ny*nz) {
00621 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00622 }
00623 #endif // BOUNDS_CHECKING
00624 return *(get_data() + pos);
00625 }
00626
00627 inline float& operator()(const int ix, const int iy) const {
00628 ptrdiff_t pos = (ix - xoff) + (iy-yoff)*nx;
00629 #ifdef BOUNDS_CHECKING
00630 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00631 {
00632 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00633 }
00634 #endif // BOUNDS_CHECKING
00635 return *(get_data() + pos);
00636 }
00637
00638
00639 inline float& operator()(const size_t ix) const {
00640 ptrdiff_t pos = ix - xoff;
00641 #ifdef BOUNDS_CHECKING
00642 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00643 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00644 #endif // BOUNDS_CHECKING
00645 return *(get_data() + pos);
00646 }
00647
00648
00650 void set_array_offsets(const int xoff_=0, const int yoff_=0,
00651 const int zoff_=0) {
00652 xoff=xoff_; yoff=yoff_; zoff=zoff_;
00653 }
00654
00655
00656 void set_array_offsets(vector<int> offsets) {
00657 set_array_offsets(offsets[0],offsets[1],offsets[2]);
00658 }
00659
00660
00661 vector<int> get_array_offsets() {
00662 vector<int> offsets;
00663 offsets.push_back(xoff);
00664 offsets.push_back(yoff);
00665 offsets.push_back(zoff);
00666 return offsets;
00667 }
00668
00669
00671 std::complex<float>& cmplx(const int ix, const int iy, const int iz) {
00672 ptrdiff_t pos = 2*(ix-xoff)+((iy-yoff)+(iz-zoff)*ny)*(size_t)nx;
00673 #ifdef BOUNDS_CHECKING
00674 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00675 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00676 #endif // BOUNDS_CHECKING
00677 float* begin = get_data() + pos;
00678 return *(reinterpret_cast<std::complex<float>* >(begin));
00679 }
00680
00681
00682 std::complex<float>& cmplx(const int ix, const int iy) {
00683 ptrdiff_t pos = 2*(ix-xoff)+(iy-yoff)*nx;
00684 #ifdef BOUNDS_CHECKING
00685 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00686 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00687 #endif // BOUNDS_CHECKING
00688 float* begin = get_data() + pos;
00689 return *(reinterpret_cast<std::complex<float>* >(begin));
00690 }
00691
00692
00693 std::complex<float>& cmplx(const int ix) {
00694 ptrdiff_t pos = 2*(ix-xoff);
00695 #ifdef BOUNDS_CHECKING
00696 if (pos < 0 || pos >= (size_t)nx*ny*nz)
00697 throw OutofRangeException(0, (size_t)nx*ny*nz-1, pos, "EMData");
00698 #endif // BOUNDS_CHECKING
00699 float* begin = get_data() + pos;
00700 return *(reinterpret_cast<std::complex<float>* >(begin));
00701 }
00702
00703
00709 EMData * power(int n) const;
00710
00715 EMData * sqrt() const;
00716
00717
00723 EMData * log() const;
00724
00725
00731 EMData * log10() const;
00732
00733
00738 EMData * real() const;
00739
00740
00746 EMData * imag() const;
00747
00753 EMData * absi() const;
00754
00755
00760 EMData * amplitude() const;
00761
00762
00767 EMData * phase() const;
00768
00769
00775 EMData * real2complex(float img = 0.0f) const;
00776
00777 #endif //emdata__core_h__