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