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
00036 #include "emdata.h"
00037 #include "ctf.h"
00038 #include "emfft.h"
00039 #include "cmp.h"
00040
00041 using namespace EMAN;
00042
00043
00044 #include <iostream>
00045 #include <cstring>
00046
00047 using std::cout;
00048 using std::endl;
00049
00050 #ifdef EMAN2_USING_CUDA
00051 #include "cuda/cuda_processor.h"
00052 #include "cuda/cuda_cmp.h"
00053 #endif // EMAN2_USING_CUDA
00054
00055 void EMData::free_memory()
00056 {
00057 ENTERFUNC;
00058 if (rdata) {
00059 EMUtil::em_free(rdata);
00060 rdata = 0;
00061 }
00062
00063 if (supp) {
00064 EMUtil::em_free(supp);
00065 supp = 0;
00066 }
00067
00068 if (rot_fp != 0)
00069 {
00070 delete rot_fp;
00071 rot_fp = 0;
00072 }
00073
00074
00075
00076
00077
00078
00079
00080 EXITFUNC;
00081 }
00082
00083 void EMData::free_rdata()
00084 {
00085 ENTERFUNC;
00086 if (rdata) {
00087 EMUtil::em_free(rdata);
00088 rdata = 0;
00089 }
00090 EXITFUNC;
00091 }
00092
00093 EMData * EMData::copy() const
00094 {
00095 ENTERFUNC;
00096
00097 EMData *ret = new EMData(*this);
00098
00099 EXITFUNC;
00100 return ret;
00101 }
00102
00103
00104 EMData *EMData::copy_head() const
00105 {
00106 ENTERFUNC;
00107 EMData *ret = new EMData();
00108 ret->attr_dict = attr_dict;
00109
00110 ret->set_size(nx, ny, nz);
00111 ret->flags = flags;
00112
00113 ret->all_translation = all_translation;
00114
00115 ret->path = path;
00116 ret->pathnum = pathnum;
00117
00118
00119
00120
00121
00122
00123
00124 ret->update();
00125
00126 EXITFUNC;
00127 return ret;
00128 }
00129
00130 std::complex<float> EMData::get_complex_at(const int &x,const int &y) const{
00131 if (abs(x)>=nx/2 || abs(y)>ny/2) return std::complex<float>(0,0);
00132 if (x>=0 && y>=0) return std::complex<float>(rdata[ x*2+y*nx], rdata[x*2+y*nx+1]);
00133 if (x>0 && y<0) return std::complex<float>( rdata[ x*2+(ny+y)*nx], rdata[x*2+(ny+y)*nx+1]);
00134 if (x<0 && y>0) return std::complex<float>( rdata[-x*2+(ny-y)*nx],-rdata[-x*2+(ny-y)*nx+1]);
00135 return std::complex<float>(rdata[-x*2-y*nx],-rdata[-x*2+-y*nx+1]);
00136 }
00137
00138 std::complex<float> EMData::get_complex_at(const int &x,const int &y,const int &z) const{
00139 if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return std::complex<float>(0,0);
00140
00141 if (x<0) {
00142 int idx=-x*2+(y<=0?-y:ny-y)*nx+(z<=0?-z:nz-z)*nxy;
00143 return std::complex<float>(rdata[idx],rdata[idx+1]);
00144 }
00145
00146 int idx=x*2+(y<0?ny+y:y)*nx+(z<0?nz+z:z)*nxy;
00147 return std::complex<float>(rdata[idx],rdata[idx+1]);
00148 }
00149
00150 size_t EMData::get_complex_index(const int &x,const int &y,const int &z) const {
00151 if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return nxyz;
00152 if (x<0) {
00153 return -x*2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
00154 }
00155 return x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
00156 }
00157
00158 size_t EMData::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 {
00159 if (abs(x)>=fullnx/2 || abs(y)>fullny/2 || abs(z)>fullnz/2) return nxyz;
00160
00161 if (x<0) {
00162 x*=-1;
00163 y*=-1;
00164 z*=-1;
00165 }
00166 if (y<0) y=fullny+y;
00167 if (z<0) z=fullnz+z;
00168
00169 if (x<subx0||y<suby0||z<subz0||x>=subx0+nx||y>=suby0+ny||z>=subz0+nz) return nxyz;
00170
00171 return (x-subx0)*2+(y-suby0)*(size_t)nx+(z-subz0)*(size_t)nx*(size_t)ny;
00172 }
00173
00174
00175 void EMData::set_complex_at(const int &x,const int &y,const std::complex<float> &val) {
00176 if (abs(x)>=nx/2 || abs(y)>ny/2) return;
00177 if (x>=0 && y>=0) { rdata[ x*2+y*nx]=val.real(); rdata[x*2+y*nx+1]=val.imag(); }
00178 else if (x>0 && y<0) { rdata[ x*2+(ny+y)*nx]=val.real(); rdata[x*2+(ny+y)*nx+1]=val.imag(); }
00179 else if (x<0 && y>0) { rdata[-x*2+(ny-y)*nx]=val.real(); rdata[-x*2+(ny-y)*nx+1]=-val.imag(); }
00180 else { rdata[-x*2-y*nx]=val.real(); rdata[-x*2+-y*nx+1]=-val.imag(); }
00181 return;
00182 }
00183
00184 void EMData::set_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val) {
00185 if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return;
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 size_t idx;
00197 if (x<0) {
00198 idx=-x*2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
00199 rdata[idx]=(float)val.real();
00200 rdata[idx+1]=-(float)val.imag();
00201 return;
00202 }
00203
00204 idx=x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
00205 rdata[idx]=(float)val.real();
00206 rdata[idx+1]=(float)val.imag();
00207
00208 return;
00209 }
00210
00211 size_t EMData::add_complex_at(const int &x,const int &y,const int &z,const std::complex<float> &val) {
00212
00213 if (x>=nx/2 || y>ny/2 || z>nz/2 || x<=-nx/2 || y<-ny/2 || z<-nz/2) return nxyz;
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223 size_t idx;
00224 if (x<0) {
00225 idx=-x*2+(y<=0?-y:ny-y)*(size_t)nx+(z<=0?-z:nz-z)*(size_t)nxy;
00226 rdata[idx]+=(float)val.real();
00227 rdata[idx+1]-=(float)val.imag();
00228 return idx;
00229 }
00230
00231 idx=x*2+(y<0?ny+y:y)*(size_t)nx+(z<0?nz+z:z)*(size_t)nxy;
00232 rdata[idx]+=(float)val.real();
00233 rdata[idx+1]+=(float)val.imag();
00234
00235 return idx;
00236 }
00237
00238 size_t EMData::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) {
00239 if (abs(x)>=fullnx/2 || abs(y)>fullny/2 || abs(z)>fullnz/2) return nxyz;
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 float cc=1.0;
00253 if (x<0) {
00254 x*=-1;
00255 y*=-1;
00256 z*=-1;
00257 cc=-1.0;
00258 }
00259 if (y<0) y=fullny+y;
00260 if (z<0) z=fullnz+z;
00261
00262 if (x<subx0||y<suby0||z<subz0||x>=subx0+nx||y>=suby0+ny||z>=subz0+nz) return nxyz;
00263
00264 size_t idx=(x-subx0)*2+(y-suby0)*(size_t)nx+(z-subz0)*(size_t)nx*ny;
00265 rdata[idx]+=(float)val.real();
00266 rdata[idx+1]+=cc*(float)val.imag();
00267 return idx;
00268 }
00269
00270
00271 void EMData::add(float f,int keepzero)
00272 {
00273 ENTERFUNC;
00274
00275 float * data = get_data();
00276 if( is_real() )
00277 {
00278 if (f != 0) {
00279
00280 size_t size = nxyz;
00281 if (keepzero) {
00282 for (size_t i = 0; i < size; i++) {
00283 if (data[i]) data[i] += f;
00284 }
00285 }
00286 else {
00287 for (size_t i = 0; i < size; i++) {
00288 data[i] += f;
00289 }
00290 }
00291 update();
00292 }
00293 }
00294 else if( is_complex() )
00295 {
00296 if( f!=0 )
00297 {
00298 update();
00299 size_t size = (size_t)nx*ny*nz;
00300 if( keepzero )
00301 {
00302 for(size_t i=0; i<size; i+=2)
00303 {
00304 if (data[i]) data[i] += f;
00305 }
00306 }
00307 else
00308 {
00309 for(size_t i=0; i<size; i+=2)
00310 {
00311 data[i] += f;
00312 }
00313 }
00314 }
00315 }
00316 else
00317 {
00318 throw ImageFormatException("This image is neither a real nor a complex image.");
00319 }
00320 update();
00321 EXITFUNC;
00322 }
00323
00324
00325
00326 void EMData::add(const EMData & image)
00327 {
00328 ENTERFUNC;
00329 if (nx != image.get_xsize() || ny != image.get_ysize() || nz != image.get_zsize()) {
00330 throw ImageFormatException( "images not same sizes");
00331 }
00332 else if( (is_real()^image.is_real()) == true )
00333 {
00334 throw ImageFormatException( "not support add between real image and complex image");
00335 }
00336 else {
00337
00338 const float *src_data = image.get_data();
00339 size_t size = nxyz;
00340 float* data = get_data();
00341
00342 for (size_t i = 0; i < size; i++) {
00343 data[i] += src_data[i];
00344 }
00345 update();
00346 }
00347 EXITFUNC;
00348 }
00349
00350
00351 void EMData::addsquare(const EMData & image)
00352 {
00353 ENTERFUNC;
00354 if (nx != image.get_xsize() || ny != image.get_ysize() || nz != image.get_zsize()) {
00355 throw ImageFormatException( "images not same sizes");
00356 }
00357 else if( this->is_complex() || image.is_complex() )
00358 {
00359 throw ImageFormatException( "Cannot addsquare() with complex images");
00360 }
00361 else {
00362
00363 const float *src_data = image.get_data();
00364 size_t size = nxyz;
00365 float* data = get_data();
00366
00367 for (size_t i = 0; i < size; i++) {
00368 data[i] += src_data[i]*src_data[i];
00369 }
00370 update();
00371 }
00372 EXITFUNC;
00373 }
00374
00375
00376 void EMData::subsquare(const EMData & image)
00377 {
00378 ENTERFUNC;
00379 if (nx != image.get_xsize() || ny != image.get_ysize() || nz != image.get_zsize()) {
00380 throw ImageFormatException( "images not same sizes");
00381 }
00382 else if( this->is_complex() || image.is_complex() )
00383 {
00384 throw ImageFormatException( "Cannot addsquare() with complex images");
00385 }
00386 else {
00387
00388 const float *src_data = image.get_data();
00389 size_t size = nxyz;
00390 float* data = get_data();
00391
00392 for (size_t i = 0; i < size; i++) {
00393 data[i] -= src_data[i]*src_data[i];
00394 }
00395 update();
00396 }
00397 EXITFUNC;
00398 }
00399
00400
00401 void EMData::sub(float f)
00402 {
00403 ENTERFUNC;
00404
00405 #ifdef EMAN2_USING_CUDA
00406 if (EMData::usecuda == 1 && cudarwdata) {
00407 if(f != 0){
00408 subtract_cuda(cudarwdata, f, nx, ny, nz);
00409 }
00410 EXITFUNC;
00411 return;
00412 }
00413 #endif // EMAN2_USING_CUDA
00414
00415 float* data = get_data();
00416 if( is_real() )
00417 {
00418 if (f != 0) {
00419 size_t size = nxyz;
00420 for (size_t i = 0; i < size; i++) {
00421 data[i] -= f;
00422 }
00423 }
00424 update();
00425 }
00426 else if( is_complex() )
00427 {
00428 if( f != 0 )
00429 {
00430 size_t size = nxyz;
00431 for( size_t i=0; i<size; i+=2 )
00432 {
00433 data[i] -= f;
00434 }
00435 }
00436 update();
00437 }
00438 else
00439 {
00440 throw ImageFormatException("This image is neither a real nor a complex image.");
00441 }
00442
00443 EXITFUNC;
00444 }
00445
00446
00447
00448 void EMData::sub(const EMData & em)
00449 {
00450 ENTERFUNC;
00451
00452 if (nx != em.get_xsize() || ny != em.get_ysize() || nz != em.get_zsize()) {
00453 throw ImageFormatException("images not same sizes");
00454 }
00455 else if( (is_real()^em.is_real()) == true )
00456 {
00457 throw ImageFormatException( "not support sub between real image and complex image");
00458 }
00459 else {
00460 const float *src_data = em.get_data();
00461 size_t size = nxyz;
00462 float* data = get_data();
00463
00464 for (size_t i = 0; i < size; i++) {
00465 data[i] -= src_data[i];
00466 }
00467 update();
00468 }
00469 EXITFUNC;
00470 }
00471
00472
00473 void EMData::mult(float f)
00474 {
00475 ENTERFUNC;
00476
00477
00478 if (is_complex()) {
00479 ap2ri();
00480 }
00481 if (f != 1.0) {
00482 #ifdef EMAN2_USING_CUDA
00483 if (EMData::usecuda == 1 && cudarwdata) {
00484
00485 emdata_processor_mult(cudarwdata,f,nx,ny,nz);
00486 EXITFUNC;
00487 return;
00488 }
00489 #endif // EMAN2_USING_CUDA
00490 float* data = get_data();
00491 size_t size = nxyz;
00492 for (size_t i = 0; i < size; i++) {
00493 data[i] *= f;
00494 }
00495 update();
00496 }
00497 EXITFUNC;
00498 }
00499
00500
00501 void EMData::mult(const EMData & em, bool prevent_complex_multiplication)
00502 {
00503 ENTERFUNC;
00504
00505 if (nx != em.get_xsize() || ny != em.get_ysize() || nz != em.get_zsize()) {
00506 throw ImageFormatException( "can not multiply images that are not the same size");
00507 }
00508 else if( (is_real()^em.is_real()) == true )
00509 {
00510 throw ImageFormatException( "can not multiply real and complex images.");
00511 }
00512 else
00513 {
00514 const float *src_data = em.get_data();
00515 size_t size = nxyz;
00516 float* data = get_data();
00517 if( is_real() || prevent_complex_multiplication )
00518 {
00519 for (size_t i = 0; i < size; i++) {
00520 data[i] *= src_data[i];
00521 }
00522 }
00523 else
00524 {
00525 typedef std::complex<float> comp;
00526 for( size_t i = 0; i < size; i+=2 )
00527 {
00528 comp c_src( src_data[i], src_data[i+1] );
00529 comp c_rdat( data[i], data[i+1] );
00530 comp c_result = c_src * c_rdat;
00531 data[i] = c_result.real();
00532 data[i+1] = c_result.imag();
00533 }
00534 }
00535 update();
00536 }
00537
00538 EXITFUNC;
00539 }
00540
00541 void EMData::mult_complex_efficient(const EMData & em, const int radius)
00542 {
00543 ENTERFUNC;
00544
00545 if( is_real() || em.is_real() )throw ImageFormatException( "can call mult_complex_efficient unless both images are complex");
00546
00547
00548 const float *src_data = em.get_data();
00549
00550 size_t i_radius = radius;
00551 size_t k_radius = 1;
00552 size_t j_radius = 1;
00553 int ndim = get_ndim();
00554
00555 if (ndim != em.get_ndim()) throw ImageDimensionException("Can't do that");
00556
00557 if ( ndim == 3 ) {
00558 k_radius = radius;
00559 j_radius = radius;
00560 } else if ( ndim == 2 ) {
00561 j_radius = radius;
00562 }
00563
00564
00565 int s_nx = em.get_xsize();
00566 int s_nxy = s_nx*em.get_ysize();
00567
00568 size_t r_size = nxyz;
00569 int s_size = s_nxy*em.get_zsize();
00570 float* data = get_data();
00571
00572 for (size_t k = 0; k < k_radius; ++k ) {
00573 for (size_t j = 0; j < j_radius; j++) {
00574 for (size_t i = 0; i < i_radius; i++) {
00575 int r_idx = k*nxy + j*nx + i;
00576 int s_idx = k*s_nxy + j*s_nx + i;
00577 data[r_idx] *= src_data[s_idx];
00578 data[r_size-r_idx-1] *= src_data[s_size-s_idx-1];
00579 }
00580 }
00581 }
00582
00583 update();
00584
00585 EXITFUNC;
00586 }
00587
00588
00589 void EMData::div(float f)
00590 {
00591 ENTERFUNC;
00592 if ( f == 0 ) {
00593 throw InvalidValueException(f,"Can not divide by zero");
00594 }
00595 mult(1.0f/f);
00596 EXITFUNC;
00597 }
00598
00599
00600 void EMData::div(const EMData & em)
00601 {
00602 ENTERFUNC;
00603
00604 if (nx != em.get_xsize() || ny != em.get_ysize() || nz != em.get_zsize()) {
00605 throw ImageFormatException( "images not same sizes");
00606 }
00607 else if( (is_real()^em.is_real()) == true )
00608 {
00609 throw ImageFormatException( "not support division between real image and complex image");
00610 }
00611 else {
00612 const float *src_data = em.get_data();
00613 size_t size = nxyz;
00614 float* data = get_data();
00615
00616 if( is_real() )
00617 {
00618 for (size_t i = 0; i < size; i++) {
00619 if(src_data[i] != 0) {
00620 data[i] /= src_data[i];
00621 }
00622 else {
00623 if (data[i]==0) continue;
00624 throw InvalidValueException(src_data[i], "divide by zero");
00625 }
00626 }
00627 }
00628 else
00629 {
00630 typedef std::complex<float> comp;
00631 for( size_t i = 0; i < size; i+=2 )
00632 {
00633 comp c_src( src_data[i], src_data[i+1] );
00634 comp c_rdat( data[i], data[i+1] );
00635 comp c_result = c_rdat / c_src;
00636 data[i] = c_result.real();
00637 data[i+1] = c_result.imag();
00638 }
00639 }
00640 update();
00641 }
00642
00643 EXITFUNC;
00644 }
00645
00646
00647
00648 float EMData::dot(EMData * with)
00649 {
00650 ENTERFUNC;
00651 if (!with) {
00652 throw NullPointerException("Null EMData Image");
00653 }
00654 DotCmp dot_cmp;
00655 float r = -dot_cmp.cmp(this, with);
00656 EXITFUNC;
00657 return r;
00658 }
00659
00660
00661 EMData *EMData::get_row(int row_index) const
00662 {
00663 ENTERFUNC;
00664
00665 if (get_ndim() > 2) {
00666 throw ImageDimensionException("1D/2D image only");
00667 }
00668
00669 EMData *ret = new EMData();
00670 ret->set_size(nx, 1, 1);
00671 memcpy(ret->get_data(), get_data() + nx * row_index, nx * sizeof(float));
00672 ret->update();
00673 EXITFUNC;
00674 return ret;
00675 }
00676
00677
00678 void EMData::set_row(const EMData * d, int row_index)
00679 {
00680 ENTERFUNC;
00681
00682 if (get_ndim() > 2) {
00683 throw ImageDimensionException("1D/2D image only");
00684 }
00685 if (d->get_ndim() != 1) {
00686 throw ImageDimensionException("1D image only");
00687 }
00688
00689 float *dst = get_data();
00690 float *src = d->get_data();
00691 memcpy(dst + nx * row_index, src, nx * sizeof(float));
00692 update();
00693 EXITFUNC;
00694 }
00695
00696 EMData *EMData::get_col(int col_index) const
00697 {
00698 ENTERFUNC;
00699
00700 if (get_ndim() != 2) {
00701 throw ImageDimensionException("2D image only");
00702 }
00703
00704 EMData *ret = new EMData();
00705 ret->set_size(ny, 1, 1);
00706 float *dst = ret->get_data();
00707 float *src = get_data();
00708
00709 for (int i = 0; i < ny; i++) {
00710 dst[i] = src[i * nx + col_index];
00711 }
00712
00713 ret->update();
00714 EXITFUNC;
00715 return ret;
00716 }
00717
00718
00719 void EMData::set_col(const EMData * d, int n)
00720 {
00721 ENTERFUNC;
00722
00723 if (get_ndim() != 2) {
00724 throw ImageDimensionException("2D image only");
00725 }
00726 if (d->get_ndim() != 1) {
00727 throw ImageDimensionException("1D image only");
00728 }
00729
00730 float *dst = get_data();
00731 float *src = d->get_data();
00732
00733 for (int i = 0; i < ny; i++) {
00734 dst[i * nx + n] = src[i];
00735 }
00736
00737 update();
00738 EXITFUNC;
00739 }
00740
00741 float& EMData::get_value_at_wrap(int x)
00742 {
00743 if (x < 0) x = nx + x;
00744 return get_data()[x];
00745 }
00746
00747 float& EMData::get_value_at_wrap(int x, int y)
00748 {
00749 if (x < 0) x = nx + x;
00750 if (y < 0) y = ny + y;
00751
00752 return get_data()[x + y * nx];
00753 }
00754
00755 float& EMData::get_value_at_wrap(int x, int y, int z)
00756 {
00757
00758 #ifdef EMAN2_USING_CUDA
00759 if(EMData::usecuda == 1 && cudarwdata){
00760 float result = get_value_at_wrap_cuda(cudarwdata, x, y, z, nx, ny, nz);
00761 return result;
00762 }
00763 #endif
00764 int lx = x;
00765 int ly = y;
00766 int lz = z;
00767
00768 if (lx < 0) lx = nx + lx;
00769 if (ly < 0) ly = ny + ly;
00770 if (lz < 0) lz = nz + lz;
00771
00772 return get_data()[lx + ly * (size_t)nx + lz * (size_t)nxy];
00773 }
00774
00775 float EMData::get_value_at_wrap(int x) const
00776 {
00777 if (x < 0) x = nx - x;
00778 return get_data()[x];
00779 }
00780
00781 float EMData::get_value_at_wrap(int x, int y) const
00782 {
00783 if (x < 0) x = nx - x;
00784 if (y < 0) y = ny - y;
00785
00786 return get_data()[x + y * nx];
00787 }
00788
00789 float EMData::get_value_at_wrap(int x, int y, int z) const
00790 {
00791 ptrdiff_t lx = x;
00792 ptrdiff_t ly = y;
00793 ptrdiff_t lz = z;
00794 if (lx < 0) lx = nx + lx;
00795 if (ly < 0) ly = ny + ly;
00796 if (lz < 0) lz = nz + lz;
00797
00798 return get_data()[lx + ly * nx + lz * nxy];
00799 }
00800
00801 float EMData::sget_value_at(int x, int y, int z) const
00802 {
00803 if (x < 0 || x >= nx || y < 0 || y >= ny || z < 0 || z >= nz) {
00804 return 0;
00805 }
00806 return get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy];
00807 }
00808
00809
00810 float EMData::sget_value_at(int x, int y) const
00811 {
00812 if (x < 0 || x >= nx || y < 0 || y >= ny) {
00813 return 0;
00814 }
00815 return get_data()[x + y * nx];
00816 }
00817
00818
00819 float EMData::sget_value_at(size_t i) const
00820 {
00821 size_t size = nx*ny;
00822 size *= nz;
00823 if (i >= size) {
00824 return 0;
00825 }
00826 return get_data()[i];
00827 }
00828
00829
00830 float EMData::sget_value_at_interp(float xx, float yy) const
00831 {
00832 int x = static_cast < int >(Util::fast_floor(xx));
00833 int y = static_cast < int >(Util::fast_floor(yy));
00834
00835 float p1 = sget_value_at(x, y);
00836 float p2 = sget_value_at(x + 1, y);
00837 float p3 = sget_value_at(x, y + 1);
00838 float p4 = sget_value_at(x + 1, y + 1);
00839
00840 float result = Util::bilinear_interpolate(p1, p2, p3, p4, xx - x, yy - y);
00841 return result;
00842 }
00843
00844
00845 float EMData::sget_value_at_interp(float xx, float yy, float zz) const
00846 {
00847 int x = (int) Util::fast_floor(xx);
00848 int y = (int) Util::fast_floor(yy);
00849 int z = (int) Util::fast_floor(zz);
00850
00851 float p1 = sget_value_at(x, y, z);
00852 float p2 = sget_value_at(x + 1, y, z);
00853 float p3 = sget_value_at(x, y + 1, z);
00854 float p4 = sget_value_at(x + 1, y + 1, z);
00855
00856 float p5 = sget_value_at(x, y, z + 1);
00857 float p6 = sget_value_at(x + 1, y, z + 1);
00858 float p7 = sget_value_at(x, y + 1, z + 1);
00859 float p8 = sget_value_at(x + 1, y + 1, z + 1);
00860
00861 float result = Util::trilinear_interpolate(p1, p2, p3, p4, p5, p6, p7, p8,
00862 xx - x, yy - y, zz - z);
00863
00864 return result;
00865 }
00866
00867
00868 EMData & EMData::operator+=(float n)
00869 {
00870 add(n);
00871 update();
00872 return *this;
00873 }
00874
00875
00876 EMData & EMData::operator-=(float n)
00877 {
00878 *this += (-n);
00879 return *this;
00880 }
00881
00882
00883 EMData & EMData::operator*=(float n)
00884 {
00885 mult(n);
00886 update();
00887 return *this;
00888 }
00889
00890
00891 EMData & EMData::operator/=(float n)
00892 {
00893 if (n == 0) {
00894 LOGERR("divided by zero");
00895 return *this;
00896 }
00897 *this *= (1.0f / n);
00898 return *this;
00899 }
00900
00901
00902 EMData & EMData::operator+=(const EMData & em)
00903 {
00904 add(em);
00905 update();
00906 return *this;
00907 }
00908
00909
00910 EMData & EMData::operator-=(const EMData & em)
00911 {
00912 sub(em);
00913 update();
00914 return *this;
00915 }
00916
00917
00918 EMData & EMData::operator*=(const EMData & em)
00919 {
00920 mult(em);
00921 update();
00922 return *this;
00923 }
00924
00925
00926 EMData & EMData::operator/=(const EMData & em)
00927 {
00928 div(em);
00929 update();
00930 return *this;
00931 }
00932
00933
00934 EMData * EMData::power(int n) const
00935 {
00936 ENTERFUNC;
00937
00938 if( n<0 ) {
00939 throw InvalidValueException(n, "the power of negative integer not supported.");
00940 }
00941
00942 EMData * r = this->copy();
00943 if( n == 0 ) {
00944 r->to_one();
00945 }
00946 else if( n>1 ) {
00947 for( int i=1; i<n; i++ ) {
00948 *r *= *this;
00949 }
00950 }
00951
00952 r->update();
00953 return r;
00954
00955 EXITFUNC;
00956 }
00957
00958
00959 EMData * EMData::sqrt() const
00960 {
00961 ENTERFUNC;
00962
00963 if (is_complex()) {
00964 throw ImageFormatException("real image only");
00965 }
00966
00967 EMData * r = this->copy();
00968 float * new_data = r->get_data();
00969 float * data = get_data();
00970 size_t size = nxyz;
00971 for (size_t i = 0; i < size; ++i) {
00972 if(data[i] < 0) {
00973 throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
00974 }
00975 else {
00976 if(data[i]) {
00977 new_data[i] = std::sqrt(data[i]);
00978 }
00979 }
00980 }
00981
00982 r->update();
00983 return r;
00984
00985 EXITFUNC;
00986 }
00987
00988
00989 EMData * EMData::log() const
00990 {
00991 ENTERFUNC;
00992
00993 if (is_complex()) {
00994 throw ImageFormatException("real image only");
00995 }
00996
00997 EMData * r = this->copy();
00998 float * new_data = r->get_data();
00999 float * data = get_data();
01000 size_t size = nxyz;
01001 for (size_t i = 0; i < size; ++i) {
01002 if(data[i] < 0) {
01003 throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
01004 }
01005 else {
01006 if(data[i]) {
01007 new_data[i] = std::log(data[i]);
01008 }
01009 }
01010 }
01011
01012 r->update();
01013 return r;
01014
01015 EXITFUNC;
01016 }
01017
01018
01019 EMData * EMData::log10() const
01020 {
01021 ENTERFUNC;
01022
01023 if (is_complex()) {
01024 throw ImageFormatException("real image only");
01025 }
01026
01027 EMData * r = this->copy();
01028 float * new_data = r->get_data();
01029 float * data = get_data();
01030 size_t size = nxyz;
01031 for (size_t i = 0; i < size; ++i) {
01032 if(data[i] < 0) {
01033 throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
01034 }
01035 else {
01036 if(data[i]) {
01037 new_data[i] = std::log10(data[i]);
01038 }
01039 }
01040 }
01041
01042 r->update();
01043 return r;
01044
01045 EXITFUNC;
01046 }
01047
01048
01049 EMData * EMData::real() const
01050 {
01051 ENTERFUNC;
01052
01053 EMData * e = new EMData();
01054
01055 if( is_real() )
01056 {
01057 e = this->copy();
01058 }
01059 else
01060 {
01061 if( !is_ri() )
01062 {
01063 delete e;
01064 throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01065 }
01066 int nx = get_xsize();
01067 int ny = get_ysize();
01068 int nz = get_zsize();
01069 e->set_size(nx/2, ny, nz);
01070 float * edata = e->get_data();
01071 float * data = get_data();
01072 size_t idx1, idx2;
01073 for( int i=0; i<nx; ++i )
01074 {
01075 for( int j=0; j<ny; ++j )
01076 {
01077 for( int k=0; k<nz; ++k )
01078 {
01079 if( i%2 == 0 )
01080 {
01081
01082 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01083 idx2 = i+j*nx+k*nx*ny;
01084 edata[idx1] = data[idx2];
01085 }
01086 }
01087 }
01088 }
01089 }
01090
01091 e->set_complex(false);
01092 if(e->get_ysize()==1 && e->get_zsize()==1) {
01093 e->set_complex_x(false);
01094 }
01095 e->update();
01096 return e;
01097
01098 EXITFUNC;
01099 }
01100
01101
01102 EMData * EMData::imag() const
01103 {
01104 ENTERFUNC;
01105
01106 EMData * e = new EMData();
01107
01108 if( is_real() ) {
01109 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01110 }
01111 else {
01112 if( !is_ri() ) {
01113 throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01114 }
01115 int nx = get_xsize();
01116 int ny = get_ysize();
01117 int nz = get_zsize();
01118 e->set_size(nx/2, ny, nz);
01119 float * edata = e->get_data();
01120 float * data = get_data();
01121 for( int i=0; i<nx; i++ ) {
01122 for( int j=0; j<ny; j++ ) {
01123 for( int k=0; k<nz; k++ ) {
01124 if( i%2 == 1 ) {
01125
01126 edata[i/2+j*(nx/2)+k*(nx/2)*ny] = data[i+j*nx+k*nx*ny];
01127 }
01128 }
01129 }
01130 }
01131 }
01132
01133 e->set_complex(false);
01134 if(e->get_ysize()==1 && e->get_zsize()==1) {
01135 e->set_complex_x(false);
01136 }
01137 e->update();
01138 return e;
01139
01140 EXITFUNC;
01141 }
01142
01143 EMData * EMData::absi() const
01144 {
01145 ENTERFUNC;
01146
01147 EMData * e = new EMData();
01148
01149 if( is_real() )
01150 {
01151 e = this->copy();
01152 int nx = get_xsize();
01153 int ny = get_ysize();
01154 int nz = get_zsize();
01155 float *edata = e->get_data();
01156 float * data = get_data();
01157 size_t idx;
01158 for( int i=0; i<nx; ++i ) {
01159 for( int j=0; j<ny; ++j ) {
01160 for( int k=0; k<nz; ++k ) {
01161 idx = i+j*nx+k*nx*ny;
01162 edata[idx] = std::abs(data[idx]);
01163 }
01164 }
01165 }
01166 }
01167 else
01168 {
01169 if( !is_ri() )
01170 {
01171 throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01172 }
01173 int nx = get_xsize();
01174 int ny = get_ysize();
01175 int nz = get_zsize();
01176 e->set_size(nx/2, ny, nz);
01177 float * edata = e->get_data();
01178 float * data = get_data();
01179 size_t idx1, idx2;
01180 for( int i=0; i<nx; ++i )
01181 {
01182 for( int j=0; j<ny; ++j )
01183 {
01184 for( int k=0; k<nz; ++k )
01185 {
01186 if( i%2 == 0 )
01187 {
01188 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01189 idx2 = i+j*nx+k*nx*ny;
01190
01191 edata[idx1] =
01192 std::sqrt(data[idx2]*data[idx2]+data[idx2+1]*data[idx2+1]);
01193 }
01194 }
01195 }
01196 }
01197 }
01198
01199 e->set_complex(false);
01200 if(e->get_ysize()==1 && e->get_zsize()==1) {
01201 e->set_complex_x(false);
01202 }
01203 e->update();
01204 return e;
01205
01206 EXITFUNC;
01207 }
01208
01209
01210 EMData * EMData::amplitude() const
01211 {
01212 ENTERFUNC;
01213
01214 EMData * e = new EMData();
01215
01216 if( is_real() ) {
01217 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01218 }
01219 else {
01220 if(is_ri()) {
01221 throw InvalidCallException("This image is in real/imaginary format, this function call require a complex image in amplitude/phase format.");
01222 }
01223
01224 int nx = get_xsize();
01225 int ny = get_ysize();
01226 int nz = get_zsize();
01227 e->set_size(nx/2, ny, nz);
01228 float * edata = e->get_data();
01229 float * data = get_data();
01230 size_t idx1, idx2;
01231 for( int i=0; i<nx; ++i )
01232 {
01233 for( int j=0; j<ny; ++j )
01234 {
01235 for( int k=0; k<nz; ++k )
01236 {
01237 if( i%2 == 0 )
01238 {
01239 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01240 idx2 = i+j*nx+k*nx*ny;
01241
01242 edata[idx1] = data[idx2];
01243 }
01244 }
01245 }
01246 }
01247 }
01248
01249 e->set_complex(false);
01250 if(e->get_ysize()==1 && e->get_zsize()==1) {
01251 e->set_complex_x(false);
01252 }
01253 e->update();
01254 return e;
01255
01256 EXITFUNC;
01257 }
01258
01259 EMData * EMData::phase() const
01260 {
01261 ENTERFUNC;
01262
01263 EMData * e = new EMData();
01264
01265 if( is_real() ) {
01266 delete e;
01267 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01268 }
01269 else {
01270 if(is_ri()) {
01271 delete e;
01272 throw InvalidCallException("This image is in real/imaginary format, this function call require a complex image in amplitude/phase format.");
01273 }
01274
01275 int nx = get_xsize();
01276 int ny = get_ysize();
01277 int nz = get_zsize();
01278 e->set_size(nx/2, ny, nz);
01279 float * edata = e->get_data();
01280 float * data = get_data();
01281 size_t idx1, idx2;
01282 for( int i=0; i<nx; ++i ) {
01283 for( int j=0; j<ny; ++j ) {
01284 for( int k=0; k<nz; ++k ) {
01285 if( i%2 == 1 ) {
01286 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01287 idx2 = i+j*nx+k*nx*ny;
01288
01289 edata[idx1] = data[idx2];
01290 }
01291 }
01292 }
01293 }
01294 }
01295
01296 e->set_complex(false);
01297 if(e->get_ysize()==1 && e->get_zsize()==1) {
01298 e->set_complex_x(false);
01299 }
01300 e->update();
01301 return e;
01302
01303 EXITFUNC;
01304 }
01305
01306 EMData * EMData::real2complex(const float img) const
01307 {
01308 ENTERFUNC;
01309
01310 if( is_complex() ) {
01311 throw InvalidCallException("This function call only apply to real image");
01312 }
01313
01314 EMData * e = new EMData();
01315 int nx = get_xsize();
01316 int ny = get_ysize();
01317 int nz = get_zsize();
01318 e->set_size(nx*2, ny, nz);
01319
01320 for( int k=0; k<nz; ++k ) {
01321 for( int j=0; j<ny; ++j ) {
01322 for( int i=0; i<nx; ++i ) {
01323 (*e)(i*2,j,k) = (*this)(i,j,k);
01324 (*e)(i*2+1,j,k) = img;
01325 }
01326 }
01327 }
01328
01329 e->set_complex(true);
01330 if(e->get_ysize()==1 && e->get_zsize()==1) {
01331 e->set_complex_x(true);
01332 }
01333 e->set_ri(true);
01334 e->update();
01335 return e;
01336
01337 EXITFUNC;
01338 }
01339
01340 void EMData::to_zero()
01341 {
01342 ENTERFUNC;
01343
01344 if (is_complex()) {
01345 set_ri(true);
01346 }
01347 else {
01348 set_ri(false);
01349 }
01350
01351
01352 to_value(0.0);
01353 update();
01354 EXITFUNC;
01355 }
01356
01357 void EMData::to_one()
01358 {
01359 ENTERFUNC;
01360
01361 if (is_complex()) {
01362 set_ri(true);
01363 }
01364 else {
01365 set_ri(false);
01366 }
01367 to_value(1.0);
01368
01369 update();
01370 EXITFUNC;
01371 }
01372
01373 void EMData::to_value(const float& value)
01374 {
01375 ENTERFUNC;
01376
01377 #ifdef EMAN2_USING_CUDA
01378 if(EMData::usecuda == 1 && cudarwdata){
01379 to_value_cuda(cudarwdata, value, nx, ny, nz);
01380 return;
01381 }
01382 #endif // EMAN2_USING_CUDA
01383 float* data = get_data();
01384
01385
01386
01387
01388
01389 std::fill(data,data+get_size(),value);
01390
01391 update();
01392 EXITFUNC;
01393 }
01394
01395