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