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 if (width != nx_fft || height != ny_fft ) {
01615 f1.set_size(width,height);
01616 f2.set_size(width,height);
01617 rslt.set_size(nx,height);
01618 nx_fft = width;
01619 ny_fft = height;
01620 }
01621
01622 #ifdef EMAN2_USING_CUDA
01623 if (EMData::usecuda == 1 && cudarwdata && with->cudarwdata) {
01624
01625 if(!f1.cudarwdata) f1.rw_alloc();
01626 if(!f2.cudarwdata) f2.rw_alloc();
01627 if(!rslt.cudarwdata) rslt.rw_alloc();
01628 cuda_dd_fft_real_to_complex_nd(cudarwdata, f1.cudarwdata, nx, 1, 1, height);
01629 cuda_dd_fft_real_to_complex_nd(with->cudarwdata, f2.cudarwdata, nx, 1, 1, height);
01630 calc_ccf_cuda(f1.cudarwdata, f2.cudarwdata, nx, ny, nz);
01631 cuda_dd_fft_complex_to_real_nd(f1.cudarwdata, rslt.cudarwdata, nx, 1, 1, height);
01632 if(no_sum){
01633 EMData* result = new EMData(rslt);
01634 return result;
01635 }
01636 EMData* cf = new EMData(0,0,nx,1,1);
01637 cf->runcuda(emdata_column_sum(rslt.cudarwdata, nx, ny));
01638 cf->update();
01639
01640 EXITFUNC;
01641 return cf;
01642 }
01643 #endif
01644
01645 float *d1 = get_data();
01646 float *d2 = with->get_data();
01647 float *f1d = f1.get_data();
01648 float *f2d = f2.get_data();
01649 for (int j = 0; j < height; j++) {
01650 EMfft::real_to_complex_1d(d1 + j * nx, f1d+j*width, nx);
01651 EMfft::real_to_complex_1d(d2 + j * nx, f2d+j*width, nx);
01652 }
01653
01654 if(flip == false) {
01655 for (int j = 0; j < height; j++) {
01656 float *f1a = f1d + j * width;
01657 float *f2a = f2d + j * width;
01658
01659 for (int i = 0; i < width / 2; i++) {
01660 float re1 = f1a[2*i];
01661 float re2 = f2a[2*i];
01662 float im1 = f1a[2*i+1];
01663 float im2 = f2a[2*i+1];
01664
01665 f1d[j*width+i*2] = re1 * re2 + im1 * im2;
01666 f1d[j*width+i*2+1] = im1 * re2 - re1 * im2;
01667 }
01668 }
01669 } else {
01670 for (int j = 0; j < height; j++) {
01671 float *f1a = f1d + j * width;
01672 float *f2a = f2d + j * width;
01673
01674 for (int i = 0; i < width / 2; i++) {
01675 float re1 = f1a[2*i];
01676 float re2 = f2a[2*i];
01677 float im1 = f1a[2*i+1];
01678 float im2 = f2a[2*i+1];
01679
01680 f1d[j*width+i*2] = re1 * re2 - im1 * im2;
01681 f1d[j*width+i*2+1] = im1 * re2 + re1 * im2;
01682 }
01683 }
01684 }
01685
01686 float* rd = rslt.get_data();
01687 for (int j = y0; j < y1; j++) {
01688 EMfft::complex_to_real_1d(f1d+j*width, rd+j*nx, nx);
01689 }
01690
01691 if (no_sum) {
01692 rslt.update();
01693 EXITFUNC;
01694 return new EMData(rslt);
01695 } else {
01696 EMData *cf = new EMData(nx,1,1);
01697 cf->to_zero();
01698 float *c = cf->get_data();
01699 for (int j = 0; j < height; j++) {
01700 for(int i = 0; i < nx; ++i) {
01701 c[i] += rd[i+j*nx];
01702 }
01703 }
01704 cf->update();
01705 EXITFUNC;
01706 return cf;
01707 }
01708 }
01709
01710 EMData *EMData::make_rotational_footprint_cmc( bool unwrap) {
01711 ENTERFUNC;
01712 update_stat();
01713
01714
01715
01716
01717
01718
01719 if ( rot_fp != 0 && unwrap == true) {
01720 return new EMData(*rot_fp);
01721 }
01722
01723 static EMData obj_filt;
01724 EMData* filt = &obj_filt;
01725 filt->set_complex(true);
01726
01727
01728
01729
01730
01731
01732 if (filt->get_xsize() != nx+2-(nx%2) || filt->get_ysize() != ny ||
01733 filt->get_zsize() != nz ) {
01734 filt->set_size(nx+2-(nx%2), ny, nz);
01735 filt->to_one();
01736
01737 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01738 }
01739
01740 EMData *ccf = this->calc_mutual_correlation(this, true,filt);
01741 ccf->sub(ccf->get_edge_mean());
01742 EMData *result = ccf->unwrap();
01743 delete ccf; ccf = 0;
01744
01745 EXITFUNC;
01746 if ( unwrap == true)
01747 {
01748
01749
01750
01751
01752
01753
01754
01755
01756 rot_fp = result;
01757 return new EMData(*rot_fp);
01758 }
01759 else return result;
01760 }
01761
01762 EMData *EMData::make_rotational_footprint( bool unwrap) {
01763 ENTERFUNC;
01764 update_stat();
01765
01766
01767
01768
01769
01770
01771 if ( rot_fp != 0 && unwrap == true) {
01772 return new EMData(*rot_fp);
01773 }
01774
01775 EMData* ccf = this->calc_ccf(this,CIRCULANT,true);
01776 ccf->sub(ccf->get_edge_mean());
01777
01778 EMData *result = ccf->unwrap();
01779 delete ccf; ccf = 0;
01780
01781 EXITFUNC;
01782 if ( unwrap == true)
01783 {
01784
01785
01786
01787
01788
01789
01790
01791 rot_fp = result;
01792 return new EMData(*rot_fp);
01793 }
01794 else return result;
01795 }
01796
01797 EMData *EMData::make_rotational_footprint_e1( bool unwrap)
01798 {
01799 ENTERFUNC;
01800
01801 update_stat();
01802
01803
01804
01805
01806
01807
01808 if ( rot_fp != 0 && unwrap == true) {
01809 return new EMData(*rot_fp);
01810 }
01811
01812 static EMData obj_filt;
01813 EMData* filt = &obj_filt;
01814 filt->set_complex(true);
01815
01816
01817
01818
01819
01820
01821 int cs = (((nx * 7 / 4) & 0xfffff8) - nx) / 2;
01822
01823 static EMData big_clip;
01824 int big_x = nx+2*cs;
01825 int big_y = ny+2*cs;
01826 int big_z = 1;
01827 if ( nz != 1 ) {
01828 big_z = nz+2*cs;
01829 }
01830
01831
01832 if ( big_clip.get_xsize() != big_x || big_clip.get_ysize() != big_y || big_clip.get_zsize() != big_z ) {
01833 big_clip.set_size(big_x,big_y,big_z);
01834 }
01835
01836
01837
01838
01839 big_clip.to_value(get_edge_mean());
01840
01841 if (nz != 1) {
01842 big_clip.insert_clip(this,IntPoint(cs,cs,cs));
01843 } else {
01844 big_clip.insert_clip(this,IntPoint(cs,cs,0));
01845 }
01846
01847
01848
01849
01850
01851 if (filt->get_xsize() != big_clip.get_xsize() +2-(big_clip.get_xsize()%2) || filt->get_ysize() != big_clip.get_ysize() ||
01852 filt->get_zsize() != big_clip.get_zsize()) {
01853 filt->set_size(big_clip.get_xsize() + 2-(big_clip.get_xsize()%2), big_clip.get_ysize(), big_clip.get_zsize());
01854 filt->to_one();
01855 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01856 #ifdef EMAN2_USING_CUDA
01857
01858
01859
01860
01861
01862
01863 #endif
01864 }
01865 #ifdef EMAN2_USING_CUDA
01866
01867
01868
01869
01870
01871
01872 #endif
01873
01874 EMData *mc = big_clip.calc_mutual_correlation(&big_clip, true,filt);
01875 mc->sub(mc->get_edge_mean());
01876
01877 static EMData sml_clip;
01878 int sml_x = nx * 3 / 2;
01879 int sml_y = ny * 3 / 2;
01880 int sml_z = 1;
01881 if ( nz != 1 ) {
01882 sml_z = nz * 3 / 2;
01883 }
01884
01885 if ( sml_clip.get_xsize() != sml_x || sml_clip.get_ysize() != sml_y || sml_clip.get_zsize() != sml_z ) {
01886 sml_clip.set_size(sml_x,sml_y,sml_z); }
01887 if (nz != 1) {
01888 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,-cs+nz/4));
01889 } else {
01890 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,0));
01891 }
01892
01893 delete mc; mc = 0;
01894 EMData * result = NULL;
01895
01896 if (nz == 1) {
01897 if (!unwrap) {
01898 #ifdef EMAN2_USING_CUDA
01899
01900 #endif
01901 result = sml_clip.process("mask.sharp", Dict("outer_radius", -1, "value", 0));
01902
01903 }
01904 else {
01905 result = sml_clip.unwrap();
01906 }
01907 }
01908 else {
01909
01910
01911
01912 result = new EMData(sml_clip);
01913 }
01914
01915 #ifdef EMAN2_USING_CUDA
01916
01917 #endif
01918 EXITFUNC;
01919 if ( unwrap == true)
01920 {
01921
01922
01923
01924 if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
01925
01926
01927
01928 rot_fp = result;
01929 return new EMData(*rot_fp);
01930 }
01931 else return result;
01932 }
01933
01934 EMData *EMData::make_footprint(int type)
01935 {
01936
01937 if (type==0) {
01938 EMData *un=make_rotational_footprint_e1();
01939 if (un->get_ysize() <= 6) {
01940 throw UnexpectedBehaviorException("In EMData::make_footprint. The rotational footprint is too small");
01941 }
01942 EMData *tmp=un->get_clip(Region(0,4,un->get_xsize(),un->get_ysize()-6));
01943 EMData *cx=tmp->calc_ccfx(tmp,0,-1,1);
01944 EMData *fp=cx->get_clip(Region(0,0,cx->get_xsize()/2,cx->get_ysize()));
01945 delete un;
01946 delete tmp;
01947 delete cx;
01948 return fp;
01949 }
01950 else if (type==1 || type==2 ||type==5 || type==6) {
01951 int i,j,kx,ky,lx,ly;
01952
01953 EMData *fft=do_fft();
01954
01955
01956 int rmax=(get_xsize()+1)/2;
01957 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
01958 for (i=0; i<rmax; i++) {
01959 for (j=0; j<rmax; j++) {
01960 #ifdef _WIN32
01961 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
01962 #else
01963 rmap[i+j*rmax]=hypot((float)i,(float)j);
01964 #endif //_WIN32
01965
01966 }
01967 }
01968
01969 EMData *fp=new EMData(rmax*2+2,rmax*2,1);
01970 fp->set_complex(1);
01971 fp->to_zero();
01972
01973
01974
01975
01976 for (kx=-rmax+1; kx<rmax; kx++) {
01977 for (ky=-rmax+1; ky<rmax; ky++) {
01978 for (lx=-rmax+1; lx<rmax; lx++) {
01979 for (ly=-rmax+1; ly<rmax; ly++) {
01980 int ax=kx+lx;
01981 int ay=ky+ly;
01982 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
01983 int r1=(int)floor(.5+rmap[abs(kx)+rmax*abs(ky)]);
01984 int r2=(int)floor(.5+rmap[abs(lx)+rmax*abs(ly)]);
01985
01986
01987 if (r1+r2>=rmax) continue;
01988
01989 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
01990 fp->set_value_at(r1*2,r2,p.real()+fp->get_value_at(r1*2,r2));
01991
01992
01993 fp->set_value_at(r1*2+1,r2,fp->get_value_at(r1*2+1,r2)+1);
01994 }
01995 }
01996 }
01997 }
01998
01999
02000 if (type==5 || type==6) {
02001 for (i=0; i<rmax*2; i+=2) {
02002 for (j=0; j<rmax; j++) {
02003 float norm=fp->get_value_at(i+1,j);
02004 #ifdef _WIN32
02005 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));
02006 fp->set_value_at(i,j,pow(fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm), 1.0f/3.0f));
02007 #else
02008 fp->set_value_at(i,rmax*2-j-1,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
02009 fp->set_value_at(i,j,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
02010 #endif //_WIN32
02011 fp->set_value_at(i+1,j,0.0);
02012 }
02013 }
02014 }
02015 else {
02016 for (i=0; i<rmax*2; i+=2) {
02017 for (j=0; j<rmax; j++) {
02018 float norm=fp->get_value_at(i+1,j);
02019 fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
02020 fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
02021 fp->set_value_at(i+1,j,0.0);
02022 }
02023 }
02024 }
02025
02026 free(rmap);
02027 if (type==2||type==6) {
02028 EMData *f2=fp->do_ift();
02029 if (f2->get_value_at(0,0)<0) f2->mult(-1.0f);
02030 f2->process_inplace("xform.phaseorigin.tocorner");
02031 delete fp;
02032 return f2;
02033 }
02034 return fp;
02035 }
02036 else if (type==3 || type==4) {
02037 int h,i,j,kx,ky,lx,ly;
02038
02039 EMData *fft=do_fft();
02040
02041
02042 int rmax=(get_xsize()+1)/2;
02043 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
02044 for (i=0; i<rmax; i++) {
02045 for (j=0; j<rmax; j++) {
02046 #ifdef _WIN32
02047 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
02048 #else
02049 rmap[i+j*rmax]=hypot((float)i,(float)j);
02050 #endif //_WIN32
02051
02052 }
02053 }
02054
02055 EMData *fp=new EMData(rmax*2+2,rmax*2,16);
02056
02057 fp->set_complex(1);
02058 fp->to_zero();
02059
02060
02061
02062
02063 for (kx=-rmax+1; kx<rmax; kx++) {
02064 for (ky=-rmax+1; ky<rmax; ky++) {
02065 for (lx=-rmax+1; lx<rmax; lx++) {
02066 for (ly=-rmax+1; ly<rmax; ly++) {
02067 int ax=kx+lx;
02068 int ay=ky+ly;
02069 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
02070 float rr1=rmap[abs(kx)+rmax*abs(ky)];
02071 float rr2=rmap[abs(lx)+rmax*abs(ly)];
02072 int r1=(int)floor(.5+rr1);
02073 int r2=(int)floor(.5+rr2);
02074
02075
02076 if (r1+r2>=rmax || rr1==0 ||rr2==0) continue;
02077
02078 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
02079 int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5);
02080 if (dot<0) dot=16+dot;
02081
02082 fp->set_value_at(r1*2,r2,dot,p.real()+fp->get_value_at(r1*2,r2,dot));
02083
02084
02085 fp->set_value_at(r1*2+1,r2,dot,fp->get_value_at(r1*2+1,r2,dot)+1);
02086 }
02087 }
02088 }
02089 }
02090
02091
02092 for (i=0; i<rmax*2; i+=2) {
02093 for (j=0; j<rmax; j++) {
02094 for (h=0; h<16; h++) {
02095 float norm=fp->get_value_at(i+1,j,h);
02096
02097
02098 fp->set_value_at(i,rmax*2-j-1,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02099 fp->set_value_at(i,j,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02100
02101
02102 fp->set_value_at(i+1,j,h,0.0);
02103 }
02104 }
02105 }
02106
02107 free(rmap);
02108 if (type==4) {
02109 EMData *f2=fp->do_ift();
02110 if (f2->get_value_at(0,0,0)<0) f2->mult(-1.0f);
02111 f2->process_inplace("xform.phaseorigin.tocorner");
02112 delete fp;
02113 return f2;
02114 }
02115 return fp;
02116 }
02117 throw UnexpectedBehaviorException("There is not implementation for the parameters you specified");
02118 }
02119
02120
02121 EMData *EMData::calc_mutual_correlation(EMData * with, bool tocenter, EMData * filter)
02122 {
02123 ENTERFUNC;
02124
02125 if (with && !EMUtil::is_same_size(this, with)) {
02126 LOGERR("images not same size");
02127 throw ImageFormatException( "images not same size");
02128 }
02129
02130 #ifdef EMAN2_USING_CUDA
02131 if(EMData::usecuda == 1 && cudarwdata && with->cudarwdata)
02132 {
02133
02134 EMData* this_fft = do_fft_cuda();
02135
02136 EMData *cf = 0;
02137 if (with && with != this) {
02138 cf = with->do_fft_cuda();
02139 }else{
02140 cf = this_fft->copy();
02141 }
02142
02143 if (filter) {
02144 if (!EMUtil::is_same_size(filter, cf)) {
02145 LOGERR("improperly sized filter");
02146 throw ImageFormatException("improperly sized filter");
02147 }
02148 mult_complex_efficient_cuda(cf->cudarwdata, filter->cudarwdata, cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), 1);
02149 mult_complex_efficient_cuda(this_fft->cudarwdata, filter->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize(), 1);
02150 }
02151
02152 mcf_cuda(this_fft->cudarwdata, cf->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize());
02153
02154 EMData *f2 = cf->do_ift_cuda();
02155
02156 if (tocenter) {
02157 f2->process_inplace("xform.phaseorigin.tocenter");
02158 }
02159
02160 if( cf )
02161 {
02162 delete cf;
02163 cf = 0;
02164 }
02165
02166 if( this_fft )
02167 {
02168 delete this_fft;
02169 this_fft = 0;
02170 }
02171
02172 f2->set_attr("label", "MCF");
02173 f2->set_path("/tmp/eman.mcf");
02174 f2->update();
02175
02176 EXITFUNC;
02177 return f2;
02178 }
02179 #endif
02180
02181 EMData *this_fft = 0;
02182 this_fft = do_fft();
02183
02184 if (!this_fft) {
02185
02186 LOGERR("FFT returns NULL image");
02187 throw NullPointerException("FFT returns NULL image");
02188 }
02189
02190 this_fft->ap2ri();
02191 EMData *cf = 0;
02192
02193 if (with && with != this) {
02194 cf = with->do_fft();
02195 if (!cf) {
02196 LOGERR("FFT returns NULL image");
02197 throw NullPointerException("FFT returns NULL image");
02198 }
02199 cf->ap2ri();
02200 }
02201 else {
02202 cf = this_fft->copy();
02203 }
02204
02205 if (filter) {
02206 if (!EMUtil::is_same_size(filter, cf)) {
02207 LOGERR("improperly sized filter");
02208 throw ImageFormatException("improperly sized filter");
02209 }
02210
02211 cf->mult_complex_efficient(*filter,true);
02212 this_fft->mult(*filter,true);
02213
02214
02215
02216
02217 }
02218
02219 float *rdata1 = this_fft->get_data();
02220 float *rdata2 = cf->get_data();
02221 size_t this_fft_size = (size_t)this_fft->get_xsize() * this_fft->get_ysize() * this_fft->get_zsize();
02222
02223 if (with == this) {
02224 for (size_t i = 0; i < this_fft_size; i += 2) {
02225 rdata2[i] = std::sqrt(rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02226 rdata2[i + 1] = 0;
02227 }
02228
02229 this_fft->update();
02230 cf->update();
02231 }
02232 else {
02233 for (size_t i = 0; i < this_fft_size; i += 2) {
02234 rdata2[i] = (rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02235 rdata2[i + 1] = (rdata1[i + 1] * rdata2[i] - rdata1[i] * rdata2[i + 1]);
02236 }
02237
02238
02239 for (size_t i = 0; i < this_fft_size; i += 2) {
02240 float t = Util::square(rdata2[i]) + Util::square(rdata2[i + 1]);
02241 if (t != 0) {
02242 t = pow(t, 0.25f);
02243 rdata2[i] /= t;
02244 rdata2[i + 1] /= t;
02245 }
02246 }
02247 this_fft->update();
02248 cf->update();
02249 }
02250
02251 EMData *f2 = cf->do_ift();
02252
02253 if (tocenter) {
02254 f2->process_inplace("xform.phaseorigin.tocenter");
02255 }
02256
02257 if( cf )
02258 {
02259 delete cf;
02260 cf = 0;
02261 }
02262
02263 if( this_fft )
02264 {
02265 delete this_fft;
02266 this_fft = 0;
02267 }
02268
02269 f2->set_attr("label", "MCF");
02270 f2->set_path("/tmp/eman.mcf");
02271
02272 EXITFUNC;
02273 return f2;
02274 }
02275
02276
02277 vector < float > EMData::calc_hist(int hist_size, float histmin, float histmax,const float& brt, const float& cont)
02278 {
02279 ENTERFUNC;
02280
02281 static size_t prime[] = { 1, 3, 7, 11, 17, 23, 37, 59, 127, 253, 511 };
02282
02283 if (histmin == histmax) {
02284 histmin = get_attr("minimum");
02285 histmax = get_attr("maximum");
02286 }
02287
02288 vector <float> hist(hist_size, 0.0);
02289
02290 int p0 = 0;
02291 int p1 = 0;
02292 size_t size = (size_t)nx * ny * nz;
02293 if (size < 300000) {
02294 p0 = 0;
02295 p1 = 0;
02296 }
02297 else if (size < 2000000) {
02298 p0 = 2;
02299 p1 = 3;
02300 }
02301 else if (size < 8000000) {
02302 p0 = 4;
02303 p1 = 6;
02304 }
02305 else {
02306 p0 = 7;
02307 p1 = 9;
02308 }
02309
02310 if (is_complex() && p0 > 0) {
02311 p0++;
02312 p1++;
02313 }
02314
02315 size_t di = 0;
02316
02317 size_t n = hist.size();
02318
02319 float * data = get_data();
02320 for (int k = p0; k <= p1; ++k) {
02321 if (is_complex()) {
02322 di = prime[k] * 2;
02323 }
02324 else {
02325 di = prime[k];
02326 }
02327
02328
02329 float w = (float)n / (histmax - histmin);
02330
02331 for(size_t i=0; i<=size-di; i += di) {
02332 float val;
02333 if (cont != 1.0f || brt != 0)val = cont*(data[i]+brt);
02334 else val = data[i];
02335 int j = Util::round((val - histmin) * w);
02336 if (j >= 0 && j < (int) n) {
02337 hist[j] += 1;
02338 }
02339 }
02340 }
02341
02342
02343
02344
02345
02346
02347
02348 return hist;
02349
02350 EXITFUNC;
02351 }
02352
02353
02354
02355
02356
02357 vector<float> EMData::calc_az_dist(int n, float a0, float da, float rmin, float rmax)
02358 {
02359 ENTERFUNC;
02360
02361 a0=a0*M_PI/180.0f;
02362 da=da*M_PI/180.0f;
02363
02364 if (get_ndim() > 2) {
02365 throw ImageDimensionException("no 3D image");
02366 }
02367
02368 float *yc = new float[n];
02369
02370 vector<float> vd(n);
02371 for (int i = 0; i < n; i++) {
02372 yc[i] = 0.00001f;
02373 }
02374
02375 float * data = get_data();
02376 if (is_complex()) {
02377 int c = 0;
02378 for (int y = 0; y < ny; y++) {
02379 for (int x = 0; x < nx; x += 2, c += 2) {
02380 int x1 = x / 2;
02381 int y1 = y<ny/2?y:y-ny;
02382 float r = (float)Util::hypot_fast(x1,y1);
02383
02384 if (r >= rmin && r <= rmax) {
02385 float a = 0;
02386
02387 if (y != ny / 2 || x != 0) {
02388 a = (atan2((float)y1, (float)x1) - a0) / da;
02389 }
02390
02391 int i = (int)(floor(a));
02392 a -= i;
02393
02394 if (i == 0) {
02395 vd[0] += data[c] * (1.0f - a);
02396 yc[0] += (1.0f - a);
02397 }
02398 else if (i == n - 1) {
02399 vd[n - 1] += data[c] * a;
02400 yc[n - 1] += a;
02401 }
02402 else if (i > 0 && i < (n - 1)) {
02403 float h = 0;
02404 if (is_ri()) {
02405 #ifdef _WIN32
02406 h = (float)_hypot(data[c], data[c + 1]);
02407 #else
02408 h = (float)hypot(data[c], data[c + 1]);
02409 #endif //_WIN32
02410 }
02411 else {
02412 h = data[c];
02413 }
02414
02415 vd[i] += h * (1.0f - a);
02416 yc[i] += (1.0f - a);
02417 vd[i + 1] += h * a;
02418 yc[i + 1] += a;
02419 }
02420 }
02421 }
02422 }
02423 }
02424 else {
02425 int c = 0;
02426 float half_nx = (nx - 1) / 2.0f;
02427 float half_ny = (ny - 1) / 2.0f;
02428
02429 for (int y = 0; y < ny; y++) {
02430 for (int x = 0; x < nx; x++, c++) {
02431 float y1 = y - half_ny;
02432 float x1 = x - half_nx;
02433 #ifdef _WIN32
02434 float r = (float)_hypot(x1, y1);
02435 #else
02436 float r = (float)hypot(x1, y1);
02437 #endif
02438
02439 if (r >= rmin && r <= rmax) {
02440 float a = 0;
02441 if (x1 != 0 || y1 != 0) {
02442 a = atan2(y1, x1);
02443 if (a < 0) {
02444 a += M_PI * 2;
02445 }
02446 }
02447
02448 a = (a - a0) / da;
02449 int i = static_cast < int >(floor(a));
02450 a -= i;
02451
02452 if (i == 0) {
02453 vd[0] += data[c] * (1.0f - a);
02454 yc[0] += (1.0f - a);
02455 }
02456 else if (i == n - 1) {
02457 vd[n - 1] += data[c] * a;
02458 yc[n - 1] += (a);
02459 }
02460 else if (i > 0 && i < (n - 1)) {
02461 vd[i] += data[c] * (1.0f - a);
02462 yc[i] += (1.0f - a);
02463 vd[i + 1] += data[c] * a;
02464 yc[i + 1] += a;
02465 }
02466 }
02467 }
02468 }
02469 }
02470
02471
02472 for (int i = 0; i < n; i++) {
02473 vd[i] /= yc[i];
02474 }
02475
02476 if( yc )
02477 {
02478 delete[]yc;
02479 yc = 0;
02480 }
02481
02482 return vd;
02483
02484 EXITFUNC;
02485 }
02486
02487
02488 EMData *EMData::unwrap(int r1, int r2, int xs, int dx, int dy, bool do360, bool weight_radial) const
02489 {
02490 ENTERFUNC;
02491
02492 if (get_ndim() != 2) {
02493 throw ImageDimensionException("2D image only");
02494 }
02495
02496 int p = 1;
02497 if (do360) {
02498 p = 2;
02499 }
02500
02501 if (xs < 1) {
02502 xs = (int) Util::fast_floor(p * M_PI * ny / 4);
02503 xs -= xs % 8;
02504 if (xs<=8) xs=16;
02505 }
02506
02507 if (r1 < 0) {
02508 r1 = 4;
02509 }
02510
02511 #ifdef _WIN32
02512 int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(_hypot(dx, dy)));
02513 #else
02514 int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(hypot(dx, dy)));
02515 #endif //_WIN32
02516 rr-=rr%2;
02517 if (r2 <= r1 || r2 > rr) {
02518 r2 = rr;
02519 }
02520
02521 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");
02522
02523 #ifdef EMAN2_USING_CUDA
02524 if (EMData::usecuda == 1 && isrodataongpu()){
02525
02526 EMData* result = new EMData(0,0,xs,(r2-r1),1);
02527 if(!result->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
02528 bindcudaarrayA(true);
02529 emdata_unwrap(result->cudarwdata, r1, r2, xs, p, dx, dy, weight_radial, nx, ny);
02530 unbindcudaarryA();
02531 result->update();
02532 return result;
02533 }
02534 #endif
02535
02536 EMData *ret = new EMData();
02537 ret->set_size(xs, r2 - r1, 1);
02538 const float *const d = get_const_data();
02539 float *dd = ret->get_data();
02540 float pfac = (float)p/(float)xs;
02541
02542 int nxon2 = nx/2;
02543 int nyon2 = ny/2;
02544 for (int x = 0; x < xs; x++) {
02545 float ang = x * M_PI * pfac;
02546 float si = sin(ang);
02547 float co = cos(ang);
02548
02549 for (int y = 0; y < r2 - r1; y++) {
02550 float ypr1 = (float)y + r1;
02551 float xx = ypr1 * co + nxon2 + dx;
02552 float yy = ypr1 * si + nyon2 + dy;
02553
02554
02555 float t = xx - (int)xx;
02556 float u = yy - (int)yy;
02557
02558 int k = (int) xx + ((int) yy) * nx;
02559 float val = Util::bilinear_interpolate(d[k], d[k + 1], d[k + nx], d[k + nx+1], t,u);
02560 if (weight_radial) val *= ypr1;
02561 dd[x + y * xs] = val;
02562 }
02563
02564 }
02565 ret->update();
02566
02567 EXITFUNC;
02568 return ret;
02569 }
02570
02571
02572
02573 void EMData::apply_radial_func(float x0, float step, vector < float >array, bool interp)
02574 {
02575 ENTERFUNC;
02576
02577 if (!is_complex()) throw ImageFormatException("apply_radial_func requires a complex image");
02578
02579 int n = static_cast < int >(array.size());
02580
02581 if (n*step>2.0) printf("Warning, apply_radial_func takes x0 and step with respect to Nyquist (0.5)\n");
02582
02583
02584
02585 ap2ri();
02586
02587 size_t ndims = get_ndim();
02588 float * data = get_data();
02589 if (ndims == 2) {
02590 int k = 0;
02591 for (int j = 0; j < ny; j++) {
02592 for (int i = 0; i < nx; i += 2, k += 2) {
02593 float r;
02594 #ifdef _WIN32
02595 if (j<ny/2) r = (float)_hypot(i/(float)(nx*2), j/(float)ny);
02596 else r = (float)_hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02597 #else
02598 if (j<ny/2) r = (float)hypot(i/(float)(nx*2), j/(float)ny);
02599 else r = (float)hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02600 #endif //_WIN32
02601 r = (r - x0) / step;
02602
02603 int l = 0;
02604 if (interp) {
02605 l = (int) floor(r);
02606 }
02607 else {
02608 l = (int) floor(r + 1);
02609 }
02610
02611
02612 float f = 0;
02613 if (l >= n - 2) {
02614 f = array[n - 1];
02615 }
02616 else {
02617 if (interp) {
02618 r -= l;
02619 f = (array[l] * (1.0f - r) + array[l + 1] * r);
02620 }
02621 else {
02622 f = array[l];
02623 }
02624 }
02625
02626 data[k] *= f;
02627 data[k + 1] *= f;
02628 }
02629 }
02630 }
02631 else if (ndims == 3) {
02632 int k = 0;
02633 for (int m = 0; m < nz; m++) {
02634 float mnz;
02635 if (m<nz/2) mnz=m*m/(float)(nz*nz);
02636 else { mnz=(nz-m)/(float)nz; mnz*=mnz; }
02637
02638 for (int j = 0; j < ny; j++) {
02639 float jny;
02640 if (j<ny/2) jny= j*j/(float)(ny*ny);
02641 else { jny=(ny-j)/(float)ny; jny*=jny; }
02642
02643 for (int i = 0; i < nx; i += 2, k += 2) {
02644 float r = std::sqrt((i * i / (nx*nx*4.0f)) + jny + mnz);
02645 r = (r - x0) / step;
02646
02647 int l = 0;
02648 if (interp) {
02649 l = (int) floor(r);
02650 }
02651 else {
02652 l = (int) floor(r + 1);
02653 }
02654
02655
02656 float f = 0;
02657 if (l >= n - 2) {
02658 f = array[n - 1];
02659 }
02660 else {
02661 if (interp) {
02662 r -= l;
02663 f = (array[l] * (1.0f - r) + array[l + 1] * r);
02664 }
02665 else {
02666 f = array[l];
02667 }
02668 }
02669
02670 data[k] *= f;
02671 data[k + 1] *= f;
02672 }
02673 }
02674 }
02675
02676 }
02677
02678 update();
02679 EXITFUNC;
02680 }
02681
02682 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, bool inten)
02683 {
02684 ENTERFUNC;
02685
02686 vector<float>ret(n);
02687 vector<float>norm(n);
02688
02689 int x,y,z,i;
02690 int step=is_complex()?2:1;
02691 int isinten=get_attr_default("is_intensity",0);
02692
02693 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
02694
02695 for (i=0; i<n; i++) ret[i]=norm[i]=0.0;
02696 float * data = get_data();
02697
02698
02699 if (nz==1) {
02700 for (y=i=0; y<ny; y++) {
02701 for (x=0; x<nx; x+=step,i+=step) {
02702 float r,v;
02703 if (step==2) {
02704 if (x==0 && y>ny/2) continue;
02705 r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y));
02706 if (!inten) {
02707 #ifdef _WIN32
02708 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02709 #else
02710 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02711 #endif
02712 else v=data[i];
02713 } else {
02714 if (isinten) v=data[i];
02715 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02716 else v=data[i]*data[i];
02717 }
02718 }
02719 else {
02720 r=(float)(Util::hypot_fast(x-nx/2,y-ny/2));
02721 if (inten) v=data[i]*data[i];
02722 else v=data[i];
02723 }
02724 r=(r-x0)/dx;
02725 int f=int(r);
02726 r-=float(f);
02727
02728 if (f>=0 && f<n) {
02729 ret[f]+=v*(1.0f-r);
02730 norm[f]+=(1.0f-r);
02731 if (f<n-1) {
02732 ret[f+1]+=v*r;
02733 norm[f+1]+=r;
02734 }
02735 }
02736 }
02737 }
02738 }
02739 else {
02740 size_t i;
02741 for (z=i=0; z<nz; ++z) {
02742 for (y=0; y<ny; ++y) {
02743 for (x=0; x<nx; x+=step,i+=step) {
02744 float r,v;
02745 if (step==2) {
02746 if (x==0 && z<nz/2) continue;
02747 if (x==0 && z==nz/2 && y<ny/2) continue;
02748 r=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z);
02749 if (!inten) {
02750 #ifdef _WIN32
02751 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02752 #else
02753 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02754 #endif //_WIN32
02755 else v=data[i];
02756 } else {
02757 if (isinten) v=data[i];
02758 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02759 else v=data[i]*data[i];
02760 }
02761 }
02762 else {
02763 r=Util::hypot3(x-nx/2,y-ny/2,z-nz/2);
02764 if (inten) v=data[i]*data[i];
02765 else v=data[i];
02766 }
02767 r=(r-x0)/dx;
02768 int f=int(r);
02769 r-=float(f);
02770 if (f>=0 && f<n) {
02771 ret[f]+=v*(1.0f-r);
02772 norm[f]+=(1.0f-r);
02773 if (f<n-1) {
02774 ret[f+1]+=v*r;
02775 norm[f+1]+=r;
02776 }
02777 }
02778 }
02779 }
02780 }
02781 }
02782
02783 for (i=0; i<n; i++) ret[i]/=norm[i]?norm[i]:1.0f;
02784
02785 EXITFUNC;
02786
02787 return ret;
02788 }
02789
02790 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, int nwedge, float offset, bool inten)
02791 {
02792 ENTERFUNC;
02793
02794 if (nz > 1) {
02795 LOGERR("2D images only.");
02796 throw ImageDimensionException("2D images only");
02797 }
02798 int isinten=get_attr_default("is_intensity",0);
02799
02800 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
02801
02802
02803 vector<float>ret(n*nwedge);
02804 vector<float>norm(n*nwedge);
02805
02806 int x,y,i;
02807 int step=is_complex()?2:1;
02808 float astep=static_cast<float>(M_PI*2.0/nwedge);
02809 if (is_complex()) astep/=2;
02810 float* data = get_data();
02811 for (i=0; i<n*nwedge; i++) ret[i]=norm[i]=0.0;
02812
02813
02814 for (y=i=0; y<ny; y++) {
02815 for (x=0; x<nx; x+=step,i+=step) {
02816 float r,v,a;
02817 int bin;
02818 if (is_complex()) {
02819 #ifdef _WIN32
02820 r=static_cast<float>(_hypot(x/2.0,y<ny/2?y:ny-y));
02821 #else
02822 r=static_cast<float>(hypot(x/2.0,y<ny/2?y:ny-y));
02823 #endif
02824 a=atan2(float(y<ny/2?y:y-ny),x/2.0f);
02825 if (!inten) {
02826 #ifdef _WIN32
02827 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02828 #else
02829 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02830 #endif //_WIN32
02831 else v=data[i];
02832 } else {
02833 if (isinten) v=data[i];
02834 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02835 else v=data[i]*data[i];
02836 }
02837 bin=n*int(floor((a+M_PI/2.0f+offset)/astep));
02838 }
02839 else {
02840 #ifdef _WIN32
02841 r=static_cast<float>(_hypot(x-nx/2,y-ny/2));
02842 #else
02843 r=static_cast<float>(hypot(x-nx/2,y-ny/2));
02844 #endif //_WIN32
02845 a=atan2(float(y-ny/2),float(x-nx/2));
02846 if (inten) v=data[i]*data[i];
02847 else v=data[i];
02848 bin=n*int(floor((a+M_PI+offset)/astep));
02849 }
02850 if (bin>=nwedge*n) bin-=nwedge*n;
02851 if (bin<0) bin+=nwedge*n;
02852 r=(r-x0)/dx;
02853 int f=int(r);
02854 r-=float(f);
02855
02856 if (f>=0 && f<n) {
02857 ret[f+bin]+=v*(1.0f-r);
02858 norm[f+bin]+=(1.0f-r);
02859 if (f<n-1) {
02860 ret[f+1+bin]+=v*r;
02861 norm[f+1+bin]+=r;
02862 }
02863 }
02864 }
02865 }
02866
02867 for (i=0; i<n*nwedge; i++) ret[i]/=norm[i]?norm[i]:1.0f;
02868 EXITFUNC;
02869
02870 return ret;
02871 }
02872
02873 void EMData::cconj() {
02874 ENTERFUNC;
02875 if (!is_complex() || !is_ri())
02876 throw ImageFormatException("EMData::conj requires a complex, ri image");
02877 int nxreal = nx -2 + int(is_fftodd());
02878 int nxhalf = nxreal/2;
02879 for (int iz = 0; iz < nz; iz++)
02880 for (int iy = 0; iy < ny; iy++)
02881 for (int ix = 0; ix <= nxhalf; ix++)
02882 cmplx(ix,iy,iz) = conj(cmplx(ix,iy,iz));
02883 EXITFUNC;
02884 }
02885
02886 void EMData::update_stat() const
02887 {
02888 ENTERFUNC;
02889
02890 if (!(flags & EMDATA_NEEDUPD))
02891 {
02892 EXITFUNC;
02893 return;
02894 }
02895 if (rdata==0) return;
02896
02897 float* data = get_data();
02898 float max = -FLT_MAX;
02899 float min = -max;
02900
02901 double sum = 0;
02902 double square_sum = 0;
02903
02904 int step = 1;
02905 if (is_complex() && !is_ri()) {
02906 step = 2;
02907 }
02908
02909 int n_nonzero = 0;
02910
02911 size_t size = (size_t)nx*ny*nz;
02912 for (size_t i = 0; i < size; i += step) {
02913 float v = data[i];
02914 #ifdef _WIN32
02915 max = _cpp_max(max,v);
02916 min = _cpp_min(min,v);
02917 #else
02918 max=std::max<float>(max,v);
02919 min=std::min<float>(min,v);
02920 #endif //_WIN32
02921 sum += v;
02922 square_sum += v * (double)(v);
02923 if (v != 0) n_nonzero++;
02924 }
02925
02926 size_t n = size / step;
02927 double mean = sum / n;
02928
02929 #ifdef _WIN32
02930 float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / n)/(n-1)));
02931 n_nonzero = _cpp_max(1,n_nonzero);
02932 double sigma_nonzero = std::sqrt( _cpp_max(0,(square_sum - sum*sum/n_nonzero)/(n_nonzero-1)));
02933 #else
02934 float sigma = (float)std::sqrt(std::max<double>(0.0,(square_sum - sum*sum / n)/(n-1)));
02935 n_nonzero = std::max<int>(1,n_nonzero);
02936 double sigma_nonzero = std::sqrt(std::max<double>(0,(square_sum - sum*sum/n_nonzero)/(n_nonzero-1)));
02937 #endif //_WIN32
02938 double mean_nonzero = sum / n_nonzero;
02939
02940 attr_dict["minimum"] = min;
02941 attr_dict["maximum"] = max;
02942 attr_dict["mean"] = (float)(mean);
02943 attr_dict["sigma"] = (float)(sigma);
02944 attr_dict["square_sum"] = (float)(square_sum);
02945 attr_dict["mean_nonzero"] = (float)(mean_nonzero);
02946 attr_dict["sigma_nonzero"] = (float)(sigma_nonzero);
02947 attr_dict["is_complex"] = (int) is_complex();
02948 attr_dict["is_complex_ri"] = (int) is_ri();
02949
02950 flags &= ~EMDATA_NEEDUPD;
02951
02952 if (rot_fp != 0)
02953 {
02954 delete rot_fp; rot_fp = 0;
02955 }
02956
02957 EXITFUNC;
02958
02959 }
02960
02965 bool EMData::operator==(const EMData& that) const {
02966 if(this != &that) {
02967 return false;
02968 }
02969 else {
02970 return true;
02971 }
02972 }
02973
02974 bool EMData::equal(const EMData& that) const {
02975 if (that.get_xsize() != nx || that.get_ysize() != ny || that.get_zsize() != nz ) return false;
02976
02977 const float* d1 = that.get_const_data();
02978 float* d2 = get_data();
02979
02980 for(size_t i =0; i < get_size(); ++i,++d1,++d2) {
02981 if ((*d1) != (*d2)) return false;
02982 }
02983
02984
02985
02986
02987
02988 return true;
02989 }
02990
02991 EMData * EMAN::operator+(const EMData & em, float n)
02992 {
02993 EMData * r = em.copy();
02994 r->add(n);
02995 return r;
02996 }
02997
02998 EMData * EMAN::operator-(const EMData & em, float n)
02999 {
03000 EMData* r = em.copy();
03001 r->sub(n);
03002 return r;
03003 }
03004
03005 EMData * EMAN::operator*(const EMData & em, float n)
03006 {
03007 EMData* r = em.copy();
03008 r ->mult(n);
03009 return r;
03010 }
03011
03012 EMData * EMAN::operator/(const EMData & em, float n)
03013 {
03014 EMData * r = em.copy();
03015 r->div(n);
03016 return r;
03017 }
03018
03019
03020 EMData * EMAN::operator+(float n, const EMData & em)
03021 {
03022 EMData * r = em.copy();
03023 r->add(n);
03024 return r;
03025 }
03026
03027 EMData * EMAN::operator-(float n, const EMData & em)
03028 {
03029 EMData * r = em.copy();
03030 r->mult(-1.0f);
03031 r->add(n);
03032 return r;
03033 }
03034
03035 EMData * EMAN::operator*(float n, const EMData & em)
03036 {
03037 EMData * r = em.copy();
03038 r->mult(n);
03039 return r;
03040 }
03041
03042 EMData * EMAN::operator/(float n, const EMData & em)
03043 {
03044 EMData * r = em.copy();
03045 r->to_one();
03046 r->mult(n);
03047 r->div(em);
03048
03049 return r;
03050 }
03051
03052 EMData * EMAN::rsub(const EMData & em, float n)
03053 {
03054 return EMAN::operator-(n, em);
03055 }
03056
03057 EMData * EMAN::rdiv(const EMData & em, float n)
03058 {
03059 return EMAN::operator/(n, em);
03060 }
03061
03062 EMData * EMAN::operator+(const EMData & a, const EMData & b)
03063 {
03064 EMData * r = a.copy();
03065 r->add(b);
03066 return r;
03067 }
03068
03069 EMData * EMAN::operator-(const EMData & a, const EMData & b)
03070 {
03071 EMData * r = a.copy();
03072 r->sub(b);
03073 return r;
03074 }
03075
03076 EMData * EMAN::operator*(const EMData & a, const EMData & b)
03077 {
03078 EMData * r = a.copy();
03079 r->mult(b);
03080 return r;
03081 }
03082
03083 EMData * EMAN::operator/(const EMData & a, const EMData & b)
03084 {
03085 EMData * r = a.copy();
03086 r->div(b);
03087 return r;
03088 }
03089
03090 void EMData::set_xyz_origin(float origin_x, float origin_y, float origin_z)
03091 {
03092 attr_dict["origin_x"] = origin_x;
03093 attr_dict["origin_y"] = origin_y;
03094 attr_dict["origin_z"] = origin_z;
03095 }
03096
03097 #if 0
03098 void EMData::calc_rcf(EMData * with, vector < float >&sum_array)
03099 {
03100 ENTERFUNC;
03101
03102 int array_size = sum_array.size();
03103 float da = 2 * M_PI / array_size;
03104 float *dat = new float[array_size + 2];
03105 float *dat2 = new float[array_size + 2];
03106 int nx2 = nx * 9 / 20;
03107
03108 float lim = 0;
03109 if (fabs(translation[0]) < fabs(translation[1])) {
03110 lim = fabs(translation[1]);
03111 }
03112 else {
03113 lim = fabs(translation[0]);
03114 }
03115
03116 nx2 -= static_cast < int >(floor(lim));
03117
03118 for (int i = 0; i < array_size; i++) {
03119 sum_array[i] = 0;
03120 }
03121
03122 float sigma = attr_dict["sigma"];
03123 float with_sigma = with->get_attr_dict().get("sigma");
03124
03125 vector<float> vdata, vdata2;
03126 for (int i = 8; i < nx2; i += 6) {
03127 vdata = calc_az_dist(array_size, 0, da, i, i + 6);
03128 vdata2 = with->calc_az_dist(array_size, 0, da, i, i + 6);
03129 Assert(vdata.size() <= array_size + 2);
03130 Assert(cdata2.size() <= array_size + 2);
03131 std::copy(vdata.begin(), vdata.end(), dat);
03132 std::copy(vdata2.begin(), vdata2.end(), dat2);
03133
03134 EMfft::real_to_complex_1d(dat, dat, array_size);
03135 EMfft::real_to_complex_1d(dat2, dat2, array_size);
03136
03137 for (int j = 0; j < array_size + 2; j += 2) {
03138 float max = dat[j] * dat2[j] + dat[j + 1] * dat2[j + 1];
03139 float max2 = dat[j + 1] * dat2[j] - dat2[j + 1] * dat[j];
03140 dat[j] = max;
03141 dat[j + 1] = max2;
03142 }
03143
03144 EMfft::complex_to_real_1d(dat, dat, array_size);
03145 float norm = array_size * array_size * (4.0f * sigma) * (4.0f * with_sigma);
03146
03147 for (int j = 0; j < array_size; j++) {
03148 sum_array[j] += dat[j] * (float) i / norm;
03149 }
03150 }
03151
03152 if( dat )
03153 {
03154 delete[]dat;
03155 dat = 0;
03156 }
03157
03158 if( dat2 )
03159 {
03160 delete[]dat2;
03161 dat2 = 0;
03162 }
03163 EXITFUNC;
03164 }
03165
03166 #endif
03167
03168 void EMData::add_incoherent(EMData * obj)
03169 {
03170 ENTERFUNC;
03171
03172 if (!obj) {
03173 LOGERR("NULL image");
03174 throw NullPointerException("NULL image");
03175 }
03176
03177 if (!obj->is_complex() || !is_complex()) {
03178 throw ImageFormatException("complex images only");
03179 }
03180
03181 if (!EMUtil::is_same_size(this, obj)) {
03182 throw ImageFormatException("images not same size");
03183 }
03184
03185 ri2ap();
03186 obj->ri2ap();
03187
03188 float *dest = get_data();
03189 float *src = obj->get_data();
03190 size_t size = (size_t)nx * ny * nz;
03191 for (size_t j = 0; j < size; j += 2) {
03192 #ifdef _WIN32
03193 dest[j] = (float) _hypot(src[j], dest[j]);
03194 #else
03195 dest[j] = (float) hypot(src[j], dest[j]);
03196 #endif //_WIN32
03197 dest[j + 1] = 0;
03198 }
03199
03200 obj->update();
03201 update();
03202 EXITFUNC;
03203 }
03204
03205
03206 float EMData::calc_dist(EMData * second_img, int y_index) const
03207 {
03208 ENTERFUNC;
03209
03210 if (get_ndim() != 1) {
03211 throw ImageDimensionException("'this' image is 1D only");
03212 }
03213
03214 if (second_img->get_xsize() != nx || ny != 1) {
03215 throw ImageFormatException("image xsize not same");
03216 }
03217
03218 if (y_index > second_img->get_ysize() || y_index < 0) {
03219 return -1;
03220 }
03221
03222 float ret = 0;
03223 float *d1 = get_data();
03224 float *d2 = second_img->get_data() + second_img->get_xsize() * y_index;
03225
03226 for (int i = 0; i < nx; i++) {
03227 ret += Util::square(d1[i] - d2[i]);
03228 }
03229 EXITFUNC;
03230 return std::sqrt(ret);
03231 }
03232
03233
03234 EMData * EMData::calc_fast_sigma_image( EMData* mask)
03235 {
03236 ENTERFUNC;
03237
03238 bool maskflag = false;
03239 if (mask == 0) {
03240 mask = new EMData(nx,ny,nz);
03241 mask->process_inplace("testimage.circlesphere");
03242 maskflag = true;
03243 }
03244
03245 if (get_ndim() != mask->get_ndim() ) throw ImageDimensionException("The dimensions do not match");
03246
03247 int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
03248
03249 if ( mnx > nx || mny > ny || mnz > nz)
03250 throw ImageDimensionException("Can not calculate variance map using an image that is larger than this image");
03251
03252 size_t P = 0;
03253 for(size_t i = 0; i < mask->get_size(); ++i){
03254 if (mask->get_value_at(i) != 0){
03255 ++P;
03256 }
03257 }
03258 float normfac = 1.0f/(float)P;
03259
03260
03261
03262 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03263
03264 Region r;
03265 if (ny == 1) r = Region((mnx-nxc)/2,nxc);
03266 else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03267 else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03268 mask->clip_inplace(r,0.0);
03269
03270
03271
03272
03273
03274
03275 Region r2;
03276 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03277 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03278 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03279 EMData* squared = get_clip(r2,get_edge_mean());
03280
03281 EMData* tmp = squared->copy();
03282 Dict pow;
03283 pow["pow"] = 2.0f;
03284 squared->process_inplace("math.pow",pow);
03285 EMData* s = mask->convolute(squared);
03286 squared->mult(normfac);
03287
03288 EMData* m = mask->convolute(tmp);
03289 m->mult(normfac);
03290 m->process_inplace("math.pow",pow);
03291 delete tmp; tmp = 0;
03292 s->sub(*m);
03293
03294 s->process_inplace("math.sqrt");
03295
03296
03297
03298
03299
03300
03301 if (maskflag) {
03302 delete mask;
03303 mask = 0;
03304 } else {
03305 Region r;
03306 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
03307 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
03308 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
03309 mask->clip_inplace(r);
03310 }
03311
03312 delete squared;
03313 delete m;
03314
03315 s->process_inplace("xform.phaseorigin.tocenter");
03316 Region r3;
03317 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03318 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03319 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03320 s->clip_inplace(r3);
03321 EXITFUNC;
03322 return s;
03323 }
03324
03325
03326
03327 EMData *EMData::calc_flcf(EMData * with)
03328 {
03329 ENTERFUNC;
03330 EMData *this_copy=this;
03331 this_copy=copy();
03332
03333 int mnx = with->get_xsize(); int mny = with->get_ysize(); int mnz = with->get_zsize();
03334 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03335
03336
03337 EMData* ones = new EMData(mnx,mny,mnz);
03338 ones->process_inplace("testimage.circlesphere");
03339
03340
03341 EMData* with_resized = with->copy();
03342 with_resized->process_inplace("normalize");
03343 with_resized->mult(*ones);
03344
03345 EMData* s = calc_fast_sigma_image(ones);
03346
03347 Region r1;
03348 if (ny == 1) r1 = Region((mnx-nxc)/2,nxc);
03349 else if (nz == 1) r1 = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03350 else r1 = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03351 with_resized->clip_inplace(r1,0.0);
03352
03353 Region r2;
03354 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03355 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03356 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03357 this_copy->clip_inplace(r2,0.0);
03358
03359 EMData* corr = this_copy->calc_ccf(with_resized);
03360
03361 corr->process_inplace("xform.phaseorigin.tocenter");
03362 Region r3;
03363 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03364 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03365 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03366 corr->clip_inplace(r3);
03367
03368 corr->div(*s);
03369
03370 delete with_resized; delete ones; delete this_copy; delete s;
03371 EXITFUNC;
03372 return corr;
03373 }
03374
03375 EMData *EMData::convolute(EMData * with)
03376 {
03377 ENTERFUNC;
03378
03379 EMData *f1 = do_fft();
03380 if (!f1) {
03381 LOGERR("FFT returns NULL image");
03382 throw NullPointerException("FFT returns NULL image");
03383 }
03384
03385 f1->ap2ri();
03386
03387 EMData *cf = 0;
03388 if (with) {
03389 cf = with->do_fft();
03390 if (!cf) {
03391 LOGERR("FFT returns NULL image");
03392 throw NullPointerException("FFT returns NULL image");
03393 }
03394 cf->ap2ri();
03395 }
03396 else {
03397 cf = f1->copy();
03398 }
03399
03400 if (with && !EMUtil::is_same_size(f1, cf)) {
03401 LOGERR("images not same size");
03402 throw ImageFormatException("images not same size");
03403 }
03404
03405 float *rdata1 = f1->get_data();
03406 float *rdata2 = cf->get_data();
03407 size_t cf_size = (size_t)cf->get_xsize() * cf->get_ysize() * cf->get_zsize();
03408
03409 float re,im;
03410
03411 for (size_t i = 0; i < cf_size; i += 2) {
03412 re = rdata1[i] * rdata2[i] - rdata1[i + 1] * rdata2[i + 1];
03413 im = rdata1[i + 1] * rdata2[i] + rdata1[i] * rdata2[i + 1];
03414 rdata2[i]=re;
03415 rdata2[i+1]=im;
03416 }
03417 cf->update();
03418 EMData *f2 = cf->do_ift();
03419
03420 if( cf )
03421 {
03422 delete cf;
03423 cf = 0;
03424 }
03425
03426 if( f1 )
03427 {
03428 delete f1;
03429 f1=0;
03430 }
03431
03432 EXITFUNC;
03433 return f2;
03434 }
03435
03436
03437 void EMData::common_lines(EMData * image1, EMData * image2,
03438 int mode, int steps, bool horizontal)
03439 {
03440 ENTERFUNC;
03441
03442 if (!image1 || !image2) {
03443 throw NullPointerException("NULL image");
03444 }
03445
03446 if (mode < 0 || mode > 2) {
03447 throw OutofRangeException(0, 2, mode, "invalid mode");
03448 }
03449
03450 if (!image1->is_complex()) {
03451 image1 = image1->do_fft();
03452 }
03453 if (!image2->is_complex()) {
03454 image2 = image2->do_fft();
03455 }
03456
03457 image1->ap2ri();
03458 image2->ap2ri();
03459
03460 if (!EMUtil::is_same_size(image1, image2)) {
03461 throw ImageFormatException("images not same sizes");
03462 }
03463
03464 int image2_nx = image2->get_xsize();
03465 int image2_ny = image2->get_ysize();
03466
03467 int rmax = image2_ny / 4 - 1;
03468 int array_size = steps * rmax * 2;
03469 float *im1 = new float[array_size];
03470 float *im2 = new float[array_size];
03471 for (int i = 0; i < array_size; i++) {
03472 im1[i] = 0;
03473 im2[i] = 0;
03474 }
03475
03476 set_size(steps * 2, steps * 2, 1);
03477
03478 float *image1_data = image1->get_data();
03479 float *image2_data = image2->get_data();
03480
03481 float da = M_PI / steps;
03482 float a = -M_PI / 2.0f + da / 2.0f;
03483 int jmax = 0;
03484
03485 for (int i = 0; i < steps * 2; i += 2, a += da) {
03486 float s1 = 0;
03487 float s2 = 0;
03488 int i2 = i * rmax;
03489 int j = 0;
03490
03491 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
03492 float x = r * cos(a);
03493 float y = r * sin(a);
03494
03495 if (x < 0) {
03496 x = -x;
03497 y = -y;
03498 LOGERR("CCL ERROR %d, %f !\n", i, -x);
03499 }
03500
03501 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
03502 int l = i2 + j;
03503 float x2 = x - floor(x);
03504 float y2 = y - floor(y);
03505
03506 im1[l] = Util::bilinear_interpolate(image1_data[k],
03507 image1_data[k + 2],
03508 image1_data[k + image2_nx],
03509 image1_data[k + 2 + image2_nx], x2, y2);
03510
03511 im2[l] = Util::bilinear_interpolate(image2_data[k],
03512 image2_data[k + 2],
03513 image2_data[k + image2_nx],
03514 image2_data[k + 2 + image2_nx], x2, y2);
03515
03516 k++;
03517
03518 im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
03519 image1_data[k + 2],
03520 image1_data[k + image2_nx],
03521 image1_data[k + 2 + image2_nx], x2, y2);
03522
03523 im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
03524 image2_data[k + 2],
03525 image2_data[k + image2_nx],
03526 image2_data[k + 2 + image2_nx], x2, y2);
03527
03528 s1 += Util::square_sum(im1[l], im1[l + 1]);
03529 s2 += Util::square_sum(im2[l], im2[l + 1]);
03530 }
03531
03532 jmax = j - 1;
03533 float sqrt_s1 = std::sqrt(s1);
03534 float sqrt_s2 = std::sqrt(s2);
03535
03536 int l = 0;
03537 for (float r = 1; r < rmax; r += 1.0f) {
03538 int i3 = i2 + l;
03539 im1[i3] /= sqrt_s1;
03540 im1[i3 + 1] /= sqrt_s1;
03541 im2[i3] /= sqrt_s2;
03542 im2[i3 + 1] /= sqrt_s2;
03543 l += 2;
03544 }
03545 }
03546 float * data = get_data();
03547
03548 switch (mode) {
03549 case 0:
03550 for (int m1 = 0; m1 < 2; m1++) {
03551 for (int m2 = 0; m2 < 2; m2++) {
03552
03553 if (m1 == 0 && m2 == 0) {
03554 for (int i = 0; i < steps; i++) {
03555 int i2 = i * rmax * 2;
03556 for (int j = 0; j < steps; j++) {
03557 int l = i + j * steps * 2;
03558 int j2 = j * rmax * 2;
03559 data[l] = 0;
03560 for (int k = 0; k < jmax; k++) {
03561 data[l] += im1[i2 + k] * im2[j2 + k];
03562 }
03563 }
03564 }
03565 }
03566 else {
03567 int steps2 = steps * m2 + steps * steps * 2 * m1;
03568
03569 for (int i = 0; i < steps; i++) {
03570 int i2 = i * rmax * 2;
03571 for (int j = 0; j < steps; j++) {
03572 int j2 = j * rmax * 2;
03573 int l = i + j * steps * 2 + steps2;
03574 data[l] = 0;
03575
03576 for (int k = 0; k < jmax; k += 2) {
03577 i2 += k;
03578 j2 += k;
03579 data[l] += im1[i2] * im2[j2];
03580 data[l] += -im1[i2 + 1] * im2[j2 + 1];
03581 }
03582 }
03583 }
03584 }
03585 }
03586 }
03587
03588 break;
03589 case 1:
03590 for (int m1 = 0; m1 < 2; m1++) {
03591 for (int m2 = 0; m2 < 2; m2++) {
03592 int steps2 = steps * m2 + steps * steps * 2 * m1;
03593 int p1_sign = 1;
03594 if (m1 != m2) {
03595 p1_sign = -1;
03596 }
03597
03598 for (int i = 0; i < steps; i++) {
03599 int i2 = i * rmax * 2;
03600
03601 for (int j = 0; j < steps; j++) {
03602 int j2 = j * rmax * 2;
03603
03604 int l = i + j * steps * 2 + steps2;
03605 data[l] = 0;
03606 float a = 0;
03607
03608 for (int k = 0; k < jmax; k += 2) {
03609 i2 += k;
03610 j2 += k;
03611
03612 #ifdef _WIN32
03613 float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
03614 #else
03615 float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
03616 #endif //_WIN32
03617 float p1 = atan2(im1[i2 + 1], im1[i2]);
03618 float p2 = atan2(im2[j2 + 1], im2[j2]);
03619
03620 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
03621 a += a1;
03622 }
03623
03624 data[l] /= (float)(a * M_PI / 180.0f);
03625 }
03626 }
03627 }
03628 }
03629
03630 break;
03631 case 2:
03632 for (int m1 = 0; m1 < 2; m1++) {
03633 for (int m2 = 0; m2 < 2; m2++) {
03634 int steps2 = steps * m2 + steps * steps * 2 * m1;
03635
03636 for (int i = 0; i < steps; i++) {
03637 int i2 = i * rmax * 2;
03638
03639 for (int j = 0; j < steps; j++) {
03640 int j2 = j * rmax * 2;
03641 int l = i + j * steps * 2 + steps2;
03642 data[l] = 0;
03643
03644 for (int k = 0; k < jmax; k += 2) {
03645 i2 += k;
03646 j2 += k;
03647 #ifdef _WIN32
03648 data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
03649 #else
03650 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
03651 #endif //_WIN32
03652 }
03653 }
03654 }
03655 }
03656 }
03657
03658 break;
03659 default:
03660 break;
03661 }
03662
03663 if (horizontal) {
03664 float *tmp_array = new float[ny];
03665 for (int i = 1; i < nx; i++) {
03666 for (int j = 0; j < ny; j++) {
03667 tmp_array[j] = get_value_at(i, j);
03668 }
03669 for (int j = 0; j < ny; j++) {
03670 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
03671 }
03672 }
03673 if( tmp_array )
03674 {
03675 delete[]tmp_array;
03676 tmp_array = 0;
03677 }
03678 }
03679
03680 if( im1 )
03681 {
03682 delete[]im1;
03683 im1 = 0;
03684 }
03685
03686 if( im2 )
03687 {
03688 delete im2;
03689 im2 = 0;
03690 }
03691
03692
03693 image1->update();
03694 image2->update();
03695 if( image1 )
03696 {
03697 delete image1;
03698 image1 = 0;
03699 }
03700 if( image2 )
03701 {
03702 delete image2;
03703 image2 = 0;
03704 }
03705 update();
03706 EXITFUNC;
03707 }
03708
03709
03710
03711 void EMData::common_lines_real(EMData * image1, EMData * image2,
03712 int steps, bool horiz)
03713 {
03714 ENTERFUNC;
03715
03716 if (!image1 || !image2) {
03717 throw NullPointerException("NULL image");
03718 }
03719
03720 if (!EMUtil::is_same_size(image1, image2)) {
03721 throw ImageFormatException("images not same size");
03722 }
03723
03724 int steps2 = steps * 2;
03725 int image_ny = image1->get_ysize();
03726 EMData *image1_copy = image1->copy();
03727 EMData *image2_copy = image2->copy();
03728
03729 float *im1 = new float[steps2 * image_ny];
03730 float *im2 = new float[steps2 * image_ny];
03731
03732 EMData *images[] = { image1_copy, image2_copy };
03733 float *ims[] = { im1, im2 };
03734
03735 for (int m = 0; m < 2; m++) {
03736 float *im = ims[m];
03737 float a = M_PI / steps2;
03738 Transform t(Dict("type","2d","alpha",-a));
03739 for (int i = 0; i < steps2; i++) {
03740 images[i]->transform(t);
03741 float *data = images[i]->get_data();
03742
03743 for (int j = 0; j < image_ny; j++) {
03744 float sum = 0;
03745 for (int k = 0; k < image_ny; k++) {
03746 sum += data[j * image_ny + k];
03747 }
03748 im[i * image_ny + j] = sum;
03749 }
03750
03751 float sum1 = 0;
03752 float sum2 = 0;
03753 for (int j = 0; j < image_ny; j++) {
03754 int l = i * image_ny + j;
03755 sum1 += im[l];
03756 sum2 += im[l] * im[l];
03757 }
03758
03759 float mean = sum1 / image_ny;
03760 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03761
03762 for (int j = 0; j < image_ny; j++) {
03763 int l = i * image_ny + j;
03764 im[l] = (im[l] - mean) / sigma;
03765 }
03766
03767 images[i]->update();
03768 a += M_PI / steps;
03769 }
03770 }
03771
03772 set_size(steps2, steps2, 1);
03773 float *data1 = get_data();
03774
03775 if (horiz) {
03776 for (int i = 0; i < steps2; i++) {
03777 for (int j = 0, l = i; j < steps2; j++, l++) {
03778 if (l == steps2) {
03779 l = 0;
03780 }
03781
03782 float sum = 0;
03783 for (int k = 0; k < image_ny; k++) {
03784 sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03785 }
03786 data1[i + j * steps2] = sum;
03787 }
03788 }
03789 }
03790 else {
03791 for (int i = 0; i < steps2; i++) {
03792 for (int j = 0; j < steps2; j++) {
03793 float sum = 0;
03794 for (int k = 0; k < image_ny; k++) {
03795 sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03796 }
03797 data1[i + j * steps2] = sum;
03798 }
03799 }
03800 }
03801
03802 update();
03803
03804 if( image1_copy )
03805 {
03806 delete image1_copy;
03807 image1_copy = 0;
03808 }
03809
03810 if( image2_copy )
03811 {
03812 delete image2_copy;
03813 image2_copy = 0;
03814 }
03815
03816 if( im1 )
03817 {
03818 delete[]im1;
03819 im1 = 0;
03820 }
03821
03822 if( im2 )
03823 {
03824 delete[]im2;
03825 im2 = 0;
03826 }
03827 EXITFUNC;
03828 }
03829
03830
03831 void EMData::cut_slice(const EMData *const map, const Transform& transform, bool interpolate)
03832 {
03833 ENTERFUNC;
03834
03835 if (!map) throw NullPointerException("NULL image");
03836
03837 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
03838 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
03839
03840 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
03841 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
03842
03843
03844 float *sdata = map->get_data();
03845 float *ddata = get_data();
03846
03847 int map_nx = map->get_xsize();
03848 int map_ny = map->get_ysize();
03849 int map_nz = map->get_zsize();
03850 int map_nxy = map_nx * map_ny;
03851
03852 int ymax = ny/2;
03853 if ( ny % 2 == 1 ) ymax += 1;
03854 int xmax = nx/2;
03855 if ( nx % 2 == 1 ) xmax += 1;
03856 for (int y = -ny/2; y < ymax; y++) {
03857 for (int x = -nx/2; x < xmax; x++) {
03858 Vec3f coord(x,y,0);
03859 Vec3f soln = transform*coord;
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870 float xx = soln[0]+map_nx/2;
03871 float yy = soln[1]+map_ny/2;
03872 float zz = soln[2]+map_nz/2;
03873
03874 int l = (x+nx/2) + (y+ny/2) * nx;
03875
03876 float t = xx - floor(xx);
03877 float u = yy - floor(yy);
03878 float v = zz - floor(zz);
03879
03880 if (xx < 0 || yy < 0 || zz < 0 ) {
03881 ddata[l] = 0;
03882 continue;
03883 }
03884 if (interpolate) {
03885 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03886 ddata[l] = 0;
03887 continue;
03888 }
03889 int k = (int) (Util::fast_floor(xx) + Util::fast_floor(yy) * map_nx + Util::fast_floor(zz) * map_nxy);
03890
03891
03892 if (xx < (map_nx - 1) && yy < (map_ny - 1) && zz < (map_nz - 1)) {
03893 ddata[l] = Util::trilinear_interpolate(sdata[k],
03894 sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
03895 sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
03896 sdata[k + map_nx + map_nxy + 1],t, u, v);
03897 }
03898 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03899 ddata[l] += sdata[k];
03900 }
03901 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) ) {
03902 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nxy],v);
03903 }
03904 else if ( xx == (map_nx - 1) && zz == (map_nz - 1) ) {
03905 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nx],u);
03906 }
03907 else if ( yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03908 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + 1],t);
03909 }
03910 else if ( xx == (map_nx - 1) ) {
03911 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + map_nx], sdata[k + map_nxy], sdata[k + map_nxy + map_nx],u,v);
03912 }
03913 else if ( yy == (map_ny - 1) ) {
03914 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nxy], sdata[k + map_nxy + 1],t,v);
03915 }
03916 else if ( zz == (map_nz - 1) ) {
03917 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nx], sdata[k + map_nx + 1],t,u);
03918 }
03919
03920
03921
03922
03923
03924
03925
03926
03927
03928
03929
03930
03931 }
03932 else {
03933 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03934 ddata[l] = 0;
03935 continue;
03936 }
03937 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
03938 ddata[l] = sdata[k];
03939 }
03940
03941 }
03942 }
03943
03944 update();
03945
03946 EXITFUNC;
03947 }
03948
03949 EMData *EMData::unwrap_largerR(int r1,int r2,int xs, float rmax_f) {
03950 float *d,*dd;
03951 int do360=2;
03952 int rmax = (int)(rmax_f+0.5f);
03953 unsigned long i;
03954 unsigned int nvox=get_xsize()*get_ysize();
03955 float maxmap=0.0f, minmap=0.0f;
03956 float temp=0.0f, diff_den=0.0f, norm=0.0f;
03957 float cut_off_va =0.0f;
03958
03959 d=get_data();
03960 maxmap=-1000000.0f;
03961 minmap=1000000.0f;
03962 for (i=0;i<nvox;i++){
03963 if(d[i]>maxmap) maxmap=d[i];
03964 if(d[i]<minmap) minmap=d[i];
03965 }
03966 diff_den = maxmap-minmap;
03967 for(i=0;i<nvox;i++) {
03968 temp = (d[i]-minmap)/diff_den;
03969 if(cut_off_va != 0.0) {
03970 if(temp < cut_off_va)
03971 d[i] = 0.0f;
03972 else
03973 d[i] = temp-cut_off_va;
03974 }
03975 else d[i] = temp;
03976 }
03977
03978 for(i=0;i<nvox;i++) {
03979 temp=d[i];
03980 norm += temp*temp;
03981 }
03982 for(i=0;i<nvox;i++) d[i] /= norm;
03983
03984 if (xs<1) {
03985 xs = (int) floor(do360*M_PI*get_ysize()/4);
03986 xs=Util::calc_best_fft_size(xs);
03987 }
03988 if (r1<0) r1=0;
03989 float maxext=ceil(0.6f*std::sqrt((float)(get_xsize()*get_xsize()+get_ysize()*get_ysize())));
03990
03991 if (r2<r1) r2=(int)maxext;
03992 EMData *ret = new EMData;
03993
03994 ret->set_size(xs,r2+1,1);
03995
03996 dd=ret->get_data();
03997
03998 for (int i=0; i<xs; i++) {
03999 float si=sin(i*M_PI*2/xs);
04000 float co=cos(i*M_PI*2/xs);
04001 for (int r=0; r<=maxext; r++) {
04002 float x=(r+r1)*co+get_xsize()/2;
04003 float y=(r+r1)*si+get_ysize()/2;
04004 if(x<0.0 || x>=get_xsize()-1.0 || y<0.0 || y>=get_ysize()-1.0 || r>rmax){
04005 for(;r<=r2;r++)
04006 dd[i+r*xs]=0.0;
04007 break;
04008 }
04009 int x1=(int)floor(x);
04010 int y1=(int)floor(y);
04011 float t=x-x1;
04012 float u=y-y1;
04013 float f11= d[x1+y1*get_xsize()];
04014 float f21= d[(x1+1)+y1*get_xsize()];
04015 float f12= d[x1+(y1+1)*get_xsize()];
04016 float f22= d[(x1+1)+(y1+1)*get_xsize()];
04017 dd[i+r*xs] = (1-t)*(1-u)*f11+t*(1-u)*f21+t*u*f22+(1-t)*u*f12;
04018 }
04019 }
04020 update();
04021 ret->update();
04022 return ret;
04023 }
04024
04025
04026 EMData *EMData::oneDfftPolar(int size, float rmax, float MAXR){
04027 float *pcs=get_data();
04028 EMData *imagepcsfft = new EMData;
04029 imagepcsfft->set_size((size+2), (int)MAXR+1, 1);
04030 float *d=imagepcsfft->get_data();
04031
04032 EMData *data_in=new EMData;
04033 data_in->set_size(size,1,1);
04034 float *in=data_in->get_data();
04035
04036 for(int row=0; row<=(int)MAXR; ++row){
04037 if(row<=(int)rmax) {
04038 for(int i=0; i<size;++i) in[i] = pcs[i+row*size];
04039 data_in->set_complex(false);
04040 data_in->do_fft_inplace();
04041 for(int j=0;j<size+2;j++) d[j+row*(size+2)]=in[j];
04042 }
04043 else for(int j=0;j<size+2;j++) d[j+row*(size+2)]=0.0;
04044 }
04045 imagepcsfft->update();
04046 delete data_in;
04047 return imagepcsfft;
04048 }
04049
04050 void EMData::uncut_slice(EMData * const map, const Transform& transform) const
04051 {
04052 ENTERFUNC;
04053
04054 if (!map) throw NullPointerException("NULL image");
04055
04056 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
04057 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
04058
04059 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
04060 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
04061
04062
04063
04064
04065
04066
04067 float *ddata = map->get_data();
04068 float *sdata = get_data();
04069
04070 int map_nx = map->get_xsize();
04071 int map_ny = map->get_ysize();
04072 int map_nz = map->get_zsize();
04073 int map_nxy = map_nx * map_ny;
04074 float map_nz_round_limit = (float) map_nz-0.5f;
04075 float map_ny_round_limit = (float) map_ny-0.5f;
04076 float map_nx_round_limit = (float) map_nx-0.5f;
04077
04078
04079
04080
04081 int ymax = ny/2;
04082 if ( ny % 2 == 1 ) ymax += 1;
04083 int xmax = nx/2;
04084 if ( nx % 2 == 1 ) xmax += 1;
04085 for (int y = -ny/2; y < ymax; y++) {
04086 for (int x = -nx/2; x < xmax; x++) {
04087 Vec3f coord(x,y,0);
04088 Vec3f soln = transform*coord;
04089
04090
04091
04092
04093
04094
04095
04096
04097 float xx = soln[0]+map_nx/2;
04098 float yy = soln[1]+map_ny/2;
04099 float zz = soln[2]+map_nz/2;
04100
04101
04102 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;
04103
04104 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
04105 int l = (x+nx/2) + (y+ny/2) * nx;
04106 ddata[k] = sdata[l];
04107 }
04108 }
04109
04110 map->update();
04111 EXITFUNC;
04112 }
04113
04114 EMData *EMData::extract_box(const Transform& cs, const Region& r)
04115 {
04116 vector<float> cs_matrix = cs.get_matrix();
04117
04118 EMData* box = new EMData();
04119 box->set_size((r.get_width()-r.x_origin()), (r.get_height()-r.y_origin()), (r.get_depth()-r.z_origin()));
04120 int box_nx = box->get_xsize();
04121 int box_ny = box->get_ysize();
04122 int box_nxy = box_nx*box_ny;
04123 float* bdata = box->get_data();
04124 float* ddata = get_data();
04125
04126 for (int x = r.x_origin(); x < r.get_width(); x++) {
04127 for (int y = r.y_origin(); y < r.get_height(); y++) {
04128 for (int z = r.z_origin(); z < r.get_depth(); z++) {
04129
04130
04131
04132 float xb = cs_matrix[0]*x + y*cs_matrix[4] + z*cs_matrix[8] + cs_matrix[3];
04133 float yb = cs_matrix[1]*x + y*cs_matrix[5] + z*cs_matrix[9] + cs_matrix[7];
04134 float zb = cs_matrix[2]*x + y*cs_matrix[6] + z*cs_matrix[10] + cs_matrix[11];
04135 float t = xb - Util::fast_floor(xb);
04136 float u = yb - Util::fast_floor(yb);
04137 float v = zb - Util::fast_floor(zb);
04138
04139
04140 int l = (x - r.x_origin()) + (y - r.y_origin())*box_nx + (z - r.z_origin())*box_nxy;
04141 int k = (int) (Util::fast_floor(xb) + Util::fast_floor(yb) * nx + Util::fast_floor(zb) * nxy);
04142
04143 if ( xb > nx - 1 || yb > ny - 1 || zb > nz - 1) {
04144 bdata[l] = 0;
04145 continue;
04146 }
04147 if (xb < 0 || yb < 0 || zb < 0){
04148 bdata[l] = 0;
04149 continue;
04150 }
04151
04152 if (xb < (nx - 1) && yb < (ny - 1) && zb < (nz - 1)) {
04153 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);
04154 }
04155 }
04156 }
04157 }
04158
04159 return box;
04160 }
04161
04162 void EMData::save_byteorder_to_dict(ImageIO * imageio)
04163 {
04164 string image_endian = "ImageEndian";
04165 string host_endian = "HostEndian";
04166
04167 if (imageio->is_image_big_endian()) {
04168 attr_dict[image_endian] = "big";
04169 }
04170 else {
04171 attr_dict[image_endian] = "little";
04172 }
04173
04174 if (ByteOrder::is_host_big_endian()) {
04175 attr_dict[host_endian] = "big";
04176 }
04177 else {
04178 attr_dict[host_endian] = "little";
04179 }
04180 }
04181
04182 EMData* EMData::compute_missingwedge(float wedgeangle, float start, float stop)
04183 {
04184 EMData* test = new EMData();
04185 test->set_size(nx,ny,nz);
04186
04187 float ratio = tan((90.0f-wedgeangle)*M_PI/180.0f);
04188
04189 int offset_i = 2*int(start*nz/2);
04190 int offset_f = int(stop*nz/2);
04191
04192 int step = 0;
04193 float sum = 0.0;
04194 double square_sum = 0.0;
04195 for (int j = 0; j < offset_f; j++){
04196 for (int k = offset_i; k < offset_f; k++) {
04197 for (int i = 0; i < nx; i+=2) {
04198 if (i < int(k*ratio)) {
04199 test->set_value_at(i, j, k, 1.0);
04200 float v = std::sqrt(pow(get_value_at_wrap(i, j, k),2) + pow(get_value_at_wrap(i+1, j, k),2));
04201 sum += v;
04202 square_sum += v * (double)(v);
04203 step++;
04204 }
04205 }
04206 }
04207 }
04208
04209 float mean = sum / step;
04210
04211 #ifdef _WIN32
04212 float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*mean)/(step-1)));
04213 #else
04214 float sigma = (float)std::sqrt(std::max<double>(0.0,(square_sum - sum*mean)/(step-1)));
04215 #endif //_WIN32
04216
04217 cout << "Mean sqr wedge amp " << mean << " Sigma Squ wedge Amp " << sigma << endl;
04218 set_attr("spt_wedge_mean", mean);
04219 set_attr("spt_wedge_sigma", sigma);
04220
04221 return test;
04222 }
04223
04224
04225
04226
04227
04228