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