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

emdata.cpp

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

Generated on Fri Apr 30 15:38:50 2010 for EMAN2 by  doxygen 1.3.9.1