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