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),
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),
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),
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 (num_bytes != 0 && that.cudarwdata != 0) {
00150
00151 rw_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 cout << "That copy constructor #2" << endl;
00199 if (num_bytes != 0 && that.cudarwdata != 0) {
00200 rw_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 }
00205 #endif //EMAN2_USING_CUDA
00206
00207 changecount = that.changecount;
00208
00209 if (that.rot_fp != 0) rot_fp = new EMData(*(that.rot_fp));
00210 else rot_fp = 0;
00211 }
00212 EXITFUNC;
00213 return *this;
00214 }
00215
00216 EMData::EMData(int nx, int ny, int nz, bool is_real) :
00217 #ifdef EMAN2_USING_CUDA
00218 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0),
00219 #endif
00220 attr_dict(), rdata(0), supp(0), flags(0), changecount(0), nx(0), ny(0), nz(0), nxy(0), nxyz(0), xoff(0), yoff(0), zoff(0),
00221 all_translation(), path(""), pathnum(0), rot_fp(0)
00222 {
00223 ENTERFUNC;
00224
00225
00226 attr_dict["apix_x"] = 1.0f;
00227 attr_dict["apix_y"] = 1.0f;
00228 attr_dict["apix_z"] = 1.0f;
00229
00230 if(is_real) {
00231 attr_dict["is_complex"] = int(0);
00232 attr_dict["is_complex_x"] = int(0);
00233 attr_dict["is_complex_ri"] = int(1);
00234 set_size(nx, ny, nz);
00235 }
00236 else {
00237 int new_nx = nx + 2 - nx%2;
00238 set_size(new_nx, ny, nz);
00239
00240 attr_dict["is_complex"] = int(1);
00241
00242 if(ny==1 && nz ==1) {
00243 attr_dict["is_complex_x"] = int(1);
00244 }
00245 else {
00246 attr_dict["is_complex_x"] = int(0);
00247 }
00248
00249 attr_dict["is_complex_ri"] = int(1);
00250 attr_dict["is_fftpad"] = int(1);
00251
00252 if(nx%2 == 1) {
00253 attr_dict["is_fftodd"] = 1;
00254 }
00255 }
00256
00257 this->to_zero();
00258 update();
00259 EMData::totalalloc++;
00260 #ifdef MEMDEBUG2
00261 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
00262 #endif
00263
00264 EXITFUNC;
00265 }
00266
00267
00268 EMData::EMData(float* data, const int x, const int y, const int z, const Dict& attr_dict) :
00269 #ifdef EMAN2_USING_CUDA
00270 cudarwdata(0), cudarodata(0), num_bytes(0), nextlistitem(0), prevlistitem(0), roneedsupdate(0),
00271 #endif
00272 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
00273 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
00274 {
00275 ENTERFUNC;
00276
00277 attr_dict["apix_x"] = 1.0f;
00278 attr_dict["apix_y"] = 1.0f;
00279 attr_dict["apix_z"] = 1.0f;
00280
00281 EMData::totalalloc++;
00282 #ifdef MEMDEBUG2
00283 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
00284 #endif
00285
00286 update();
00287 EXITFUNC;
00288 }
00289
00290 #ifdef EMAN2_USING_CUDA
00291
00292 EMData::EMData(float* data, float* cudadata, const int x, const int y, const int z, const Dict& attr_dict) :
00293 cudarwdata(cudadata), cudarodata(0), num_bytes(x*y*z*sizeof(float)), nextlistitem(0), prevlistitem(0), roneedsupdate(0),
00294 attr_dict(attr_dict), rdata(data), supp(0), flags(0), changecount(0), nx(x), ny(y), nz(z), nxy(x*y), nxyz((size_t)x*y*z), xoff(0),
00295 yoff(0), zoff(0), all_translation(), path(""), pathnum(0), rot_fp(0)
00296 {
00297 ENTERFUNC;
00298
00299
00300 attr_dict["apix_x"] = 1.0f;
00301 attr_dict["apix_y"] = 1.0f;
00302 attr_dict["apix_z"] = 1.0f;
00303
00304 EMData::totalalloc++;
00305 #ifdef MEMDEBUG2
00306 printf("EMDATA+ %4d %p\n",EMData::totalalloc,this);
00307 #endif
00308
00309 update();
00310 EXITFUNC;
00311 }
00312
00313 #endif //EMAN2_USING_CUDA
00314
00315
00316 using std::cout;
00317 using std::endl;
00318 EMData::~EMData()
00319 {
00320 ENTERFUNC;
00321 free_memory();
00322
00323 #ifdef EMAN2_USING_CUDA
00324 if(cudarwdata){rw_free();}
00325 if(cudarodata){ro_free();}
00326 #endif // EMAN2_USING_CUDA
00327 EMData::totalalloc--;
00328 #ifdef MEMDEBUG2
00329 printf("EMDATA- %4d %p\n",EMData::totalalloc,this);
00330 #endif
00331 EXITFUNC;
00332 }
00333
00334 void EMData::clip_inplace(const Region & area,const float& fill_value)
00335 {
00336
00337 ENTERFUNC;
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("math.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 Vec3f v = Vec3f(r0*cos((float)i), r0*sin((float)i), 0);
00953 Vec3f d = t*v-v;
00954 float dd = d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
00955 if (dd > ddmax) ddmax = dd;
00956 }
00957 return std::sqrt(ddmax);
00958 }
00959
00960 void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy, float dz)
00961 {
00962 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
00963
00964 Transform t;
00965 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
00966 t.set_trans(dx, dy, dz);
00967 rotate_translate(t);
00968 }
00969
00970
00971 void EMData::rotate_translate(float az, float alt, float phi, float dx, float dy,
00972 float dz, float pdx, float pdy, float pdz)
00973 {
00974 cout << "Deprecation warning in EMData::rotate_translate. Please consider using EMData::transform() instead " << endl;
00975
00976
00977
00978 Transform t;
00979 t.set_pre_trans(Vec3f(dx, dy, dz));
00980 t.set_rotation(Dict("type", "eman", "az", az, "alt", alt, "phi", phi));
00981 t.set_trans(pdx, pdy, pdz);
00982 rotate_translate(t);
00983 }
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00999
01001
01002
01003
01004
01006
01007
01008
01009
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
01071
01072
01073
01074
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 void EMData::rotate_x(int dx)
01175 {
01176 ENTERFUNC;
01177
01178 if (get_ndim() > 2) {
01179 throw ImageDimensionException("no 3D image");
01180 }
01181
01182
01183 size_t row_size = nx * sizeof(float);
01184 float *tmp = (float*)EMUtil::em_malloc(row_size);
01185 float * data = get_data();
01186
01187 for (int y = 0; y < ny; y++) {
01188 int y_nx = y * nx;
01189 for (int x = 0; x < nx; x++) {
01190 tmp[x] = data[y_nx + (x + dx) % nx];
01191 }
01192 EMUtil::em_memcpy(&data[y_nx], tmp, row_size);
01193 }
01194
01195 update();
01196 if( tmp )
01197 {
01198 delete[]tmp;
01199 tmp = 0;
01200 }
01201 EXITFUNC;
01202 }
01203
01204 double EMData::dot_rotate_translate(EMData * with, float dx, float dy, float da, const bool mirror)
01205 {
01206 ENTERFUNC;
01207
01208 if (!EMUtil::is_same_size(this, with)) {
01209 LOGERR("images not same size");
01210 throw ImageFormatException("images not same size");
01211 }
01212
01213 if (get_ndim() == 3) {
01214 LOGERR("1D/2D Images only");
01215 throw ImageDimensionException("1D/2D only");
01216 }
01217
01218 float *this_data = 0;
01219
01220 this_data = get_data();
01221
01222 float da_rad = da*(float)M_PI/180.0f;
01223
01224 float *with_data = with->get_data();
01225 float mx0 = cos(da_rad);
01226 float mx1 = sin(da_rad);
01227 float y = -ny / 2.0f;
01228 float my0 = mx0 * (-nx / 2.0f - 1.0f) + nx / 2.0f - dx;
01229 float my1 = -mx1 * (-nx / 2.0f - 1.0f) + ny / 2.0f - dy;
01230 double result = 0;
01231
01232 for (int j = 0; j < ny; j++) {
01233 float x2 = my0 + mx1 * y;
01234 float y2 = my1 + mx0 * y;
01235
01236 int ii = Util::fast_floor(x2);
01237 int jj = Util::fast_floor(y2);
01238 float t = x2 - ii;
01239 float u = y2 - jj;
01240
01241 for (int i = 0; i < nx; i++) {
01242 t += mx0;
01243 u -= mx1;
01244
01245 if (t >= 1.0f) {
01246 ii++;
01247 t -= 1.0f;
01248 }
01249
01250 if (u >= 1.0f) {
01251 jj++;
01252 u -= 1.0f;
01253 }
01254
01255 if (t < 0) {
01256 ii--;
01257 t += 1.0f;
01258 }
01259
01260 if (u < 0) {
01261 jj--;
01262 u += 1.0f;
01263 }
01264
01265 if (ii >= 0 && ii <= nx - 2 && jj >= 0 && jj <= ny - 2) {
01266 int k0 = ii + jj * nx;
01267 int k1 = k0 + 1;
01268 int k2 = k0 + nx + 1;
01269 int k3 = k0 + nx;
01270
01271 float tt = 1 - t;
01272 float uu = 1 - u;
01273 int idx = i + j * nx;
01274 if (mirror) idx = nx-1-i+j*nx;
01275 result += (this_data[k0] * tt * uu + this_data[k1] * t * uu +
01276 this_data[k2] * t * u + this_data[k3] * tt * u) * with_data[idx];
01277 }
01278 }
01279 y += 1.0f;
01280 }
01281
01282 EXITFUNC;
01283 return result;
01284 }
01285
01286
01287 EMData *EMData::little_big_dot(EMData * with, bool do_sigma)
01288 {
01289 ENTERFUNC;
01290
01291 if (get_ndim() > 2) {
01292 throw ImageDimensionException("1D/2D only");
01293 }
01294
01295 EMData *ret = copy_head();
01296 ret->set_size(nx,ny,nz);
01297 ret->to_zero();
01298
01299 int nx2 = with->get_xsize();
01300 int ny2 = with->get_ysize();
01301 float em = with->get_edge_mean();
01302
01303 float *data = get_data();
01304 float *with_data = with->get_data();
01305 float *ret_data = ret->get_data();
01306
01307 float sum2 = (Util::square((float)with->get_attr("sigma")) +
01308 Util::square((float)with->get_attr("mean")));
01309 if (do_sigma) {
01310 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
01311 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
01312 float sum = 0;
01313 float sum1 = 0;
01314 float summ = 0;
01315 int k = 0;
01316
01317 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01318 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01319 int l = ii + jj * nx;
01320 sum1 += Util::square(data[l]);
01321 summ += data[l];
01322 sum += data[l] * with_data[k];
01323 k++;
01324 }
01325 }
01326 float tmp_f1 = (sum1 / 2.0f - sum) / (nx2 * ny2);
01327 float tmp_f2 = Util::square((float)with->get_attr("mean") -
01328 summ / (nx2 * ny2));
01329 ret_data[i + j * nx] = sum2 + tmp_f1 - tmp_f2;
01330 }
01331 }
01332 }
01333 else {
01334 for (int j = ny2 / 2; j < ny - ny2 / 2; j++) {
01335 for (int i = nx2 / 2; i < nx - nx2 / 2; i++) {
01336 float eml = 0;
01337 float dot = 0;
01338 float dot2 = 0;
01339
01340 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01341 eml += data[ii + (j - ny2 / 2) * nx] + data[ii + (j + ny2 / 2 - 1) * nx];
01342 }
01343
01344 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01345 eml += data[i - nx2 / 2 + jj * nx] + data[i + nx2 / 2 - 1 + jj * nx];
01346 }
01347
01348 eml /= (nx2 + ny2) * 2.0f;
01349 int k = 0;
01350
01351 for (int jj = j - ny2 / 2; jj < j + ny2 / 2; jj++) {
01352 for (int ii = i - nx2 / 2; ii < i + nx2 / 2; ii++) {
01353 dot += (data[ii + jj * nx] - eml) * (with_data[k] - em);
01354 dot2 += Util::square(data[ii + jj * nx] - eml);
01355 k++;
01356 }
01357 }
01358
01359 dot2 = std::sqrt(dot2);
01360
01361 if (dot2 == 0) {
01362 ret_data[i + j * nx] = 0;
01363 }
01364 else {
01365 ret_data[i + j * nx] = dot / (nx2 * ny2 * dot2 * (float)with->get_attr("sigma"));
01366 }
01367 }
01368 }
01369 }
01370
01371 ret->update();
01372
01373 EXITFUNC;
01374 return ret;
01375 }
01376
01377
01378 EMData *EMData::do_radon()
01379 {
01380 ENTERFUNC;
01381
01382 if (get_ndim() != 2) {
01383 throw ImageDimensionException("2D only");
01384 }
01385
01386 if (nx != ny) {
01387 throw ImageFormatException("square image only");
01388 }
01389
01390 EMData *result = new EMData();
01391 result->set_size(nx, ny, 1);
01392 result->to_zero();
01393 float *result_data = result->get_data();
01394
01395 EMData *this_copy = this;
01396 this_copy = copy();
01397
01398 for (int i = 0; i < nx; i++) {
01399 Transform t(Dict("type","2d","alpha",(float) M_PI * 2.0f * i / nx));
01400 this_copy->transform(t);
01401
01402 float *copy_data = this_copy->get_data();
01403
01404 for (int y = 0; y < nx; y++) {
01405 for (int x = 0; x < nx; x++) {
01406 if (Util::square(x - nx / 2) + Util::square(y - nx / 2) <= nx * nx / 4) {
01407 result_data[i + y * nx] += copy_data[x + y * nx];
01408 }
01409 }
01410 }
01411
01412 this_copy->update();
01413 }
01414
01415 result->update();
01416
01417 if( this_copy )
01418 {
01419 delete this_copy;
01420 this_copy = 0;
01421 }
01422
01423 EXITFUNC;
01424 return result;
01425 }
01426
01427 void EMData::zero_corner_circulant(const int radius)
01428 {
01429 if ( nz > 1 && nz < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nz is too small");
01430 if ( ny > 1 && ny < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - ny is too small");
01431 if ( nx > 1 && nx < (2*radius+1) ) throw ImageDimensionException("Error: cannot zero corner - nx is too small");
01432
01433 int it_z = radius;
01434 int it_y = radius;
01435 int it_x = radius;
01436
01437 if ( nz == 1 ) it_z = 0;
01438 if ( ny == 1 ) it_y = 0;
01439 if ( nx == 1 ) it_z = 0;
01440
01441 if ( nz == 1 && ny == 1 )
01442 {
01443 for ( int x = -it_x; x <= it_x; ++x )
01444 get_value_at_wrap(x) = 0;
01445
01446 }
01447 else if ( nz == 1 )
01448 {
01449 for ( int y = -it_y; y <= it_y; ++y)
01450 for ( int x = -it_x; x <= it_x; ++x )
01451 get_value_at_wrap(x,y) = 0;
01452 }
01453 else
01454 {
01455 for( int z = -it_z; z <= it_z; ++z )
01456 for ( int y = -it_y; y <= it_y; ++y)
01457 for ( int x = -it_x; x < it_x; ++x )
01458 get_value_at_wrap(x,y,z) = 0;
01459
01460 }
01461
01462 }
01463
01464 EMData *EMData::calc_ccf(EMData * with, fp_flag fpflag,bool center)
01465 {
01466 ENTERFUNC;
01467 if( with == 0 ) {
01468 EXITFUNC;
01469 return convolution(this,this,fpflag, center);
01470 }
01471 else if ( with == this ){
01472 EXITFUNC;
01473 return correlation(this, this, fpflag,center);
01474 }
01475 else {
01476
01477 #ifdef EMAN2_USING_CUDA //CUDA
01478
01479
01480
01481 if(cudarwdata && with->cudarwdata) {
01482
01483 EMData* afft = 0;
01484 EMData* bfft = 0;
01485 bool delafft = false, delbfft = false;
01486 int offset = 0;
01487
01488
01489 if(!is_complex()){
01490 afft = do_fft_cuda();
01491 delafft = true;
01492 offset = 2 - nx%2;
01493
01494 }else{
01495 afft = this;
01496 }
01497 if(!with->is_complex()){
01498 bfft = with->do_fft_cuda();
01499
01500 delbfft = true;
01501 }else{
01502 bfft = with;
01503 }
01504
01505 calc_ccf_cuda(afft->cudarwdata,bfft->cudarwdata,nx + offset, ny, nz);
01506
01507 if(delbfft) delete bfft;
01508
01509 EMData * corr = afft->do_ift_cuda();
01510 if(delafft) delete afft;
01511
01512
01513 return corr;
01514 }
01515 #endif
01516
01517
01518 bool undoresize = false;
01519 int wnx = with->get_xsize(); int wny = with->get_ysize(); int wnz = with->get_zsize();
01520 if (!(is_complex() ^ with->is_complex()) && (wnx != nx || wny != ny || wnz != nz)) {
01521 cout << "Warning!!! resizing image before CCF calculation" << endl;
01522 Region r((wnx-nx)/2, (wny-ny)/2, (wnz-nz)/2,nx,ny,nz);
01523 with->clip_inplace(r);
01524 undoresize = true;
01525 }
01526
01527 EMData* cor = correlation(this, with, fpflag, center);
01528
01529
01530 if ( undoresize ) {
01531 Region r((nx-wnx)/2, (ny-wny)/2,(nz-wnz)/2,wnx,wny,wnz);
01532 with->clip_inplace(r);
01533 }
01534
01535 EXITFUNC;
01536 return cor;
01537 }
01538 }
01539
01540 EMData *EMData::calc_ccfx( EMData * const with, int y0, int y1, bool no_sum, bool flip)
01541 {
01542 ENTERFUNC;
01543
01544 if (!with) {
01545 LOGERR("NULL 'with' image. ");
01546 throw NullPointerException("NULL input image");
01547 }
01548
01549 if (!EMUtil::is_same_size(this, with)) {
01550 LOGERR("images not same size: (%d,%d,%d) != (%d,%d,%d)",
01551 nx, ny, nz,
01552 with->get_xsize(), with->get_ysize(), with->get_zsize());
01553 throw ImageFormatException("images not same size");
01554 }
01555 if (get_ndim() > 2) {
01556 LOGERR("2D images only");
01557 throw ImageDimensionException("2D images only");
01558 }
01559
01560 if (y1 <= y0) {
01561 y1 = ny;
01562 }
01563
01564 if (y0 >= y1) {
01565 y0 = 0;
01566 }
01567
01568 if (y0 < 0) {
01569 y0 = 0;
01570 }
01571
01572 if (y1 > ny) {
01573 y1 = ny;
01574 }
01575 if (is_complex_x() || with->is_complex_x() ) throw;
01576
01577 static int nx_fft = 0;
01578 static int ny_fft = 0;
01579 static EMData f1;
01580 static EMData f2;
01581 static EMData rslt;
01582
01583 int height = y1-y0;
01584 int width = (nx+2-(nx%2));
01585 if (width != nx_fft || height != ny_fft ) {
01586 f1.set_size(width,height);
01587 f2.set_size(width,height);
01588 rslt.set_size(nx,height);
01589 nx_fft = width;
01590 ny_fft = height;
01591 }
01592
01593 #ifdef EMAN2_USING_CUDA
01594 if (cudarwdata && with->cudarwdata) {
01595
01596 if(!f1.cudarwdata) f1.rw_alloc();
01597 if(!f2.cudarwdata) f2.rw_alloc();
01598 if(!rslt.cudarwdata) rslt.rw_alloc();
01599 cuda_dd_fft_real_to_complex_nd(cudarwdata, f1.cudarwdata, nx, 1, 1, height);
01600 cuda_dd_fft_real_to_complex_nd(with->cudarwdata, f2.cudarwdata, nx, 1, 1, height);
01601 calc_ccf_cuda(f1.cudarwdata, f2.cudarwdata, nx, ny, nz);
01602 cuda_dd_fft_complex_to_real_nd(f1.cudarwdata, rslt.cudarwdata, nx, 1, 1, height);
01603 if(no_sum){
01604 EMData* result = new EMData(rslt);
01605 return result;
01606 }
01607 EMData* cf = new EMData(0,0,nx,1,1);
01608 cf->runcuda(emdata_column_sum(rslt.cudarwdata, nx, ny));
01609
01610 EXITFUNC;
01611 return cf;
01612 }
01613 #endif
01614
01615 float *d1 = get_data();
01616 float *d2 = with->get_data();
01617 float *f1d = f1.get_data();
01618 float *f2d = f2.get_data();
01619 for (int j = 0; j < height; j++) {
01620 EMfft::real_to_complex_1d(d1 + j * nx, f1d+j*width, nx);
01621 EMfft::real_to_complex_1d(d2 + j * nx, f2d+j*width, nx);
01622 }
01623
01624 if(flip == false) {
01625 for (int j = 0; j < height; j++) {
01626 float *f1a = f1d + j * width;
01627 float *f2a = f2d + j * width;
01628
01629 for (int i = 0; i < width / 2; i++) {
01630 float re1 = f1a[2*i];
01631 float re2 = f2a[2*i];
01632 float im1 = f1a[2*i+1];
01633 float im2 = f2a[2*i+1];
01634
01635 f1d[j*width+i*2] = re1 * re2 + im1 * im2;
01636 f1d[j*width+i*2+1] = im1 * re2 - re1 * im2;
01637 }
01638 }
01639 } else {
01640 for (int j = 0; j < height; j++) {
01641 float *f1a = f1d + j * width;
01642 float *f2a = f2d + j * width;
01643
01644 for (int i = 0; i < width / 2; i++) {
01645 float re1 = f1a[2*i];
01646 float re2 = f2a[2*i];
01647 float im1 = f1a[2*i+1];
01648 float im2 = f2a[2*i+1];
01649
01650 f1d[j*width+i*2] = re1 * re2 - im1 * im2;
01651 f1d[j*width+i*2+1] = im1 * re2 + re1 * im2;
01652 }
01653 }
01654 }
01655
01656 float* rd = rslt.get_data();
01657 for (int j = y0; j < y1; j++) {
01658 EMfft::complex_to_real_1d(f1d+j*width, rd+j*nx, nx);
01659 }
01660
01661 if (no_sum) {
01662 rslt.update();
01663 EXITFUNC;
01664 return new EMData(rslt);
01665 } else {
01666 EMData *cf = new EMData(nx,1,1);
01667 cf->to_zero();
01668 float *c = cf->get_data();
01669 for (int j = 0; j < height; j++) {
01670 for(int i = 0; i < nx; ++i) {
01671 c[i] += rd[i+j*nx];
01672 }
01673 }
01674 cf->update();
01675 EXITFUNC;
01676 return cf;
01677 }
01678 }
01679
01680 EMData *EMData::make_rotational_footprint_cmc( bool unwrap) {
01681 ENTERFUNC;
01682 update_stat();
01683
01684
01685
01686
01687
01688
01689 if ( rot_fp != 0 && unwrap == true) {
01690 return new EMData(*rot_fp);
01691 }
01692
01693 static EMData obj_filt;
01694 EMData* filt = &obj_filt;
01695 filt->set_complex(true);
01696
01697
01698
01699
01700
01701
01702 if (filt->get_xsize() != nx+2-(nx%2) || filt->get_ysize() != ny ||
01703 filt->get_zsize() != nz ) {
01704 filt->set_size(nx+2-(nx%2), ny, nz);
01705 filt->to_one();
01706
01707 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01708 }
01709
01710 EMData *ccf = this->calc_mutual_correlation(this, true,filt);
01711 ccf->sub(ccf->get_edge_mean());
01712 EMData *result = ccf->unwrap();
01713 delete ccf; ccf = 0;
01714
01715 EXITFUNC;
01716 if ( unwrap == true)
01717 {
01718
01719
01720
01721
01722
01723
01724
01725
01726 rot_fp = result;
01727 return new EMData(*rot_fp);
01728 }
01729 else return result;
01730 }
01731
01732 EMData *EMData::make_rotational_footprint( bool unwrap) {
01733 ENTERFUNC;
01734 update_stat();
01735
01736
01737
01738
01739
01740
01741 if ( rot_fp != 0 && unwrap == true) {
01742 return new EMData(*rot_fp);
01743 }
01744
01745 EMData* ccf = this->calc_ccf(this,CIRCULANT,true);
01746 ccf->sub(ccf->get_edge_mean());
01747
01748 EMData *result = ccf->unwrap();
01749 delete ccf; ccf = 0;
01750
01751 EXITFUNC;
01752 if ( unwrap == true)
01753 {
01754
01755
01756
01757
01758
01759
01760
01761 rot_fp = result;
01762 return new EMData(*rot_fp);
01763 }
01764 else return result;
01765 }
01766
01767 EMData *EMData::make_rotational_footprint_e1( bool unwrap)
01768 {
01769 ENTERFUNC;
01770
01771 update_stat();
01772
01773
01774
01775
01776
01777
01778 if ( rot_fp != 0 && unwrap == true) {
01779 return new EMData(*rot_fp);
01780 }
01781
01782 static EMData obj_filt;
01783 EMData* filt = &obj_filt;
01784 filt->set_complex(true);
01785
01786
01787
01788
01789
01790
01791 int cs = (((nx * 7 / 4) & 0xfffff8) - nx) / 2;
01792
01793 static EMData big_clip;
01794 int big_x = nx+2*cs;
01795 int big_y = ny+2*cs;
01796 int big_z = 1;
01797 if ( nz != 1 ) {
01798 big_z = nz+2*cs;
01799 }
01800
01801
01802 if ( big_clip.get_xsize() != big_x || big_clip.get_ysize() != big_y || big_clip.get_zsize() != big_z ) {
01803 big_clip.set_size(big_x,big_y,big_z);
01804 }
01805
01806
01807
01808
01809 big_clip.to_value(get_edge_mean());
01810
01811 if (nz != 1) {
01812 big_clip.insert_clip(this,IntPoint(cs,cs,cs));
01813 } else {
01814 big_clip.insert_clip(this,IntPoint(cs,cs,0));
01815 }
01816
01817
01818
01819
01820
01821 if (filt->get_xsize() != big_clip.get_xsize() +2-(big_clip.get_xsize()%2) || filt->get_ysize() != big_clip.get_ysize() ||
01822 filt->get_zsize() != big_clip.get_zsize()) {
01823 filt->set_size(big_clip.get_xsize() + 2-(big_clip.get_xsize()%2), big_clip.get_ysize(), big_clip.get_zsize());
01824 filt->to_one();
01825 filt->process_inplace("filter.highpass.gauss", Dict("cutoff_abs", 1.5f/nx));
01826 #ifdef EMAN2_USING_CUDA
01827 if(big_clip.cudarwdata)
01828 {
01829 filt->copy_to_cuda();
01830 }
01831 #endif
01832 }
01833 #ifdef EMAN2_USING_CUDA
01834 if(big_clip.cudarwdata && !filt->cudarwdata)
01835 {
01836 filt->copy_to_cuda();
01837 }
01838 #endif
01839
01840 EMData *mc = big_clip.calc_mutual_correlation(&big_clip, true,filt);
01841 mc->sub(mc->get_edge_mean());
01842
01843 static EMData sml_clip;
01844 int sml_x = nx * 3 / 2;
01845 int sml_y = ny * 3 / 2;
01846 int sml_z = 1;
01847 if ( nz != 1 ) {
01848 sml_z = nz * 3 / 2;
01849 }
01850
01851 if ( sml_clip.get_xsize() != sml_x || sml_clip.get_ysize() != sml_y || sml_clip.get_zsize() != sml_z ) {
01852 sml_clip.set_size(sml_x,sml_y,sml_z); }
01853 if (nz != 1) {
01854 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,-cs+nz/4));
01855 } else {
01856 sml_clip.insert_clip(mc,IntPoint(-cs+nx/4,-cs+ny/4,0));
01857 }
01858
01859 delete mc; mc = 0;
01860 EMData * result = NULL;
01861
01862 if (nz == 1) {
01863 if (!unwrap) {
01864 #ifdef EMAN2_USING_CUDA
01865 if(sml_clip.cudarwdata) throw UnexpectedBehaviorException("shap masking is not yet supported by CUDA");
01866 #endif
01867 result = sml_clip.process("mask.sharp", Dict("outer_radius", -1, "value", 0));
01868
01869 }
01870 else {
01871 result = sml_clip.unwrap();
01872 }
01873 }
01874 else {
01875
01876
01877
01878 result = new EMData(sml_clip);
01879 }
01880
01881 #ifdef EMAN2_USING_CUDA
01882 sml_clip.roneedsanupdate();
01883 #endif
01884 EXITFUNC;
01885 if ( unwrap == true)
01886 {
01887
01888
01889
01890 if ( rot_fp != 0 ) throw UnexpectedBehaviorException("The rotational foot print is only expected to be cached if it is not NULL");
01891
01892
01893
01894 rot_fp = result;
01895 return new EMData(*rot_fp);
01896 }
01897 else return result;
01898 }
01899
01900 EMData *EMData::make_footprint(int type)
01901 {
01902
01903 if (type==0) {
01904 EMData *un=make_rotational_footprint_e1();
01905 if (un->get_ysize() <= 6) {
01906 throw UnexpectedBehaviorException("In EMData::make_footprint. The rotational footprint is too small");
01907 }
01908 EMData *tmp=un->get_clip(Region(0,4,un->get_xsize(),un->get_ysize()-6));
01909 EMData *cx=tmp->calc_ccfx(tmp,0,-1,1);
01910 EMData *fp=cx->get_clip(Region(0,0,cx->get_xsize()/2,cx->get_ysize()));
01911 delete un;
01912 delete tmp;
01913 delete cx;
01914 return fp;
01915 }
01916 else if (type==1 || type==2 ||type==5 || type==6) {
01917 int i,j,kx,ky,lx,ly;
01918
01919 EMData *fft=do_fft();
01920
01921
01922 int rmax=(get_xsize()+1)/2;
01923 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
01924 for (i=0; i<rmax; i++) {
01925 for (j=0; j<rmax; j++) {
01926 #ifdef _WIN32
01927 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
01928 #else
01929 rmap[i+j*rmax]=hypot((float)i,(float)j);
01930 #endif //_WIN32
01931
01932 }
01933 }
01934
01935 EMData *fp=new EMData(rmax*2+2,rmax*2,1);
01936 fp->set_complex(1);
01937 fp->to_zero();
01938
01939
01940
01941
01942 for (kx=-rmax+1; kx<rmax; kx++) {
01943 for (ky=-rmax+1; ky<rmax; ky++) {
01944 for (lx=-rmax+1; lx<rmax; lx++) {
01945 for (ly=-rmax+1; ly<rmax; ly++) {
01946 int ax=kx+lx;
01947 int ay=ky+ly;
01948 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
01949 int r1=(int)floor(.5+rmap[abs(kx)+rmax*abs(ky)]);
01950 int r2=(int)floor(.5+rmap[abs(lx)+rmax*abs(ly)]);
01951
01952
01953 if (r1+r2>=rmax) continue;
01954
01955 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
01956 fp->set_value_at(r1*2,r2,p.real()+fp->get_value_at(r1*2,r2));
01957
01958
01959 fp->set_value_at(r1*2+1,r2,fp->get_value_at(r1*2+1,r2)+1);
01960 }
01961 }
01962 }
01963 }
01964
01965
01966 if (type==5 || type==6) {
01967 for (i=0; i<rmax*2; i+=2) {
01968 for (j=0; j<rmax; j++) {
01969 float norm=fp->get_value_at(i+1,j);
01970 #ifdef _WIN32
01971 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));
01972 fp->set_value_at(i,j,pow(fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm), 1.0f/3.0f));
01973 #else
01974 fp->set_value_at(i,rmax*2-j-1,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
01975 fp->set_value_at(i,j,cbrt(fp->get_value_at(i,j)/(norm==0?1.0:norm)));
01976 #endif //_WIN32
01977 fp->set_value_at(i+1,j,0.0);
01978 }
01979 }
01980 }
01981 else {
01982 for (i=0; i<rmax*2; i+=2) {
01983 for (j=0; j<rmax; j++) {
01984 float norm=fp->get_value_at(i+1,j);
01985 fp->set_value_at(i,rmax*2-j-1,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
01986 fp->set_value_at(i,j,fp->get_value_at(i,j)/(norm==0.0f?1.0f:norm));
01987 fp->set_value_at(i+1,j,0.0);
01988 }
01989 }
01990 }
01991
01992 free(rmap);
01993 if (type==2||type==6) {
01994 EMData *f2=fp->do_ift();
01995 if (f2->get_value_at(0,0)<0) f2->mult(-1.0f);
01996 f2->process_inplace("xform.phaseorigin.tocorner");
01997 delete fp;
01998 return f2;
01999 }
02000 return fp;
02001 }
02002 else if (type==3 || type==4) {
02003 int h,i,j,kx,ky,lx,ly;
02004
02005 EMData *fft=do_fft();
02006
02007
02008 int rmax=(get_xsize()+1)/2;
02009 float *rmap=(float *)malloc(rmax*rmax*sizeof(float));
02010 for (i=0; i<rmax; i++) {
02011 for (j=0; j<rmax; j++) {
02012 #ifdef _WIN32
02013 rmap[i+j*rmax]=_hypotf((float)i,(float)j);
02014 #else
02015 rmap[i+j*rmax]=hypot((float)i,(float)j);
02016 #endif //_WIN32
02017
02018 }
02019 }
02020
02021 EMData *fp=new EMData(rmax*2+2,rmax*2,16);
02022
02023 fp->set_complex(1);
02024 fp->to_zero();
02025
02026
02027
02028
02029 for (kx=-rmax+1; kx<rmax; kx++) {
02030 for (ky=-rmax+1; ky<rmax; ky++) {
02031 for (lx=-rmax+1; lx<rmax; lx++) {
02032 for (ly=-rmax+1; ly<rmax; ly++) {
02033 int ax=kx+lx;
02034 int ay=ky+ly;
02035 if (abs(ax)>=rmax || abs(ay)>=rmax) continue;
02036 float rr1=rmap[abs(kx)+rmax*abs(ky)];
02037 float rr2=rmap[abs(lx)+rmax*abs(ly)];
02038 int r1=(int)floor(.5+rr1);
02039 int r2=(int)floor(.5+rr2);
02040
02041
02042 if (r1+r2>=rmax || rr1==0 ||rr2==0) continue;
02043
02044 std::complex<float> p=fft->get_complex_at(kx,ky)*fft->get_complex_at(lx,ly)*conj(fft->get_complex_at(ax,ay));
02045 int dot=(int)floor((kx*lx+ky*ly)/(rr1*rr2)*7.5);
02046 if (dot<0) dot=16+dot;
02047
02048 fp->set_value_at(r1*2,r2,dot,p.real()+fp->get_value_at(r1*2,r2,dot));
02049
02050
02051 fp->set_value_at(r1*2+1,r2,dot,fp->get_value_at(r1*2+1,r2,dot)+1);
02052 }
02053 }
02054 }
02055 }
02056
02057
02058 for (i=0; i<rmax*2; i+=2) {
02059 for (j=0; j<rmax; j++) {
02060 for (h=0; h<16; h++) {
02061 float norm=fp->get_value_at(i+1,j,h);
02062
02063
02064 fp->set_value_at(i,rmax*2-j-1,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02065 fp->set_value_at(i,j,h,(fp->get_value_at(i,j,h)/(norm==0.0f?1.0f:norm)));
02066
02067
02068 fp->set_value_at(i+1,j,h,0.0);
02069 }
02070 }
02071 }
02072
02073 free(rmap);
02074 if (type==4) {
02075 EMData *f2=fp->do_ift();
02076 if (f2->get_value_at(0,0,0)<0) f2->mult(-1.0f);
02077 f2->process_inplace("xform.phaseorigin.tocorner");
02078 delete fp;
02079 return f2;
02080 }
02081 return fp;
02082 }
02083 throw UnexpectedBehaviorException("There is not implementation for the parameters you specified");
02084 }
02085
02086
02087 EMData *EMData::calc_mutual_correlation(EMData * with, bool tocenter, EMData * filter)
02088 {
02089 ENTERFUNC;
02090
02091 if (with && !EMUtil::is_same_size(this, with)) {
02092 LOGERR("images not same size");
02093 throw ImageFormatException( "images not same size");
02094 }
02095
02096 #ifdef EMAN2_USING_CUDA
02097 if(cudarwdata && with->cudarwdata)
02098 {
02099
02100 EMData* this_fft = do_fft_cuda();
02101
02102 EMData *cf = 0;
02103 if (with && with != this) {
02104 cf = with->do_fft_cuda();
02105 }else{
02106 cf = this_fft->copy();
02107 }
02108
02109 if (filter) {
02110 if (!EMUtil::is_same_size(filter, cf)) {
02111 LOGERR("improperly sized filter");
02112 throw ImageFormatException("improperly sized filter");
02113 }
02114 mult_complex_efficient_cuda(cf->cudarwdata, filter->cudarwdata, cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), 1);
02115 mult_complex_efficient_cuda(this_fft->cudarwdata, filter->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize(), 1);
02116 }
02117
02118 mcf_cuda(this_fft->cudarwdata, cf->cudarwdata, this_fft->get_xsize(), this_fft->get_ysize(), this_fft->get_zsize());
02119
02120 EMData *f2 = cf->do_ift_cuda();
02121
02122 if (tocenter) {
02123 f2->process_inplace("xform.phaseorigin.tocenter");
02124 }
02125
02126 if( cf )
02127 {
02128 delete cf;
02129 cf = 0;
02130 }
02131
02132 if( this_fft )
02133 {
02134 delete this_fft;
02135 this_fft = 0;
02136 }
02137
02138 f2->set_attr("label", "MCF");
02139 f2->set_path("/tmp/eman.mcf");
02140
02141 EXITFUNC;
02142 return f2;
02143 }
02144 #endif
02145
02146 EMData *this_fft = 0;
02147 this_fft = do_fft();
02148
02149 if (!this_fft) {
02150
02151 LOGERR("FFT returns NULL image");
02152 throw NullPointerException("FFT returns NULL image");
02153 }
02154
02155 this_fft->ap2ri();
02156 EMData *cf = 0;
02157
02158 if (with && with != this) {
02159 cf = with->do_fft();
02160 if (!cf) {
02161 LOGERR("FFT returns NULL image");
02162 throw NullPointerException("FFT returns NULL image");
02163 }
02164 cf->ap2ri();
02165 }
02166 else {
02167 cf = this_fft->copy();
02168 }
02169
02170 if (filter) {
02171 if (!EMUtil::is_same_size(filter, cf)) {
02172 LOGERR("improperly sized filter");
02173 throw ImageFormatException("improperly sized filter");
02174 }
02175
02176 cf->mult_complex_efficient(*filter,true);
02177 this_fft->mult(*filter,true);
02178
02179
02180
02181
02182 }
02183
02184 float *rdata1 = this_fft->get_data();
02185 float *rdata2 = cf->get_data();
02186 size_t this_fft_size = (size_t)this_fft->get_xsize() * this_fft->get_ysize() * this_fft->get_zsize();
02187
02188 if (with == this) {
02189 for (size_t i = 0; i < this_fft_size; i += 2) {
02190 rdata2[i] = std::sqrt(rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02191 rdata2[i + 1] = 0;
02192 }
02193
02194 this_fft->update();
02195 cf->update();
02196 }
02197 else {
02198 for (size_t i = 0; i < this_fft_size; i += 2) {
02199 rdata2[i] = (rdata1[i] * rdata2[i] + rdata1[i + 1] * rdata2[i + 1]);
02200 rdata2[i + 1] = (rdata1[i + 1] * rdata2[i] - rdata1[i] * rdata2[i + 1]);
02201 }
02202
02203
02204 for (size_t i = 0; i < this_fft_size; i += 2) {
02205 float t = Util::square(rdata2[i]) + Util::square(rdata2[i + 1]);
02206 if (t != 0) {
02207 t = pow(t, 0.25f);
02208 rdata2[i] /= t;
02209 rdata2[i + 1] /= t;
02210 }
02211 }
02212 this_fft->update();
02213 cf->update();
02214 }
02215
02216 EMData *f2 = cf->do_ift();
02217
02218 if (tocenter) {
02219 f2->process_inplace("xform.phaseorigin.tocenter");
02220 }
02221
02222 if( cf )
02223 {
02224 delete cf;
02225 cf = 0;
02226 }
02227
02228 if( this_fft )
02229 {
02230 delete this_fft;
02231 this_fft = 0;
02232 }
02233
02234 f2->set_attr("label", "MCF");
02235 f2->set_path("/tmp/eman.mcf");
02236
02237 EXITFUNC;
02238 return f2;
02239 }
02240
02241
02242 vector < float > EMData::calc_hist(int hist_size, float histmin, float histmax,const float& brt, const float& cont)
02243 {
02244 ENTERFUNC;
02245
02246 static size_t prime[] = { 1, 3, 7, 11, 17, 23, 37, 59, 127, 253, 511 };
02247
02248 if (histmin == histmax) {
02249 histmin = get_attr("minimum");
02250 histmax = get_attr("maximum");
02251 }
02252
02253 vector <float> hist(hist_size, 0.0);
02254
02255 int p0 = 0;
02256 int p1 = 0;
02257 size_t size = (size_t)nx * ny * nz;
02258 if (size < 300000) {
02259 p0 = 0;
02260 p1 = 0;
02261 }
02262 else if (size < 2000000) {
02263 p0 = 2;
02264 p1 = 3;
02265 }
02266 else if (size < 8000000) {
02267 p0 = 4;
02268 p1 = 6;
02269 }
02270 else {
02271 p0 = 7;
02272 p1 = 9;
02273 }
02274
02275 if (is_complex() && p0 > 0) {
02276 p0++;
02277 p1++;
02278 }
02279
02280 size_t di = 0;
02281
02282 size_t n = hist.size();
02283
02284 float * data = get_data();
02285 for (int k = p0; k <= p1; ++k) {
02286 if (is_complex()) {
02287 di = prime[k] * 2;
02288 }
02289 else {
02290 di = prime[k];
02291 }
02292
02293
02294 float w = (float)n / (histmax - histmin);
02295
02296 for(size_t i=0; i<=size-di; i += di) {
02297 float val;
02298 if (cont != 1.0f || brt != 0)val = cont*(data[i]+brt);
02299 else val = data[i];
02300 int j = Util::round((val - histmin) * w);
02301 if (j >= 0 && j < (int) n) {
02302 hist[j] += 1;
02303 }
02304 }
02305 }
02306
02307
02308
02309
02310
02311
02312
02313 return hist;
02314
02315 EXITFUNC;
02316 }
02317
02318
02319
02320
02321
02322 vector<float> EMData::calc_az_dist(int n, float a0, float da, float rmin, float rmax)
02323 {
02324 ENTERFUNC;
02325
02326 a0=a0*M_PI/180.0f;
02327 da=da*M_PI/180.0f;
02328
02329 if (get_ndim() > 2) {
02330 throw ImageDimensionException("no 3D image");
02331 }
02332
02333 float *yc = new float[n];
02334
02335 vector<float> vd(n);
02336 for (int i = 0; i < n; i++) {
02337 yc[i] = 0.00001f;
02338 }
02339
02340 float * data = get_data();
02341 if (is_complex()) {
02342 int c = 0;
02343 for (int y = 0; y < ny; y++) {
02344 for (int x = 0; x < nx; x += 2, c += 2) {
02345 int x1 = x / 2;
02346 int y1 = y<ny/2?y:y-ny;
02347 float r = (float)Util::hypot_fast(x1,y1);
02348
02349 if (r >= rmin && r <= rmax) {
02350 float a = 0;
02351
02352 if (y != ny / 2 || x != 0) {
02353 a = (atan2((float)y1, (float)x1) - a0) / da;
02354 }
02355
02356 int i = (int)(floor(a));
02357 a -= i;
02358
02359 if (i == 0) {
02360 vd[0] += data[c] * (1.0f - a);
02361 yc[0] += (1.0f - a);
02362 }
02363 else if (i == n - 1) {
02364 vd[n - 1] += data[c] * a;
02365 yc[n - 1] += a;
02366 }
02367 else if (i > 0 && i < (n - 1)) {
02368 float h = 0;
02369 if (is_ri()) {
02370 #ifdef _WIN32
02371 h = (float)_hypot(data[c], data[c + 1]);
02372 #else
02373 h = (float)hypot(data[c], data[c + 1]);
02374 #endif //_WIN32
02375 }
02376 else {
02377 h = data[c];
02378 }
02379
02380 vd[i] += h * (1.0f - a);
02381 yc[i] += (1.0f - a);
02382 vd[i + 1] += h * a;
02383 yc[i + 1] += a;
02384 }
02385 }
02386 }
02387 }
02388 }
02389 else {
02390 int c = 0;
02391 float half_nx = (nx - 1) / 2.0f;
02392 float half_ny = (ny - 1) / 2.0f;
02393
02394 for (int y = 0; y < ny; y++) {
02395 for (int x = 0; x < nx; x++, c++) {
02396 float y1 = y - half_ny;
02397 float x1 = x - half_nx;
02398 #ifdef _WIN32
02399 float r = (float)_hypot(x1, y1);
02400 #else
02401 float r = (float)hypot(x1, y1);
02402 #endif
02403
02404 if (r >= rmin && r <= rmax) {
02405 float a = 0;
02406 if (x1 != 0 || y1 != 0) {
02407 a = atan2(y1, x1);
02408 if (a < 0) {
02409 a += M_PI * 2;
02410 }
02411 }
02412
02413 a = (a - a0) / da;
02414 int i = static_cast < int >(floor(a));
02415 a -= i;
02416
02417 if (i == 0) {
02418 vd[0] += data[c] * (1.0f - a);
02419 yc[0] += (1.0f - a);
02420 }
02421 else if (i == n - 1) {
02422 vd[n - 1] += data[c] * a;
02423 yc[n - 1] += (a);
02424 }
02425 else if (i > 0 && i < (n - 1)) {
02426 vd[i] += data[c] * (1.0f - a);
02427 yc[i] += (1.0f - a);
02428 vd[i + 1] += data[c] * a;
02429 yc[i + 1] += a;
02430 }
02431 }
02432 }
02433 }
02434 }
02435
02436
02437 for (int i = 0; i < n; i++) {
02438 vd[i] /= yc[i];
02439 }
02440
02441 if( yc )
02442 {
02443 delete[]yc;
02444 yc = 0;
02445 }
02446
02447 return vd;
02448
02449 EXITFUNC;
02450 }
02451
02452
02453 EMData *EMData::unwrap(int r1, int r2, int xs, int dx, int dy, bool do360, bool weight_radial) const
02454 {
02455 ENTERFUNC;
02456
02457 if (get_ndim() != 2) {
02458 throw ImageDimensionException("2D image only");
02459 }
02460
02461 int p = 1;
02462 if (do360) {
02463 p = 2;
02464 }
02465
02466 if (xs < 1) {
02467 xs = (int) Util::fast_floor(p * M_PI * ny / 4);
02468 xs -= xs % 8;
02469 if (xs<=8) xs=16;
02470 }
02471
02472 if (r1 < 0) {
02473 r1 = 4;
02474 }
02475
02476 #ifdef _WIN32
02477 int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(_hypot(dx, dy)));
02478 #else
02479 int rr = ny / 2 - 2 - (int) Util::fast_floor(static_cast<float>(hypot(dx, dy)));
02480 #endif //_WIN32
02481 rr-=rr%2;
02482 if (r2 <= r1 || r2 > rr) {
02483 r2 = rr;
02484 }
02485
02486 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");
02487
02488 #ifdef EMAN2_USING_CUDA
02489 if (isrodataongpu()){
02490
02491 EMData* result = new EMData(0,0,xs,(r2-r1),1);
02492 result->rw_alloc();
02493 bindcudaarrayA(true);
02494 emdata_unwrap(result->cudarwdata, r1, r2, xs, p, dx, dy, weight_radial, nx, ny);
02495 unbindcudaarryA();
02496 return result;
02497 }
02498 #endif
02499
02500 EMData *ret = new EMData();
02501 ret->set_size(xs, r2 - r1, 1);
02502 const float *const d = get_const_data();
02503 float *dd = ret->get_data();
02504 float pfac = (float)p/(float)xs;
02505
02506 int nxon2 = nx/2;
02507 int nyon2 = ny/2;
02508 for (int x = 0; x < xs; x++) {
02509 float ang = x * M_PI * pfac;
02510 float si = sin(ang);
02511 float co = cos(ang);
02512
02513 for (int y = 0; y < r2 - r1; y++) {
02514 float ypr1 = (float)y + r1;
02515 float xx = ypr1 * co + nxon2 + dx;
02516 float yy = ypr1 * si + nyon2 + dy;
02517
02518
02519 float t = xx - (int)xx;
02520 float u = yy - (int)yy;
02521
02522 int k = (int) xx + ((int) yy) * nx;
02523 float val = Util::bilinear_interpolate(d[k], d[k + 1], d[k + nx], d[k + nx+1], t,u);
02524 if (weight_radial) val *= ypr1;
02525 dd[x + y * xs] = val;
02526 }
02527
02528 }
02529 ret->update();
02530
02531 EXITFUNC;
02532 return ret;
02533 }
02534
02535
02536
02537 void EMData::apply_radial_func(float x0, float step, vector < float >array, bool interp)
02538 {
02539 ENTERFUNC;
02540
02541 if (!is_complex()) throw ImageFormatException("apply_radial_func requires a complex image");
02542
02543 int n = static_cast < int >(array.size());
02544
02545 if (n*step>2.0) printf("Warning, apply_radial_func takes x0 and step with respect to Nyquist (0.5)\n");
02546
02547
02548
02549 ap2ri();
02550
02551 size_t ndims = get_ndim();
02552 float * data = get_data();
02553 if (ndims == 2) {
02554 int k = 0;
02555 for (int j = 0; j < ny; j++) {
02556 for (int i = 0; i < nx; i += 2, k += 2) {
02557 float r;
02558 #ifdef _WIN32
02559 if (j<ny/2) r = (float)_hypot(i/(float)(nx*2), j/(float)ny);
02560 else r = (float)_hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02561 #else
02562 if (j<ny/2) r = (float)hypot(i/(float)(nx*2), j/(float)ny);
02563 else r = (float)hypot(i/(float)(nx*2), (ny-j)/(float)ny);
02564 #endif //_WIN32
02565 r = (r - x0) / step;
02566
02567 int l = 0;
02568 if (interp) {
02569 l = (int) floor(r);
02570 }
02571 else {
02572 l = (int) floor(r + 1);
02573 }
02574
02575
02576 float f = 0;
02577 if (l >= n - 2) {
02578 f = array[n - 1];
02579 }
02580 else {
02581 if (interp) {
02582 r -= l;
02583 f = (array[l] * (1.0f - r) + array[l + 1] * r);
02584 }
02585 else {
02586 f = array[l];
02587 }
02588 }
02589
02590 data[k] *= f;
02591 data[k + 1] *= f;
02592 }
02593 }
02594 }
02595 else if (ndims == 3) {
02596 int k = 0;
02597 for (int m = 0; m < nz; m++) {
02598 float mnz;
02599 if (m<nz/2) mnz=m*m/(float)(nz*nz);
02600 else { mnz=(nz-m)/(float)nz; mnz*=mnz; }
02601
02602 for (int j = 0; j < ny; j++) {
02603 float jny;
02604 if (j<ny/2) jny= j*j/(float)(ny*ny);
02605 else { jny=(ny-j)/(float)ny; jny*=jny; }
02606
02607 for (int i = 0; i < nx; i += 2, k += 2) {
02608 float r = std::sqrt((i * i / (nx*nx*4.0f)) + jny + mnz);
02609 r = (r - x0) / step;
02610
02611 int l = 0;
02612 if (interp) {
02613 l = (int) floor(r);
02614 }
02615 else {
02616 l = (int) floor(r + 1);
02617 }
02618
02619
02620 float f = 0;
02621 if (l >= n - 2) {
02622 f = array[n - 1];
02623 }
02624 else {
02625 if (interp) {
02626 r -= l;
02627 f = (array[l] * (1.0f - r) + array[l + 1] * r);
02628 }
02629 else {
02630 f = array[l];
02631 }
02632 }
02633
02634 data[k] *= f;
02635 data[k + 1] *= f;
02636 }
02637 }
02638 }
02639
02640 }
02641
02642 update();
02643 EXITFUNC;
02644 }
02645
02646 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, bool inten)
02647 {
02648 ENTERFUNC;
02649
02650 vector<float>ret(n);
02651 vector<float>norm(n);
02652
02653 int x,y,z,i;
02654 int step=is_complex()?2:1;
02655 int isinten=get_attr_default("is_intensity",0);
02656
02657 if (isinten&&!inten) { throw InvalidParameterException("Must set inten for calc_radial_dist with intensity image"); }
02658
02659 for (i=0; i<n; i++) ret[i]=norm[i]=0.0;
02660 float * data = get_data();
02661
02662
02663 if (nz==1) {
02664 for (y=i=0; y<ny; y++) {
02665 for (x=0; x<nx; x+=step,i+=step) {
02666 float r,v;
02667 if (step==2) {
02668 if (x==0 && y>ny/2) continue;
02669 r=(float)(Util::hypot_fast(x/2,y<ny/2?y:ny-y));
02670 if (!inten) {
02671 #ifdef _WIN32
02672 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02673 #else
02674 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02675 #endif
02676 else v=data[i];
02677 } else {
02678 if (isinten) v=data[i];
02679 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02680 else v=data[i]*data[i];
02681 }
02682 }
02683 else {
02684 r=(float)(Util::hypot_fast(x-nx/2,y-ny/2));
02685 if (inten) v=data[i]*data[i];
02686 else v=data[i];
02687 }
02688 r=(r-x0)/dx;
02689 int f=int(r);
02690 r-=float(f);
02691
02692 if (f>=0 && f<n) {
02693 ret[f]+=v*(1.0f-r);
02694 norm[f]+=(1.0f-r);
02695 if (f<n-1) {
02696 ret[f+1]+=v*r;
02697 norm[f+1]+=r;
02698 }
02699 }
02700 }
02701 }
02702 }
02703 else {
02704 size_t i;
02705 for (z=i=0; z<nz; ++z) {
02706 for (y=0; y<ny; ++y) {
02707 for (x=0; x<nx; x+=step,i+=step) {
02708 float r,v;
02709 if (step==2) {
02710 if (x==0 && z<nz/2) continue;
02711 if (x==0 && z==nz/2 && y<ny/2) continue;
02712 r=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z);
02713 if (!inten) {
02714 #ifdef _WIN32
02715 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02716 #else
02717 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02718 #endif //_WIN32
02719 else v=data[i];
02720 } else {
02721 if (isinten) v=data[i];
02722 else if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02723 else v=data[i]*data[i];
02724 }
02725 }
02726 else {
02727 r=Util::hypot3(x-nx/2,y-ny/2,z-nz/2);
02728 if (inten) v=data[i]*data[i];
02729 else v=data[i];
02730 }
02731 r=(r-x0)/dx;
02732 int f=int(r);
02733 r-=float(f);
02734 if (f>=0 && f<n) {
02735 ret[f]+=v*(1.0f-r);
02736 norm[f]+=(1.0f-r);
02737 if (f<n-1) {
02738 ret[f+1]+=v*r;
02739 norm[f+1]+=r;
02740 }
02741 }
02742 }
02743 }
02744 }
02745 }
02746
02747 for (i=0; i<n; i++) ret[i]/=norm[i]?norm[i]:1.0f;
02748
02749 EXITFUNC;
02750
02751 return ret;
02752 }
02753
02754 vector<float> EMData::calc_radial_dist(int n, float x0, float dx, int nwedge, bool inten)
02755 {
02756 ENTERFUNC;
02757
02758 if (nz > 1) {
02759 LOGERR("2D images only.");
02760 throw ImageDimensionException("2D images only");
02761 }
02762
02763 vector<float>ret(n*nwedge);
02764 vector<float>norm(n*nwedge);
02765
02766 int x,y,i;
02767 int step=is_complex()?2:1;
02768 float astep=static_cast<float>(M_PI*2.0/nwedge);
02769 float* data = get_data();
02770 for (i=0; i<n*nwedge; i++) ret[i]=norm[i]=0.0;
02771
02772
02773 for (y=i=0; y<ny; y++) {
02774 for (x=0; x<nx; x+=step,i+=step) {
02775 float r,v,a;
02776 if (is_complex()) {
02777 #ifdef _WIN32
02778 r=static_cast<float>(_hypot(x/2.0,y<ny/2?y:ny-y));
02779 #else
02780 r=static_cast<float>(hypot(x/2.0,y<ny/2?y:ny-y));
02781 #endif
02782 a=atan2(float(y<ny/2?y:ny-y),x/2.0f);
02783 if (!inten) {
02784 #ifdef _WIN32
02785 if (is_ri()) v=static_cast<float>(_hypot(data[i],data[i+1]));
02786 #else
02787 if (is_ri()) v=static_cast<float>(hypot(data[i],data[i+1]));
02788 #endif //_WIN32
02789 else v=data[i];
02790 } else {
02791 if (is_ri()) v=data[i]*data[i]+data[i+1]*data[i+1];
02792 else v=data[i]*data[i];
02793 }
02794 }
02795 else {
02796 #ifdef _WIN32
02797 r=static_cast<float>(_hypot(x-nx/2,y-ny/2));
02798 #else
02799 r=static_cast<float>(hypot(x-nx/2,y-ny/2));
02800 #endif //_WIN32
02801 a=atan2(float(y-ny/2),float(x-nx/2));
02802 if (inten) v=data[i]*data[i];
02803 else v=data[i];
02804 }
02805 int bin=n*int((a+M_PI)/astep);
02806 if (bin>=nwedge) bin=nwedge-1;
02807 r=(r-x0)/dx;
02808 int f=int(r);
02809 r-=float(f);
02810 if (f>=0 && f<n) {
02811 ret[f+bin]+=v*(1.0f-r);
02812 norm[f+bin]+=(1.0f-r);
02813 if (f<n-1) {
02814 ret[f+1+bin]+=v*r;
02815 norm[f+1+bin]+=r;
02816 }
02817 }
02818 }
02819 }
02820
02821 for (i=0; i<n*nwedge; i++) ret[i]/=norm[i]?norm[i]:1.0f;
02822 EXITFUNC;
02823
02824 return ret;
02825 }
02826
02827 void EMData::cconj() {
02828 ENTERFUNC;
02829 if (!is_complex() || !is_ri())
02830 throw ImageFormatException("EMData::conj requires a complex, ri image");
02831 int nxreal = nx -2 + int(is_fftodd());
02832 int nxhalf = nxreal/2;
02833 for (int iz = 0; iz < nz; iz++)
02834 for (int iy = 0; iy < ny; iy++)
02835 for (int ix = 0; ix <= nxhalf; ix++)
02836 cmplx(ix,iy,iz) = conj(cmplx(ix,iy,iz));
02837 EXITFUNC;
02838 }
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851
02852 void EMData::update_stat() const
02853 {
02854 ENTERFUNC;
02855
02856 if (!(flags & EMDATA_NEEDUPD))
02857 {
02858 EXITFUNC;
02859 return;
02860 }
02861
02862 float* data = get_data();
02863 float max = -FLT_MAX;
02864 float min = -max;
02865
02866 double sum = 0;
02867 double square_sum = 0;
02868
02869 int step = 1;
02870 if (is_complex() && !is_ri()) {
02871 step = 2;
02872 }
02873
02874 int n_nonzero = 0;
02875
02876 size_t size = (size_t)nx*ny*nz;
02877 for (size_t i = 0; i < size; i += step) {
02878 float v = data[i];
02879 #ifdef _WIN32
02880 max = _cpp_max(max,v);
02881 min = _cpp_min(min,v);
02882 #else
02883 max=std::max<float>(max,v);
02884 min=std::min<float>(min,v);
02885 #endif //_WIN32
02886 sum += v;
02887 square_sum += v * (double)(v);
02888 if (v != 0) n_nonzero++;
02889 }
02890
02891 size_t n = size / step;
02892 double mean = sum / n;
02893
02894 #ifdef _WIN32
02895 float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / n)/(n-1)));
02896 n_nonzero = _cpp_max(1,n_nonzero);
02897 double sigma_nonzero = std::sqrt( _cpp_max(0,(square_sum - sum*sum/n_nonzero)/(n_nonzero-1)));
02898 #else
02899 float sigma = (float)std::sqrt(std::max<double>(0.0,(square_sum - sum*sum / n)/(n-1)));
02900 n_nonzero = std::max<int>(1,n_nonzero);
02901 double sigma_nonzero = std::sqrt(std::max<double>(0,(square_sum - sum*sum/n_nonzero)/(n_nonzero-1)));
02902 #endif //_WIN32
02903 double mean_nonzero = sum / n_nonzero;
02904
02905 attr_dict["minimum"] = min;
02906 attr_dict["maximum"] = max;
02907 attr_dict["mean"] = (float)(mean);
02908 attr_dict["sigma"] = (float)(sigma);
02909 attr_dict["square_sum"] = (float)(square_sum);
02910 attr_dict["mean_nonzero"] = (float)(mean_nonzero);
02911 attr_dict["sigma_nonzero"] = (float)(sigma_nonzero);
02912 attr_dict["is_complex"] = (int) is_complex();
02913 attr_dict["is_complex_ri"] = (int) is_ri();
02914
02915 flags &= ~EMDATA_NEEDUPD;
02916
02917 if (rot_fp != 0)
02918 {
02919 delete rot_fp; rot_fp = 0;
02920 }
02921
02922 EXITFUNC;
02923
02924 }
02925
02926 bool EMData::operator==(const EMData& that) const {
02927 if (that.get_xsize() != nx || that.get_ysize() != ny || that.get_zsize() != nz ) return false;
02928
02929 const float* d1 = that.get_const_data();
02930 float* d2 = get_data();
02931
02932 for(size_t i =0; i < get_size(); ++i,++d1,++d2) {
02933 if ((*d1) != (*d2)) return false;
02934 }
02935 return true;
02936
02937 }
02938
02939 EMData * EMAN::operator+(const EMData & em, float n)
02940 {
02941 EMData * r = em.copy();
02942 r->add(n);
02943 return r;
02944 }
02945
02946 EMData * EMAN::operator-(const EMData & em, float n)
02947 {
02948 EMData* r = em.copy();
02949 r->sub(n);
02950 return r;
02951 }
02952
02953 EMData * EMAN::operator*(const EMData & em, float n)
02954 {
02955 EMData* r = em.copy();
02956 r ->mult(n);
02957 return r;
02958 }
02959
02960 EMData * EMAN::operator/(const EMData & em, float n)
02961 {
02962 EMData * r = em.copy();
02963 r->div(n);
02964 return r;
02965 }
02966
02967
02968 EMData * EMAN::operator+(float n, const EMData & em)
02969 {
02970 EMData * r = em.copy();
02971 r->add(n);
02972 return r;
02973 }
02974
02975 EMData * EMAN::operator-(float n, const EMData & em)
02976 {
02977 EMData * r = em.copy();
02978 r->mult(-1.0f);
02979 r->add(n);
02980 return r;
02981 }
02982
02983 EMData * EMAN::operator*(float n, const EMData & em)
02984 {
02985 EMData * r = em.copy();
02986 r->mult(n);
02987 return r;
02988 }
02989
02990 EMData * EMAN::operator/(float n, const EMData & em)
02991 {
02992 EMData * r = em.copy();
02993 r->to_one();
02994 r->mult(n);
02995 r->div(em);
02996
02997 return r;
02998 }
02999
03000 EMData * EMAN::rsub(const EMData & em, float n)
03001 {
03002 return EMAN::operator-(n, em);
03003 }
03004
03005 EMData * EMAN::rdiv(const EMData & em, float n)
03006 {
03007 return EMAN::operator/(n, em);
03008 }
03009
03010 EMData * EMAN::operator+(const EMData & a, const EMData & b)
03011 {
03012 EMData * r = a.copy();
03013 r->add(b);
03014 return r;
03015 }
03016
03017 EMData * EMAN::operator-(const EMData & a, const EMData & b)
03018 {
03019 EMData * r = a.copy();
03020 r->sub(b);
03021 return r;
03022 }
03023
03024 EMData * EMAN::operator*(const EMData & a, const EMData & b)
03025 {
03026 EMData * r = a.copy();
03027 r->mult(b);
03028 return r;
03029 }
03030
03031 EMData * EMAN::operator/(const EMData & a, const EMData & b)
03032 {
03033 EMData * r = a.copy();
03034 r->div(b);
03035 return r;
03036 }
03037
03038 void EMData::set_xyz_origin(float origin_x, float origin_y, float origin_z)
03039 {
03040 attr_dict["origin_x"] = origin_x;
03041 attr_dict["origin_y"] = origin_y;
03042 attr_dict["origin_z"] = origin_z;
03043 }
03044
03045 #if 0
03046 void EMData::calc_rcf(EMData * with, vector < float >&sum_array)
03047 {
03048 ENTERFUNC;
03049
03050 int array_size = sum_array.size();
03051 float da = 2 * M_PI / array_size;
03052 float *dat = new float[array_size + 2];
03053 float *dat2 = new float[array_size + 2];
03054 int nx2 = nx * 9 / 20;
03055
03056 float lim = 0;
03057 if (fabs(translation[0]) < fabs(translation[1])) {
03058 lim = fabs(translation[1]);
03059 }
03060 else {
03061 lim = fabs(translation[0]);
03062 }
03063
03064 nx2 -= static_cast < int >(floor(lim));
03065
03066 for (int i = 0; i < array_size; i++) {
03067 sum_array[i] = 0;
03068 }
03069
03070 float sigma = attr_dict["sigma"];
03071 float with_sigma = with->get_attr_dict().get("sigma");
03072
03073 vector<float> vdata, vdata2;
03074 for (int i = 8; i < nx2; i += 6) {
03075 vdata = calc_az_dist(array_size, 0, da, i, i + 6);
03076 vdata2 = with->calc_az_dist(array_size, 0, da, i, i + 6);
03077 Assert(vdata.size() <= array_size + 2);
03078 Assert(cdata2.size() <= array_size + 2);
03079 std::copy(vdata.begin(), vdata.end(), dat);
03080 std::copy(vdata2.begin(), vdata2.end(), dat2);
03081
03082 EMfft::real_to_complex_1d(dat, dat, array_size);
03083 EMfft::real_to_complex_1d(dat2, dat2, array_size);
03084
03085 for (int j = 0; j < array_size + 2; j += 2) {
03086 float max = dat[j] * dat2[j] + dat[j + 1] * dat2[j + 1];
03087 float max2 = dat[j + 1] * dat2[j] - dat2[j + 1] * dat[j];
03088 dat[j] = max;
03089 dat[j + 1] = max2;
03090 }
03091
03092 EMfft::complex_to_real_1d(dat, dat, array_size);
03093 float norm = array_size * array_size * (4.0f * sigma) * (4.0f * with_sigma);
03094
03095 for (int j = 0; j < array_size; j++) {
03096 sum_array[j] += dat[j] * (float) i / norm;
03097 }
03098 }
03099
03100 if( dat )
03101 {
03102 delete[]dat;
03103 dat = 0;
03104 }
03105
03106 if( dat2 )
03107 {
03108 delete[]dat2;
03109 dat2 = 0;
03110 }
03111 EXITFUNC;
03112 }
03113
03114 #endif
03115
03116 void EMData::add_incoherent(EMData * obj)
03117 {
03118 ENTERFUNC;
03119
03120 if (!obj) {
03121 LOGERR("NULL image");
03122 throw NullPointerException("NULL image");
03123 }
03124
03125 if (!obj->is_complex() || !is_complex()) {
03126 throw ImageFormatException("complex images only");
03127 }
03128
03129 if (!EMUtil::is_same_size(this, obj)) {
03130 throw ImageFormatException("images not same size");
03131 }
03132
03133 ri2ap();
03134 obj->ri2ap();
03135
03136 float *dest = get_data();
03137 float *src = obj->get_data();
03138 size_t size = (size_t)nx * ny * nz;
03139 for (size_t j = 0; j < size; j += 2) {
03140 #ifdef _WIN32
03141 dest[j] = (float) _hypot(src[j], dest[j]);
03142 #else
03143 dest[j] = (float) hypot(src[j], dest[j]);
03144 #endif //_WIN32
03145 dest[j + 1] = 0;
03146 }
03147
03148 obj->update();
03149 update();
03150 EXITFUNC;
03151 }
03152
03153
03154 float EMData::calc_dist(EMData * second_img, int y_index) const
03155 {
03156 ENTERFUNC;
03157
03158 if (get_ndim() != 1) {
03159 throw ImageDimensionException("'this' image is 1D only");
03160 }
03161
03162 if (second_img->get_xsize() != nx || ny != 1) {
03163 throw ImageFormatException("image xsize not same");
03164 }
03165
03166 if (y_index > second_img->get_ysize() || y_index < 0) {
03167 return -1;
03168 }
03169
03170 float ret = 0;
03171 float *d1 = get_data();
03172 float *d2 = second_img->get_data() + second_img->get_xsize() * y_index;
03173
03174 for (int i = 0; i < nx; i++) {
03175 ret += Util::square(d1[i] - d2[i]);
03176 }
03177 EXITFUNC;
03178 return std::sqrt(ret);
03179 }
03180
03181
03182 EMData * EMData::calc_fast_sigma_image( EMData* mask)
03183 {
03184 ENTERFUNC;
03185
03186 bool maskflag = false;
03187 if (mask == 0) {
03188 mask = new EMData(nx,ny,nz);
03189 mask->process_inplace("testimage.circlesphere");
03190 maskflag = true;
03191 }
03192
03193 if (get_ndim() != mask->get_ndim() ) throw ImageDimensionException("The dimensions do not match");
03194
03195 int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
03196
03197 if ( mnx > nx || mny > ny || mnz > nz)
03198 throw ImageDimensionException("Can not calculate variance map using an image that is larger than this image");
03199
03200 size_t P = 0;
03201 for(size_t i = 0; i < mask->get_size(); ++i){
03202 if (mask->get_value_at(i) != 0){
03203 ++P;
03204 }
03205 }
03206 float normfac = 1.0f/(float)P;
03207
03208
03209
03210 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03211
03212 Region r;
03213 if (ny == 1) r = Region((mnx-nxc)/2,nxc);
03214 else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03215 else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03216 mask->clip_inplace(r,0.0);
03217
03218
03219
03220
03221
03222
03223 Region r2;
03224 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03225 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03226 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03227 EMData* squared = get_clip(r2,get_edge_mean());
03228
03229 EMData* tmp = squared->copy();
03230 Dict pow;
03231 pow["pow"] = 2.0f;
03232 squared->process_inplace("math.pow",pow);
03233 EMData* s = mask->convolute(squared);
03234 squared->mult(normfac);
03235
03236 EMData* m = mask->convolute(tmp);
03237 m->mult(normfac);
03238 m->process_inplace("math.pow",pow);
03239 delete tmp; tmp = 0;
03240 s->sub(*m);
03241
03242 s->process_inplace("math.sqrt");
03243
03244
03245
03246
03247
03248
03249 if (maskflag) {
03250 delete mask;
03251 mask = 0;
03252 } else {
03253 Region r;
03254 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
03255 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
03256 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
03257 mask->clip_inplace(r);
03258 }
03259
03260 delete squared;
03261 delete m;
03262
03263 s->process_inplace("xform.phaseorigin.tocenter");
03264 Region r3;
03265 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03266 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03267 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03268 s->clip_inplace(r3);
03269 EXITFUNC;
03270 return s;
03271 }
03272
03273
03274
03275 EMData *EMData::calc_flcf(EMData * with)
03276 {
03277 ENTERFUNC;
03278 EMData *this_copy=this;
03279 this_copy=copy();
03280
03281 int mnx = with->get_xsize(); int mny = with->get_ysize(); int mnz = with->get_zsize();
03282 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
03283
03284
03285 EMData* ones = new EMData(mnx,mny,mnz);
03286 ones->process_inplace("testimage.circlesphere");
03287
03288
03289 EMData* with_resized = with->copy();
03290 with_resized->process_inplace("normalize");
03291 with_resized->mult(*ones);
03292
03293 EMData* s = calc_fast_sigma_image(ones);
03294
03295 Region r1;
03296 if (ny == 1) r1 = Region((mnx-nxc)/2,nxc);
03297 else if (nz == 1) r1 = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
03298 else r1 = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
03299 with_resized->clip_inplace(r1,0.0);
03300
03301 Region r2;
03302 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
03303 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
03304 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
03305 this_copy->clip_inplace(r2,0.0);
03306
03307 EMData* corr = this_copy->calc_ccf(with_resized);
03308
03309 corr->process_inplace("xform.phaseorigin.tocenter");
03310 Region r3;
03311 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
03312 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
03313 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
03314 corr->clip_inplace(r3);
03315
03316 corr->div(*s);
03317
03318 delete with_resized; delete ones; delete this_copy; delete s;
03319 EXITFUNC;
03320 return corr;
03321 }
03322
03323 EMData *EMData::convolute(EMData * with)
03324 {
03325 ENTERFUNC;
03326
03327 EMData *f1 = do_fft();
03328 if (!f1) {
03329 LOGERR("FFT returns NULL image");
03330 throw NullPointerException("FFT returns NULL image");
03331 }
03332
03333 f1->ap2ri();
03334
03335 EMData *cf = 0;
03336 if (with) {
03337 cf = with->do_fft();
03338 if (!cf) {
03339 LOGERR("FFT returns NULL image");
03340 throw NullPointerException("FFT returns NULL image");
03341 }
03342 cf->ap2ri();
03343 }
03344 else {
03345 cf = f1->copy();
03346 }
03347
03348 if (with && !EMUtil::is_same_size(f1, cf)) {
03349 LOGERR("images not same size");
03350 throw ImageFormatException("images not same size");
03351 }
03352
03353 float *rdata1 = f1->get_data();
03354 float *rdata2 = cf->get_data();
03355 size_t cf_size = (size_t)cf->get_xsize() * cf->get_ysize() * cf->get_zsize();
03356
03357 float re,im;
03358
03359 for (size_t i = 0; i < cf_size; i += 2) {
03360 re = rdata1[i] * rdata2[i] - rdata1[i + 1] * rdata2[i + 1];
03361 im = rdata1[i + 1] * rdata2[i] + rdata1[i] * rdata2[i + 1];
03362 rdata2[i]=re;
03363 rdata2[i+1]=im;
03364 }
03365 cf->update();
03366 EMData *f2 = cf->do_ift();
03367
03368 if( cf )
03369 {
03370 delete cf;
03371 cf = 0;
03372 }
03373
03374 if( f1 )
03375 {
03376 delete f1;
03377 f1=0;
03378 }
03379
03380 EXITFUNC;
03381 return f2;
03382 }
03383
03384
03385 void EMData::common_lines(EMData * image1, EMData * image2,
03386 int mode, int steps, bool horizontal)
03387 {
03388 ENTERFUNC;
03389
03390 if (!image1 || !image2) {
03391 throw NullPointerException("NULL image");
03392 }
03393
03394 if (mode < 0 || mode > 2) {
03395 throw OutofRangeException(0, 2, mode, "invalid mode");
03396 }
03397
03398 if (!image1->is_complex()) {
03399 image1 = image1->do_fft();
03400 }
03401 if (!image2->is_complex()) {
03402 image2 = image2->do_fft();
03403 }
03404
03405 image1->ap2ri();
03406 image2->ap2ri();
03407
03408 if (!EMUtil::is_same_size(image1, image2)) {
03409 throw ImageFormatException("images not same sizes");
03410 }
03411
03412 int image2_nx = image2->get_xsize();
03413 int image2_ny = image2->get_ysize();
03414
03415 int rmax = image2_ny / 4 - 1;
03416 int array_size = steps * rmax * 2;
03417 float *im1 = new float[array_size];
03418 float *im2 = new float[array_size];
03419 for (int i = 0; i < array_size; i++) {
03420 im1[i] = 0;
03421 im2[i] = 0;
03422 }
03423
03424 set_size(steps * 2, steps * 2, 1);
03425
03426 float *image1_data = image1->get_data();
03427 float *image2_data = image2->get_data();
03428
03429 float da = M_PI / steps;
03430 float a = -M_PI / 2.0f + da / 2.0f;
03431 int jmax = 0;
03432
03433 for (int i = 0; i < steps * 2; i += 2, a += da) {
03434 float s1 = 0;
03435 float s2 = 0;
03436 int i2 = i * rmax;
03437 int j = 0;
03438
03439 for (float r = 3.0f; r < rmax - 3.0f; j += 2, r += 1.0f) {
03440 float x = r * cos(a);
03441 float y = r * sin(a);
03442
03443 if (x < 0) {
03444 x = -x;
03445 y = -y;
03446 LOGERR("CCL ERROR %d, %f !\n", i, -x);
03447 }
03448
03449 int k = (int) (floor(x) * 2 + floor(y + image2_ny / 2) * image2_nx);
03450 int l = i2 + j;
03451 float x2 = x - floor(x);
03452 float y2 = y - floor(y);
03453
03454 im1[l] = Util::bilinear_interpolate(image1_data[k],
03455 image1_data[k + 2],
03456 image1_data[k + image2_nx],
03457 image1_data[k + 2 + image2_nx], x2, y2);
03458
03459 im2[l] = Util::bilinear_interpolate(image2_data[k],
03460 image2_data[k + 2],
03461 image2_data[k + image2_nx],
03462 image2_data[k + 2 + image2_nx], x2, y2);
03463
03464 k++;
03465
03466 im1[l + 1] = Util::bilinear_interpolate(image1_data[k],
03467 image1_data[k + 2],
03468 image1_data[k + image2_nx],
03469 image1_data[k + 2 + image2_nx], x2, y2);
03470
03471 im2[l + 1] = Util::bilinear_interpolate(image2_data[k],
03472 image2_data[k + 2],
03473 image2_data[k + image2_nx],
03474 image2_data[k + 2 + image2_nx], x2, y2);
03475
03476 s1 += Util::square_sum(im1[l], im1[l + 1]);
03477 s2 += Util::square_sum(im2[l], im2[l + 1]);
03478 }
03479
03480 jmax = j - 1;
03481 float sqrt_s1 = std::sqrt(s1);
03482 float sqrt_s2 = std::sqrt(s2);
03483
03484 int l = 0;
03485 for (float r = 1; r < rmax; r += 1.0f) {
03486 int i3 = i2 + l;
03487 im1[i3] /= sqrt_s1;
03488 im1[i3 + 1] /= sqrt_s1;
03489 im2[i3] /= sqrt_s2;
03490 im2[i3 + 1] /= sqrt_s2;
03491 l += 2;
03492 }
03493 }
03494 float * data = get_data();
03495
03496 switch (mode) {
03497 case 0:
03498 for (int m1 = 0; m1 < 2; m1++) {
03499 for (int m2 = 0; m2 < 2; m2++) {
03500
03501 if (m1 == 0 && m2 == 0) {
03502 for (int i = 0; i < steps; i++) {
03503 int i2 = i * rmax * 2;
03504 for (int j = 0; j < steps; j++) {
03505 int l = i + j * steps * 2;
03506 int j2 = j * rmax * 2;
03507 data[l] = 0;
03508 for (int k = 0; k < jmax; k++) {
03509 data[l] += im1[i2 + k] * im2[j2 + k];
03510 }
03511 }
03512 }
03513 }
03514 else {
03515 int steps2 = steps * m2 + steps * steps * 2 * m1;
03516
03517 for (int i = 0; i < steps; i++) {
03518 int i2 = i * rmax * 2;
03519 for (int j = 0; j < steps; j++) {
03520 int j2 = j * rmax * 2;
03521 int l = i + j * steps * 2 + steps2;
03522 data[l] = 0;
03523
03524 for (int k = 0; k < jmax; k += 2) {
03525 i2 += k;
03526 j2 += k;
03527 data[l] += im1[i2] * im2[j2];
03528 data[l] += -im1[i2 + 1] * im2[j2 + 1];
03529 }
03530 }
03531 }
03532 }
03533 }
03534 }
03535
03536 break;
03537 case 1:
03538 for (int m1 = 0; m1 < 2; m1++) {
03539 for (int m2 = 0; m2 < 2; m2++) {
03540 int steps2 = steps * m2 + steps * steps * 2 * m1;
03541 int p1_sign = 1;
03542 if (m1 != m2) {
03543 p1_sign = -1;
03544 }
03545
03546 for (int i = 0; i < steps; i++) {
03547 int i2 = i * rmax * 2;
03548
03549 for (int j = 0; j < steps; j++) {
03550 int j2 = j * rmax * 2;
03551
03552 int l = i + j * steps * 2 + steps2;
03553 data[l] = 0;
03554 float a = 0;
03555
03556 for (int k = 0; k < jmax; k += 2) {
03557 i2 += k;
03558 j2 += k;
03559
03560 #ifdef _WIN32
03561 float a1 = (float) _hypot(im1[i2], im1[i2 + 1]);
03562 #else
03563 float a1 = (float) hypot(im1[i2], im1[i2 + 1]);
03564 #endif //_WIN32
03565 float p1 = atan2(im1[i2 + 1], im1[i2]);
03566 float p2 = atan2(im2[j2 + 1], im2[j2]);
03567
03568 data[l] += Util::angle_sub_2pi(p1_sign * p1, p2) * a1;
03569 a += a1;
03570 }
03571
03572 data[l] /= (float)(a * M_PI / 180.0f);
03573 }
03574 }
03575 }
03576 }
03577
03578 break;
03579 case 2:
03580 for (int m1 = 0; m1 < 2; m1++) {
03581 for (int m2 = 0; m2 < 2; m2++) {
03582 int steps2 = steps * m2 + steps * steps * 2 * m1;
03583
03584 for (int i = 0; i < steps; i++) {
03585 int i2 = i * rmax * 2;
03586
03587 for (int j = 0; j < steps; j++) {
03588 int j2 = j * rmax * 2;
03589 int l = i + j * steps * 2 + steps2;
03590 data[l] = 0;
03591
03592 for (int k = 0; k < jmax; k += 2) {
03593 i2 += k;
03594 j2 += k;
03595 #ifdef _WIN32
03596 data[l] += (float) (_hypot(im1[i2], im1[i2 + 1]) * _hypot(im2[j2], im2[j2 + 1]));
03597 #else
03598 data[l] += (float) (hypot(im1[i2], im1[i2 + 1]) * hypot(im2[j2], im2[j2 + 1]));
03599 #endif //_WIN32
03600 }
03601 }
03602 }
03603 }
03604 }
03605
03606 break;
03607 default:
03608 break;
03609 }
03610
03611 if (horizontal) {
03612 float *tmp_array = new float[ny];
03613 for (int i = 1; i < nx; i++) {
03614 for (int j = 0; j < ny; j++) {
03615 tmp_array[j] = get_value_at(i, j);
03616 }
03617 for (int j = 0; j < ny; j++) {
03618 set_value_at(i, j, 0, tmp_array[(j + i) % ny]);
03619 }
03620 }
03621 if( tmp_array )
03622 {
03623 delete[]tmp_array;
03624 tmp_array = 0;
03625 }
03626 }
03627
03628 if( im1 )
03629 {
03630 delete[]im1;
03631 im1 = 0;
03632 }
03633
03634 if( im2 )
03635 {
03636 delete im2;
03637 im2 = 0;
03638 }
03639
03640
03641 image1->update();
03642 image2->update();
03643 if( image1 )
03644 {
03645 delete image1;
03646 image1 = 0;
03647 }
03648 if( image2 )
03649 {
03650 delete image2;
03651 image2 = 0;
03652 }
03653 update();
03654 EXITFUNC;
03655 }
03656
03657
03658
03659 void EMData::common_lines_real(EMData * image1, EMData * image2,
03660 int steps, bool horiz)
03661 {
03662 ENTERFUNC;
03663
03664 if (!image1 || !image2) {
03665 throw NullPointerException("NULL image");
03666 }
03667
03668 if (!EMUtil::is_same_size(image1, image2)) {
03669 throw ImageFormatException("images not same size");
03670 }
03671
03672 int steps2 = steps * 2;
03673 int image_ny = image1->get_ysize();
03674 EMData *image1_copy = image1->copy();
03675 EMData *image2_copy = image2->copy();
03676
03677 float *im1 = new float[steps2 * image_ny];
03678 float *im2 = new float[steps2 * image_ny];
03679
03680 EMData *images[] = { image1_copy, image2_copy };
03681 float *ims[] = { im1, im2 };
03682
03683 for (int m = 0; m < 2; m++) {
03684 float *im = ims[m];
03685 float a = M_PI / steps2;
03686 Transform t(Dict("type","2d","alpha",-a));
03687 for (int i = 0; i < steps2; i++) {
03688 images[i]->transform(t);
03689 float *data = images[i]->get_data();
03690
03691 for (int j = 0; j < image_ny; j++) {
03692 float sum = 0;
03693 for (int k = 0; k < image_ny; k++) {
03694 sum += data[j * image_ny + k];
03695 }
03696 im[i * image_ny + j] = sum;
03697 }
03698
03699 float sum1 = 0;
03700 float sum2 = 0;
03701 for (int j = 0; j < image_ny; j++) {
03702 int l = i * image_ny + j;
03703 sum1 += im[l];
03704 sum2 += im[l] * im[l];
03705 }
03706
03707 float mean = sum1 / image_ny;
03708 float sigma = std::sqrt(sum2 / image_ny - sum1 * sum1);
03709
03710 for (int j = 0; j < image_ny; j++) {
03711 int l = i * image_ny + j;
03712 im[l] = (im[l] - mean) / sigma;
03713 }
03714
03715 images[i]->update();
03716 a += M_PI / steps;
03717 }
03718 }
03719
03720 set_size(steps2, steps2, 1);
03721 float *data1 = get_data();
03722
03723 if (horiz) {
03724 for (int i = 0; i < steps2; i++) {
03725 for (int j = 0, l = i; j < steps2; j++, l++) {
03726 if (l == steps2) {
03727 l = 0;
03728 }
03729
03730 float sum = 0;
03731 for (int k = 0; k < image_ny; k++) {
03732 sum += im1[i * image_ny + k] * im2[l * image_ny + k];
03733 }
03734 data1[i + j * steps2] = sum;
03735 }
03736 }
03737 }
03738 else {
03739 for (int i = 0; i < steps2; i++) {
03740 for (int j = 0; j < steps2; j++) {
03741 float sum = 0;
03742 for (int k = 0; k < image_ny; k++) {
03743 sum += im1[i * image_ny + k] * im2[j * image_ny + k];
03744 }
03745 data1[i + j * steps2] = sum;
03746 }
03747 }
03748 }
03749
03750 update();
03751
03752 if( image1_copy )
03753 {
03754 delete image1_copy;
03755 image1_copy = 0;
03756 }
03757
03758 if( image2_copy )
03759 {
03760 delete image2_copy;
03761 image2_copy = 0;
03762 }
03763
03764 if( im1 )
03765 {
03766 delete[]im1;
03767 im1 = 0;
03768 }
03769
03770 if( im2 )
03771 {
03772 delete[]im2;
03773 im2 = 0;
03774 }
03775 EXITFUNC;
03776 }
03777
03778
03779 void EMData::cut_slice(const EMData *const map, const Transform& transform, bool interpolate)
03780 {
03781 ENTERFUNC;
03782
03783 if (!map) throw NullPointerException("NULL image");
03784
03785 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
03786 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
03787
03788 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
03789 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
03790
03791
03792 float *sdata = map->get_data();
03793 float *ddata = get_data();
03794
03795 int map_nx = map->get_xsize();
03796 int map_ny = map->get_ysize();
03797 int map_nz = map->get_zsize();
03798 int map_nxy = map_nx * map_ny;
03799
03800 int ymax = ny/2;
03801 if ( ny % 2 == 1 ) ymax += 1;
03802 int xmax = nx/2;
03803 if ( nx % 2 == 1 ) xmax += 1;
03804 for (int y = -ny/2; y < ymax; y++) {
03805 for (int x = -nx/2; x < xmax; x++) {
03806 Vec3f coord(x,y,0);
03807 Vec3f soln = transform*coord;
03808
03809
03810
03811
03812
03813
03814
03815
03816
03817
03818 float xx = soln[0]+map_nx/2;
03819 float yy = soln[1]+map_ny/2;
03820 float zz = soln[2]+map_nz/2;
03821
03822 int l = (x+nx/2) + (y+ny/2) * nx;
03823
03824 float t = xx - floor(xx);
03825 float u = yy - floor(yy);
03826 float v = zz - floor(zz);
03827
03828 if (xx < 0 || yy < 0 || zz < 0 ) {
03829 ddata[l] = 0;
03830 continue;
03831 }
03832 if (interpolate) {
03833 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03834 ddata[l] = 0;
03835 continue;
03836 }
03837 int k = (int) (Util::fast_floor(xx) + Util::fast_floor(yy) * map_nx + Util::fast_floor(zz) * map_nxy);
03838
03839
03840 if (xx < (map_nx - 1) && yy < (map_ny - 1) && zz < (map_nz - 1)) {
03841 ddata[l] = Util::trilinear_interpolate(sdata[k],
03842 sdata[k + 1], sdata[k + map_nx],sdata[k + map_nx + 1],
03843 sdata[k + map_nxy], sdata[k + map_nxy + 1], sdata[k + map_nx + map_nxy],
03844 sdata[k + map_nx + map_nxy + 1],t, u, v);
03845 }
03846 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03847 ddata[l] += sdata[k];
03848 }
03849 else if ( xx == (map_nx - 1) && yy == (map_ny - 1) ) {
03850 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nxy],v);
03851 }
03852 else if ( xx == (map_nx - 1) && zz == (map_nz - 1) ) {
03853 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + map_nx],u);
03854 }
03855 else if ( yy == (map_ny - 1) && zz == (map_nz - 1) ) {
03856 ddata[l] += Util::linear_interpolate(sdata[k], sdata[k + 1],t);
03857 }
03858 else if ( xx == (map_nx - 1) ) {
03859 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + map_nx], sdata[k + map_nxy], sdata[k + map_nxy + map_nx],u,v);
03860 }
03861 else if ( yy == (map_ny - 1) ) {
03862 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nxy], sdata[k + map_nxy + 1],t,v);
03863 }
03864 else if ( zz == (map_nz - 1) ) {
03865 ddata[l] += Util::bilinear_interpolate(sdata[k], sdata[k + 1], sdata[k + map_nx], sdata[k + map_nx + 1],t,u);
03866 }
03867
03868
03869
03870
03871
03872
03873
03874
03875
03876
03877
03878
03879 }
03880 else {
03881 if ( xx > map_nx - 1 || yy > map_ny - 1 || zz > map_nz - 1) {
03882 ddata[l] = 0;
03883 continue;
03884 }
03885 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
03886 ddata[l] = sdata[k];
03887 }
03888
03889 }
03890 }
03891
03892 update();
03893
03894 EXITFUNC;
03895 }
03896
03897 EMData *EMData::unwrap_largerR(int r1,int r2,int xs, float rmax_f) {
03898 float *d,*dd;
03899 int do360=2;
03900 int rmax = (int)(rmax_f+0.5f);
03901 unsigned long i;
03902 unsigned int nvox=get_xsize()*get_ysize();
03903 float maxmap=0.0f, minmap=0.0f;
03904 float temp=0.0f, diff_den=0.0f, norm=0.0f;
03905 float cut_off_va =0.0f;
03906
03907 d=get_data();
03908 maxmap=-1000000.0f;
03909 minmap=1000000.0f;
03910 for (i=0;i<nvox;i++){
03911 if(d[i]>maxmap) maxmap=d[i];
03912 if(d[i]<minmap) minmap=d[i];
03913 }
03914 diff_den = maxmap-minmap;
03915 for(i=0;i<nvox;i++) {
03916 temp = (d[i]-minmap)/diff_den;
03917 if(cut_off_va != 0.0) {
03918 if(temp < cut_off_va)
03919 d[i] = 0.0f;
03920 else
03921 d[i] = temp-cut_off_va;
03922 }
03923 else d[i] = temp;
03924 }
03925
03926 for(i=0;i<nvox;i++) {
03927 temp=d[i];
03928 norm += temp*temp;
03929 }
03930 for(i=0;i<nvox;i++) d[i] /= norm;
03931
03932 if (xs<1) {
03933 xs = (int) floor(do360*M_PI*get_ysize()/4);
03934 xs=Util::calc_best_fft_size(xs);
03935 }
03936 if (r1<0) r1=0;
03937 float maxext=ceil(0.6f*std::sqrt((float)(get_xsize()*get_xsize()+get_ysize()*get_ysize())));
03938
03939 if (r2<r1) r2=(int)maxext;
03940 EMData *ret = new EMData;
03941
03942 ret->set_size(xs,r2+1,1);
03943
03944 dd=ret->get_data();
03945
03946 for (int i=0; i<xs; i++) {
03947 float si=sin(i*M_PI*2/xs);
03948 float co=cos(i*M_PI*2/xs);
03949 for (int r=0; r<=maxext; r++) {
03950 float x=(r+r1)*co+get_xsize()/2;
03951 float y=(r+r1)*si+get_ysize()/2;
03952 if(x<0.0 || x>=get_xsize()-1.0 || y<0.0 || y>=get_ysize()-1.0 || r>rmax){
03953 for(;r<=r2;r++)
03954 dd[i+r*xs]=0.0;
03955 break;
03956 }
03957 int x1=(int)floor(x);
03958 int y1=(int)floor(y);
03959 float t=x-x1;
03960 float u=y-y1;
03961 float f11= d[x1+y1*get_xsize()];
03962 float f21= d[(x1+1)+y1*get_xsize()];
03963 float f12= d[x1+(y1+1)*get_xsize()];
03964 float f22= d[(x1+1)+(y1+1)*get_xsize()];
03965 dd[i+r*xs] = (1-t)*(1-u)*f11+t*(1-u)*f21+t*u*f22+(1-t)*u*f12;
03966 }
03967 }
03968 update();
03969 ret->update();
03970 return ret;
03971 }
03972
03973
03974 EMData *EMData::oneDfftPolar(int size, float rmax, float MAXR){
03975 float *pcs=get_data();
03976 EMData *imagepcsfft = new EMData;
03977 imagepcsfft->set_size((size+2), (int)MAXR+1, 1);
03978 float *d=imagepcsfft->get_data();
03979
03980 EMData *data_in=new EMData;
03981 data_in->set_size(size,1,1);
03982 float *in=data_in->get_data();
03983
03984 for(int row=0; row<=(int)MAXR; ++row){
03985 if(row<=(int)rmax) {
03986 for(int i=0; i<size;++i) in[i] = pcs[i+row*size];
03987 data_in->set_complex(false);
03988 data_in->do_fft_inplace();
03989 for(int j=0;j<size+2;j++) d[j+row*(size+2)]=in[j];
03990 }
03991 else for(int j=0;j<size+2;j++) d[j+row*(size+2)]=0.0;
03992 }
03993 imagepcsfft->update();
03994 delete data_in;
03995 return imagepcsfft;
03996 }
03997
03998
03999
04000
04001
04002
04003
04004
04005
04006
04007
04008
04009
04010
04011
04012
04013
04014
04015
04016
04017
04018
04019
04020
04021
04022
04023
04024
04025
04026
04027
04028
04029
04030
04031 void EMData::uncut_slice(EMData * const map, const Transform& transform) const
04032 {
04033 ENTERFUNC;
04034
04035 if (!map) throw NullPointerException("NULL image");
04036
04037 if ( get_ndim() != 2 ) throw ImageDimensionException("Can not call cut slice on an image that is not 2D");
04038 if ( map->get_ndim() != 3 ) throw ImageDimensionException("Can not cut slice from an image that is not 3D");
04039
04040 if ( is_complex() ) throw ImageFormatException("Can not call cut slice on an image that is complex");
04041 if ( map->is_complex() ) throw ImageFormatException("Can not cut slice from a complex image");
04042
04043
04044
04045
04046
04047
04048 float *ddata = map->get_data();
04049 float *sdata = get_data();
04050
04051 int map_nx = map->get_xsize();
04052 int map_ny = map->get_ysize();
04053 int map_nz = map->get_zsize();
04054 int map_nxy = map_nx * map_ny;
04055 float map_nz_round_limit = (float) map_nz-0.5f;
04056 float map_ny_round_limit = (float) map_ny-0.5f;
04057 float map_nx_round_limit = (float) map_nx-0.5f;
04058
04059
04060
04061
04062 int ymax = ny/2;
04063 if ( ny % 2 == 1 ) ymax += 1;
04064 int xmax = nx/2;
04065 if ( nx % 2 == 1 ) xmax += 1;
04066 for (int y = -ny/2; y < ymax; y++) {
04067 for (int x = -nx/2; x < xmax; x++) {
04068 Vec3f coord(x,y,0);
04069 Vec3f soln = transform*coord;
04070
04071
04072
04073
04074
04075
04076
04077
04078 float xx = soln[0]+map_nx/2;
04079 float yy = soln[1]+map_ny/2;
04080 float zz = soln[2]+map_nz/2;
04081
04082
04083 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;
04084
04085 size_t k = Util::round(xx) + Util::round(yy) * map_nx + Util::round(zz) * (size_t)map_nxy;
04086 int l = (x+nx/2) + (y+ny/2) * nx;
04087 ddata[k] = sdata[l];
04088 }
04089 }
04090
04091 map->update();
04092 EXITFUNC;
04093 }
04094
04095
04096 void EMData::save_byteorder_to_dict(ImageIO * imageio)
04097 {
04098 string image_endian = "ImageEndian";
04099 string host_endian = "HostEndian";
04100
04101 if (imageio->is_image_big_endian()) {
04102 attr_dict[image_endian] = "big";
04103 }
04104 else {
04105 attr_dict[image_endian] = "little";
04106 }
04107
04108 if (ByteOrder::is_host_big_endian()) {
04109 attr_dict[host_endian] = "big";
04110 }
04111 else {
04112 attr_dict[host_endian] = "little";
04113 }
04114 }
04115