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

Generated on Mon Mar 7 18:15:25 2011 for EMAN2 by  doxygen 1.3.9.1