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

emdata.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 "all_imageio.h"
00038 #include "ctf.h"
00039 #include "processor.h"
00040 #include "aligner.h"
00041 #include "cmp.h"
00042 #include "emfft.h"
00043 #include "projector.h"
00044 #include "geometry.h"
00045 
00046 #include <gsl/gsl_sf_bessel.h>
00047 #include <gsl/gsl_errno.h>
00048 
00049 #include <iomanip>
00050 #include <complex>
00051 
00052 #include <algorithm> // fill
00053 #include <cmath>
00054 
00055 #ifdef WIN32
00056         #define M_PI 3.14159265358979323846f
00057 #endif  //WIN32
00058 
00059 #define EMDATA_EMAN2_DEBUG 0
00060 
00061 #ifdef EMAN2_USING_CUDA
00062 //#include <driver_functions.h>
00063 #include "cuda/cuda_processor.h"
00064 #include "cuda/cuda_emfft.h"
00065 #endif // EMAN2_USING_CUDA
00066 
00067 using namespace EMAN;
00068 using namespace std;
00069 using namespace boost;
00070 
00071 int EMData::totalalloc=0;               // mainly used for debugging/memory leak purposes
00072 
00073 EMData::EMData() :
00074 #ifdef EMAN2_USING_CUDA
00075                 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
00076 #endif //EMAN2_USING_CUDA
00077                 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0),
00078                 zoff(0), all_translation(),     path(""), pathnum(0), rot_fp(0)
00079 
00080 {
00081         ENTERFUNC;
00082 
00083         attr_dict["apix_x"] = 1.0f;
00084         attr_dict["apix_y"] = 1.0f;
00085         attr_dict["apix_z"] = 1.0f;
00086 
00087         attr_dict["is_complex"] = int(0);
00088         attr_dict["is_complex_x"] = int(0);
00089         attr_dict["is_complex_ri"] = int(1);
00090 
00091         attr_dict["datatype"] = (int)EMUtil::EM_FLOAT;
00092 
00093         EMData::totalalloc++;
00094 #ifdef MEMDEBUG2
00095         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00096 #endif
00097 
00098         EXITFUNC;
00099 }
00100 
00101 EMData::EMData(const string& filename, int image_index) :
00102 #ifdef EMAN2_USING_CUDA
00103                 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
00104 #endif //EMAN2_USING_CUDA
00105                 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0), zoff(0),
00106                 all_translation(),      path(filename), pathnum(image_index), rot_fp(0)
00107 {
00108         ENTERFUNC;
00109 
00110         attr_dict["apix_x"] = 1.0f;
00111         attr_dict["apix_y"] = 1.0f;
00112         attr_dict["apix_z"] = 1.0f;
00113 
00114         attr_dict["is_complex"] = int(0);
00115         attr_dict["is_complex_x"] = int(0);
00116         attr_dict["is_complex_ri"] = int(1);
00117 
00118         attr_dict["datatype"] = (int)EMUtil::EM_FLOAT;
00119 
00120         this->read_image(filename, image_index);
00121 
00122         update();
00123         EMData::totalalloc++;
00124 #ifdef MEMDEBUG2
00125         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00126 #endif
00127 
00128         EXITFUNC;
00129 }
00130 
00131 EMData::EMData(const EMData& that) :
00132 #ifdef EMAN2_USING_CUDA
00133                 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
00134 #endif //EMAN2_USING_CUDA
00135                 attr_dict(that.attr_dict), rdata(0), supp(0), flags(that.flags), changecount(that.changecount), nx(that.nx), ny(that.ny), nz(that.nz),
00136                 nxy(that.nx*that.ny), nxyz((size_t)that.nx*that.ny*that.nz), xoff(that.xoff), yoff(that.yoff), zoff(that.zoff),all_translation(that.all_translation),   path(that.path),
00137                 pathnum(that.pathnum), rot_fp(0)
00138 {
00139         ENTERFUNC;
00140         
00141         float* data = that.rdata;
00142         size_t num_bytes = (size_t)nx*ny*nz*sizeof(float);
00143         if (data && num_bytes != 0)
00144         {
00145                 rdata = (float*)EMUtil::em_malloc(num_bytes);
00146                 EMUtil::em_memcpy(rdata, data, num_bytes);
00147         }
00148 #ifdef EMAN2_USING_CUDA
00149         if (EMData::usecuda == 1 && num_bytes != 0 && that.cudarwdata != 0) {
00150                 //cout << "That copy constructor" << endl;
00151                 if(!rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
00152                 cudaError_t error = cudaMemcpy(cudarwdata,that.cudarwdata,num_bytes,cudaMemcpyDeviceToDevice);
00153                 if ( error != cudaSuccess ) throw UnexpectedBehaviorException("cudaMemcpy failed in EMData copy construction with error: " + string(cudaGetErrorString(error)));
00154         }
00155 #endif //EMAN2_USING_CUDA
00156 
00157         if (that.rot_fp != 0) rot_fp = new EMData(*(that.rot_fp));
00158 
00159         EMData::totalalloc++;
00160 #ifdef MEMDEBUG2
00161         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00162 #endif
00163 
00164         ENTERFUNC;
00165 }
00166 
00167 EMData& EMData::operator=(const EMData& that)
00168 {
00169         ENTERFUNC;
00170 
00171         if ( this != &that )
00172         {
00173                 free_memory(); // Free memory sets nx,ny and nz to 0
00174 
00175                 // Only copy the rdata if it exists, we could be in a scenario where only the header has been read
00176                 float* data = that.rdata;
00177                 size_t num_bytes = that.nx*that.ny*that.nz*sizeof(float);
00178                 if (data && num_bytes != 0)
00179                 {
00180                         nx = 1; // This prevents a memset in set_size
00181                         set_size(that.nx,that.ny,that.nz);
00182                         EMUtil::em_memcpy(rdata, data, num_bytes);
00183                 }
00184 
00185                 flags = that.flags;
00186 
00187                 all_translation = that.all_translation;
00188 
00189                 path = that.path;
00190                 pathnum = that.pathnum;
00191                 attr_dict = that.attr_dict;
00192 
00193                 xoff = that.xoff;
00194                 yoff = that.yoff;
00195                 zoff = that.zoff;
00196 
00197 #ifdef EMAN2_USING_CUDA
00198                 if (EMData::usecuda == 1 && num_bytes != 0 && that.cudarwdata != 0) {
00199                         if(cudarwdata){rw_free();}
00200                         if(!rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");;
00201                         cudaError_t error = cudaMemcpy(cudarwdata,that.cudarwdata,num_bytes,cudaMemcpyDeviceToDevice);
00202                         if ( error != cudaSuccess ) throw UnexpectedBehaviorException("cudaMemcpy failed in EMData copy construction with error: " + string(cudaGetErrorString(error)));
00203                 }
00204 #endif //EMAN2_USING_CUDA
00205 
00206                 changecount = that.changecount;
00207 
00208                 if (that.rot_fp != 0) rot_fp = new EMData(*(that.rot_fp));
00209                 else rot_fp = 0;
00210         }
00211         EXITFUNC;
00212         return *this;
00213 }
00214 
00215 EMData::EMData(int nx, int ny, int nz, bool is_real) :
00216 #ifdef EMAN2_USING_CUDA
00217                 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
00218 #endif //EMAN2_USING_CUDA
00219                 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0), zoff(0),
00220                 all_translation(),      path(""), pathnum(0), rot_fp(0)
00221 {
00222         ENTERFUNC;
00223 
00224         // used to replace cube 'pixel'
00225         attr_dict["apix_x"] = 1.0f;
00226         attr_dict["apix_y"] = 1.0f;
00227         attr_dict["apix_z"] = 1.0f;
00228 
00229         if(is_real) {   // create a real image [nx, ny, nz]
00230                 attr_dict["is_complex"] = int(0);
00231                 attr_dict["is_complex_x"] = int(0);
00232                 attr_dict["is_complex_ri"] = int(1);
00233                 set_size(nx, ny, nz);
00234         }
00235         else {  //create a complex image which real dimension is [ny, ny, nz]
00236                 int new_nx = nx + 2 - nx%2;
00237                 set_size(new_nx, ny, nz);
00238 
00239                 attr_dict["is_complex"] = int(1);
00240 
00241                 if(ny==1 && nz ==1)     {
00242                         attr_dict["is_complex_x"] = int(1);
00243                 }
00244                 else {
00245                         attr_dict["is_complex_x"] = int(0);
00246                 }
00247 
00248                 attr_dict["is_complex_ri"] = int(1);
00249                 attr_dict["is_fftpad"] = int(1);
00250 
00251                 if(nx%2 == 1) {
00252                         attr_dict["is_fftodd"] = 1;
00253                 }
00254         }
00255 
00256         this->to_zero();
00257         update();
00258         EMData::totalalloc++;
00259 #ifdef MEMDEBUG2
00260         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00261 #endif
00262 
00263         EXITFUNC;
00264 }
00265 
00266 
00267 EMData::EMData(float* data, const int x, const int y, const int z, const Dict& attr_dict) :
00268 #ifdef EMAN2_USING_CUDA
00269                 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(0),
00270 #endif //EMAN2_USING_CUDA
00271                 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
00272                 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
00273 {
00274         ENTERFUNC;
00275         // used to replace cube 'pixel'
00276         attr_dict["apix_x"] = 1.0f;
00277         attr_dict["apix_y"] = 1.0f;
00278         attr_dict["apix_z"] = 1.0f;
00279 
00280         EMData::totalalloc++;
00281 #ifdef MEMDEBUG2
00282         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00283 #endif
00284 
00285         update();
00286         EXITFUNC;
00287 }
00288 
00289 #ifdef EMAN2_USING_CUDA
00290 
00291 EMData::EMData(float* data, float* cudadata, const int x, const int y, const int z, const Dict& attr_dict) :
00292                 cudarwdata(cudadata), cudarodata(0), num_bytes(x*y*z*sizeof(float)), nextlistitem(0), prevlistitem(0), roneedsupdate(0), cudadirtybit(1),
00293                 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
00294                 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
00295 {
00296         ENTERFUNC;
00297 
00298         // used to replace cube 'pixel'
00299         attr_dict["apix_x"] = 1.0f;
00300         attr_dict["apix_y"] = 1.0f;
00301         attr_dict["apix_z"] = 1.0f;
00302 
00303         EMData::totalalloc++;
00304 #ifdef MEMDEBUG2
00305         printf("EMDATA+  %4d    %p\n",EMData::totalalloc,this);
00306 #endif
00307 
00308         update();
00309         EXITFUNC;
00310 }
00311 
00312 #endif //EMAN2_USING_CUDA
00313 
00314 //debug
00315 using std::cout;
00316 using std::endl;
00317 EMData::~EMData()
00318 {
00319         ENTERFUNC;
00320         free_memory();
00321 
00322 #ifdef EMAN2_USING_CUDA
00323         if(cudarwdata){rw_free();}
00324         if(cudarodata){ro_free();}
00325 #endif // EMAN2_USING_CUDA
00326         EMData::totalalloc--;
00327 #ifdef MEMDEBUG2
00328         printf("EMDATA-  %4d    %p\n",EMData::totalalloc,this);
00329 #endif
00330         EXITFUNC;
00331 }
00332 
00333 void EMData::clip_inplace(const Region & area,const float& fill_value)
00334 {
00335         // Added by d.woolford
00336         ENTERFUNC;
00337 
00338 //      printf("cip %d %d %d %d %d %d %f %d %d %d\n",area.origin[0],area.origin[1],area.origin[2],area.size[0],area.size[1],area.size[2],fill_value,nx,ny,nz);
00339         // Store the current dimension values
00340         int prev_nx = nx, prev_ny = ny, prev_nz = nz;
00341         size_t prev_size = (size_t)nx*ny*nz;
00342 
00343         // Get the zsize, ysize and xsize of the final area, these are the new dimension sizes of the pixel data
00344         int new_nz = ( area.size[2]==0 ? 1 : (int)area.size[2]);
00345         int new_ny = ( area.size[1]==0 ? 1 : (int)area.size[1]);
00346         int new_nx = (int)area.size[0];
00347 
00348         if ( new_nz < 0 || new_ny < 0 || new_nx < 0 )
00349         {
00350                 // Negative image dimensions were never tested nor considered when creating this implementation
00351                 throw ImageDimensionException("New image dimensions are negative - this is not supported in the clip_inplace operation");
00352         }
00353 
00354         size_t new_size = (size_t)new_nz*new_ny*new_nx;
00355 
00356         // Get the translation values, they are used to construct the ClipInplaceVariables object
00357         int x0 = (int) area.origin[0];
00358         int y0 = (int) area.origin[1];
00359         int z0 = (int) area.origin[2];
00360 
00361         // Get a object that calculates all the interesting variables associated with the clip inplace operation
00362         ClipInplaceVariables civ(prev_nx, prev_ny, prev_nz, new_nx, new_ny, new_nz, x0, y0, z0);
00363 
00364         get_data(); // Do this here to make sure rdata is up to date, applicable if GPU stuff is occuring
00365         // Double check to see if any memory shifting even has to occur
00366         if ( x0 > prev_nx || y0 > prev_ny || z0 > prev_nz || civ.x_iter == 0 || civ.y_iter == 0 || civ.z_iter == 0)
00367         {
00368                 // In this case the volume has been shifted beyond the location of the pixel rdata and
00369                 // the client should expect to see a volume with nothing in it.
00370 
00371                 // Set size calls realloc,
00372                 set_size(new_nx, new_ny, new_nz);
00373 
00374                 // Set pixel memory to zero - the client should expect to see nothing
00375                 EMUtil::em_memset(rdata, 0, (size_t)new_nx*new_ny*new_nz);
00376 
00377                 return;
00378         }
00379 
00380         // Resize the volume before memory shifting occurs if the new volume is larger than the previous volume
00381         // All of the pixel rdata is guaranteed to be at the start of the new volume because realloc (called in set size)
00382         // guarantees this.
00383         if ( new_size > prev_size )
00384                 set_size(new_nx, new_ny, new_nz);
00385 
00386         // Store the clipped row size.
00387         size_t clipped_row_size = (civ.x_iter) * sizeof(float);
00388 
00389         // Get the new sector sizes to save multiplication later.
00390         size_t new_sec_size = new_nx * new_ny;
00391         size_t prev_sec_size = prev_nx * prev_ny;
00392 
00393         // Determine the memory locations of the source and destination pixels - at the point nearest
00394         // to the beginning of the volume (rdata)
00395         size_t src_it_begin = civ.prv_z_bottom*prev_sec_size + civ.prv_y_front*prev_nx + civ.prv_x_left;
00396         size_t dst_it_begin = civ.new_z_bottom*new_sec_size + civ.new_y_front*new_nx + civ.new_x_left;
00397 
00398         // This loop is in the forward direction (starting at points nearest to the beginning of the volume)
00399         // it copies memory only when the destination pointer is less the source pointer - therefore
00400         // ensuring that no memory "copied to" is yet to be "copied from"
00401         for (int i = 0; i < civ.z_iter; ++i) {
00402                 for (int j = 0; j < civ.y_iter; ++j) {
00403 
00404                         // Determine the memory increments as dependent on i and j
00405                         // This could be optimized so that not so many multiplications are occurring...
00406                         size_t dst_inc = dst_it_begin + j*new_nx + i*new_sec_size;
00407                         size_t src_inc = src_it_begin + j*prev_nx + i*prev_sec_size;
00408                         float* local_dst = rdata + dst_inc;
00409                         float* local_src = rdata + src_inc;
00410 
00411                         if ( dst_inc >= src_inc )
00412                         {
00413                                 // this is fine, it will happen now and then and it will be necessary to continue.
00414                                 // the tempatation is to break, but you can't do that (because the point where memory intersects
00415                                 // could be in this slice - and yes, this aspect could be optimized).
00416                                 continue;
00417                         }
00418 
00419                         // Asserts are compiled only in debug mode
00420                         // This situation not encountered in testing thus far
00421                         Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
00422 
00423                         // Finally copy the memory
00424                         EMUtil::em_memcpy(local_dst, local_src, clipped_row_size);
00425                 }
00426         }
00427 
00428         // Determine the memory locations of the source and destination pixels - at the point nearest
00429         // to the end of the volume (rdata+new_size)
00430         size_t src_it_end = prev_size - civ.prv_z_top*prev_sec_size - civ.prv_y_back*prev_nx - prev_nx + civ.prv_x_left;
00431         size_t dst_it_end = new_size - civ.new_z_top*new_sec_size - civ.new_y_back*new_nx - new_nx + civ.new_x_left;
00432 
00433         // This loop is in the reverse direction (starting at points nearest to the end of the volume).
00434         // It copies memory only when the destination pointer is greater than  the source pointer therefore
00435         // ensuring that no memory "copied to" is yet to be "copied from"
00436         for (int i = 0; i < civ.z_iter; ++i) {
00437                 for (int j = 0; j < civ.y_iter; ++j) {
00438 
00439                         // Determine the memory increments as dependent on i and j
00440                         size_t dst_inc = dst_it_end - j*new_nx - i*new_sec_size;
00441                         size_t src_inc = src_it_end - j*prev_nx - i*prev_sec_size;
00442                         float* local_dst = rdata + dst_inc;
00443                         float* local_src = rdata + src_inc;
00444 
00445                         if (dst_inc <= (src_inc + civ.x_iter ))
00446                         {
00447                                 // Overlap
00448                                 if ( dst_inc > src_inc )
00449                                 {
00450                                         // Because the memcpy operation is the forward direction, and this "reverse
00451                                         // direction" loop is proceeding in a backwards direction, it is possible
00452                                         // that memory copied to is yet to be copied from (because memcpy goes forward).
00453                                         // In this scenario pixel memory "copied to" is yet to be "copied from"
00454                                         // i.e. there is overlap
00455 
00456                                         // memmove handles overlapping cases.
00457                                         // memmove could use a temporary buffer, or could go just go backwards
00458                                         // the specification doesn't say how the function behaves...
00459                                         // If memmove creates a temporary buffer is clip_inplace no longer inplace?
00460                                         memmove(local_dst, local_src, clipped_row_size);
00461                                 }
00462                                 continue;
00463                         }
00464 
00465                         // This situation not encountered in testing thus far
00466                         Assert( dst_inc < new_size && src_inc < prev_size && dst_inc >= 0 && src_inc >= 0 );
00467 
00468                         // Perform the memory copy
00469                         EMUtil::em_memcpy(local_dst, local_src, clipped_row_size);
00470                 }
00471         }
00472 
00473         // Resize the volume after memory shifting occurs if the new volume is smaller than the previous volume
00474         // set_size calls realloc, guaranteeing that the pixel rdata is in the right location.
00475         if ( new_size < prev_size )
00476                 set_size(new_nx, new_ny, new_nz);
00477 
00478         // Now set all the edges to zero
00479 
00480         // Set the extra bottom z slices to the fill_value
00481         if (  z0 < 0 )
00482         {
00483                 //EMUtil::em_memset(rdata, 0, (-z0)*new_sec_size*sizeof(float));
00484                 size_t inc = (-z0)*new_sec_size;
00485                 std::fill(rdata,rdata+inc,fill_value);
00486         }
00487 
00488         // Set the extra top z slices to the fill_value
00489         if (  civ.new_z_top > 0 )
00490         {
00491                 float* begin_pointer = rdata + (new_nz-civ.new_z_top)*new_sec_size;
00492                 //EMUtil::em_memset(begin_pointer, 0, (civ.new_z_top)*new_sec_size*sizeof(float));
00493                 float* end_pointer = begin_pointer+(civ.new_z_top)*new_sec_size;
00494                 std::fill(begin_pointer,end_pointer,fill_value);
00495         }
00496 
00497         // Next deal with x and y edges by iterating through each slice
00498         for ( int i = civ.new_z_bottom; i < civ.new_z_bottom + civ.z_iter; ++i )
00499         {
00500                 // Set the extra front y components to the fill_value
00501                 if ( y0 < 0 )
00502                 {
00503                         float* begin_pointer = rdata + i*new_sec_size;
00504                         //EMUtil::em_memset(begin_pointer, 0, (-y0)*new_nx*sizeof(float));
00505                         float* end_pointer = begin_pointer+(-y0)*new_nx;
00506                         std::fill(begin_pointer,end_pointer,fill_value);
00507                 }
00508 
00509                 // Set the extra back y components to the fill_value
00510                 if ( civ.new_y_back > 0 )
00511                 {
00512                         float* begin_pointer = rdata + i*new_sec_size + (new_ny-civ.new_y_back)*new_nx;
00513                         //EMUtil::em_memset(begin_pointer, 0, (civ.new_y_back)*new_nx*sizeof(float));
00514                         float* end_pointer = begin_pointer+(civ.new_y_back)*new_nx;
00515                         std::fill(begin_pointer,end_pointer,fill_value);
00516                 }
00517 
00518                 // Iterate through the y to set each correct x component to the fill_value
00519                 for (int j = civ.new_y_front; j <civ.new_y_front + civ.y_iter; ++j)
00520                 {
00521                         // Set the extra left x components to the fill_value
00522                         if ( x0 < 0 )
00523                         {
00524                                 float* begin_pointer = rdata + i*new_sec_size + j*new_nx;
00525                                 //EMUtil::em_memset(begin_pointer, 0, (-x0)*sizeof(float));
00526                                 float* end_pointer = begin_pointer+(-x0);
00527                                 std::fill(begin_pointer,end_pointer,fill_value);
00528                         }
00529 
00530                         // Set the extra right x components to the fill_value
00531                         if ( civ.new_x_right > 0 )
00532                         {
00533                                 float* begin_pointer = rdata + i*new_sec_size + j*new_nx + (new_nx - civ.new_x_right);
00534                                 //EMUtil::em_memset(begin_pointer, 0, (civ.new_x_right)*sizeof(float));
00535                                 float* end_pointer = begin_pointer+(civ.new_x_right);
00536                                 std::fill(begin_pointer,end_pointer,fill_value);
00537                         }
00538 
00539                 }
00540         }
00541 
00542 // These couts may be useful
00543 //      cout << "start starts " << civ.prv_x_left << " " << civ.prv_y_front << " " << civ.prv_z_bottom << endl;
00544 //      cout << "start ends " << civ.prv_x_right << " " << civ.prv_y_back << " " << civ.prv_z_top << endl;
00545 //      cout << "dst starts " << civ.new_x_left << " " << civ.new_y_front << " " << civ.new_z_bottom << endl;
00546 //      cout << "dst ends " << civ.new_x_right << " " << civ.new_y_back << " " << civ.new_z_top << endl;
00547 //      cout << "total iter z - " << civ.z_iter << " y - " << civ.y_iter << " x - " << civ.x_iter << endl;
00548 //      cout << "=====" << endl;
00549 //      cout << "dst_end is " << dst_it_end << " src end is " << src_it_end << endl;
00550 //      cout << "dst_begin is " << dst_it_begin << " src begin is " << src_it_begin << endl;
00551 
00552         // Update appropriate attributes (Copied and pasted from get_clip)
00553         if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
00554         attr_dict.has_key("origin_z") )
00555         {
00556                 float xorigin = attr_dict["origin_x"];
00557                 float yorigin = attr_dict["origin_y"];
00558                 float zorigin = attr_dict["origin_z"];
00559 
00560                 float apix_x = attr_dict["apix_x"];
00561                 float apix_y = attr_dict["apix_y"];
00562                 float apix_z = attr_dict["apix_z"];
00563 
00564                 set_xyz_origin(xorigin + apix_x * area.origin[0],
00565                         yorigin + apix_y * area.origin[1],
00566                         zorigin + apix_z * area.origin[2]);
00567         }
00568 
00569         // Set the update flag because the size of the image has changed and stats should probably be recalculated if requested.
00570         update();
00571 
00572         EXITFUNC;
00573 }
00574 
00575 EMData *EMData::get_clip(const Region & area, const float fill) const
00576 {
00577         ENTERFUNC;
00578         if (get_ndim() != area.get_ndim()) {
00579                 LOGERR("cannot get %dD clip out of %dD image", area.get_ndim(),get_ndim());
00580                 return 0;
00581         }
00582 
00583         EMData *result = new EMData();
00584 
00585         // Ensure that all of the metadata of this is stored in the new object
00586         // Originally added to ensure that euler angles were retained when preprocessing (zero padding) images
00587         // prior to insertion into the 3D for volume in the reconstruction phase (see reconstructor.cpp/h).
00588         result->attr_dict = this->attr_dict;
00589         int zsize = (int)area.size[2];
00590         if (zsize == 0 && nz <= 1) {
00591                 zsize = 1;
00592         }
00593         int ysize = (ny<=1 && (int)area.size[1]==0 ? 1 : (int)area.size[1]);
00594 
00595         if ( (int)area.size[0] < 0 || ysize < 0 || zsize < 0 )
00596         {
00597                 // Negative image dimensions not supported - added retrospectively by d.woolford (who didn't write get_clip but wrote clip_inplace)
00598                 throw ImageDimensionException("New image dimensions are negative - this is not supported in the the get_clip operation");
00599         }
00600 
00601 //#ifdef EMAN2_USING_CUDA
00602         // Strategy is always to prefer using the GPU if possible
00603 //      bool use_gpu = false;
00604 //      if ( gpu_operation_preferred() ) {
00605 //              result->set_size_cuda((int)area.size[0], ysize, zsize);
00606                 //CudaDataLock lock(this); // Just so we never have to recopy this data to and from the GPU
00607 //              result->get_cuda_data(); // Force the allocation - set_size_cuda is lazy
00608                 // Setting the value is necessary seeing as cuda data is not automatically zeroed
00609 //              result->to_value(fill); // This will automatically use the GPU.
00610 //              use_gpu = true;
00611 //      } else { // cpu == True
00612 //              result->set_size((int)area.size[0], ysize, zsize);
00613 //              if (fill != 0.0) { result->to_value(fill); };
00614 //      }
00615 //#else
00616         result->set_size((int)area.size[0], ysize, zsize);
00617         if (fill != 0.0) { result->to_value(fill); };
00618 //#endif //EMAN2_USING_CUDA
00619 
00620         int x0 = (int) area.origin[0];
00621         x0 = x0 < 0 ? 0 : x0;
00622 
00623         int y0 = (int) area.origin[1];
00624         y0 = y0 < 0 ? 0 : y0;
00625 
00626         int z0 = (int) area.origin[2];
00627         z0 = z0 < 0 ? 0 : z0;
00628 
00629         int x1 = (int) (area.origin[0] + area.size[0]);
00630         x1 = x1 > nx ? nx : x1;
00631 
00632         int y1 = (int) (area.origin[1] + area.size[1]);
00633         y1 = y1 > ny ? ny : y1;
00634 
00635         int z1 = (int) (area.origin[2] + area.size[2]);
00636         z1 = z1 > nz ? nz : z1;
00637         if (z1 <= 0) {
00638                 z1 = 1;
00639         }
00640 
00641         result->insert_clip(this,-((IntPoint)area.origin));
00642 
00643         if( attr_dict.has_key("apix_x") && attr_dict.has_key("apix_y") &&
00644                 attr_dict.has_key("apix_z") )
00645         {
00646                 if( attr_dict.has_key("origin_x") && attr_dict.has_key("origin_y") &&
00647                     attr_dict.has_key("origin_z") )
00648                 {
00649                         float xorigin = attr_dict["origin_x"];
00650                         float yorigin = attr_dict["origin_y"];
00651                         float zorigin = attr_dict["origin_z"];
00652 
00653                         float apix_x = attr_dict["apix_x"];
00654                         float apix_y = attr_dict["apix_y"];
00655                         float apix_z = attr_dict["apix_z"];
00656 
00657                         result->set_xyz_origin(xorigin + apix_x * area.origin[0],
00658                                                                    yorigin + apix_y * area.origin[1],
00659                                                                zorigin + apix_z * area.origin[2]);
00660                 }
00661         }
00662 
00663 //#ifdef EMAN2_USING_CUDA
00664 //      if (use_gpu) result->gpu_update();
00665 //      else result->update();
00666 //#else
00667         result->update();
00668 //#endif // EMAN2_USING_CUDA
00669 
00670 
00671         result->set_path(path);
00672         result->set_pathnum(pathnum);
00673 
00674         EXITFUNC;
00675         return result;
00676 }
00677 
00678 
00679 EMData *EMData::get_top_half() const
00680 {
00681         ENTERFUNC;
00682 
00683         if (get_ndim() != 3) {
00684                 throw ImageDimensionException("3D only");
00685         }
00686 
00687         EMData *half = new EMData();
00688         half->attr_dict = attr_dict;
00689         half->set_size(nx, ny, nz / 2);
00690 
00691         float *half_data = half->get_data();
00692         EMUtil::em_memcpy(half_data, &(get_data()[(size_t)nz / 2 * (size_t)nx * (size_t)ny]), sizeof(float) * (size_t)nx * (size_t)ny * (size_t)nz / 2lu);
00693 
00694         float apix_z = attr_dict["apix_z"];
00695         float origin_z = attr_dict["origin_z"];
00696         origin_z += apix_z * nz / 2;
00697         half->attr_dict["origin_z"] = origin_z;
00698         half->update();
00699 
00700         EXITFUNC;
00701         return half;
00702 }
00703 
00704 
00705 EMData *EMData::get_rotated_clip(const Transform &xform,
00706                                                                  const IntSize &size, float)
00707 {
00708         EMData *result = new EMData();
00709         result->set_size(size[0],size[1],size[2]);
00710 
00711         if (nz==1) {
00712                 for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
00713                         for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
00714                                 Vec3f xv=xform.transform(Vec3f((float)x,(float)y,0.0f));
00715                                 float v = 0;
00716 
00717                                 if (xv[0]<0||xv[1]<0||xv[0]>nx-2||xv[1]>ny-2) v=0.;
00718                                 else v=sget_value_at_interp(xv[0],xv[1]);
00719                                 result->set_value_at(x+size[0]/2,y+size[1]/2,v);
00720                         }
00721                 }
00722         }
00723         else {
00724                 for (int z=-size[2]/2; z<(size[2]+1)/2; z++) {
00725                         for (int y=-size[1]/2; y<(size[1]+1)/2; y++) {
00726                                 for (int x=-size[0]/2; x<(size[0]+1)/2; x++) {
00727                                         Vec3f xv=xform.transform(Vec3f((float)x,(float)y,0.0f));
00728                                         float v = 0;
00729 
00730                                         if (xv[0]<0||xv[1]<0||xv[2]<0||xv[0]>nx-2||xv[1]>ny-2||xv[2]>nz-2) v=0.;
00731                                         else v=sget_value_at_interp(xv[0],xv[1],xv[2]);
00732                                         result->set_value_at(x+size[0]/2,y+size[1]/2,z+size[2]/2,v);
00733                                 }
00734                         }
00735                 }
00736         }
00737         result->update();
00738 
00739         return result;
00740 }
00741 
00742 
00743 EMData* EMData::window_center(int l) {
00744         ENTERFUNC;
00745         // sanity checks
00746         int n = nx;
00747         if (is_complex()) {
00748                 LOGERR("Need real-space data for window_center()");
00749                 throw ImageFormatException(
00750                         "Complex input image; real-space expected.");
00751         }
00752         if (is_fftpadded()) {
00753                 // image has been fft-padded, compute the real-space size
00754                 n -= (2 - int(is_fftodd()));
00755         }
00756         int corner = n/2 - l/2;
00757         int ndim = get_ndim();
00758         EMData* ret;
00759         switch (ndim) {
00760                 case 3:
00761                         if ((n != ny) || (n != nz)) {
00762                                 LOGERR("Need the real-space image to be cubic.");
00763                                 throw ImageFormatException(
00764                                                 "Need cubic real-space image.");
00765                         }
00766                         ret = get_clip(Region(corner, corner, corner, l, l, l));
00767                         break;
00768                 case 2:
00769                         if (n != ny) {
00770                                 LOGERR("Need the real-space image to be square.");
00771                                 throw ImageFormatException(
00772                                                 "Need square real-space image.");
00773                         }
00774                         //cout << "Using corner " << corner << endl;
00775                         ret = get_clip(Region(corner, corner, l, l));
00776                         break;
00777                 case 1:
00778                         ret = get_clip(Region(corner, l));
00779                         break;
00780                 default:
00781                         throw ImageDimensionException(
00782                                         "window_center only supports 1-d, 2-d, and 3-d images");
00783         }
00784         return ret;
00785         EXITFUNC;
00786 }
00787 
00788 
00789 float *EMData::setup4slice(bool redo)
00790 {
00791         ENTERFUNC;
00792 
00793         if (!is_complex()) {
00794                 throw ImageFormatException("complex image only");
00795         }
00796 
00797         if (get_ndim() != 3) {
00798                 throw ImageDimensionException("3D only");
00799         }
00800 
00801         if (supp) {
00802                 if (redo) {
00803                         EMUtil::em_free(supp);
00804                         supp = 0;
00805                 }
00806                 else {
00807                         EXITFUNC;
00808                         return supp;
00809                 }
00810         }
00811 
00812         const int SUPP_ROW_SIZE = 8;
00813         const int SUPP_ROW_OFFSET = 4;
00814         const int supp_size = SUPP_ROW_SIZE + SUPP_ROW_OFFSET;
00815 
00816         supp = (float *) EMUtil::em_calloc(supp_size * ny * nz, sizeof(float));
00817         int nxy = nx * ny;
00818         int supp_xy = supp_size * ny;
00819         float * data = get_data();
00820 
00821         for (int z = 0; z < nz; z++) {
00822                 size_t cur_z1 = z * nxy;
00823                 size_t cur_z2 = z * supp_xy;
00824 
00825                 for (int y = 0; y < ny; y++) {
00826                         size_t cur_y1 = y * nx;
00827                         size_t cur_y2 = y * supp_size;
00828 
00829                         for (int x = 0; x < SUPP_ROW_SIZE; x++) {
00830                                 size_t k = (x + SUPP_ROW_OFFSET) + cur_y2 + cur_z2;
00831                                 supp[k] = data[x + cur_y1 + cur_z1];
00832                         }
00833                 }
00834         }
00835 
00836         for (int z = 1, zz = nz - 1; z < nz; z++, zz--) {
00837                 size_t cur_z1 = zz * nxy;
00838                 size_t cur_z2 = z * supp_xy;
00839 
00840                 for (int y = 1, yy = ny - 1; y < ny; y++, yy--) {
00841                         supp[y * 12 + cur_z2] = data[4 + yy * nx + cur_z1];
00842                         supp[1 + y * 12 + cur_z2] = -data[5 + yy * nx + cur_z1];
00843                         supp[2 + y * 12 + cur_z2] = data[2 + yy * nx + cur_z1];
00844                         supp[3 + y * 12 + cur_z2] = -data[3 + yy * nx + cur_z1];
00845                 }
00846         }
00847 
00848         EXITFUNC;
00849         return supp;
00850 }
00851 
00852 
00853 void EMData::scale(float s)
00854 {
00855         ENTERFUNC;
00856         Transform t;
00857         t.set_scale(s);
00858         transform(t);
00859         EXITFUNC;
00860 }
00861 
00862 
00863 void EMData::translate(int dx, int dy, int dz)
00864 {
00865         ENTERFUNC;
00866         translate(Vec3i(dx, dy, dz));
00867         EXITFUNC;
00868 }
00869 
00870 
00871 void EMData::translate(float dx, float dy, float dz)
00872 {
00873         ENTERFUNC;
00874         int dx_ = Util::round(dx);
00875         int dy_ = Util::round(dy);
00876         int dz_ = Util::round(dz);
00877         if( ( (dx-dx_) == 0 ) && ( (dy-dy_) == 0 ) && ( (dz-dz_) == 0 )) {
00878                 translate(dx_, dy_, dz_);
00879         }
00880         else {
00881                 translate(Vec3f(dx, dy, dz));
00882         }
00883         EXITFUNC;
00884 }
00885 
00886 
00887 void EMData::translate(const Vec3i &translation)
00888 {
00889         ENTERFUNC;
00890 
00891         //if traslation is 0, do nothing
00892         if( translation[0] == 0 && translation[1] == 0 && translation[2] == 0) {
00893                 EXITFUNC;
00894                 return;
00895         }
00896 
00897         Dict params("trans",static_cast< vector<int> >(translation));
00898         process_inplace("xform.translate.int",params);
00899 
00900         // update() - clip_inplace does the update
00901         all_translation += translation;
00902 
00903         EXITFUNC;
00904 }
00905 
00906 
00907 void EMData::translate(const Vec3f &translation)
00908 {
00909         ENTERFUNC;
00910 
00911         if( translation[0] == 0.0f && translation[1] == 0.0f && translation[2] == 0.0f ) {
00912                 EXITFUNC;
00913                 return;
00914         }
00915 
00916         Transform* t = new Transform();
00917         t->set_trans(translation);
00918         process_inplace("xform",Dict("transform",t));
00919         delete t;
00920 
00921         all_translation += translation;
00922         EXITFUNC;
00923 }
00924 
00925 
00926 void EMData::rotate(float az, float alt, float phi)
00927 {
00928         Dict d("type","eman");
00929         d["az"] = az;
00930         d["alt"] = alt;
00931         d["phi"] = phi;
00932         Transform t(d);
00933         transform(t);
00934 }
00935 
00936 
00937 
00938 void EMData::rotate(const Transform & t)
00939 {
00940         cout << "Deprecation warning in EMData::rotate. Please consider using EMData::transform() instead " << endl;
00941         transform(t);
00942 }
00943 
00944 float EMData::max_3D_pixel_error(const Transform &t1, const Transform & t2, float r) {
00945         
00946         Transform t;
00947         int r0 = (int)r;
00948         float ddmax = 0.0f;
00949 
00950         t = t2*t1.inverse();
00951         for (int i=0; i<int(2*M_PI*r0+0.5); i++) {
00952                 float ang = (float)i/r;
00953                 Vec3f v = Vec3f(r0*cos(ang), r0*sin(ang), 0.0f);
00954                 Vec3f d = t*v-v;
00955 #ifdef _WIN32
00956                 ddmax = _cpp_max(ddmax,d[0]*d[0]+d[1]*d[1]+d[2]*d[2]);
00957 #else
00958                 ddmax = std::max(ddmax,d[0]*d[0]+d[1]*d[1]+d[2]*d[2]);
00959 #endif  //_WIN32
00960         }
00961         return std::sqrt(ddmax);
00962 }
00963 
00964 void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy, float dz)
00965 {
00966         cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
00967 //      Transform3D t( az, alt, phi,Vec3f(dx, dy, dz));
00968         Transform t;
00969         t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
00970         t.set_trans(dx, dy, dz);
00971         rotate_translate(t);
00972 }
00973 
00974 
00975 void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy,
00976                                                           float dz, float pdx, float pdy, float pdz)
00977 {
00978         cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
00979 //      Transform3D t(Vec3f(dx, dy, dz), az, alt, phi, Vec3f(pdx,pdy,pdz));
00980 //      rotate_translate(t);
00981 
00982         Transform t;
00983         t.set_pre_trans(Vec3f(dx, dy, dz));
00984         t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
00985         t.set_trans(pdx, pdy, pdz);
00986         rotate_translate(t);
00987 }
00988 
00989 //void EMData::rotate_translate(const Transform3D & RA)
00990 //{
00991 //      cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
00992 //      ENTERFUNC;
00993 //
00994 //#if EMDATA_EMAN2_DEBUG
00995 //      std::cout << "start rotate_translate..." << std::endl;
00996 //#endif
00997 //
00998 //      float scale       = RA.get_scale();
00999 //      Vec3f dcenter     = RA.get_center();
01000 //      Vec3f translation = RA.get_posttrans();
01001 //      Dict rotation      = RA.get_rotation(Transform3D::EMAN);
01003 //      Transform3D RAInv = RA.inverse(); // We're rotating the coordinate system, not the data
01005 //#if EMDATA_EMAN2_DEBUG
01006 //      vector<string> keys = rotation.keys();
01007 //      vector<string>::const_iterator it;
01008 //      for(it=keys.begin(); it!=keys.end(); ++it) {
01010 //              std::cout << *it << " : " << (float)rotation.get(*it) << std::endl;
01011 //      }
01012 //#endif
01013 //      float inv_scale = 1.0f;
01014 //
01015 //      if (scale != 0) {
01016 //              inv_scale = 1.0f / scale;
01017 //      }
01018 //
01019 //      float *src_data = 0;
01020 //      float *des_data = 0;
01021 //
01022 //      src_data = get_data();
01023 //      des_data = (float *) EMUtil::em_malloc(nx * ny * nz * sizeof(float));
01024 //
01025 //      if (nz == 1) {
01026 //              float x2c =  nx / 2 - dcenter[0] + RAInv[0][3];
01027 //              float y2c =  ny / 2 - dcenter[1] + RAInv[1][3];
01028 //              float y   = -ny / 2 + dcenter[1]; // changed 0 to 1 in dcenter and below
01029 //              for (int j = 0; j < ny; j++, y += 1.0f) {
01030 //                      float x = -nx / 2 + dcenter[0];
01031 //                      for (int i = 0; i < nx; i++, x += 1.0f) {
01032 //                              float x2 = RAInv[0][0]*x + RAInv[0][1]*y + x2c;
01033 //                              float y2 = RAInv[1][0]*x + RAInv[1][1]*y + y2c;
01034 //
01035 //                              if (x2 < 0 || x2 >= nx || y2 < 0 || y2 >= ny ) {
01036 //                                      des_data[i + j * nx] = 0; // It may be tempting to set this value to the
01037 //                                      // mean but in fact this is not a good thing to do. Talk to S.Ludtke about it.
01038 //                              }
01039 //                              else {
01040 //                                      int ii = Util::fast_floor(x2);
01041 //                                      int jj = Util::fast_floor(y2);
01042 //                                      int k0 = ii + jj * nx;
01043 //                                      int k1 = k0 + 1;
01044 //                                      int k2 = k0 + nx;
01045 //                                      int k3 = k0 + nx + 1;
01046 //
01047 //                                      if (ii == nx - 1) {
01048 //                                              k1--;
01049 //                                              k3--;
01050 //                                      }
01051 //                                      if (jj == ny - 1) {
01052 //                                              k2 -= nx;
01053 //                                              k3 -= nx;
01054 //                                      }
01055 //
01056 //                                      float t = x2 - ii;
01057 //                                      float u = y2 - jj;
01058 //
01059 //                                      des_data[i + j * nx] = Util::bilinear_interpolate(src_data[k0],src_data[k1], src_data[k2], src_data[k3],t,u); // This is essentially basic interpolation
01060 //                              }
01061 //                      }
01062 //              }
01063 //      }
01064 //      else {
01065 //#if EMDATA_EMAN2_DEBUG
01066 //              std::cout << "This is the 3D case." << std::endl    ;
01067 //#endif
01068 //
01069 //              Transform3D mx = RA;
01070 //              mx.set_scale(inv_scale);
01071 //
01072 //#if EMDATA_EMAN2_DEBUG
01075 //#endif
01076 //
01077 //              int nxy = nx * ny;
01078 //              int l = 0;
01079 //
01080 //              float x2c =  nx / 2 - dcenter[0] + RAInv[0][3];;
01081 //              float y2c =  ny / 2 - dcenter[1] + RAInv[1][3];;
01082 //              float z2c =  nz / 2 - dcenter[2] + RAInv[2][3];;
01083 //
01084 //              float z   = -nz / 2 + dcenter[2]; //
01085 //
01086 //              size_t ii, k0, k1, k2, k3, k4, k5, k6, k7;
01087 //              for (int k = 0; k < nz; k++, z += 1.0f) {
01088 //                      float y   = -ny / 2 + dcenter[1]; //
01089 //                      for (int j = 0; j < ny; j++,   y += 1.0f) {
01090 //                              float x = -nx / 2 + dcenter[0];
01091 //                              for (int i = 0; i < nx; i++, l++ ,  x += 1.0f) {
01092 //                                      float x2 = RAInv[0][0] * x + RAInv[0][1] * y + RAInv[0][2] * z + x2c;
01093 //                                      float y2 = RAInv[1][0] * x + RAInv[1][1] * y + RAInv[1][2] * z + y2c;
01094 //                                      float z2 = RAInv[2][0] * x + RAInv[2][1] * y + RAInv[2][2] * z + z2c;
01095 //
01096 //
01097 //                                      if (x2 < 0 || y2 < 0 || z2 < 0 ||
01098 //                                              x2 >= nx  || y2 >= ny  || z2>= nz ) {
01099 //                                              des_data[l] = 0;
01100 //                                      }
01101 //                                      else {
01102 //                                              int ix = Util::fast_floor(x2);
01103 //                                              int iy = Util::fast_floor(y2);
01104 //                                              int iz = Util::fast_floor(z2);
01105 //                                              float tuvx = x2-ix;
01106 //                                              float tuvy = y2-iy;
01107 //                                              float tuvz = z2-iz;
01108 //                                              ii = ix + iy * nx + iz * nxy;
01109 //
01110 //                                              k0 = ii;
01111 //                                              k1 = k0 + 1;
01112 //                                              k2 = k0 + nx;
01113 //                                              k3 = k0 + nx+1;
01114 //                                              k4 = k0 + nxy;
01115 //                                              k5 = k1 + nxy;
01116 //                                              k6 = k2 + nxy;
01117 //                                              k7 = k3 + nxy;
01118 //
01119 //                                              if (ix == nx - 1) {
01120 //                                                      k1--;
01121 //                                                      k3--;
01122 //                                                      k5--;
01123 //                                                      k7--;
01124 //                                              }
01125 //                                              if (iy == ny - 1) {
01126 //                                                      k2 -= nx;
01127 //                                                      k3 -= nx;
01128 //                                                      k6 -= nx;
01129 //                                                      k7 -= nx;
01130 //                                              }
01131 //                                              if (iz == nz - 1) {
01132 //                                                      k4 -= nxy;
01133 //                                                      k5 -= nxy;
01134 //                                                      k6 -= nxy;
01135 //                                                      k7 -= nxy;
01136 //                                              }
01137 //
01138 //                                              des_data[l] = Util::trilinear_interpolate(src_data[k0],
01139 //                                                        src_data[k1],
01140 //                                                        src_data[k2],
01141 //                                                        src_data[k3],
01142 //                                                        src_data[k4],
01143 //                                                        src_data[k5],
01144 //                                                        src_data[k6],
01145 //                                                        src_data[k7],
01146 //                                                        tuvx, tuvy, tuvz);
01147 //#if EMDATA_EMAN2_DEBUG
01148 //                                              printf(" ix=%d \t iy=%d \t iz=%d \t value=%f \n", ix ,iy, iz, des_data[l] );
01149 //                                              std::cout << src_data[ii] << std::endl;
01150 //#endif
01151 //                                      }
01152 //                              }
01153 //                      }
01154 //              }
01155 //      }
01156 //
01157 //      if( rdata )
01158 //      {
01159 //              EMUtil::em_free(rdata);
01160 //              rdata = 0;
01161 //      }
01162 //      rdata = des_data;
01163 //
01164 //      scale_pixel(inv_scale);
01165 //
01166 //      attr_dict["origin_x"] = (float) attr_dict["origin_x"] * inv_scale;
01167 //      attr_dict["origin_y"] = (float) attr_dict["origin_y"] * inv_scale;
01168 //      attr_dict["origin_z"] = (float) attr_dict["origin_z"] * inv_scale;
01169 //
01170 //      update();
01171 //      all_translation += translation;
01172 //      EXITFUNC;
01173 //}
01174 
01175 
01176 
01177 
01178 void EMData::rotate_x(int dx)
01179 {
01180         ENTERFUNC;
01181 
01182         if (get_ndim() > 2) {
01183                 throw ImageDimensionException("no 3D image");
01184         }
01185 
01186 
01187         size_t row_size = nx * sizeof(float);
01188         float *tmp = (float*)EMUtil::em_malloc(row_size);
01189         float * data = get_data();
01190 
01191         for (int y = 0; y < ny; y++) {
01192                 int y_nx = y * nx;
01193                 for (int x = 0; x < nx; x++) {
01194                         tmp[x] = data[y_nx + (x + dx) % nx];
01195                 }
01196                 EMUtil::em_memcpy(&data[y_nx], tmp, row_size);
01197         }
01198 
01199         update();
01200         if( tmp )
01201         {
01202                 delete[]tmp;
01203                 tmp = 0;
01204         }
01205         EXITFUNC;
01206 }
01207 
01208 double EMData::dot_rotate_translate(EMData * with, float dx, float dy, float da, const bool mirror)
01209 {
01210         ENTERFUNC;
01211 
01212         if (!EMUtil::is_same_size(this, with)) {
01213                 LOGERR("images not same size");
01214                 throw ImageFormatException("images not same size");
01215         }
01216 
01217         if (get_ndim() == 3) {
01218                 LOGERR("1D/2D Images only");
01219                 throw ImageDimensionException("1D/2D only");
01220         }
01221 
01222         float *this_data = 0;
01223 
01224         this_data = get_data();
01225 
01226         float da_rad = da*(float)M_PI/180.0f;
01227 
01228         float *with_data = with->get_data();
01229         float mx0 = cos(da_rad);
01230         float mx1 = sin(da_rad);
01231         float y = -ny / 2.0f;
01232         float my0 = mx0 * (-nx / 2.0f - 1.0f) + nx / 2.0f - dx;
01233         float my1 = -mx1 * (-nx / 2.0f - 1.0f) + ny / 2.0f - dy;
01234         double result = 0;
01235 
01236         for (int j = 0; j < ny; j++) {
01237                 float x2 = my0 + mx1 * y;
01238                 float y2 = my1 + mx0 * y;
01239 
01240                 int ii = Util::fast_floor(x2);
01241                 int jj = Util::fast_floor(y2);
01242                 float t = x2 - ii;
01243                 float u = y2 - jj;
01244 
01245                 for (int i = 0; i < nx; i++) {
01246                         t += mx0;
01247                         u -= mx1;
01248 
01249                         if (t >= 1.0f) {
01250                                 ii++;
01251                                 t -= 1.0f;
01252                         }
01253 
01254                         if (u >= 1.0f) {
01255                                 jj++;
01256                                 u -= 1.0f;
01257                         }
01258 
01259                         if (t < 0) {
01260                                 ii--;
01261                                 t += 1.0f;
01262                         }
01263 
01264                         if (u < 0) {
01265                                 jj--;
01266                                 u += 1.0f;
01267                         }
01268 
01269                         if (ii >= 0 && ii <= nx - 2 && jj >= 0 && jj <= ny - 2) {
01270                                 int k0 = ii + jj * nx;
01271                                 int k1 = k0 + 1;
01272                                 int k2 = k0 + nx + 1;
01273                                 int k3 = k0 + nx;
01274 
01275                                 float tt = 1 - t;
01276                                 float uu = 1 - u;
01277                                 int idx = i + j * nx;
01278                                 if (mirror) idx = nx-1-i+j*nx; // mirroring of Transforms is always about the y axis
01279                                 result += (this_data[k0] * tt * uu + this_data[k1] * t * uu +
01280                                                    this_data[k2] * t * u + this_data[k3] * tt * u) * with_data[idx];
01281                         }
01282                 }
01283                 y += 1.0f;
01284         }
01285 
01286         EXITFUNC;
01287         return result;
01288 }
01289 
01290 
01291 EMData *EMData::little_big_dot(EMData * with, bool do_sigma)
01292 {
01293         ENTERFUNC;
01294 
01295         if (get_ndim() > 2) {
01296                 throw ImageDimensionException("1D/2D only");
01297         }
01298 
01299         EMData *ret = copy_head();
01300         ret->set_size(nx,ny,nz);
01301         ret->to_zero();
01302 
01303         int nx2 = with->get_xsize();
01304         int ny2 = with->get_ysize();
01305         float em = with->get_edge_mean();
01306 
01307         float *data = get_data();
01308         float *with_data = with->get_data();
01309         float *ret_data = ret->get_data();
01310 
01311         float sum2 = (Util::square((float)with->get_attr("sigma")) +
01312                                   Util::square((float)with->get_attr("mean")));
01313         if (do_sigma) {
01314                 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
01315                         for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
01316                                 float sum = 0;
01317                                 float sum1 = 0;
01318                                 float summ = 0;
01319                                 int k = 0;
01320 
01321                                 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01322                                         for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01323                                                 int l = ii + jj * nx;
01324                                                 sum1 += Util::square(data[l]);
01325                                                 summ += data[l];
01326                                                 sum += data[l] * with_data[k];
01327                                                 k++;
01328                                         }
01329                                 }
01330                                 float tmp_f1 = (sum1 / 2.0f - sum) / (nx2 * ny2);
01331                                 float tmp_f2 = Util::square((float)with->get_attr("mean") -
01332                                                                                         summ / (nx2 * ny2));
01333                                 ret_data[i + j * nx] = sum2 + tmp_f1 - tmp_f2;
01334                         }
01335                 }
01336         }
01337         else {
01338                 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
01339                         for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
01340                                 float eml = 0;
01341                                 float dot = 0;
01342                                 float dot2 = 0;
01343 
01344                                 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01345                                         eml += data[ii + (j - ny2 / 2) * nx] + data[ii + (j + ny2 / 2 - 1) * nx];
01346                                 }
01347 
01348                                 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01349                                         eml += data[i - nx2 / 2 + jj * nx] + data[i + nx2 / 2 - 1 + jj * nx];
01350                                 }
01351 
01352                                 eml /= (nx2 + ny2) * 2.0f;
01353                                 int k = 0;
01354 
01355                                 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01356                                         for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01357                                                 dot += (data[ii + jj * nx] - eml) * (with_data[k] - em);
01358                                                 dot2 += Util::square(data[ii + jj * nx] - eml);
01359                                                 k++;
01360                                         }
01361                                 }
01362 
01363                                 dot2 = std::sqrt(dot2);
01364 
01365                                 if (dot2 == 0) {
01366                                         ret_data[i + j * nx] = 0;
01367                                 }
01368                                 else {
01369                                         ret_data[i + j * nx] = dot / (nx2 * ny2 * dot2 * (float)with->get_attr("sigma"));
01370                                 }
01371                         }
01372                 }
01373         }
01374 
01375         ret->update();
01376 
01377         EXITFUNC;
01378         return ret;
01379 }
01380 
01381 
01382 EMData *EMData::do_radon()
01383 {
01384         ENTERFUNC;
01385 
01386         if (get_ndim() != 2) {
01387                 throw ImageDimensionException("2D only");
01388         }
01389 
01390         if (nx != ny) {
01391                 throw ImageFormatException("square image only");
01392         }
01393 
01394         EMData *result = new EMData();
01395         result->set_size(nx, ny, 1);
01396         result->to_zero();
01397         float *result_data = result->get_data();
01398 
01399         EMData *this_copy = this;
01400         this_copy = copy();
01401 
01402         for (int i = 0; i < nx; i++) {
01403                 Transform t(Dict("type","2d","alpha",(float) M_PI * 2.0f * i / nx));
01404                 this_copy->transform(t);
01405 
01406                 float *copy_data = this_copy->get_data();
01407 
01408                 for (int y = 0; y < nx; y++) {
01409                         for (int x = 0; x < nx; x++) {
01410                                 if (Util::square(x - nx / 2) + Util::square(y - nx / 2) <= nx * nx / 4) {
01411                                         result_data[i + y * nx] += copy_data[x + y * nx];
01412                                 }
01413                         }
01414                 }
01415 
01416                 this_copy->update();
01417         }
01418 
01419         result->update();
01420 
01421         if( this_copy )
01422         {
01423                 delete this_copy;
01424                 this_copy = 0;
01425         }
01426 
01427         EXITFUNC;
01428         return result;
01429 }
01430 
01431 void EMData::zero_corner_circulant(const int radius)
01432 {
01433         if ( nz > 1 && nz < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nz is too small");
01434         if ( ny > 1 && ny < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - ny is too small");
01435         if ( nx > 1 && nx < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nx is too small");
01436 
01437         int it_z = radius;
01438         int it_y = radius;
01439         int it_x = radius;
01440 
01441         if ( nz == 1 ) it_z = 0;
01442         if ( ny == 1 ) it_y = 0;
01443         if ( nx == 1 ) it_z = 0;
01444 
01445         if ( nz == 1 && ny == 1 )
01446         {
01447                 for ( int x = -it_x; x <= it_x; ++x )
01448                         get_value_at_wrap(x) = 0;
01449 
01450         }
01451         else if ( nz == 1 )
01452         {
01453                 for ( int y = -it_y; y <= it_y; ++y)
01454                         for ( int x = -it_x; x <= it_x; ++x )
01455                                 get_value_at_wrap(x,y) = 0;
01456         }
01457         else
01458         {
01459                 for( int z = -it_z; z <= it_z; ++z )
01460                         for ( int y = -it_y; y <= it_y; ++y)
01461                                 for ( int x = -it_x; x < it_x; ++x )
01462                                         get_value_at_wrap(x,y,z) = 0;
01463 
01464         }
01465 
01466 }
01467 
01468 EMData *EMData::calc_ccf(EMData * with, fp_flag fpflag,bool center)
01469 {
01470         ENTERFUNC;
01471 
01472         if( with == 0 ) {
01473 #ifdef EMAN2_USING_CUDA //CUDA 
01474         if(EMData::usecuda == 1 && cudarwdata) {
01475                 //cout << "calc ccf" << endl;
01476                 EMData* ifft = 0;
01477                 bool delifft = false;
01478                 int offset = 0;
01479                 
01480                 //do fft if not alreay done
01481                 if(!is_complex()){
01482                         ifft = do_fft_cuda();
01483                         delifft = true;
01484                         offset = 2 - nx%2;
01485                 }else{
01486                         ifft = this;
01487                 }
01488                 calc_conv_cuda(ifft->cudarwdata,ifft->cudarwdata,nx + offset, ny, nz); //this is the business end, results go in afft
01489                 
01490                 EMData * conv = ifft->do_ift_cuda();
01491                 if(delifft) delete ifft;
01492                 conv->update();
01493                         
01494                 return conv;
01495         }
01496 #endif
01497                 EXITFUNC;
01498                 return convolution(this,this,fpflag, center);
01499         }
01500         else if ( with == this ){ // this if statement is not necessary, the correlation function tests to see if with == this
01501                 EXITFUNC;
01502                 return correlation(this, this, fpflag,center);
01503         }
01504         else {
01505 
01506 #ifdef EMAN2_USING_CUDA //CUDA 
01507                 // assume always get rw data (makes life a lot easier!!! 
01508                 // also assume that both images are the same size. When using CUDA we are only interested in speed, not flexibility!!
01509                 // P.S. (I feel like I am pounding square pegs into a round holes with CUDA)
01510                 if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata) {
01511                         //cout << "using CUDA for ccf" << endl;
01512                         EMData* afft = 0;
01513                         EMData* bfft = 0;
01514                         bool delafft = false, delbfft = false;
01515                         int offset = 0;
01516                         
01517                         //do ffts if not alreay done
01518                         if(!is_complex()){
01519                                 afft = do_fft_cuda();
01520                                 delafft = true;
01521                                 offset = 2 - nx%2;
01522                                 //cout << "Do cuda FFT A" << endl;
01523                         }else{
01524                                 afft = this;
01525                         }
01526                         if(!with->is_complex()){
01527                                 bfft = with->do_fft_cuda();
01528                                 //cout << "Do cuda FFT B" << endl;
01529                                 delbfft = true;
01530                         }else{
01531                                 bfft = with;
01532                         }
01533 
01534                         calc_ccf_cuda(afft->cudarwdata,bfft->cudarwdata,nx + offset, ny, nz); //this is the business end, results go in afft
01535                         
01536                         if(delbfft) delete bfft;
01537                         
01538                         EMData * corr = afft->do_ift_cuda();
01539                         if(delafft) delete afft;
01540                         //cor->do_ift_inplace_cuda();//a bit faster, but I'll alos need to rearrnage the mem structure for it to work, BUT this is very SLOW.
01541                         corr->update();
01542                         
01543                         return corr;
01544                 }
01545 #endif
01546                 
01547                 // If the argument EMData pointer is not the same size we automatically resize it
01548                 bool undoresize = false;
01549                 int wnx = with->get_xsize(); int wny = with->get_ysize(); int wnz = with->get_zsize();
01550                 if (!(is_complex()^with->is_complex()) && (wnx != nx || wny != ny || wnz != nz) ) {
01551                         Region r((wnx-nx)/2, (wny-ny)/2, (wnz-nz)/2,nx,ny,nz);
01552                         with->clip_inplace(r);
01553                         undoresize = true;
01554                 }
01555 
01556                 EMData* cor = correlation(this, with, fpflag, center);
01557 
01558                 // If the argument EMData pointer was resized, it is returned to its original dimensions
01559                 if ( undoresize ) {
01560                         Region r((nx-wnx)/2, (ny-wny)/2,(nz-wnz)/2,wnx,wny,wnz);
01561                         with->clip_inplace(r);
01562                 }
01563 
01564                 EXITFUNC;
01565                 return cor;
01566         }
01567 }
01568 
01569 EMData *EMData::calc_ccfx( EMData * const with, int y0, int y1, bool no_sum, bool flip)
01570 {
01571         ENTERFUNC;
01572 
01573         if (!with) {
01574                 LOGERR("NULL 'with' image. ");
01575                 throw NullPointerException("NULL input image");
01576         }
01577 
01578         if (!EMUtil::is_same_size(this, with)) {
01579                 LOGERR("images not same size: (%d,%d,%d) != (%d,%d,%d)",
01580                            nx, ny, nz,
01581                            with->get_xsize(), with->get_ysize(), with->get_zsize());
01582                 throw ImageFormatException("images not same size");
01583         }
01584         if (get_ndim() > 2) {
01585                 LOGERR("2D images only");
01586                 throw ImageDimensionException("2D images only");
01587         }
01588 
01589         if (y1 <= y0) {
01590                 y1 = ny;
01591         }
01592 
01593         if (y0 >= y1) {
01594                 y0 = 0;
01595         }
01596 
01597         if (y0 < 0) {
01598                 y0 = 0;
01599         }
01600 
01601         if (y1 > ny) {
01602                 y1 = ny;
01603         }
01604         if (is_complex_x() || with->is_complex_x() ) throw; // Woops don't support this anymore!
01605 
01606         static int nx_fft = 0;
01607         static int ny_fft = 0;
01608         static EMData f1;
01609         static EMData f2;
01610         static EMData rslt;
01611 
01612         int height = y1-y0;
01613         int width = (nx+2-(nx%2));
01614         int wpad = ((width+3)/4)*4;                     // This is for 128 bit alignment of rows to prevent SSE crashes
01615         if (wpad != nx_fft || height != ny_fft ) {      // Seems meaningless, but due to static definitions above. f1,f2 are cached to prevent multiple reallocations
01616                 f1.set_size(wpad,height);
01617                 f2.set_size(wpad,height);
01618                 rslt.set_size(nx,height);
01619                 nx_fft = wpad;
01620                 ny_fft = height;
01621         }
01622 
01623 #ifdef EMAN2_USING_CUDA
01624         // FIXME : Not tested with new wpad change
01625         if (EMData::usecuda == 1 && cudarwdata && with->cudarwdata) {
01626                 //cout << "calc_ccfx CUDA" << endl;
01627                 if(!f1.cudarwdata) f1.rw_alloc();
01628                 if(!f2.cudarwdata) f2.rw_alloc();
01629                 if(!rslt.cudarwdata) rslt.rw_alloc();
01630                 cuda_dd_fft_real_to_complex_nd(cudarwdata, f1.cudarwdata, nx, 1, 1, height);
01631                 cuda_dd_fft_real_to_complex_nd(with->cudarwdata, f2.cudarwdata, nx, 1, 1, height);
01632                 calc_ccf_cuda(f1.cudarwdata, f2.cudarwdata, nx, ny, nz);
01633                 cuda_dd_fft_complex_to_real_nd(f1.cudarwdata, rslt.cudarwdata, nx, 1, 1, height);
01634                 if(no_sum){
01635                         EMData* result = new EMData(rslt);
01636                         return result;
01637                 }
01638                 EMData* cf = new EMData(0,0,nx,1,1); //cuda constructor
01639                 cf->runcuda(emdata_column_sum(rslt.cudarwdata, nx, ny));
01640                 cf->update();
01641                 
01642                 EXITFUNC;
01643                 return cf;
01644         }
01645 #endif
01646 
01647 //      printf("%d %d %d\n",(int)get_attr("nx"),(int)f2.get_attr("nx"),width);
01648 
01649         float *d1 = get_data();
01650         float *d2 = with->get_data();
01651         float *f1d = f1.get_data();
01652         float *f2d = f2.get_data();
01653         for (int j = 0; j < height; j++) {
01654                 EMfft::real_to_complex_1d(d1 + j * nx, f1d+j*wpad, nx);
01655                 EMfft::real_to_complex_1d(d2 + j * nx, f2d+j*wpad, nx);
01656         }
01657 
01658         if(flip == false) {
01659                 for (int j = 0; j < height; j++) {
01660                         float *f1a = f1d + j * wpad;
01661                         float *f2a = f2d + j * wpad;
01662 
01663                         for (int i = 0; i < width / 2; i++) {
01664                                 float re1 = f1a[2*i];
01665                                 float re2 = f2a[2*i];
01666                                 float im1 = f1a[2*i+1];
01667                                 float im2 = f2a[2*i+1];
01668 
01669                                 f1d[j*wpad+i*2] = re1 * re2 + im1 * im2;
01670                                 f1d[j*wpad+i*2+1] = im1 * re2 - re1 * im2;
01671                         }
01672                 }
01673         } else {
01674                 for (int j = 0; j < height; j++) {
01675                         float *f1a = f1d + j * wpad;
01676                         float *f2a = f2d + j * wpad;
01677 
01678                         for (int i = 0; i < width / 2; i++) {
01679                                 float re1 = f1a[2*i];
01680                                 float re2 = f2a[2*i];
01681                                 float im1 = f1a[2*i+1];
01682                                 float im2 = f2a[2*i+1];
01683 
01684                                 f1d[j*wpad+i*2] = re1 * re2 - im1 * im2;
01685                                 f1d[j*wpad+i*2+1] = im1 * re2 + re1 * im2;
01686                         }
01687                 }
01688         }
01689 
01690         float* rd = rslt.get_data();
01691         for (int j = y0; j < y1; j++) {
01692                 EMfft::complex_to_real_1d(f1d+j*wpad, rd+j*nx, nx);
01693         }
01694 
01695         if (no_sum) {
01696                 rslt.update(); // This is important in terms of the copy - the returned object won't have the correct flags unless we do this
01697                 EXITFUNC;
01698                 return new EMData(rslt);
01699         } else {
01700                 EMData *cf = new EMData(nx,1,1);
01701                 cf->to_zero();
01702                 float *c = cf->get_data();
01703                 for (int j = 0; j < height; j++) {
01704                         for(int i = 0; i < nx; ++i) {
01705                                 c[i] += rd[i+j*nx];
01706                         }
01707                 }
01708                 cf->update();
01709                 EXITFUNC;
01710                 return cf;
01711         }
01712 }
01713 
01714 EMData *EMData::make_rotational_footprint_cmc( bool unwrap) {
01715         ENTERFUNC;
01716         update_stat();
01717         // Note that rotational_footprint caching saves a large amount of time
01718         // but this is at the expense of memory. Note that a policy is hardcoded here,
01719         // that is that caching is only employed when premasked is false and unwrap
01720         // is true - this is probably going to be what is used in most scenarios
01721         // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
01722         // generated by e2speedtest.
01723         if ( rot_fp != 0 && unwrap == true) {
01724                 return new EMData(*rot_fp);
01725         }
01726 
01727         static EMData obj_filt;
01728         EMData* filt = &obj_filt;
01729         filt->set_complex(true);
01730 
01731 
01732         // The filter object is nothing more than a cached high pass filter
01733         // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
01734         // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
01735         // set to true, which is used for speed reasons.
01736         if (filt->get_xsize() != nx+2-(nx%2) || filt->get_ysize() != ny ||
01737                    filt->get_zsize() != nz ) {
01738                 filt->set_size(nx+2-(nx%2), ny, nz);
01739                 filt->to_one();
01740 
01741                 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01742         }
01743 
01744         EMData *ccf = this->calc_mutual_correlation(this, true,filt);
01745         ccf->sub(ccf->get_edge_mean());
01746         EMData *result = ccf->unwrap();
01747         delete ccf; ccf = 0;
01748 
01749         EXITFUNC;
01750         if ( unwrap == true)
01751         {
01752         // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
01753 
01754 // Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
01755 // to throw any exception
01756 // if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
01757 
01758 // Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
01759 // The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
01760                 rot_fp = result;
01761                 return new EMData(*rot_fp);
01762         }
01763         else return result;
01764 }
01765 
01766 EMData *EMData::make_rotational_footprint( bool unwrap) {
01767         ENTERFUNC;
01768         update_stat();
01769         // Note that rotational_footprint caching saves a large amount of time
01770         // but this is at the expense of memory. Note that a policy is hardcoded here,
01771         // that is that caching is only employed when premasked is false and unwrap
01772         // is true - this is probably going to be what is used in most scenarios
01773         // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
01774         // generated by e2speedtest.
01775         if ( rot_fp != 0 && unwrap == true) {
01776                 return new EMData(*rot_fp);
01777         }
01778 
01779         EMData* ccf = this->calc_ccf(this,CIRCULANT,true);
01780 //      EMData* ccf = this->calc_ccf(this,PADDED,true);         # this would probably be a bit better, but takes 4x longer  :^/
01781         ccf->sub(ccf->get_edge_mean());
01782         //ccf->process_inplace("xform.phaseorigin.tocenter"); ccf did the centering
01783         EMData *result = ccf->unwrap();
01784         delete ccf; ccf = 0;
01785 
01786         EXITFUNC;
01787         if ( unwrap == true)
01788         { // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
01789 
01790 // Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
01791 // to throw any exception
01792 // if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
01793 
01794 // Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
01795 // The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
01796                 rot_fp = result;
01797                 return new EMData(*rot_fp);
01798         }
01799         else return result;
01800 }
01801 
01802 EMData *EMData::make_rotational_footprint_e1( bool unwrap)
01803 {
01804         ENTERFUNC;
01805 
01806         update_stat();
01807         // Note that rotational_footprint caching saves a large amount of time
01808         // but this is at the expense of memory. Note that a policy is hardcoded here,
01809         // that is that caching is only employed when premasked is false and unwrap
01810         // is true - this is probably going to be what is used in most scenarios
01811         // as advised by Steve Ludtke - In terms of performance this caching doubles the metric
01812         // generated by e2speedtest.
01813         if ( rot_fp != 0 && unwrap == true) {
01814                 return new EMData(*rot_fp);
01815         }
01816 
01817         static EMData obj_filt;
01818         EMData* filt = &obj_filt;
01819         filt->set_complex(true);
01820 //      Region filt_region;
01821 
01822 //      if (nx & 1) {
01823 //              LOGERR("even image xsize only");                throw ImageFormatException("even image xsize only");
01824 //      }
01825 
01826         int cs = (((nx * 7 / 4) & 0xfffff8) - nx) / 2; // this pads the image to 1 3/4 * size with result divis. by 8
01827 
01828         static EMData big_clip;
01829         int big_x = nx+2*cs;
01830         int big_y = ny+2*cs;
01831         int big_z = 1;
01832         if ( nz != 1 ) {
01833                 big_z = nz+2*cs;
01834         }
01835 
01836 
01837         if ( big_clip.get_xsize() != big_x || big_clip.get_ysize() != big_y || big_clip.get_zsize() != big_z ) {
01838                 big_clip.set_size(big_x,big_y,big_z);
01839         }
01840         // It is important to set all newly established pixels around the boundaries to the mean
01841         // If this is not done then the associated rotational alignment routine breaks, in fact
01842         // everythin just goes foo.
01843 
01844         big_clip.to_value(get_edge_mean());
01845 
01846         if (nz != 1) {
01847                 big_clip.insert_clip(this,IntPoint(cs,cs,cs));
01848         } else  {
01849                 big_clip.insert_clip(this,IntPoint(cs,cs,0));
01850         }
01851         
01852         // The filter object is nothing more than a cached high pass filter
01853         // Ultimately it is used an argument to the EMData::mult(EMData,prevent_complex_multiplication (bool))
01854         // function in calc_mutual_correlation. Note that in the function the prevent_complex_multiplication
01855         // set to true, which is used for speed reasons.
01856         if (filt->get_xsize() != big_clip.get_xsize() +2-(big_clip.get_xsize()%2) || filt->get_ysize() != big_clip.get_ysize() ||
01857                    filt->get_zsize() != big_clip.get_zsize()) {
01858                 filt->set_size(big_clip.get_xsize() + 2-(big_clip.get_xsize()%2), big_clip.get_ysize(), big_clip.get_zsize());
01859                 filt->to_one();
01860                 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01861 #ifdef EMAN2_USING_CUDA
01862                 /*
01863                 if(EMData::usecuda == 1 && big_clip.cudarwdata)
01864                 {
01865                         filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
01866                 }
01867                 */
01868 #endif
01869         }
01870 #ifdef EMAN2_USING_CUDA
01871         /*
01872         if(EMData::usecuda == 1 && big_clip.cudarwdata && !filt->cudarwdata)
01873         {
01874                 filt->copy_to_cuda(); // since this occurs just once for many images, we don't pay much of a speed pentalty here, and we avoid the hassel of messing with sparx
01875         }
01876         */
01877 #endif
01878         
01879         EMData *mc = big_clip.calc_mutual_correlation(&big_clip, true,filt);
01880         mc->sub(mc->get_edge_mean());
01881 
01882         static EMData sml_clip;
01883         int sml_x = nx * 3 / 2;
01884         int sml_y = ny * 3 / 2;
01885         int sml_z = 1;
01886         if ( nz != 1 ) {
01887                 sml_z = nz * 3 / 2;
01888         }
01889 
01890         if ( sml_clip.get_xsize() != sml_x || sml_clip.get_ysize() != sml_y || sml_clip.get_zsize() != sml_z ) {
01891                 sml_clip.set_size(sml_x,sml_y,sml_z);   }
01892         if (nz != 1) {
01893                 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,-cs+nz/4));
01894         } else {
01895                 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,0));         // same as padding change above
01896         }
01897                 
01898         delete mc; mc = 0;
01899         EMData * result = NULL;
01900         
01901         if (nz == 1) {
01902                 if (!unwrap) {
01903 #ifdef EMAN2_USING_CUDA
01904                         //if(EMData::usecuda == 1 && sml_clip.cudarwdata) throw UnexpectedBehaviorException("shap masking is not yet supported by CUDA");
01905 #endif
01906                         result = sml_clip.process("mask.sharp", Dict("outer_radius", -1, "value", 0));
01907 
01908                 }
01909                 else {
01910                         result = sml_clip.unwrap();
01911                 }
01912         }
01913         else {
01914                 // I am not sure why there is any consideration of non 2D images, but it was here
01915                 // in the first port so I kept when I cleaned this function up (d.woolford)
01916 //              result = clipped_mc;
01917                 result = new EMData(sml_clip);
01918         }
01919         
01920 #ifdef EMAN2_USING_CUDA
01921         //if (EMData::usecuda == 1) sml_clip.roneedsanupdate(); //If we didn't do this then unwrap would use data from the previous call of this function, happens b/c sml_clip is static
01922 #endif
01923         EXITFUNC;
01924         if ( unwrap == true)
01925         { // this if statement reflects a strict policy of caching in only one scenario see comments at beginning of function block
01926 
01927                 // Note that the if statement at the beginning of this function ensures that rot_fp is not zero, so there is no need
01928                 // to throw any exception
01929                 if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
01930 
01931                 // Here is where the caching occurs - the rot_fp takes ownsherhip of the pointer, and a deep copied EMData object is returned.
01932                 // The deep copy invokes a cost in terms of CPU cycles and memory, but prevents the need for complicated memory management (reference counting)
01933                 rot_fp = result;
01934                 return new EMData(*rot_fp);
01935         }
01936         else return result;
01937 }
01938 
01939 EMData *EMData::make_footprint(int type)
01940 {
01941 //      printf("Make fp %d\n",type);
01942         if (type==0) {
01943                 EMData *un=make_rotational_footprint_e1(); // Use EMAN1's footprint strategy
01944                 if (un->get_ysize() <= 6) {
01945                         throw UnexpectedBehaviorException("In EMData::make_footprint. The rotational footprint is too small");
01946                 }
01947                 EMData *tmp=un->get_clip(Region(0,4,un->get_xsize(),un->get_ysize()-6));        // 4 and 6 are empirical
01948                 EMData *cx=tmp->calc_ccfx(tmp,0,-1,1);
01949                 EMData *fp=cx->get_clip(Region(0,0,cx->get_xsize()/2,cx->get_ysize()));
01950                 delete un;
01951                 delete tmp;
01952                 delete cx;
01953                 return fp;
01954         }
01955         else if (type==1 || type==2 ||type==5 || type==6) {
01956                 int i,j,kx,ky,lx,ly;
01957 
01958                 EMData *fft=do_fft();
01959 
01960                 // map for x,y -> radius for speed
01961                 int rmax=(get_xsize()+1)/2;
01962                 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
01963                 for (i=0; i<rmax; i++) {
01964                         for (j=0; j<rmax; j++) {
01965 #ifdef _WIN32
01966                                 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
01967 #else
01968                                 rmap[i+j*rmax]=hypot((float)i,(float)j);
01969 #endif  //_WIN32
01970 //                              printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
01971                         }
01972                 }
01973 
01974                 EMData *fp=new EMData(rmax*2+2,rmax*2,1);
01975                 fp->set_complex(1);
01976                 fp->to_zero();
01977 
01978                 // Two vectors in to complex space (kx,ky) and (lx,ly)
01979                 // We are computing the bispectrum, f(k).f(l).f*(k+l)
01980                 // but integrating out two dimensions, leaving |k|,|l|
01981                 for (kx=-rmax+1; kx<rmax; kx++) {
01982                         for (ky=-rmax+1; ky<rmax; ky++) {
01983                                 for (lx=-rmax+1; lx<rmax; lx++) {
01984                                         for (ly=-rmax+1; ly<rmax; ly++) {
01985                                                 int ax=kx+lx;
01986                                                 int ay=ky+ly;
01987                                                 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
01988                                                 int r1=(int)floor(.5+rmap[abs(kx)+rmax*abs(ky)]);
01989                                                 int r2=(int)floor(.5+rmap[abs(lx)+rmax*abs(ly)]);
01990 //                                              if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
01991 //                                              float r3=rmap[ax+rmax*ay];
01992                                                 if (r1+r2>=rmax) continue;
01993 
01994                                                 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
01995                                                 fp->set_value_at(r1*2,r2,p.real()+fp->get_value_at(r1*2,r2));           // We keep only the real component in anticipation of zero phase sum
01996 //                                              fp->set_value_at(r1*2,rmax*2-r2-1,  fp->get_value_at(r1*2,r2));         // We keep only the real component in anticipation of zero phase sum
01997 //                                              fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2));               // We keep only the real component in anticipation of zero phase sum
01998                                                 fp->set_value_at(r1*2+1,r2,fp->get_value_at(r1*2+1,r2)+1);                      // a normalization counter
01999                                         }
02000                                 }
02001                         }
02002                 }
02003 
02004                 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
02005                 if (type==5 || type==6) {
02006                         for (i=0; i<rmax*2; i+=2) {
02007                                 for (j=0; j<rmax; j++) {
02008                                         float norm=fp->get_value_at(i+1,j);
02009 #ifdef _WIN32
02010                                         fp->set_value_at(i,rmax*2-j-1,pow(fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm), 1.0f/3.0f));
02011                                         fp->set_value_at(i,j,pow(fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm), 1.0f/3.0f));
02012 #else
02013                                         fp->set_value_at(i,rmax*2-j-1,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
02014                                         fp->set_value_at(i,j,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
02015 #endif  //_WIN32
02016                                         fp->set_value_at(i+1,j,0.0);
02017                                 }
02018                         }
02019                 }
02020                 else {
02021                         for (i=0; i<rmax*2; i+=2) {
02022                                 for (j=0; j<rmax; j++) {
02023                                         float norm=fp->get_value_at(i+1,j);
02024                                         fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
02025                                         fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
02026                                         fp->set_value_at(i+1,j,0.0);
02027                                 }
02028                         }
02029                 }
02030 
02031                 free(rmap);
02032                 if (type==2||type==6) {
02033                         EMData *f2=fp->do_ift();
02034                         if (f2->get_value_at(0,0)<0) f2->mult(-1.0f);
02035                         f2->process_inplace("xform.phaseorigin.tocorner");
02036                         delete fp;
02037                         return f2;
02038                 }
02039                 return fp;
02040         }
02041         else if (type==3 || type==4) {
02042                 int h,i,j,kx,ky,lx,ly;
02043 
02044                 EMData *fft=do_fft();
02045 
02046                 // map for x,y -> radius for speed
02047                 int rmax=(get_xsize()+1)/2;
02048                 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
02049                 for (i=0; i<rmax; i++) {
02050                         for (j=0; j<rmax; j++) {
02051 #ifdef _WIN32
02052                                 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
02053 #else
02054                                 rmap[i+j*rmax]=hypot((float)i,(float)j);
02055 #endif  //_WIN32
02056 //                              printf("%d\t%d\t%f\n",i,j,rmap[i+j*rmax]);
02057                         }
02058                 }
02059 
02060                 EMData *fp=new EMData(rmax*2+2,rmax*2,16);
02061 
02062                 fp->set_complex(1);
02063                 fp->to_zero();
02064 
02065                 // Two vectors in to complex space (kx,ky) and (lx,ly)
02066                 // We are computing the bispectrum, f(k).f(l).f*(k+l)
02067                 // but integrating out two dimensions, leaving |k|,|l|
02068                 for (kx=-rmax+1; kx<rmax; kx++) {
02069                         for (ky=-rmax+1; ky<rmax; ky++) {
02070                                 for (lx=-rmax+1; lx<rmax; lx++) {
02071                                         for (ly=-rmax+1; ly<rmax; ly++) {
02072                                                 int ax=kx+lx;
02073                                                 int ay=ky+ly;
02074                                                 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
02075                                                 float rr1=rmap[abs(kx)+rmax*abs(ky)];
02076                                                 float rr2=rmap[abs(lx)+rmax*abs(ly)];
02077                                                 int r1=(int)floor(.5+rr1);
02078                                                 int r2=(int)floor(.5+rr2);
02079 //                                              if (r1>500 ||r2>500) printf("%d\t%d\t%d\t%d\t%d\t%d\n",kx,ky,lx,ly,r1,r2);
02080 //                                              float r3=rmap[ax+rmax*ay];
02081                                                 if (r1+r2>=rmax || rr1==0 ||rr2==0) continue;
02082 
02083                                                 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
02084                                                 int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5);                                        // projection of k on l 0-31
02085                                                 if (dot<0) dot=16+dot;
02086 //                                              int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5+8.0);                                    // projection of k on l 0-15
02087                                                 fp->set_value_at(r1*2,r2,dot,p.real()+fp->get_value_at(r1*2,r2,dot));           // We keep only the real component in anticipation of zero phase sum
02088 //                                              fp->set_value_at(r1*2,rmax*2-r2-1,  fp->get_value_at(r1*2,r2));         // We keep only the real component in anticipation of zero phase sum
02089 //                                              fp->set_value_at(r1*2+1,r2,p.real()+fp->get_value_at(r1*2+1,r2));               // We keep only the real component in anticipation of zero phase sum
02090                                                 fp->set_value_at(r1*2+1,r2,dot,fp->get_value_at(r1*2+1,r2,dot)+1);                      // a normalization counter
02091                                         }
02092                                 }
02093                         }
02094                 }
02095 
02096                 // Normalizes the pixels based on the accumulated counts then sets the imaginary components back to zero
02097                 for (i=0; i<rmax*2; i+=2) {
02098                         for (j=0; j<rmax; j++) {
02099                                 for (h=0; h<16; h++) {
02100                                         float norm=fp->get_value_at(i+1,j,h);
02101 //                                      fp->set_value_at(i,rmax*2-j-1,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
02102 //                                      fp->set_value_at(i,j,h,cbrt(fp->get_value_at(i,j,h)/(norm==0?1.0:norm)));
02103                                         fp->set_value_at(i,rmax*2-j-1,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02104                                         fp->set_value_at(i,j,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02105         //                              fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0?1.0:norm));
02106         //                              fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0?1.0:norm));
02107                                         fp->set_value_at(i+1,j,h,0.0);
02108                                 }
02109                         }
02110                 }
02111 
02112                 free(rmap);
02113                 if (type==4) {
02114                         EMData *f2=fp->do_ift();
02115                         if (f2->get_value_at(0,0,0)<0) f2->mult(-1.0f);
02116                         f2->process_inplace("xform.phaseorigin.tocorner");
02117                         delete fp;
02118                         return f2;
02119                 }
02120                 return fp;
02121         }
02122         throw UnexpectedBehaviorException("There is not implementation for the parameters you specified");
02123 }
02124 
02125 
02126 EMData *EMData::calc_mutual_correlation(EMData * with, bool tocenter, EMData * filter)
02127 {
02128         ENTERFUNC;
02129 
02130         if (with && !EMUtil::is_same_size(this, with)) {
02131                 LOGERR("images not same size");
02132                 throw ImageFormatException( "images not same size");
02133         }
02134 
02135 #ifdef EMAN2_USING_CUDA
02136         if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata)
02137         {       
02138 
02139                 EMData* this_fft = do_fft_cuda();
02140 
02141                 EMData *cf = 0;
02142                 if (with && with != this) {
02143                         cf = with->do_fft_cuda();
02144                 }else{
02145                         cf = this_fft->copy();
02146                 }
02147                 
02148                 if (filter) {
02149                         if (!EMUtil::is_same_size(filter, cf)) {
02150                                 LOGERR("improperly sized filter");
02151                                 throw ImageFormatException("improperly sized filter");
02152                         }
02153                         mult_complex_efficient_cuda(cf->cudarwdata, filter->cudarwdata, cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), 1);
02154                         mult_complex_efficient_cuda(this_fft->cudarwdata, filter->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize(), 1);
02155                 }
02156                 
02157                 mcf_cuda(this_fft->cudarwdata, cf->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize());
02158                 
02159                 EMData *f2 = cf->do_ift_cuda();
02160 
02161                 if (tocenter) {
02162                         f2->process_inplace("xform.phaseorigin.tocenter");
02163                 }
02164 
02165                 if( cf )
02166                 {
02167                         delete cf;
02168                         cf = 0;
02169                 }
02170 
02171                 if( this_fft )
02172                 {
02173                         delete this_fft;
02174                         this_fft = 0;
02175                 }
02176 
02177                 f2->set_attr("label", "MCF");
02178                 f2->set_path("/tmp/eman.mcf");
02179                 f2->update();
02180 
02181                 EXITFUNC;
02182                 return f2;
02183         }
02184 #endif
02185 
02186         EMData *this_fft = 0;
02187         this_fft = do_fft();
02188 
02189         if (!this_fft) {
02190 
02191                 LOGERR("FFT returns NULL image");
02192                 throw NullPointerException("FFT returns NULL image");
02193         }
02194 
02195         this_fft->ap2ri(); //this is not needed!
02196         EMData *cf = 0;
02197 
02198         if (with && with != this) {
02199                 cf = with->do_fft();
02200                 if (!cf) {
02201                         LOGERR("FFT returns NULL image");
02202                         throw NullPointerException("FFT returns NULL image");
02203                 }
02204                 cf->ap2ri(); //nor is this!
02205         }
02206         else {
02207                 cf = this_fft->copy();
02208         }
02209         
02210         if (filter) {
02211                 if (!EMUtil::is_same_size(filter, cf)) {
02212                         LOGERR("improperly sized filter");
02213                         throw ImageFormatException("improperly sized filter");
02214                 }
02215                 
02216                 cf->mult_complex_efficient(*filter,true); //insanely this is required....
02217                 this_fft->mult(*filter,true);
02218                 //cf->mult_complex_efficient(*filter,7); // takes advantage of the fact that the filter is 1 everywhere but near the origin
02219                 //this_fft->mult_complex_efficient(*filter,7);
02220                 /*cf->mult_complex_efficient(*filter,5);
02221                 this_fft->mult_complex_efficient(*filter,5);*/
02222         }
02223 
02224         float *rdata1 = this_fft->get_data();
02225         float *rdata2 = cf->get_data();
02226         size_t this_fft_size = (size_t)this_fft->get_xsize() * this_fft->get_ysize() * this_fft->get_zsize();
02227 
02228         if (with == this) {
02229                 for (size_t i = 0; i < this_fft_size; i += 2) {
02230                         rdata2[i] = std::sqrt(rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02231                         rdata2[i + 1] = 0;
02232                 }
02233 
02234                 this_fft->update();
02235                 cf->update();
02236         }
02237         else {
02238                 for (size_t i = 0; i < this_fft_size; i += 2) {
02239                         rdata2[i] = (rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02240                         rdata2[i + 1] = (rdata1[i + 1] * rdata2[i] - rdata1[i] * rdata2[i + 1]);
02241                 }
02242                 
02243                 //This seems like a bug, but it probably is never used....
02244                 for (size_t i = 0; i < this_fft_size; i += 2) {
02245                         float t = Util::square(rdata2[i]) + Util::square(rdata2[i + 1]);
02246                         if (t != 0) {
02247                                 t = pow(t, 0.25f);
02248                                 rdata2[i] /= t;
02249                                 rdata2[i + 1] /= t;
02250                         }
02251                 }
02252                 this_fft->update();
02253                 cf->update();
02254         }
02255 
02256         EMData *f2 = cf->do_ift();
02257 
02258         if (tocenter) {
02259                 f2->process_inplace("xform.phaseorigin.tocenter");
02260         }
02261 
02262         if( cf )
02263         {
02264                 delete cf;
02265                 cf = 0;
02266         }
02267 
02268         if( this_fft )
02269         {
02270                 delete this_fft;
02271                 this_fft = 0;
02272         }
02273 
02274         f2->set_attr("label", "MCF");
02275         f2->set_path("/tmp/eman.mcf");
02276 
02277         EXITFUNC;
02278         return f2;
02279 }
02280 
02281 
02282 vector < float > EMData::calc_hist(int hist_size, float histmin, float histmax,const float& brt, const float& cont)
02283 {
02284         ENTERFUNC;
02285 
02286         static size_t prime[] = { 1, 3, 7, 11, 17, 23, 37, 59, 127, 253, 511 };
02287 
02288         if (histmin == histmax) {
02289                 histmin = get_attr("minimum");
02290                 histmax = get_attr("maximum");
02291         }
02292 
02293         vector <float> hist(hist_size, 0.0);
02294 
02295         int p0 = 0;
02296         int p1 = 0;
02297         size_t size = (size_t)nx * ny * nz;
02298         if (size < 300000) {
02299                 p0 = 0;
02300                 p1 = 0;
02301         }
02302         else if (size < 2000000) {
02303                 p0 = 2;
02304                 p1 = 3;
02305         }
02306         else if (size < 8000000) {
02307                 p0 = 4;
02308                 p1 = 6;
02309         }
02310         else {
02311                 p0 = 7;
02312                 p1 = 9;
02313         }
02314 
02315         if (is_complex() && p0 > 0) {
02316                 p0++;
02317                 p1++;
02318         }
02319 
02320         size_t di = 0;
02321 //      float norm = 0;
02322         size_t n = hist.size();
02323 
02324         float * data = get_data();
02325         for (int k = p0; k <= p1; ++k) {
02326                 if (is_complex()) {
02327                         di = prime[k] * 2;
02328                 }
02329                 else {
02330                         di = prime[k];
02331                 }
02332 
02333 //              norm += (float)size / (float) di;
02334                 float w = (float)n / (histmax - histmin);
02335 
02336                 for(size_t i=0; i<=size-di; i += di) {
02337                         float val;
02338                         if (cont != 1.0f || brt != 0)val = cont*(data[i]+brt);
02339                         else val = data[i];
02340                         int j = Util::round((val - histmin) * w);
02341                         if (j >= 0 && j < (int) n) {
02342                                 hist[j] += 1;
02343                         }
02344                 }
02345         }
02346 /*
02347         for (size_t i = 0; i < hist.size(); ++i) {
02348                 if (norm != 0) {
02349                         hist[i] = hist[i] / norm;
02350                 }
02351         }
02352 */
02353         return hist;
02354 
02355         EXITFUNC;
02356 }
02357 
02358 
02359 
02360 
02361 
02362 vector<float> EMData::calc_az_dist(int n, float a0, float da, float rmin, float rmax)
02363 {
02364         ENTERFUNC;
02365 
02366         a0=a0*M_PI/180.0f;
02367         da=da*M_PI/180.0f;
02368 
02369         if (get_ndim() > 2) {
02370                 throw ImageDimensionException("no 3D image");
02371         }
02372 
02373         float *yc = new float[n];
02374 
02375         vector<float>   vd(n);
02376         for (int i = 0; i < n; i++) {
02377                 yc[i] = 0.00001f;
02378         }
02379 
02380         float * data = get_data();
02381         if (is_complex()) {
02382                 int c = 0;
02383                 for (int y = 0; y < ny; y++) {
02384                         for (int x = 0; x < nx; x += 2, c += 2) {
02385                                 int x1 = x / 2;
02386                                 int y1 = y<ny/2?y:y-ny;
02387                                 float r = (float)Util::hypot_fast(x1,y1);
02388 
02389                                 if (r >= rmin && r <= rmax) {
02390                                         float a = 0;
02391 
02392                                         if (y != ny / 2 || x != 0) {
02393                                                 a = (atan2((float)y1, (float)x1) - a0) / da;
02394                                         }
02395 
02396                                         int i = (int)(floor(a));
02397                                         a -= i;
02398 
02399                                         if (i == 0) {
02400                                                 vd[0] += data[c] * (1.0f - a);
02401                                                 yc[0] += (1.0f - a);
02402                                         }
02403                                         else if (i == n - 1) {
02404                                                 vd[n - 1] += data[c] * a;
02405                                                 yc[n - 1] += a;
02406                                         }
02407                                         else if (i > 0 && i < (n - 1)) {
02408                                                 float h = 0;
02409                                                 if (is_ri()) {
02410 #ifdef  _WIN32
02411                                                         h = (float)_hypot(data[c], data[c + 1]);
02412 #else
02413                                                         h = (float)hypot(data[c], data[c + 1]);
02414 #endif  //_WIN32
02415                                                 }
02416                                                 else {
02417                                                         h = data[c];
02418                                                 }
02419 
02420                                                 vd[i] += h * (1.0f - a);
02421                                                 yc[i] += (1.0f - a);
02422                                                 vd[i + 1] += h * a;
02423                                                 yc[i + 1] += a;
02424                                         }
02425                                 }
02426                         }
02427                 }
02428         }
02429         else {
02430                 int c = 0;
02431                 float half_nx = (nx - 1) / 2.0f;
02432                 float half_ny = (ny - 1) / 2.0f;
02433 
02434                 for (int y = 0; y < ny; y++) {
02435                         for (int x = 0; x < nx; x++, c++) {
02436                                 float y1 = y - half_ny;
02437                                 float x1 = x - half_nx;
02438 #ifdef  _WIN32
02439                                 float r = (float)_hypot(x1, y1);
02440 #else
02441                                 float r = (float)hypot(x1, y1);
02442 #endif
02443 
02444                                 if (r >= rmin && r <= rmax) {
02445                                         float a = 0;
02446                                         if (x1 != 0 || y1 != 0) {
02447                                                 a = atan2(y1, x1);
02448                                                 if (a < 0) {
02449                                                         a += M_PI * 2;
02450                                                 }
02451                                         }
02452 
02453                                         a = (a - a0) / da;
02454                                         int i = static_cast < int >(floor(a));
02455                                         a -= i;
02456 
02457                                         if (i == 0) {
02458                                                 vd[0] += data[c] * (1.0f - a);
02459                                                 yc[0] += (1.0f - a);
02460                                         }
02461                                         else if (i == n - 1) {
02462                                                 vd[n - 1] += data[c] * a;
02463                                                 yc[n - 1] += (a);
02464                                         }
02465                                         else if (i > 0 && i < (n - 1)) {
02466                                                 vd[i] += data[c] * (1.0f - a);
02467                                                 yc[i] += (1.0f - a);
02468                                                 vd[i + 1] += data[c] * a;
02469                                                 yc[i + 1] += a;
02470                                         }
02471                                 }
02472                         }
02473                 }
02474         }
02475 
02476 
02477         for (int i = 0; i < n; i++) {
02478                 vd[i] /= yc[i];
02479         }
02480 
02481         if( yc )
02482         {
02483                 delete[]yc;
02484                 yc = 0;
02485         }
02486 
02487         return vd;
02488 
02489         EXITFUNC;
02490 }
02491 
02492 
02493 EMData *EMData::unwrap(int r1, int r2, int xs, int dx, int dy, bool do360, bool weight_radial) const
02494 {
02495         ENTERFUNC;
02496 
02497         if (get_ndim() != 2) {
02498                 throw ImageDimensionException("2D image only");
02499         }
02500 
02501         int p = 1;
02502         if (do360) {
02503                 p = 2;
02504         }
02505 
02506         if (xs < 1) {
02507                 xs = (int) Util::fast_floor(p * M_PI * ny / 5.0);
02508                 xs = Util::calc_best_fft_size(xs);
02509                 if (xs<=8) xs=16;
02510         }
02511 
02512         if (r1 < 0) {
02513                 r1 = 4;
02514         }
02515 
02516 #ifdef  _WIN32
02517         int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(_hypot(dx, dy)));
02518 #else
02519         int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(hypot(dx, dy)));
02520 #endif  //_WIN32
02521         rr-=rr%2;
02522         if (r2 <= r1 || r2 > rr) {
02523                 r2 = rr;
02524         }
02525 
02526         if ( (r2-r1) < 0 ) throw UnexpectedBehaviorException("The combination of function the arguments and the image dimensions causes unexpected behavior internally. Use a larger image, or a smaller value of r1, or a combination of both");
02527 
02528 #ifdef EMAN2_USING_CUDA
02529         if (EMData::usecuda == 1 && isrodataongpu()){
02530                 //cout << " CUDA unwrap" << endl;
02531                 EMData* result = new EMData(0,0,xs,(r2-r1),1);
02532                 if(!result->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
02533                 bindcudaarrayA(true);
02534                 emdata_unwrap(result->cudarwdata, r1, r2, xs, p, dx, dy, weight_radial, nx, ny);
02535                 unbindcudaarryA();
02536                 result->update();
02537                 return result;
02538         }
02539 #endif
02540 
02541         EMData *ret = new EMData();
02542         ret->set_size(xs, r2 - r1, 1);
02543         const float *const d = get_const_data();
02544         float *dd = ret->get_data();
02545         float pfac = (float)p/(float)xs;
02546 
02547         int nxon2 = nx/2;
02548         int nyon2 = ny/2;
02549         for (int x = 0; x < xs; x++) {
02550                 float ang = x * M_PI * pfac;
02551                 float si = sin(ang);
02552                 float co = cos(ang);
02553 
02554                 for (int y = 0; y < r2 - r1; y++) {
02555                         float ypr1 = (float)y + r1;
02556                         float xx = ypr1 * co + nxon2 + dx;
02557                         float yy = ypr1 * si + nyon2 + dy;
02558 //                      float t = xx - Util::fast_floor(xx);
02559 //                      float u = yy - Util::fast_floor(yy);
02560                         float t = xx - (int)xx;
02561                         float u = yy - (int)yy;
02562 //                      int k = (int) Util::fast_floor(xx) + (int) (Util::fast_floor(yy)) * nx;
02563                         int k = (int) xx + ((int) yy) * nx;
02564                         float val = Util::bilinear_interpolate(d[k], d[k + 1], d[k + nx], d[k + nx+1], t,u);
02565                         if (weight_radial) val *=  ypr1;
02566                         dd[x + y * xs] = val;
02567                 }
02568 
02569         }
02570         ret->update();
02571 
02572         EXITFUNC;
02573         return ret;
02574 }
02575 
02576 // NOTE : x axis is from 0 to 0.5  (Nyquist), and thus properly handles non-square images
02577 // complex only
02578 void EMData::apply_radial_func(float x0, float step, vector < float >array, bool interp)
02579 {
02580         ENTERFUNC;
02581 
02582         if (!is_complex()) throw ImageFormatException("apply_radial_func requires a complex image");
02583 
02584         int n = static_cast < int >(array.size());
02585 
02586         if (n*step>2.0) printf("Warning, apply_radial_func takes x0 and step with respect to Nyquist (0.5)\n");
02587 
02588 //      printf("%f %f %f\n",array[0],array[25],array[50]);
02589 
02590         ap2ri();
02591 
02592         size_t ndims = get_ndim();
02593         float * data = get_data();
02594         if (ndims == 2) {
02595                 int k = 0;
02596                 for (int j = 0; j < ny; j++) {
02597                         for (int i = 0; i < nx; i += 2, k += 2) {
02598                                 float r;
02599 #ifdef  _WIN32
02600                                 if (j<ny/2) r = (float)_hypot(i/(float)(nx*2), j/(float)ny);
02601                                 else r = (float)_hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02602 #else
02603                                 if (j<ny/2) r = (float)hypot(i/(float)(nx*2), j/(float)ny);
02604                                 else r = (float)hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02605 #endif  //_WIN32
02606                                 r = (r - x0) / step;
02607 
02608                                 int l = 0;
02609                                 if (interp) {
02610                                         l = (int) floor(r);
02611                                 }
02612                                 else {
02613                                         l = (int) floor(r + 1);
02614                                 }
02615 
02616 
02617                                 float f = 0;
02618                                 if (l >= n - 2) {
02619                                         f = array[n - 1];
02620                                 }
02621                                 else {
02622                                         if (interp) {
02623                                                 r -= l;
02624                                                 f = (array[l] * (1.0f - r) + array[l + 1] * r);
02625                                         }
02626                                         else {
02627                                                 f = array[l];
02628                                         }
02629                                 }
02630 
02631                                 data[k] *= f;
02632                                 data[k + 1] *= f;
02633                         }
02634                 }
02635         }
02636         else if (ndims == 3) {
02637                 int k = 0;
02638                 for (int m = 0; m < nz; m++) {
02639                         float mnz;
02640                         if (m<nz/2) mnz=m*m/(float)(nz*nz);
02641                         else { mnz=(nz-m)/(float)nz; mnz*=mnz; }
02642 
02643                         for (int j = 0; j < ny; j++) {
02644                                 float jny;
02645                                 if (j<ny/2) jny= j*j/(float)(ny*ny);
02646                                 else { jny=(ny-j)/(float)ny; jny*=jny; }
02647 
02648                                 for (int i = 0; i < nx; i += 2, k += 2) {
02649                                         float r = std::sqrt((i * i / (nx*nx*4.0f)) + jny + mnz);
02650                                         r = (r - x0) / step;
02651 
02652                                         int l = 0;
02653                                         if (interp) {
02654                                                 l = (int) floor(r);
02655                                         }
02656                                         else {
02657                                                 l = (int) floor(r + 1);
02658                                         }
02659 
02660 
02661                                         float f = 0;
02662                                         if (l >= n - 2) {
02663                                                 f = array[n - 1];
02664                                         }
02665                                         else {
02666                                                 if (interp) {
02667                                                         r -= l;
02668                                                         f = (array[l] * (1.0f - r) + array[l + 1] * r);
02669                                                 }
02670                                                 else {
02671                                                         f = array[l];
02672                                                 }
02673                                         }
02674 
02675                                         data[k] *= f;
02676                                         data[k + 1] *= f;
02677                                 }
02678                         }
02679                 }
02680 
02681         }
02682 
02683         update();
02684         EXITFUNC;
02685 }
02686 
02687 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, bool inten)
02688 {
02689         ENTERFUNC;
02690 
02691         vector<float>ret(n);
02692         vector<float>norm(n);
02693 
02694         int x,y,z,i;
02695         int step=is_complex()?2:1;
02696         int isinten=get_attr_default("is_intensity",0);
02697 
02698         if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
02699 
02700         for (i=0; i<n; i++) ret[i]=norm[i]=0.0;
02701         float * data = get_data();
02702 
02703         // We do 2D separately to avoid the hypot3 call
02704         if (nz==1) {
02705                 for (y=i=0; y<ny; y++) {
02706                         for (x=0; x<nx; x+=step,i+=step) {
02707                                 float r,v;
02708                                 if (step==2) {          //complex
02709                                         if (x==0 && y>ny/2) continue;
02710                                         r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y));         // origin at 0,0; periodic
02711                                         if (!inten) {
02712 #ifdef  _WIN32
02713                                                 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));   // real/imag, compute amplitude
02714 #else
02715                                                 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));    // real/imag, compute amplitude
02716 #endif
02717                                                 else v=data[i];                                                 // amp/phase, just get amp
02718                                         } else {
02719                                                 if (isinten) v=data[i];
02720                                                 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02721                                                 else v=data[i]*data[i];
02722                                         }
02723                                 }
02724                                 else {
02725                                         r=(float)(Util::hypot_fast(x-nx/2,y-ny/2));
02726                                         if (inten) v=data[i]*data[i];
02727                                         else v=data[i];
02728                                 }
02729                                 r=(r-x0)/dx;
02730                                 int f=int(r);   // safe truncation, so floor isn't needed
02731                                 r-=float(f);    // r is now the fractional spacing between bins
02732 //                              printf("%d\t%d\t%d\t%1.3f\t%d\t%1.3f\t%1.4g\n",x,y,f,r,step,Util::hypot_fast(x/2,y<ny/2?y:ny-y),v);
02733                                 if (f>=0 && f<n) {
02734                                         ret[f]+=v*(1.0f-r);
02735                                         norm[f]+=(1.0f-r);
02736                                         if (f<n-1) {
02737                                                 ret[f+1]+=v*r;
02738                                                 norm[f+1]+=r;
02739                                         }
02740                                 }
02741                         }
02742                 }
02743         }
02744         else {
02745                 size_t i;       //3D file may have >2G size
02746                 for (z=i=0; z<nz; ++z) {
02747                         for (y=0; y<ny; ++y) {
02748                                 for (x=0; x<nx; x+=step,i+=step) {
02749                                         float r,v;
02750                                         if (step==2) {  //complex
02751                                                 if (x==0 && z<nz/2) continue;
02752                                                 if (x==0 && z==nz/2 && y<ny/2) continue;
02753                                                 r=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z);        // origin at 0,0; periodic
02754                                                 if (!inten) {
02755 #ifdef  _WIN32
02756                                                         if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));   // real/imag, compute amplitude
02757 #else
02758                                                         if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));    // real/imag, compute amplitude
02759 #endif  //_WIN32
02760                                                         else v=data[i];                                                 // amp/phase, just get amp
02761                                                 } else {
02762                                                         if (isinten) v=data[i];
02763                                                         else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02764                                                         else v=data[i]*data[i];
02765                                                 }
02766                                         }
02767                                         else {
02768                                                 r=Util::hypot3(x-nx/2,y-ny/2,z-nz/2);
02769                                                 if (inten) v=data[i]*data[i];
02770                                                 else v=data[i];
02771                                         }
02772                                         r=(r-x0)/dx;
02773                                         int f=int(r);   // safe truncation, so floor isn't needed
02774                                         r-=float(f);    // r is now the fractional spacing between bins
02775                                         if (f>=0 && f<n) {
02776                                                 ret[f]+=v*(1.0f-r);
02777                                                 norm[f]+=(1.0f-r);
02778                                                 if (f<n-1) {
02779                                                         ret[f+1]+=v*r;
02780                                                         norm[f+1]+=r;
02781                                                 }
02782                                         }
02783                                 }
02784                         }
02785                 }
02786         }
02787 
02788         for (i=0; i<n; i++) ret[i]/=norm[i]?norm[i]:1.0f;       // Normalize
02789 
02790         EXITFUNC;
02791 
02792         return ret;
02793 }
02794 
02795 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, int nwedge, float offset, bool inten)
02796 {
02797         ENTERFUNC;
02798 
02799         if (nz > 1) {
02800                 LOGERR("2D images only.");
02801                 throw ImageDimensionException("2D images only");
02802         }
02803         int isinten=get_attr_default("is_intensity",0);
02804 
02805         if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
02806 
02807 
02808         vector<float>ret(n*nwedge);
02809         vector<float>norm(n*nwedge);
02810 
02811         int x,y,i;
02812         int step=is_complex()?2:1;
02813         float astep=static_cast<float>(M_PI*2.0/nwedge);
02814         if (is_complex()) astep/=2;                                                     // Since we only have the right 1/2 of Fourier space
02815         float* data = get_data();
02816         for (i=0; i<n*nwedge; i++) ret[i]=norm[i]=0.0;
02817 
02818         // We do 2D separately to avoid the hypot3 call
02819         for (y=i=0; y<ny; y++) {
02820                 for (x=0; x<nx; x+=step,i+=step) {
02821                         float r,v,a;
02822                         int bin;
02823                         if (is_complex()) {
02824 #ifdef  _WIN32
02825                                 r=static_cast<float>(_hypot(x/2.0,y<ny/2?y:ny-y));              // origin at 0,0; periodic
02826 #else
02827                                 r=static_cast<float>(hypot(x/2.0,y<ny/2?y:ny-y));               // origin at 0,0; periodic
02828 #endif
02829                                 a=atan2(float(y<ny/2?y:y-ny),x/2.0f);
02830                                 if (!inten) {
02831 #ifdef  _WIN32
02832                                         if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));   // real/imag, compute amplitude
02833 #else
02834                                         if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));    // real/imag, compute amplitude
02835 #endif  //_WIN32
02836                                         else v=data[i];                                                 // amp/phase, just get amp
02837                                 } else {
02838                                         if (isinten) v=data[i];
02839                                         else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02840                                         else v=data[i]*data[i];
02841                                 }
02842                                 bin=n*int(floor((a+M_PI/2.0f+offset)/astep));
02843                         }
02844                         else {
02845 #ifdef  _WIN32
02846                                 r=static_cast<float>(_hypot(x-nx/2,y-ny/2));
02847 #else
02848                                 r=static_cast<float>(hypot(x-nx/2,y-ny/2));
02849 #endif  //_WIN32
02850                                 a=atan2(float(y-ny/2),float(x-nx/2));
02851                                 if (inten) v=data[i]*data[i];
02852                                 else v=data[i];
02853                                 bin=n*int(floor((a+M_PI+offset)/astep));
02854                         }
02855                         if (bin>=nwedge*n) bin-=nwedge*n;
02856                         if (bin<0) bin+=nwedge*n;
02857                         r=(r-x0)/dx;
02858                         int f=int(r);   // safe truncation, so floor isn't needed
02859                         r-=float(f);    // r is now the fractional spacing between bins
02860 //                      printf("%d %d %d %d %1.3f %1.3f\n",x,y,bin,f,r,a);
02861                         if (f>=0 && f<n) {
02862                                 ret[f+bin]+=v*(1.0f-r);
02863                                 norm[f+bin]+=(1.0f-r);
02864                                 if (f<n-1) {
02865                                         ret[f+1+bin]+=v*r;
02866                                         norm[f+1+bin]+=r;
02867                                 }
02868                         }
02869                 }
02870         }
02871 
02872         for (i=0; i<n*nwedge; i++) ret[i]/=norm[i]?norm[i]:1.0f;        // Normalize
02873         EXITFUNC;
02874 
02875         return ret;
02876 }
02877 
02878 void EMData::cconj() {
02879         ENTERFUNC;
02880         if (!is_complex() || !is_ri())
02881                 throw ImageFormatException("EMData::conj requires a complex, ri image");
02882         int nxreal = nx -2 + int(is_fftodd());
02883         int nxhalf = nxreal/2;
02884         for (int iz = 0; iz < nz; iz++)
02885                 for (int iy = 0; iy < ny; iy++)
02886                         for (int ix = 0; ix <= nxhalf; ix++)
02887                                 cmplx(ix,iy,iz) = conj(cmplx(ix,iy,iz));
02888         EXITFUNC;
02889 }
02890 
02891 void EMData::update_stat() const
02892 {
02893         ENTERFUNC;
02894 //      printf("update stat %f %d\n",(float)attr_dict["mean"],flags);
02895         if (!(flags & EMDATA_NEEDUPD))
02896         {
02897                 EXITFUNC;
02898                 return;
02899         }
02900         if (rdata==0) return;
02901 
02902         float* data = get_data();
02903         float max = -FLT_MAX;
02904         float min = -max;
02905 
02906         double sum = 0;
02907         double square_sum = 0;
02908 
02909         int step = 1;
02910         if (is_complex() && !is_ri()) {
02911                 step = 2;
02912         }
02913 
02914         int n_nonzero = 0;
02915 
02916         size_t size = (size_t)nx*ny*nz;
02917         for (size_t i = 0; i < size; i += step) {
02918                 float v = data[i];
02919         #ifdef _WIN32
02920                 max = _cpp_max(max,v);
02921                 min = _cpp_min(min,v);
02922         #else
02923                 max=std::max<float>(max,v);
02924                 min=std::min<float>(min,v);
02925         #endif  //_WIN32
02926                 sum += v;
02927                 square_sum += v * (double)(v);
02928                 if (v != 0) n_nonzero++;
02929         }
02930 
02931         size_t n = size / step;
02932         double mean = sum / n;
02933 
02934 #ifdef _WIN32
02935         float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / n)/(n-1)));
02936         n_nonzero = _cpp_max(1,n_nonzero);
02937         double sigma_nonzero = std::sqrt( _cpp_max(0,(square_sum  - sum*sum/n_nonzero)/(n_nonzero-1)));
02938 #else
02939         float sigma = (float)std::sqrt(std::max<double>(0.0,(square_sum - sum*sum / n)/(n-1)));
02940         n_nonzero = std::max<int>(1,n_nonzero);
02941         double sigma_nonzero = std::sqrt(std::max<double>(0,(square_sum  - sum*sum/n_nonzero)/(n_nonzero-1)));
02942 #endif  //_WIN32
02943         double mean_nonzero = sum / n_nonzero; // previous version overcounted! G2
02944 
02945         attr_dict["minimum"] = min;
02946         attr_dict["maximum"] = max;
02947         attr_dict["mean"] = (float)(mean);
02948         attr_dict["sigma"] = (float)(sigma);
02949         attr_dict["square_sum"] = (float)(square_sum);
02950         attr_dict["mean_nonzero"] = (float)(mean_nonzero);
02951         attr_dict["sigma_nonzero"] = (float)(sigma_nonzero);
02952         attr_dict["is_complex"] = (int) is_complex();
02953         attr_dict["is_complex_ri"] = (int) is_ri();
02954 
02955         flags &= ~EMDATA_NEEDUPD;
02956 
02957         if (rot_fp != 0)
02958         {
02959                 delete rot_fp; rot_fp = 0;
02960         }
02961 
02962         EXITFUNC;
02963 //      printf("done stat %f %f %f\n",(float)mean,(float)max,(float)sigma);
02964 }
02965 
02970 bool EMData::operator==(const EMData& that) const {
02971         if(this != &that) {
02972                 return false;
02973         }
02974         else {
02975                 return true;
02976         }
02977 }
02978 
02979 bool EMData::equal(const EMData& that) const {
02980         if (that.get_xsize() != nx || that.get_ysize() != ny || that.get_zsize() != nz ) return false;
02981 
02982         const float*  d1 = that.get_const_data();
02983         float* d2 = get_data();
02984 
02985         for(size_t i =0; i < get_size(); ++i,++d1,++d2) {
02986                 if ((*d1) != (*d2)) return false;
02987         }
02988 
02989 //      if(attr_dict != that.attr_dict) {
02990 //              return false;
02991 //      }
02992 
02993         return true;
02994 }
02995 
02996 EMData * EMAN::operator+(const EMData & em, float n)
02997 {
02998         EMData * r = em.copy();
02999         r->add(n);
03000         return r;
03001 }
03002 
03003 EMData * EMAN::operator-(const EMData & em, float n)
03004 {
03005         EMData* r = em.copy();
03006         r->sub(n);
03007         return r;
03008 }
03009 
03010 EMData * EMAN::operator*(const EMData & em, float n)
03011 {
03012         EMData* r = em.copy();
03013         r ->mult(n);
03014         return r;
03015 }
03016 
03017 EMData * EMAN::operator/(const EMData & em, float n)
03018 {
03019         EMData * r = em.copy();
03020         r->div(n);
03021         return r;
03022 }
03023 
03024 
03025 EMData * EMAN::operator+(float n, const EMData & em)
03026 {
03027         EMData * r = em.copy();
03028         r->add(n);
03029         return r;
03030 }
03031 
03032 EMData * EMAN::operator-(float n, const EMData & em)
03033 {
03034         EMData * r = em.copy();
03035         r->mult(-1.0f);
03036         r->add(n);
03037         return r;
03038 }
03039 
03040 EMData * EMAN::operator*(float n, const EMData & em)
03041 {
03042         EMData * r = em.copy();
03043         r->mult(n);
03044         return r;
03045 }
03046 
03047 EMData * EMAN::operator/(float n, const EMData & em)
03048 {
03049         EMData * r = em.copy();
03050         r->to_one();
03051         r->mult(n);
03052         r->div(em);
03053 
03054         return r;
03055 }
03056 
03057 EMData * EMAN::rsub(const EMData & em, float n)
03058 {
03059         return EMAN::operator-(n, em);
03060 }
03061 
03062 EMData * EMAN::rdiv(const EMData & em, float n)
03063 {
03064         return EMAN::operator/(n, em);
03065 }
03066 
03067 EMData * EMAN::operator+(const EMData & a, const EMData & b)
03068 {
03069         EMData * r = a.copy();
03070         r->add(b);
03071         return r;
03072 }
03073 
03074 EMData * EMAN::operator-(const EMData & a, const EMData & b)
03075 {
03076         EMData * r = a.copy();
03077         r->sub(b);
03078         return r;
03079 }
03080 
03081 EMData * EMAN::operator*(const EMData & a, const EMData & b)
03082 {
03083         EMData * r = a.copy();
03084         r->mult(b);
03085         return r;
03086 }
03087 
03088 EMData * EMAN::operator/(const EMData & a, const EMData & b)
03089 {
03090         EMData * r = a.copy();
03091         r->div(b);
03092         return r;
03093 }
03094 
03095 void EMData::set_xyz_origin(float origin_x, float origin_y, float origin_z)
03096 {
03097         attr_dict["origin_x"] = origin_x;
03098         attr_dict["origin_y"] = origin_y;
03099         attr_dict["origin_z"] = origin_z;
03100 }
03101 
03102 #if 0
03103 void EMData::calc_rcf(EMData * with, vector < float >&sum_array)
03104 {
03105         ENTERFUNC;
03106 
03107         int array_size = sum_array.size();
03108         float da = 2 * M_PI / array_size;
03109         float *dat = new float[array_size + 2];
03110         float *dat2 = new float[array_size + 2];
03111         int nx2 = nx * 9 / 20;
03112 
03113         float lim = 0;
03114         if (fabs(translation[0]) < fabs(translation[1])) {
03115                 lim = fabs(translation[1]);
03116         }
03117         else {
03118                 lim = fabs(translation[0]);
03119         }
03120 
03121         nx2 -= static_cast < int >(floor(lim));
03122 
03123         for (int i = 0; i < array_size; i++) {
03124                 sum_array[i] = 0;
03125         }
03126 
03127         float sigma = attr_dict["sigma"];
03128         float with_sigma = with->get_attr_dict().get("sigma");
03129 
03130         vector<float> vdata, vdata2;
03131         for (int i = 8; i < nx2; i += 6) {
03132                 vdata = calc_az_dist(array_size, 0, da, i, i + 6);
03133                 vdata2 = with->calc_az_dist(array_size, 0, da, i, i + 6);
03134                 Assert(vdata.size() <= array_size + 2);
03135                 Assert(cdata2.size() <= array_size + 2);
03136                 std::copy(vdata.begin(), vdata.end(), dat);
03137                 std::copy(vdata2.begin(), vdata2.end(), dat2);
03138 
03139                 EMfft::real_to_complex_1d(dat, dat, array_size);
03140                 EMfft::real_to_complex_1d(dat2, dat2, array_size);
03141 
03142                 for (int j = 0; j < array_size + 2; j += 2) {
03143                         float max = dat[j] * dat2[j] + dat[j + 1] * dat2[j + 1];
03144                         float max2 = dat[j + 1] * dat2[j] - dat2[j + 1] * dat[j];
03145                         dat[j] = max;
03146                         dat[j + 1] = max2;
03147                 }
03148 
03149                 EMfft::complex_to_real_1d(dat, dat, array_size);
03150                 float norm = array_size * array_size * (4.0f * sigma) * (4.0f * with_sigma);
03151 
03152                 for (int j = 0; j < array_size; j++) {
03153                         sum_array[j] += dat[j] * (float) i / norm;
03154                 }
03155         }
03156 
03157         if( dat )
03158         {
03159                 delete[]dat;
03160                 dat = 0;
03161         }
03162 
03163         if( dat2 )
03164         {
03165                 delete[]dat2;
03166                 dat2 = 0;
03167         }
03168         EXITFUNC;
03169 }
03170 
03171 #endif
03172 
03173 void EMData::add_incoherent(EMData * obj)
03174 {
03175         ENTERFUNC;
03176 
03177         if (!obj) {
03178                 LOGERR("NULL image");
03179                 throw NullPointerException("NULL image");
03180         }
03181 
03182         if (!obj->is_complex() || !is_complex()) {
03183                 throw ImageFormatException("complex images only");
03184         }
03185 
03186         if (!EMUtil::is_same_size(this, obj)) {
03187                 throw ImageFormatException("images not same size");
03188         }
03189 
03190         ri2ap();
03191         obj->ri2ap();
03192 
03193         float *dest = get_data();
03194         float *src = obj->get_data();
03195         size_t size = (size_t)nx * ny * nz;
03196         for (size_t j = 0; j < size; j += 2) {
03197 #ifdef  _WIN32
03198                 dest[j] = (float) _hypot(src[j], dest[j]);
03199 #else
03200                 dest[j] = (float) hypot(src[j], dest[j]);
03201 #endif  //_WIN32
03202                 dest[j + 1] = 0;
03203         }
03204 
03205         obj->update();
03206         update();
03207         EXITFUNC;
03208 }
03209 
03210 
03211 float EMData::calc_dist(EMData * second_img, int y_index) const
03212 {
03213         ENTERFUNC;
03214 
03215         if (get_ndim() != 1) {
03216                 throw ImageDimensionException("'this' image is 1D only");
03217         }
03218 
03219         if (second_img->get_xsize() != nx || ny != 1) {
03220                 throw ImageFormatException("image xsize not same");
03221         }
03222 
03223         if (y_index > second_img->get_ysize() || y_index < 0) {
03224                 return -1;
03225         }
03226 
03227         float ret = 0;
03228         float *d1 = get_data();
03229         float *d2 = second_img->get_data() + second_img->get_xsize() * y_index;
03230 
03231         for (int i = 0; i < nx; i++) {
03232                 ret += Util::square(d1[i] - d2[i]);
03233         }
03234         EXITFUNC;
03235         return std::sqrt(ret);
03236 }
03237 
03238 
03239 EMData * EMData::calc_fast_sigma_image( EMData* mask)
03240 {
03241         ENTERFUNC;
03242 
03243         bool maskflag = false;
03244         if (mask == 0) {
03245                 mask = new EMData(nx,ny,nz);
03246                 mask->process_inplace("testimage.circlesphere");
03247                 maskflag = true;
03248         }
03249 
03250         if (get_ndim() != mask->get_ndim() ) throw ImageDimensionException("The dimensions do not match");
03251 
03252         int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
03253 
03254         if ( mnx > nx || mny > ny || mnz > nz)
03255                 throw ImageDimensionException("Can not calculate variance map using an image that is larger than this image");
03256 
03257         size_t P = 0;
03258         for(size_t i = 0; i < mask->get_size(); ++i){
03259                 if (mask->get_value_at(i) != 0){
03260                         ++P;
03261                 }
03262         }
03263         float normfac = 1.0f/(float)P;
03264 
03265 //      bool undoclip = false;
03266 
03267         int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03268 //      if ( mnx < nx || mny < ny || mnz < nz) {
03269         Region r;
03270         if (ny == 1) r = Region((mnx-nxc)/2,nxc);
03271         else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03272         else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03273         mask->clip_inplace(r,0.0);
03274         //Region r((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03275         //mask->clip_inplace(r);
03276         //undoclip = true;
03277         //}
03278 
03279         // Here we generate the local average of the squares
03280         Region r2;
03281         if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03282         else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03283         else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03284         EMData* squared = get_clip(r2,get_edge_mean());
03285 
03286         EMData* tmp = squared->copy();
03287         Dict pow;
03288         pow["pow"] = 2.0f;
03289         squared->process_inplace("math.pow",pow);
03290         EMData* s = mask->convolute(squared);//ming, mask squared exchange
03291         squared->mult(normfac);
03292 
03293         EMData* m = mask->convolute(tmp);//ming, tmp mask exchange
03294         m->mult(normfac);
03295         m->process_inplace("math.pow",pow);
03296         delete tmp; tmp = 0;
03297         s->sub(*m);
03298         // Here we finally generate the standard deviation image
03299         s->process_inplace("math.sqrt");
03300 
03301 //      if ( undoclip ) {
03302 //              Region r((nx-mnx)/2, (ny-mny)/2, (nz-mnz)/2,mnx,mny,mnz);
03303 //              mask->clip_inplace(r);
03304 //      }
03305 
03306         if (maskflag) {
03307                 delete mask;
03308                 mask = 0;
03309         } else {
03310                 Region r;
03311                 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
03312                 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
03313                 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
03314                 mask->clip_inplace(r);
03315         }
03316 
03317         delete squared;
03318         delete m;
03319 
03320         s->process_inplace("xform.phaseorigin.tocenter");
03321         Region r3;
03322         if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03323         else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03324         else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03325         s->clip_inplace(r3);
03326         EXITFUNC;
03327         return s;
03328 }
03329 
03330 //  The following code looks strange - does anybody know it?  Please let me know, pawel.a.penczek@uth.tmc.edu  04/09/06.
03331 // This is just an implementation of "Roseman's" fast normalized cross-correlation (Ultramicroscopy, 2003). But the contents of this function have changed dramatically since you wrote that comment (d.woolford).
03332 EMData *EMData::calc_flcf(EMData * with)
03333 {
03334         ENTERFUNC;
03335         EMData *this_copy=this;
03336         this_copy=copy();
03337 
03338         int mnx = with->get_xsize(); int mny = with->get_ysize(); int mnz = with->get_zsize();
03339         int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03340 
03341         // Ones is a circular/spherical mask, consisting of 1s.
03342         EMData* ones = new EMData(mnx,mny,mnz);
03343         ones->process_inplace("testimage.circlesphere");
03344 
03345         // Get a copy of with, we will eventually resize it
03346         EMData* with_resized = with->copy();
03347         with_resized->process_inplace("normalize");
03348         with_resized->mult(*ones);
03349 
03350         EMData* s = calc_fast_sigma_image(ones);// Get the local sigma image
03351 
03352         Region r1;
03353         if (ny == 1) r1 = Region((mnx-nxc)/2,nxc);
03354         else if (nz == 1) r1 = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03355         else r1 = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03356         with_resized->clip_inplace(r1,0.0);
03357 
03358         Region r2;
03359         if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03360         else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03361         else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03362         this_copy->clip_inplace(r2,0.0);
03363 
03364         EMData* corr = this_copy->calc_ccf(with_resized); // the ccf results should have same size as sigma
03365 
03366         corr->process_inplace("xform.phaseorigin.tocenter");
03367         Region r3;
03368         if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03369         else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03370         else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03371         corr->clip_inplace(r3);
03372 
03373         corr->div(*s);
03374 
03375         delete with_resized; delete ones; delete this_copy; delete s;
03376         EXITFUNC;
03377         return corr;
03378 }
03379 
03380 EMData *EMData::convolute(EMData * with)
03381 {
03382         ENTERFUNC;
03383 
03384         EMData *f1 = do_fft();
03385         if (!f1) {
03386                 LOGERR("FFT returns NULL image");
03387                 throw NullPointerException("FFT returns NULL image");
03388         }
03389 
03390         f1->ap2ri();
03391 
03392         EMData *cf = 0;
03393         if (with) {
03394                 cf = with->do_fft();
03395                 if (!cf) {
03396                         LOGERR("FFT returns NULL image");
03397                         throw NullPointerException("FFT returns NULL image");
03398                 }
03399                 cf->ap2ri();
03400         }
03401         else {
03402                 cf = f1->copy();
03403         }
03404         //printf("cf_x=%d, f1y=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f1->get_ysize(),this->get_xsize(),with->get_xsize());
03405         if (with && !EMUtil::is_same_size(f1, cf)) {
03406                 LOGERR("images not same size");
03407                 throw ImageFormatException("images not same size");
03408         }
03409 
03410         float *rdata1 = f1->get_data();
03411         float *rdata2 = cf->get_data();
03412         size_t cf_size = (size_t)cf->get_xsize() * cf->get_ysize() * cf->get_zsize();
03413 
03414         float re,im;
03415 
03416         for (size_t i = 0; i < cf_size; i += 2) {
03417                 re = rdata1[i] * rdata2[i] - rdata1[i + 1] * rdata2[i + 1];
03418                 im = rdata1[i + 1] * rdata2[i] + rdata1[i] * rdata2[i + 1];
03419                 rdata2[i]=re;
03420                 rdata2[i+1]=im;
03421         }
03422         cf->update();
03423         EMData *f2 = cf->do_ift();//ming change cf to cf_temp
03424         //printf("cf_x=%d, f2x=%d, thisx=%d, withx=%d\n",cf->get_xsize(),f2->get_xsize(),this->get_xsize(),with->get_xsize());
03425         if( cf )
03426         {
03427                 delete cf;
03428                 cf = 0;
03429         }
03430 
03431         if( f1 )
03432         {
03433                 delete f1;
03434                 f1=0;
03435         }
03436 
03437         EXITFUNC;
03438         return f2;
03439 }
03440 
03441 
03442 void EMData::common_lines(EMData * image1, EMData * image2,
03443                                                   int mode, int steps, bool horizontal)
03444 {
03445         ENTERFUNC;
03446 
03447         if (!image1 || !image2) {
03448                 throw NullPointerException("NULL image");
03449         }
03450 
03451         if (mode < 0 || mode > 2) {
03452                 throw OutofRangeException(0, 2, mode, "invalid mode");
03453         }
03454 
03455         if (!image1->is_complex()) {
03456                 image1 = image1->do_fft();
03457         }
03458         if (!image2->is_complex()) {
03459                 image2 = image2->do_fft();
03460         }
03461 
03462         image1->ap2ri();
03463         image2->ap2ri();
03464 
03465         if (!EMUtil::is_same_size(image1, image2)) {
03466                 throw ImageFormatException("images not same sizes");
03467         }
03468 
03469         int image2_nx = image2->get_xsize();
03470         int image2_ny = image2->get_ysize();
03471 
03472         int rmax = image2_ny / 4 - 1;
03473         int array_size = steps * rmax * 2;
03474         float *im1 = new float[array_size];
03475         float *im2 = new float[array_size];
03476         for (int i = 0; i < array_size; i++) {
03477                 im1[i] = 0;
03478                 im2[i] = 0;
03479         }
03480 
03481         set_size(steps * 2, steps * 2, 1);
03482 
03483         float *image1_data = image1->get_data();
03484         float *image2_data = image2->get_data();
03485 
03486         float da = M_PI / steps;
03487         float a = -M_PI / 2.0f + da / 2.0f;
03488         int jmax = 0;
03489 
03490         for (int i = 0; i < steps * 2; i += 2, a += da) {
03491                 float s1 = 0;
03492                 float s2 = 0;
03493                 int i2 = i * rmax;
03494                 int j = 0;
03495 
03496                 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
03497                         float x = r * cos(a);
03498                         float y = r * sin(a);
03499 
03500                         if (x < 0) {
03501                                 x = -x;
03502                                 y = -y;
03503                                 LOGERR("CCL ERROR %d, %f !\n", i, -x);
03504                         }
03505 
03506                         int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
03507                         int l = i2 + j;
03508                         float x2 = x - floor(x);
03509                         float y2 = y - floor(y);
03510 
03511                         im1[l] = Util::bilinear_interpolate(image1_data[k],
03512                                                                                                 image1_data[k + 2],
03513                                                                                                 image1_data[k + image2_nx],
03514                                                                                                 image1_data[k + 2 + image2_nx], x2, y2);
03515 
03516                         im2[l] = Util::bilinear_interpolate(image2_data[k],
03517                                                                                                 image2_data[k + 2],
03518                                                                                                 image2_data[k + image2_nx],
03519                                                                                                 image2_data[k + 2 + image2_nx], x2, y2);
03520 
03521                         k++;
03522 
03523                         im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
03524                                                                                                         image1_data[k + 2],
03525                                                                                                         image1_data[k + image2_nx],
03526                                                                                                         image1_data[k + 2 + image2_nx], x2, y2);
03527 
03528                         im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
03529                                                                                                         image2_data[k + 2],
03530                                                                                                         image2_data[k + image2_nx],
03531                                                                                                         image2_data[k + 2 + image2_nx], x2, y2);
03532 
03533                         s1 += Util::square_sum(im1[l], im1[l + 1]);
03534                         s2 += Util::square_sum(im2[l], im2[l + 1]);
03535                 }
03536 
03537                 jmax = j - 1;
03538                 float sqrt_s1 = std::sqrt(s1);
03539                 float sqrt_s2 = std::sqrt(s2);
03540 
03541                 int l = 0;
03542                 for (float r = 1; r < rmax; r += 1.0f) {
03543                         int i3 = i2 + l;
03544                         im1[i3] /= sqrt_s1;
03545                         im1[i3 + 1] /= sqrt_s1;
03546                         im2[i3] /= sqrt_s2;
03547                         im2[i3 + 1] /= sqrt_s2;
03548                         l += 2;
03549                 }
03550         }
03551         float * data = get_data();
03552 
03553         switch (mode) {
03554         case 0:
03555                 for (int m1 = 0; m1 < 2; m1++) {
03556                         for (int m2 = 0; m2 < 2; m2++) {
03557 
03558                                 if (m1 == 0 && m2 == 0) {
03559                                         for (int i = 0; i < steps; i++) {
03560                                                 int i2 = i * rmax * 2;
03561                                                 for (int j = 0; j < steps; j++) {
03562                                                         int l = i + j * steps * 2;
03563                                                         int j2 = j * rmax * 2;
03564                                                         data[l] = 0;
03565                                                         for (int k = 0; k < jmax; k++) {
03566                                                                 data[l] += im1[i2 + k] * im2[j2 + k];
03567                                                         }
03568                                                 }
03569                                         }
03570                                 }
03571                                 else {
03572                                         int steps2 = steps * m2 + steps * steps * 2 * m1;
03573 
03574                                         for (int i = 0; i < steps; i++) {
03575                                                 int i2 = i * rmax * 2;
03576                                                 for (int j = 0; j < steps; j++) {
03577                                                         int j2 = j * rmax * 2;
03578                                                         int l = i + j * steps * 2 + steps2;
03579                                                         data[l] = 0;
03580 
03581                                                         for (int k = 0; k < jmax; k += 2) {
03582                                                                 i2 += k;
03583                                                                 j2 += k;
03584                                                                 data[l] += im1[i2] * im2[j2];
03585                                                                 data[l] += -im1[i2 + 1] * im2[j2 + 1];
03586                                                         }
03587                                                 }
03588                                         }
03589                                 }
03590                         }
03591                 }
03592 
03593                 break;
03594         case 1:
03595                 for (int m1 = 0; m1 < 2; m1++) {
03596                         for (int m2 = 0; m2 < 2; m2++) {
03597                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03598                                 int p1_sign = 1;
03599                                 if (m1 != m2) {
03600                                         p1_sign = -1;
03601                                 }
03602 
03603                                 for (int i = 0; i < steps; i++) {
03604                                         int i2 = i * rmax * 2;
03605 
03606                                         for (int j = 0; j < steps; j++) {
03607                                                 int j2 = j * rmax * 2;
03608 
03609                                                 int l = i + j * steps * 2 + steps2;
03610                                                 data[l] = 0;
03611                                                 float a = 0;
03612 
03613                                                 for (int k = 0; k < jmax; k += 2) {
03614                                                         i2 += k;
03615                                                         j2 += k;
03616 
03617 #ifdef  _WIN32
03618                                                         float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
03619 #else
03620                                                         float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
03621 #endif  //_WIN32
03622                                                         float p1 = atan2(im1[i2 + 1], im1[i2]);
03623                                                         float p2 = atan2(im2[j2 + 1], im2[j2]);
03624 
03625                                                         data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
03626                                                         a += a1;
03627                                                 }
03628 
03629                                                 data[l] /= (float)(a * M_PI / 180.0f);
03630                                         }
03631                                 }
03632                         }
03633                 }
03634 
03635                 break;
03636         case 2:
03637                 for (int m1 = 0; m1 < 2; m1++) {
03638                         for (int m2 = 0; m2 < 2; m2++) {
03639                                 int steps2 = steps * m2 + steps * steps * 2 * m1;
03640 
03641                                 for (int i = 0; i < steps; i++) {
03642                                         int i2 = i * rmax * 2;
03643 
03644                                         for (int j = 0; j < steps; j++) {
03645                                                 int j2 = j * rmax * 2;
03646                                                 int l = i + j * steps * 2 + steps2;
03647                                                 data[l] = 0;
03648 
03649                                                 for (int k = 0; k < jmax; k += 2) {
03650                                                         i2 += k;
03651                                                         j2 += k;
03652 #ifdef  _WIN32
03653                                                         data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
03654 #else
03655                                                         data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
03656 #endif  //_WIN32
03657                                                 }
03658                                         }
03659                                 }
03660                         }
03661                 }
03662 
03663                 break;
03664         default:
03665                 break;
03666         }
03667 
03668         if (horizontal) {
03669                 float *tmp_array = new float[ny];
03670                 for (int i = 1; i < nx; i++) {
03671                         for (int j = 0; j < ny; j++) {
03672                                 tmp_array[j] = get_value_at(i, j);
03673                         }
03674                         for (int j = 0; j < ny; j++) {
03675                                 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
03676                         }
03677                 }
03678                 if( tmp_array )
03679                 {
03680                         delete[]tmp_array;
03681                         tmp_array = 0;
03682                 }
03683         }
03684 
03685         if( im1 )
03686         {
03687                 delete[]im1;
03688                 im1 = 0;
03689         }
03690 
03691         if( im2 )
03692         {
03693                 delete im2;
03694                 im2 = 0;
03695         }
03696 
03697 
03698         image1->update();
03699         image2->update();
03700         if( image1 )
03701         {
03702                 delete image1;
03703                 image1 = 0;
03704         }
03705         if( image2 )
03706         {
03707                 delete image2;
03708                 image2 = 0;
03709         }
03710         update();
03711         EXITFUNC;
03712 }
03713 
03714 
03715 
03716 void EMData::common_lines_real(EMData * image1, EMData * image2,
03717                                                            int steps, bool horiz)
03718 {
03719         ENTERFUNC;
03720 
03721         if (!image1 || !image2) {
03722                 throw NullPointerException("NULL image");
03723         }
03724 
03725         if (!EMUtil::is_same_size(image1, image2)) {
03726                 throw ImageFormatException("images not same size");
03727         }
03728 
03729         int steps2 = steps * 2;
03730         int image_ny = image1->get_ysize();
03731         EMData *image1_copy = image1->copy();
03732         EMData *image2_copy = image2->copy();
03733 
03734         float *im1 = new float[steps2 * image_ny];
03735         float *im2 = new float[steps2 * image_ny];
03736 
03737         EMData *images[] = { image1_copy, image2_copy };
03738         float *ims[] = { im1, im2 };
03739 
03740         for (int m = 0; m < 2; m++) {
03741                 float *im = ims[m];
03742                 float a = M_PI / steps2;
03743                 Transform t(Dict("type","2d","alpha",-a));
03744                 for (int i = 0; i < steps2; i++) {
03745                         images[i]->transform(t);
03746                         float *data = images[i]->get_data();
03747 
03748                         for (int j = 0; j < image_ny; j++) {
03749                                 float sum = 0;
03750                                 for (int k = 0; k < image_ny; k++) {
03751                                         sum += data[j * image_ny + k];
03752                                 }
03753                                 im[i * image_ny + j] = sum;
03754                         }
03755 
03756                         float sum1 = 0;
03757                         float sum2 = 0;
03758                         for (int j = 0; j < image_ny; j++) {
03759                                 int l = i * image_ny + j;
03760                                 sum1 += im[l];
03761                                 sum2 += im[l] * im[l];
03762                         }
03763 
03764                         float mean = sum1 / image_ny;
03765                         float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03766 
03767                         for (int j = 0; j < image_ny; j++) {
03768                                 int l = i * image_ny + j;
03769                                 im[l] = (im[l] - mean) / sigma;
03770                         }
03771 
03772                         images[i]->update();
03773                         a += M_PI / steps;
03774                 }
03775         }
03776 
03777         set_size(steps2, steps2, 1);
03778         float *data1 = get_data();
03779 
03780         if (horiz) {
03781                 for (int i = 0; i < steps2; i++) {
03782                         for (int j = 0, l = i; j < steps2; j++, l++) {
03783                                 if (l == steps2) {
03784                                         l = 0;
03785                                 }
03786 
03787                                 float sum = 0;
03788                                 for (int k = 0; k < image_ny; k++) {
03789                                         sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03790                                 }
03791                                 data1[i + j * steps2] = sum;
03792                         }
03793                 }
03794         }
03795         else {
03796                 for (int i = 0; i < steps2; i++) {
03797                         for (int j = 0; j < steps2; j++) {
03798                                 float sum = 0;
03799                                 for (int k = 0; k < image_ny; k++) {
03800                                         sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03801                                 }
03802                                 data1[i + j * steps2] = sum;
03803                         }
03804                 }
03805         }
03806 
03807         update();
03808 
03809         if( image1_copy )
03810         {
03811                 delete image1_copy;
03812                 image1_copy = 0;
03813         }
03814 
03815         if( image2_copy )
03816         {
03817                 delete image2_copy;
03818                 image2_copy = 0;
03819         }
03820 
03821         if( im1 )
03822         {
03823                 delete[]im1;
03824                 im1 = 0;
03825         }
03826 
03827         if( im2 )
03828         {
03829                 delete[]im2;
03830                 im2 = 0;
03831         }
03832         EXITFUNC;
03833 }
03834 
03835 
03836 void EMData::cut_slice(const EMData *const map, const Transform& transform, bool interpolate)
03837 {
03838         ENTERFUNC;
03839 
03840         if (!map) throw NullPointerException("NULL image");
03841         // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
03842         if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
03843         if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
03844         // Now check for complex images - this is really just being thorough
03845         if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
03846         if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
03847 
03848 
03849         float *sdata = map->get_data();
03850         float *ddata = get_data();
03851 
03852         int map_nx = map->get_xsize();
03853         int map_ny = map->get_ysize();
03854         int map_nz = map->get_zsize();
03855         int map_nxy = map_nx * map_ny;
03856 
03857         int ymax = ny/2;
03858         if ( ny % 2 == 1 ) ymax += 1;
03859         int xmax = nx/2;
03860         if ( nx % 2 == 1 ) xmax += 1;
03861         for (int y = -ny/2; y < ymax; y++) {
03862                 for (int x = -nx/2; x < xmax; x++) {
03863                         Vec3f coord(x,y,0);
03864                         Vec3f soln = transform*coord;
03865 
03866 //                      float xx = (x+pretrans[0]) * (*ort)[0][0] +  (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
03867 //                      float yy = (x+pretrans[0]) * (*ort)[1][0] +  (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
03868 //                      float zz = (x+pretrans[0]) * (*ort)[2][0] +  (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
03869 
03870 
03871 //                      xx += map_nx/2;
03872 //                      yy += map_ny/2;
03873 //                      zz += map_nz/2;
03874 
03875                         float xx = soln[0]+map_nx/2;
03876                         float yy = soln[1]+map_ny/2;
03877                         float zz = soln[2]+map_nz/2;
03878 
03879                         int l = (x+nx/2) + (y+ny/2) * nx;
03880 
03881                         float t = xx - floor(xx);
03882                         float u = yy - floor(yy);
03883                         float v = zz - floor(zz);
03884 
03885                         if (xx < 0 || yy < 0 || zz < 0 ) {
03886                                 ddata[l] = 0;
03887                                 continue;
03888                         }
03889                         if (interpolate) {
03890                                 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03891                                         ddata[l] = 0;
03892                                         continue;
03893                                 }
03894                                 int k = (int) (Util::fast_floor(xx) + Util::fast_floor(yy) * map_nx + Util::fast_floor(zz) * map_nxy);
03895 
03896 
03897                                 if (xx < (map_nx - 1) && yy < (map_ny - 1) && zz < (map_nz - 1)) {
03898                                         ddata[l] = Util::trilinear_interpolate(sdata[k],
03899                                                                 sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
03900                                                                 sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
03901                                                                 sdata[k + map_nx + map_nxy + 1],t, u, v);
03902                                 }
03903                                 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03904                                         ddata[l] += sdata[k];
03905                                 }
03906                                 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) ) {
03907                                         ddata[l] +=     Util::linear_interpolate(sdata[k], sdata[k + map_nxy],v);
03908                                 }
03909                                 else if ( xx == (map_nx - 1) && zz == (map_nz - 1) ) {
03910                                         ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nx],u);
03911                                 }
03912                                 else if ( yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03913                                         ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + 1],t);
03914                                 }
03915                                 else if ( xx == (map_nx - 1) ) {
03916                                         ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + map_nx], sdata[k + map_nxy], sdata[k + map_nxy + map_nx],u,v);
03917                                 }
03918                                 else if ( yy == (map_ny - 1) ) {
03919                                         ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nxy], sdata[k + map_nxy + 1],t,v);
03920                                 }
03921                                 else if ( zz == (map_nz - 1) ) {
03922                                         ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nx], sdata[k + map_nx + 1],t,u);
03923                                 }
03924 
03925 //                              if (k >= map->get_size()) {
03926 //                                      cout << xx << " " << yy << " " <<  zz << " " << endl;
03927 //                                      cout << k << " " << get_size() << endl;
03928 //                                      cout << get_xsize() << " " << get_ysize() << " " << get_zsize() << endl;
03929 //                                      throw;
03930 //                                      }
03931 //
03932 //                              ddata[l] = Util::trilinear_interpolate(sdata[k],
03933 //                                              sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
03934 //                                              sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
03935 //                                              sdata[k + map_nx + map_nxy + 1],t, u, v);
03936                         }
03937                         else {
03938                                 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03939                                         ddata[l] = 0;
03940                                         continue;
03941                                 }
03942                                 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
03943                                 ddata[l] = sdata[k];
03944                         }
03945 
03946                 }
03947         }
03948 
03949         update();
03950 
03951         EXITFUNC;
03952 }
03953 
03954 EMData *EMData::unwrap_largerR(int r1,int r2,int xs, float rmax_f) {
03955         float *d,*dd;
03956         int do360=2;
03957         int rmax = (int)(rmax_f+0.5f);
03958         unsigned long i;
03959         unsigned int nvox=get_xsize()*get_ysize();//ming
03960         float maxmap=0.0f, minmap=0.0f;
03961         float temp=0.0f, diff_den=0.0f, norm=0.0f;
03962         float cut_off_va =0.0f;
03963 
03964         d=get_data();
03965         maxmap=-1000000.0f;
03966         minmap=1000000.0f;
03967         for (i=0;i<nvox;i++){
03968                 if(d[i]>maxmap) maxmap=d[i];
03969                 if(d[i]<minmap) minmap=d[i];
03970         }
03971         diff_den = maxmap-minmap;
03972         for(i=0;i<nvox;i++) {
03973                 temp = (d[i]-minmap)/diff_den;
03974                 if(cut_off_va != 0.0) {               // cut off the lowerset ?% noisy information
03975                         if(temp < cut_off_va)
03976                                 d[i] = 0.0f;                   // set the empty part density=0.0
03977                         else
03978                                 d[i] = temp-cut_off_va;
03979                 }
03980                 else    d[i] = temp;
03981         }
03982 
03983         for(i=0;i<nvox;i++) {
03984                 temp=d[i];
03985                 norm += temp*temp;
03986         }
03987         for(i=0;i<nvox;i++)             d[i] /= norm;                      //  y' = y/norm(y)
03988 
03989         if (xs<1) {
03990                 xs = (int) floor(do360*M_PI*get_ysize()/4); // ming
03991                 xs=Util::calc_best_fft_size(xs); // ming
03992         }
03993         if (r1<0) r1=0;
03994         float maxext=ceil(0.6f*std::sqrt((float)(get_xsize()*get_xsize()+get_ysize()*get_ysize())));// ming add std::
03995 
03996         if (r2<r1) r2=(int)maxext;
03997         EMData *ret = new EMData;
03998 
03999         ret->set_size(xs,r2+1,1);
04000 
04001         dd=ret->get_data();
04002 
04003         for (int i=0; i<xs; i++) {
04004                 float si=sin(i*M_PI*2/xs);
04005                 float co=cos(i*M_PI*2/xs);
04006                 for (int r=0; r<=maxext; r++) {
04007                         float x=(r+r1)*co+get_xsize()/2; // ming
04008                         float y=(r+r1)*si+get_ysize()/2; // ming
04009                         if(x<0.0 || x>=get_xsize()-1.0 || y<0.0 || y>=get_ysize()-1.0 || r>rmax){    //Ming , ~~~~ rmax need pass here
04010                                 for(;r<=r2;r++)                                   // here r2=MAXR
04011                                         dd[i+r*xs]=0.0;
04012                         break;
04013                     }
04014                         int x1=(int)floor(x);
04015                         int y1=(int)floor(y);
04016                         float t=x-x1;
04017                         float u=y-y1;
04018                         float f11= d[x1+y1*get_xsize()]; // ming
04019                         float f21= d[(x1+1)+y1*get_xsize()]; // ming
04020                         float f12= d[x1+(y1+1)*get_xsize()]; // ming
04021                         float f22= d[(x1+1)+(y1+1)*get_xsize()]; // ming
04022                         dd[i+r*xs] = (1-t)*(1-u)*f11+t*(1-u)*f21+t*u*f22+(1-t)*u*f12;
04023                 }
04024         }
04025         update();
04026         ret->update();
04027         return ret;
04028 }
04029 
04030 
04031 EMData *EMData::oneDfftPolar(int size, float rmax, float MAXR){         // sent MAXR value here later!!
04032         float *pcs=get_data();
04033         EMData *imagepcsfft = new EMData;
04034         imagepcsfft->set_size((size+2), (int)MAXR+1, 1);
04035         float *d=imagepcsfft->get_data();
04036 
04037         EMData *data_in=new EMData;
04038         data_in->set_size(size,1,1);
04039         float *in=data_in->get_data();
04040 
04041         for(int row=0; row<=(int)MAXR; ++row){
04042                 if(row<=(int)rmax) {
04043                         for(int i=0; i<size;++i)        in[i] = pcs[i+row*size]; // ming
04044                         data_in->set_complex(false);
04045                         data_in->do_fft_inplace();
04046                         for(int j=0;j<size+2;j++)  d[j+row*(size+2)]=in[j];
04047                 }
04048                 else for(int j=0;j<size+2;j++) d[j+row*(size+2)]=0.0;
04049         }
04050         imagepcsfft->update();
04051         delete data_in;
04052         return imagepcsfft;
04053 }
04054 
04055 void EMData::uncut_slice(EMData * const map, const Transform& transform) const
04056 {
04057         ENTERFUNC;
04058 
04059         if (!map) throw NullPointerException("NULL image");
04060         // These restrictions should be ultimately restricted so that all that matters is get_ndim() = (map->get_ndim() -1)
04061         if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
04062         if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
04063         // Now check for complex images - this is really just being thorough
04064         if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
04065         if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
04066 
04067 //      Transform3D r( 0, 0, 0); // EMAN by default
04068 //      if (!ort) {
04069 //              ort = &r;
04070 //      }
04071 
04072         float *ddata = map->get_data();
04073         float *sdata = get_data();
04074 
04075         int map_nx = map->get_xsize();
04076         int map_ny = map->get_ysize();
04077         int map_nz = map->get_zsize();
04078         int map_nxy = map_nx * map_ny;
04079         float map_nz_round_limit = (float) map_nz-0.5f;
04080         float map_ny_round_limit = (float) map_ny-0.5f;
04081         float map_nx_round_limit = (float) map_nx-0.5f;
04082 /*
04083         Vec3f posttrans = ort->get_posttrans();
04084         Vec3f pretrans = ort->get_pretrans();*/
04085 
04086         int ymax = ny/2;
04087         if ( ny % 2 == 1 ) ymax += 1;
04088         int xmax = nx/2;
04089         if ( nx % 2 == 1 ) xmax += 1;
04090         for (int y = -ny/2; y < ymax; y++) {
04091                 for (int x = -nx/2; x < xmax; x++) {
04092                         Vec3f coord(x,y,0);
04093                         Vec3f soln = transform*coord;
04094 //                      float xx = (x+pretrans[0]) * (*ort)[0][0] +  (y+pretrans[1]) * (*ort)[0][1] + pretrans[2] * (*ort)[0][2] + posttrans[0];
04095 //                      float yy = (x+pretrans[0]) * (*ort)[1][0] +  (y+pretrans[1]) * (*ort)[1][1] + pretrans[2] * (*ort)[1][2] + posttrans[1];
04096 //                      float zz = (x+pretrans[0]) * (*ort)[2][0] +  (y+pretrans[1]) * (*ort)[2][1] + pretrans[2] * (*ort)[2][2] + posttrans[2];
04097 //
04098 //                      xx += map_nx/2;
04099 //                      yy += map_ny/2;
04100 //                      zz += map_nz/2;
04101 //
04102                         float xx = soln[0]+map_nx/2;
04103                         float yy = soln[1]+map_ny/2;
04104                         float zz = soln[2]+map_nz/2;
04105 
04106                         // These 0.5 offsets are here because the round function rounds to the nearest whole number.
04107                         if (xx < -0.5 || yy < -0.5 || zz < -0.5 || xx >= map_nx_round_limit || yy >= map_ny_round_limit || zz >= map_nz_round_limit) continue;
04108 
04109                         size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
04110                         int l = (x+nx/2) + (y+ny/2) * nx;
04111                         ddata[k] = sdata[l];
04112                 }
04113         }
04114 
04115         map->update();
04116         EXITFUNC;
04117 }
04118 
04119 EMData *EMData::extract_box(const Transform& cs, const Region& r)
04120 {
04121         vector<float> cs_matrix = cs.get_matrix();
04122         
04123         EMData* box = new EMData();
04124         box->set_size((r.get_width()-r.x_origin()), (r.get_height()-r.y_origin()), (r.get_depth()-r.z_origin()));
04125         int box_nx = box->get_xsize();
04126         int box_ny = box->get_ysize();
04127         int box_nxy = box_nx*box_ny;
04128         float* bdata = box->get_data();
04129         float* ddata = get_data();
04130         
04131         for (int x = r.x_origin(); x < r.get_width(); x++) {
04132                 for (int y = r.y_origin(); y < r.get_height(); y++) {
04133                         for (int z = r.z_origin(); z < r.get_depth(); z++) {
04134                                 //float xb = cs_matrix[0]*x + cs_matrix[1]*y + cs_matrix[2]*z + cs_matrix[3];
04135                                 //float yb = cs_matrix[4]*x + cs_matrix[5]*y + cs_matrix[6]*z + cs_matrix[7];
04136                                 //float zb = cs_matrix[8]*x + cs_matrix[9]*y + cs_matrix[10]*z + cs_matrix[11];
04137                                 float xb = cs_matrix[0]*x + y*cs_matrix[4] + z*cs_matrix[8] + cs_matrix[3];
04138                                 float yb = cs_matrix[1]*x + y*cs_matrix[5] + z*cs_matrix[9] + cs_matrix[7];
04139                                 float zb = cs_matrix[2]*x + y*cs_matrix[6] + z*cs_matrix[10] + cs_matrix[11];
04140                                 float t = xb - Util::fast_floor(xb);
04141                                 float u = yb - Util::fast_floor(yb);
04142                                 float v = zb - Util::fast_floor(zb);
04143                                 
04144                                 //cout << x << " " << y << " " << z << " Box " << xb << " " << yb << " " << zb << endl;
04145                                 int l = (x - r.x_origin()) + (y - r.y_origin())*box_nx + (z - r.z_origin())*box_nxy;
04146                                 int k = (int) (Util::fast_floor(xb) + Util::fast_floor(yb) * nx + Util::fast_floor(zb) * nxy);
04147                                 //cout << k << " " << l << endl;
04148                                 if ( xb > nx - 1 || yb > ny - 1 || zb > nz - 1) {
04149                                         bdata[l] = 0;
04150                                         continue;
04151                                 }
04152                                 if (xb < 0 || yb < 0 || zb < 0){
04153                                         bdata[l] = 0;
04154                                         continue;
04155                                 }
04156 
04157                                 if (xb < (nx - 1) && yb < (ny - 1) && zb < (nz - 1)) {
04158                                         bdata[l] = Util::trilinear_interpolate(ddata[k], ddata[k + 1], ddata[k + nx],ddata[k + nx + 1], ddata[k + nxy], ddata[k + nxy + 1], ddata[k + nx + nxy], ddata[k + nx + nxy + 1],t, u, v);
04159                                 }
04160                         }
04161                 }
04162         }
04163         
04164         return box;
04165 }
04166 
04167 void EMData::save_byteorder_to_dict(ImageIO * imageio)
04168 {
04169         string image_endian = "ImageEndian";
04170         string host_endian = "HostEndian";
04171 
04172         if (imageio->is_image_big_endian()) {
04173                 attr_dict[image_endian] = "big";
04174         }
04175         else {
04176                 attr_dict[image_endian] = "little";
04177         }
04178 
04179         if (ByteOrder::is_host_big_endian()) {
04180                 attr_dict[host_endian] = "big";
04181         }
04182         else {
04183                 attr_dict[host_endian] = "little";
04184         }
04185 }
04186 
04187 EMData* EMData::compute_missingwedge(float wedgeangle, float start, float stop)
04188 {               
04189         EMData* test = new EMData();
04190         test->set_size(nx,ny,nz);
04191         
04192         float ratio = tan((90.0f-wedgeangle)*M_PI/180.0f);
04193         
04194         int offset_i = 2*int(start*nz/2);
04195         int offset_f = int(stop*nz/2);
04196         
04197         int step = 0;
04198         float sum = 0.0;
04199         double square_sum = 0.0;
04200         for (int j = 0; j < offset_f; j++){
04201                 for (int k = offset_i; k < offset_f; k++) {
04202                         for (int i = 0; i < nx; i+=2) {
04203                                 if (i < int(k*ratio)) {
04204                                         test->set_value_at(i, j, k, 1.0);
04205                                         float v = std::sqrt(pow(get_value_at_wrap(i, j, k),2) + pow(get_value_at_wrap(i+1, j, k),2));
04206                                         sum += v;
04207                                         square_sum += v * (double)(v);
04208                                         step++;
04209                                 }
04210                         }
04211                 }
04212         }
04213         
04214         float mean = sum / step;
04215         
04216         #ifdef _WIN32
04217         float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*mean)/(step-1)));
04218         #else
04219         float sigma = (float)std::sqrt(std::max<double>(0.0,(square_sum - sum*mean)/(step-1)));
04220         #endif  //_WIN32
04221         
04222         cout << "Mean sqr wedge amp " << mean << " Sigma Squ wedge Amp " << sigma << endl;
04223         set_attr("spt_wedge_mean", mean);
04224         set_attr("spt_wedge_sigma", sigma);
04225         
04226         return test;
04227 }
04228 
04229 
04230 
04231 
04232         
04233         

Generated on Fri Aug 10 16:35:12 2012 for EMAN2 by  doxygen 1.3.9.1