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