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

Generated on Tue Jun 11 13:40:37 2013 for EMAN2 by  doxygen 1.3.9.1