EMAN2
emdata_core.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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 // debug only
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         nx = 0;
00075         ny = 0;
00076         nz = 0;
00077         nxy = 0;
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 // should these be here? d.woolford I did not comment them out, merely place them here (commented out) to draw attention
00119 //      ret->xoff = xoff;
00120 //      ret->yoff = yoff;
00121 //      ret->zoff = zoff;
00122 //      ret->changecount = changecount;
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 //if (x==0 && (y!=0 || z!=0)) set_complex_at(0,-y,-z,conj(val));
00187 
00188 // for x=0, we need to insert the value in 2 places
00189 // complex conjugate insertion. Removed due to ambiguity with returned index
00190 /*if (x==0 && (y!=0 || z!=0)) {
00191         size_t idx=(y<=0?-y:ny-y)*nx+(z<=0?-z:nz-z)*nx*ny;
00192         rdata[idx]=(float)val.real();
00193         rdata[idx+1]=(float)-val.imag();
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 //if (abs(x)>=nx/2 || abs(y)>ny/2 || abs(z)>nz/2) return nxyz;
00213 if (x>=nx/2 || y>ny/2 || z>nz/2 || x<=-nx/2 || y<-ny/2 || z<-nz/2) return nxyz;
00214 
00215 // for x=0, we need to insert the value in 2 places
00216 // complex conjugate insertion. Removed due to ambiguity with returned index
00217 /*if (x==0 && (y!=0 || z!=0)) {
00218         size_t idx=(y<=0?-y:ny-y)*nx+(z<=0?-z:nz-z)*nx*ny;
00219         rdata[idx]+=(float)val.real();
00220         rdata[idx+1]-=(float)val.imag();
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 //if (x==0 && (y!=0 || z!=0)) add_complex_at(0,-y,-z,subx0,suby0,subz0,fullnx,fullny,fullnz,conj(val));
00241 // complex conjugate insertion. Removed due to ambiguity with returned index
00242 /*if (x==0&& (y!=0 || z!=0)) {
00243         int yy=y<=0?-y:fullny-y;
00244         int zz=z<=0?-z:fullnz-z;
00245 
00246         if (yy<suby0||zz<subz0||yy>=suby0+ny||zz>=subz0+nz) return nx*ny*nz;
00247 
00248         size_t idx=(yy-suby0)*nx+(zz-subz0)*nx*ny;
00249         rdata[idx]+=(float)val.real();
00250         rdata[idx+1]+=(float)-val.imag();
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; //size of data
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 //for add operation, real and complex image is the same
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 //for add operation, real and complex image is the same
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 //for add operation, real and complex image is the same
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 //for sub operation, real and complex image is the same
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 // this will cause a crash if CUDA is used(no rdata) and a complex map is given.....
00478         if (is_complex()) {
00479                 ap2ri();
00480         }
00481         if (f != 1.0) {
00482 #ifdef EMAN2_USING_CUDA
00483                 if (EMData::usecuda == 1 && cudarwdata) { //doesn't make any sense to use RO, esp on compute devices >= 2.0
00484                         //cout << "CUDA mult" << endl;
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 // just a shortcut for cmp("dot")
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); // this should work....
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]) {   //do nothing with pixel has value zero
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]) {   //do nothing with pixel has value zero
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]) {   //do nothing with pixel has value zero
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 //real part has half of x dimension for a complex image
01050 {
01051         ENTERFUNC;
01052 
01053         EMData * e = new EMData();
01054 
01055         if( is_real() ) // a real image, return a copy of itself
01056         {
01057                 e = this->copy();
01058         }
01059         else //for a complex image
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                                                 //complex data in format [real, complex, real, complex...]
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() ) {       //a real image has no imaginary part, throw exception
01109                 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01110         }
01111         else {  //for complex image
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                                                 //complex data in format [real, complex, real, complex...]
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//abs has half of x dimension for a complex image
01144 {
01145         ENTERFUNC;
01146 
01147         EMData * e = new EMData();
01148 
01149         if( is_real() ) // a real image
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 //for a complex image
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                                                 //complex data in format [real, complex, real, complex...]
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                                                 //complex data in format [amp, phase, amp, phase...]
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                                                 //complex data in format [real, complex, real, complex...]
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         //EMUtil::em_memset(get_data(), 0, nxyz * sizeof(float));
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         //the em_memset has segfault for >8GB image, use std::fill() instead, though may be slower
01386 //      if ( value != 0 ) std::fill(data,data+get_size(),value);
01387 //      else EMUtil::em_memset(data, 0, nxyz * sizeof(float)); // This might be faster, I don't know
01388 
01389         std::fill(data,data+get_size(),value);
01390 
01391         update();
01392         EXITFUNC;
01393 }
01394 
01395