Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 
00281 //#ifdef EMAN2_USING_CUDA
00282 //                      if ( gpu_operation_preferred () && !keepzero ) {
00283 //                              EMDataForCuda tmp = get_data_struct_for_cuda();
00284 //                              emdata_processor_add(&tmp,f);
00285 //                              gpu_update();
00286 //                              EXITFUNC;
00287 //                              return;
00288 //                      }
00289 //#endif // EMAN2_USING_CUDA
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; //size of data
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 //for add operation, real and complex image is the same
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 //for add operation, real and complex image is the same
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 //for add operation, real and complex image is the same
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 (EMData::usecuda == 1 && 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 //for sub operation, real and complex image is the same
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 // this will cause a crash if CUDA is used(no rdata) and a complex map is given.....
00488         if (is_complex()) {
00489                 ap2ri();
00490         }
00491         if (f != 1.0) {
00492 #ifdef EMAN2_USING_CUDA
00493                 if (EMData::usecuda == 1 && cudarwdata) { //doesn't make any sense to use RO, esp on compute devices >= 2.0
00494                         //cout << "CUDA mult" << endl;
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 // just a shortcut for cmp("dot")
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(EMData::usecuda == 1 && cudarwdata){
00770                 float result = get_value_at_wrap_cuda(cudarwdata, x, y, z, nx, ny, nz); // this should work....
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 float EMData::get_value_at_wrap(int x) const
00786 {
00787         if (x < 0) x = nx - x;
00788         return get_data()[x];
00789 }
00790 
00791 float EMData::get_value_at_wrap(int x, int y) const
00792 {
00793         if (x < 0) x = nx - x;
00794         if (y < 0) y = ny - y;
00795 
00796         return get_data()[x + y * nx];
00797 }
00798 
00799 float EMData::get_value_at_wrap(int x, int y, int z) const
00800 {
00801         ptrdiff_t lx = x;
00802         ptrdiff_t ly = y;
00803         ptrdiff_t lz = z;
00804         if (lx < 0) lx = nx + lx;
00805         if (ly < 0) ly = ny + ly;
00806         if (lz < 0) lz = nz + lz;
00807 
00808         return get_data()[lx + ly * nx + lz * nxy];
00809 }
00810 
00811 float EMData::sget_value_at(int x, int y, int z) const
00812 {
00813         if (x < 0 || x >= nx || y < 0 || y >= ny || z < 0 || z >= nz) {
00814                 return 0;
00815         }
00816         return get_data()[(size_t)x + (size_t)y * (size_t)nx + (size_t)z * (size_t)nxy];
00817 }
00818 
00819 
00820 float EMData::sget_value_at(int x, int y) const
00821 {
00822         if (x < 0 || x >= nx || y < 0 || y >= ny) {
00823                 return 0;
00824         }
00825         return get_data()[x + y * nx];
00826 }
00827 
00828 
00829 float EMData::sget_value_at(size_t i) const
00830 {
00831         size_t size = nx*ny;
00832         size *= nz;
00833         if (i >= size) {
00834                 return 0;
00835         }
00836         return get_data()[i];
00837 }
00838 
00839 
00840 float EMData::sget_value_at_interp(float xx, float yy) const
00841 {
00842         int x = static_cast < int >(Util::fast_floor(xx));
00843         int y = static_cast < int >(Util::fast_floor(yy));
00844 
00845         float p1 = sget_value_at(x, y);
00846         float p2 = sget_value_at(x + 1, y);
00847         float p3 = sget_value_at(x, y + 1);
00848         float p4 = sget_value_at(x + 1, y + 1);
00849 
00850         float result = Util::bilinear_interpolate(p1, p2, p3, p4, xx - x, yy - y);
00851         return result;
00852 }
00853 
00854 
00855 float EMData::sget_value_at_interp(float xx, float yy, float zz) const
00856 {
00857         int x = (int) Util::fast_floor(xx);
00858         int y = (int) Util::fast_floor(yy);
00859         int z = (int) Util::fast_floor(zz);
00860 
00861         float p1 = sget_value_at(x, y, z);
00862         float p2 = sget_value_at(x + 1, y, z);
00863         float p3 = sget_value_at(x, y + 1, z);
00864         float p4 = sget_value_at(x + 1, y + 1, z);
00865 
00866         float p5 = sget_value_at(x, y, z + 1);
00867         float p6 = sget_value_at(x + 1, y, z + 1);
00868         float p7 = sget_value_at(x, y + 1, z + 1);
00869         float p8 = sget_value_at(x + 1, y + 1, z + 1);
00870 
00871         float result = Util::trilinear_interpolate(p1, p2, p3, p4, p5, p6, p7, p8,
00872                                                                                            xx - x, yy - y, zz - z);
00873 
00874         return result;
00875 }
00876 
00877 
00878 EMData & EMData::operator+=(float n)
00879 {
00880         add(n);
00881         update();
00882         return *this;
00883 }
00884 
00885 
00886 EMData & EMData::operator-=(float n)
00887 {
00888         *this += (-n);
00889         return *this;
00890 }
00891 
00892 
00893 EMData & EMData::operator*=(float n)
00894 {
00895         mult(n);
00896         update();
00897         return *this;
00898 }
00899 
00900 
00901 EMData & EMData::operator/=(float n)
00902 {
00903         if (n == 0) {
00904                 LOGERR("divided by zero");
00905                 return *this;
00906         }
00907         *this *= (1.0f / n);
00908         return *this;
00909 }
00910 
00911 
00912 EMData & EMData::operator+=(const EMData & em)
00913 {
00914         add(em);
00915         update();
00916         return *this;
00917 }
00918 
00919 
00920 EMData & EMData::operator-=(const EMData & em)
00921 {
00922         sub(em);
00923         update();
00924         return *this;
00925 }
00926 
00927 
00928 EMData & EMData::operator*=(const EMData & em)
00929 {
00930         mult(em);
00931         update();
00932         return *this;
00933 }
00934 
00935 
00936 EMData & EMData::operator/=(const EMData & em)
00937 {
00938         div(em);
00939         update();
00940         return *this;
00941 }
00942 
00943 
00944 EMData * EMData::power(int n) const
00945 {
00946         ENTERFUNC;
00947 
00948         if( n<0 ) {
00949                 throw InvalidValueException(n, "the power of negative integer not supported.");
00950         }
00951 
00952         EMData * r = this->copy();
00953         if( n == 0 ) {
00954                 r->to_one();
00955         }
00956         else if( n>1 ) {
00957                 for( int i=1; i<n; i++ ) {
00958                         *r *= *this;
00959                 }
00960         }
00961 
00962         r->update();
00963         return r;
00964 
00965         EXITFUNC;
00966 }
00967 
00968 
00969 EMData * EMData::sqrt() const
00970 {
00971         ENTERFUNC;
00972 
00973         if (is_complex()) {
00974                 throw ImageFormatException("real image only");
00975         }
00976 
00977         EMData * r = this->copy();
00978         float * new_data = r->get_data();
00979         float * data = get_data();
00980         size_t size = nxyz;
00981         for (size_t i = 0; i < size; ++i) {
00982                 if(data[i] < 0) {
00983                         throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
00984                 }
00985                 else {
00986                         if(data[i]) {   //do nothing with pixel has value zero
00987                                 new_data[i] = std::sqrt(data[i]);
00988                         }
00989                 }
00990         }
00991 
00992         r->update();
00993         return r;
00994 
00995         EXITFUNC;
00996 }
00997 
00998 
00999 EMData * EMData::log() const
01000 {
01001         ENTERFUNC;
01002 
01003         if (is_complex()) {
01004                 throw ImageFormatException("real image only");
01005         }
01006 
01007         EMData * r = this->copy();
01008         float * new_data = r->get_data();
01009         float * data = get_data();
01010         size_t size = nxyz;
01011         for (size_t i = 0; i < size; ++i) {
01012                 if(data[i] < 0) {
01013                         throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
01014                 }
01015                 else {
01016                         if(data[i]) {   //do nothing with pixel has value zero
01017                                 new_data[i] = std::log(data[i]);
01018                         }
01019                 }
01020         }
01021 
01022         r->update();
01023         return r;
01024 
01025         EXITFUNC;
01026 }
01027 
01028 
01029 EMData * EMData::log10() const
01030 {
01031         ENTERFUNC;
01032 
01033         if (is_complex()) {
01034                 throw ImageFormatException("real image only");
01035         }
01036 
01037         EMData * r = this->copy();
01038         float * new_data = r->get_data();
01039         float * data = get_data();
01040         size_t size = nxyz;
01041         for (size_t i = 0; i < size; ++i) {
01042                 if(data[i] < 0) {
01043                         throw InvalidValueException(data[i], "pixel value must be non-negative for logrithm");
01044                 }
01045                 else {
01046                         if(data[i]) {   //do nothing with pixel has value zero
01047                                 new_data[i] = std::log10(data[i]);
01048                         }
01049                 }
01050         }
01051 
01052         r->update();
01053         return r;
01054 
01055         EXITFUNC;
01056 }
01057 
01058 
01059 EMData * EMData::real() const //real part has half of x dimension for a complex image
01060 {
01061         ENTERFUNC;
01062 
01063         EMData * e = new EMData();
01064 
01065         if( is_real() ) // a real image, return a copy of itself
01066         {
01067                 e = this->copy();
01068         }
01069         else //for a complex image
01070         {
01071                 if( !is_ri() )
01072                 {
01073                         delete e;
01074                         throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01075                 }
01076                 int nx = get_xsize();
01077                 int ny = get_ysize();
01078                 int nz = get_zsize();
01079                 e->set_size(nx/2, ny, nz);
01080                 float * edata = e->get_data();
01081                 float * data = get_data();
01082                 size_t idx1, idx2;
01083                 for( int i=0; i<nx; ++i )
01084                 {
01085                         for( int j=0; j<ny; ++j )
01086                         {
01087                                 for( int k=0; k<nz; ++k )
01088                                 {
01089                                         if( i%2 == 0 )
01090                                         {
01091                                                 //complex data in format [real, complex, real, complex...]
01092                                                 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01093                                                 idx2 = i+j*nx+k*nx*ny;
01094                                                 edata[idx1] = data[idx2];
01095                                         }
01096                                 }
01097                         }
01098                 }
01099         }
01100 
01101         e->set_complex(false);
01102         if(e->get_ysize()==1 && e->get_zsize()==1) {
01103                 e->set_complex_x(false);
01104         }
01105         e->update();
01106         return e;
01107 
01108         EXITFUNC;
01109 }
01110 
01111 
01112 EMData * EMData::imag() const
01113 {
01114         ENTERFUNC;
01115 
01116         EMData * e = new EMData();
01117 
01118         if( is_real() ) {       //a real image has no imaginary part, throw exception
01119                 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01120         }
01121         else {  //for complex image
01122                 if( !is_ri() ) {
01123                         throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01124                 }
01125                 int nx = get_xsize();
01126                 int ny = get_ysize();
01127                 int nz = get_zsize();
01128                 e->set_size(nx/2, ny, nz);
01129                 float * edata = e->get_data();
01130                 float * data = get_data();
01131                 for( int i=0; i<nx; i++ ) {
01132                         for( int j=0; j<ny; j++ ) {
01133                                 for( int k=0; k<nz; k++ ) {
01134                                         if( i%2 == 1 ) {
01135                                                 //complex data in format [real, complex, real, complex...]
01136                                                 edata[i/2+j*(nx/2)+k*(nx/2)*ny] = data[i+j*nx+k*nx*ny];
01137                                         }
01138                                 }
01139                         }
01140                 }
01141         }
01142 
01143         e->set_complex(false);
01144         if(e->get_ysize()==1 && e->get_zsize()==1) {
01145                 e->set_complex_x(false);
01146         }
01147         e->update();
01148         return e;
01149 
01150         EXITFUNC;
01151 }
01152 
01153 EMData * EMData::absi() const//abs has half of x dimension for a complex image
01154 {
01155         ENTERFUNC;
01156 
01157         EMData * e = new EMData();
01158 
01159         if( is_real() ) // a real image
01160         {
01161                 e = this->copy();
01162                 int nx = get_xsize();
01163                 int ny = get_ysize();
01164                 int nz = get_zsize();
01165                 float *edata = e->get_data();
01166                 float * data = get_data();
01167                 size_t idx;
01168                 for( int i=0; i<nx; ++i ) {
01169                         for( int j=0; j<ny; ++j ) {
01170                                 for( int k=0; k<nz; ++k ) {
01171                                         idx = i+j*nx+k*nx*ny;
01172                                         edata[idx] = std::abs(data[idx]);
01173                                 }
01174                         }
01175                 }
01176         }
01177         else //for a complex image
01178         {
01179                 if( !is_ri() )
01180                 {
01181                         throw InvalidCallException("This image is in amplitude/phase format, this function call require a complex image in real/imaginary format.");
01182                 }
01183                 int nx = get_xsize();
01184                 int ny = get_ysize();
01185                 int nz = get_zsize();
01186                 e->set_size(nx/2, ny, nz);
01187                 float * edata = e->get_data();
01188                 float * data = get_data();
01189                 size_t idx1, idx2;
01190                 for( int i=0; i<nx; ++i )
01191                 {
01192                         for( int j=0; j<ny; ++j )
01193                         {
01194                                 for( int k=0; k<nz; ++k )
01195                                 {
01196                                         if( i%2 == 0 )
01197                                         {
01198                                                 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01199                                                 idx2 = i+j*nx+k*nx*ny;
01200                                                 //complex data in format [real, complex, real, complex...]
01201                                                 edata[idx1] =
01202                                                 std::sqrt(data[idx2]*data[idx2]+data[idx2+1]*data[idx2+1]);
01203                                         }
01204                                 }
01205                         }
01206                 }
01207         }
01208 
01209         e->set_complex(false);
01210         if(e->get_ysize()==1 && e->get_zsize()==1) {
01211                 e->set_complex_x(false);
01212         }
01213         e->update();
01214         return e;
01215 
01216         EXITFUNC;
01217 }
01218 
01219 
01220 EMData * EMData::amplitude() const
01221 {
01222         ENTERFUNC;
01223 
01224         EMData * e = new EMData();
01225 
01226         if( is_real() ) {
01227                 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01228         }
01229         else {
01230                 if(is_ri()) {
01231                         throw InvalidCallException("This image is in real/imaginary format, this function call require a complex image in amplitude/phase format.");
01232                 }
01233 
01234                 int nx = get_xsize();
01235                 int ny = get_ysize();
01236                 int nz = get_zsize();
01237                 e->set_size(nx/2, ny, nz);
01238                 float * edata = e->get_data();
01239                 float * data = get_data();
01240                 size_t idx1, idx2;
01241                 for( int i=0; i<nx; ++i )
01242                 {
01243                         for( int j=0; j<ny; ++j )
01244                         {
01245                                 for( int k=0; k<nz; ++k )
01246                                 {
01247                                         if( i%2 == 0 )
01248                                         {
01249                                                 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01250                                                 idx2 = i+j*nx+k*nx*ny;
01251                                                 //complex data in format [amp, phase, amp, phase...]
01252                                                 edata[idx1] = data[idx2];
01253                                         }
01254                                 }
01255                         }
01256                 }
01257         }
01258 
01259         e->set_complex(false);
01260         if(e->get_ysize()==1 && e->get_zsize()==1) {
01261                 e->set_complex_x(false);
01262         }
01263         e->update();
01264         return e;
01265 
01266         EXITFUNC;
01267 }
01268 
01269 EMData * EMData::phase() const
01270 {
01271         ENTERFUNC;
01272 
01273         EMData * e = new EMData();
01274 
01275         if( is_real() ) {
01276                 delete e;
01277                 throw InvalidCallException("No imaginary part for a real image, this function call require a complex image.");
01278         }
01279         else {
01280                 if(is_ri()) {
01281                         delete e;
01282                         throw InvalidCallException("This image is in real/imaginary format, this function call require a complex image in amplitude/phase format.");
01283                 }
01284 
01285                 int nx = get_xsize();
01286                 int ny = get_ysize();
01287                 int nz = get_zsize();
01288                 e->set_size(nx/2, ny, nz);
01289                 float * edata = e->get_data();
01290                 float * data = get_data();
01291                 size_t idx1, idx2;
01292                 for( int i=0; i<nx; ++i ) {
01293                         for( int j=0; j<ny; ++j ) {
01294                                 for( int k=0; k<nz; ++k ) {
01295                                         if( i%2 == 1 ) {
01296                                                 idx1 = i/2+j*(nx/2)+k*(nx/2)*ny;
01297                                                 idx2 = i+j*nx+k*nx*ny;
01298                                                 //complex data in format [real, complex, real, complex...]
01299                                                 edata[idx1] = data[idx2];
01300                                         }
01301                                 }
01302                         }
01303                 }
01304         }
01305 
01306         e->set_complex(false);
01307         if(e->get_ysize()==1 && e->get_zsize()==1) {
01308                 e->set_complex_x(false);
01309         }
01310         e->update();
01311         return e;
01312 
01313         EXITFUNC;
01314 }
01315 
01316 EMData * EMData::real2complex(const float img) const
01317 {
01318         ENTERFUNC;
01319 
01320         if( is_complex() ) {
01321                 throw InvalidCallException("This function call only apply to real image");
01322         }
01323 
01324         EMData * e = new EMData();
01325         int nx = get_xsize();
01326         int ny = get_ysize();
01327         int nz = get_zsize();
01328         e->set_size(nx*2, ny, nz);
01329 
01330         for( int k=0; k<nz; ++k ) {
01331                 for( int j=0; j<ny; ++j ) {
01332                         for( int i=0; i<nx; ++i ) {
01333                                 (*e)(i*2,j,k) = (*this)(i,j,k);
01334                                 (*e)(i*2+1,j,k) = img;
01335                         }
01336                 }
01337         }
01338 
01339         e->set_complex(true);
01340         if(e->get_ysize()==1 && e->get_zsize()==1) {
01341                 e->set_complex_x(true);
01342         }
01343         e->set_ri(true);
01344         e->update();
01345         return e;
01346 
01347         EXITFUNC;
01348 }
01349 
01350 void EMData::to_zero()
01351 {
01352         ENTERFUNC;
01353 
01354         if (is_complex()) {
01355                 set_ri(true);
01356         }
01357         else {
01358                 set_ri(false);
01359         }
01360 
01361         //EMUtil::em_memset(get_data(), 0, nxyz * sizeof(float));
01362         to_value(0.0);
01363         update();
01364         EXITFUNC;
01365 }
01366 
01367 void EMData::to_one()
01368 {
01369         ENTERFUNC;
01370 
01371         if (is_complex()) {
01372                 set_ri(true);
01373         }
01374         else {
01375                 set_ri(false);
01376         }
01377         to_value(1.0);
01378 
01379         update();
01380         EXITFUNC;
01381 }
01382 
01383 void EMData::to_value(const float& value)
01384 {
01385         ENTERFUNC;
01386 
01387 #ifdef EMAN2_USING_CUDA
01388         if(EMData::usecuda == 1 && cudarwdata){
01389                 to_value_cuda(cudarwdata, value, nx, ny, nz);
01390                 return;
01391         }
01392 #endif // EMAN2_USING_CUDA
01393         float* data = get_data();
01394 
01395         //the em_memset has segfault for >8GB image, use std::fill() instead, though may be slower
01396 //      if ( value != 0 ) std::fill(data,data+get_size(),value);
01397 //      else EMUtil::em_memset(data, 0, nxyz * sizeof(float)); // This might be faster, I don't know
01398 
01399         std::fill(data,data+get_size(),value);
01400 
01401         update();
01402         EXITFUNC;
01403 }
01404 
01405 

Generated on Thu Nov 17 12:43:44 2011 for EMAN2 by  doxygen 1.3.9.1