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