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 "emfft.h"
00038
00039 #include <cstring>
00040 #include <cstdio>
00041
00042 #include "gsl_sf_result.h"
00043 #include "gsl_sf_bessel.h"
00044 #include <iostream>
00045 #include <algorithm>
00046 #include <vector>
00047 #include <utility>
00048 #include <cmath>
00049 #include "util.h"
00050
00051 #ifdef EMAN2_USING_CUDA
00052 #include "cuda/cuda_processor.h"
00053 #endif
00054
00055 using namespace EMAN;
00056 using namespace std;
00057 typedef vector< pair<float,int> > vp;
00058
00059 #ifdef EMAN2_USING_CUDA
00060
00061 #include "cuda/cuda_emfft.h"
00062
00063 EMData *EMData::do_fft_cuda() const
00064 {
00065 ENTERFUNC;
00066
00067 if ( is_complex() ) {
00068 LOGERR("real image expected. Input image is complex image.");
00069 throw ImageFormatException("real image expected. Input image is complex image.");
00070 }
00071
00072 int offset;
00073 int ndim = get_ndim();
00074 EMData* dat = new EMData();
00075
00076 offset = 2 - nx%2;
00077 dat->set_size_cuda(nx+offset,ny, nz);
00078 float *d = dat->get_cuda_data();
00079 if ( ndim == 1 ) {
00080 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, nx, 1,1);
00081 } else if (ndim == 2) {
00082 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, ny, nx, 1);
00083 } else if (ndim == 3) {
00084 cuda_dd_fft_real_to_complex_nd(get_cuda_data(), d, nz, ny, nx);
00085 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00086
00087 if (offset == 1) dat->set_fftodd(true);
00088 else dat->set_fftodd(false);
00089
00090 dat->set_fftpad(true);
00091 dat->set_complex(true);
00092 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
00093 dat->set_ri(true);
00094 dat->gpu_update();
00095
00096 EXITFUNC;
00097 return dat;
00098 }
00099
00100 EMData *EMData::do_ift_cuda(bool preserve_input) const
00101 {
00102 ENTERFUNC;
00103
00104 if (!is_complex()) {
00105 LOGERR("complex image expected. Input image is real image.");
00106 throw ImageFormatException("complex image expected. Input image is real image.");
00107 }
00108
00109 if (!is_ri()) {
00110 throw ImageFormatException("complex ri expected. Got amplitude/phase.");
00111 }
00112
00113 int offset = is_fftodd() ? 1 : 2;
00114 EMData* dat = new EMData();
00115 int ndim = get_ndim();
00116 dat->set_size_cuda(nx-offset, ny, nz);
00117 float *d = dat->get_cuda_data();
00118 float *this_d;
00119 EMData* tmp = 0;
00120 if (preserve_input){
00121 tmp = new EMData(*this);
00122 this_d = tmp->get_cuda_data();
00123 } else {
00124 this_d = get_cuda_data();
00125 }
00126 if ( ndim == 1 ) {
00127 cuda_dd_fft_complex_to_real_nd(this_d,d, nx-offset,1,1);
00128 } else if (ndim == 2) {
00129 cuda_dd_fft_complex_to_real_nd(this_d,d, ny,nx-offset,1);
00130 } else if (ndim == 3) {
00131 cuda_dd_fft_complex_to_real_nd(this_d,d, nz,ny,nx-offset);
00132 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00133
00134 if (tmp != 0) delete tmp;
00135
00136
00137 float scale = 1.0f/static_cast<float>((dat->get_size()));
00138 dat->mult(scale);
00139
00140 dat->set_fftpad(false);
00141 dat->set_fftodd(false);
00142 dat->set_complex(false);
00143 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
00144 dat->set_ri(false);
00145 dat->gpu_update();
00146
00147 EXITFUNC;
00148 return dat;
00149 }
00150
00151 #endif //EMAN2_USING_CUDA
00152
00153 EMData *EMData::do_fft() const
00154 {
00155 ENTERFUNC;
00156
00157 #ifdef EMAN2_USING_CUDA
00158 if (gpu_operation_preferred()) {
00159 EXITFUNC;
00160 return do_fft_cuda();
00161 }
00162 #endif
00163
00164 if ( is_complex() ) {
00165 LOGERR("real image expected. Input image is complex image.");
00166 throw ImageFormatException("real image expected. Input image is complex image.");
00167 }
00168
00169 int nxreal = nx;
00170 int offset = 2 - nx%2;
00171 int nx2 = nx + offset;
00172 EMData* dat = copy_head();
00173 dat->set_size(nx2, ny, nz);
00174
00175 if (offset == 1) dat->set_fftodd(true);
00176 else dat->set_fftodd(false);
00177
00178 float *d = dat->get_data();
00179
00180 EMfft::real_to_complex_nd(get_data(), d, nxreal, ny, nz);
00181
00182 dat->update();
00183 dat->set_fftpad(true);
00184 dat->set_complex(true);
00185 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
00186 dat->set_ri(true);
00187
00188 EXITFUNC;
00189 return dat;
00190 }
00191
00192 EMData *EMData::do_fft_inplace()
00193 {
00194 ENTERFUNC;
00195
00196 if ( is_complex() ) {
00197 LOGERR("real image expected. Input image is complex image.");
00198 throw ImageFormatException("real image expected. Input image is complex image.");
00199 }
00200
00201 size_t offset;
00202 int nxreal;
00203 get_data();
00204 if (!is_fftpadded()) {
00205
00206
00207 nxreal = nx;
00208 offset = 2 - nx%2;
00209 if (1 == offset) set_fftodd(true);
00210 else set_fftodd(false);
00211 int nxnew = nx + offset;
00212 set_size(nxnew, ny, nz);
00213
00214 for (int iz = nz-1; iz >= 0; iz--) {
00215 for (int iy = ny-1; iy >= 0; iy--) {
00216 for (int ix = nxreal-1; ix >= 0; ix--) {
00217 size_t oldxpos = ix + (iy + iz*ny)*nxreal;
00218 size_t newxpos = ix + (iy + iz*ny)*nxnew;
00219 (*this)(newxpos) = (*this)(oldxpos);
00220 }
00221 }
00222 }
00223 set_fftpad(true);
00224 } else {
00225 offset = is_fftodd() ? 1 : 2;
00226 nxreal = nx - offset;
00227 }
00228 EMfft::real_to_complex_nd(rdata, rdata, nxreal, ny, nz);
00229
00230 set_complex(true);
00231 if(ny==1 && nz==1) set_complex_x(true);
00232 set_ri(true);
00233
00234 update();
00235
00236 EXITFUNC;
00237 return this;
00238 }
00239
00240 EMData *EMData::do_ift()
00241 {
00242 ENTERFUNC;
00243
00244 #ifdef EMAN2_USING_CUDA
00245 if (gpu_operation_preferred()) {
00246 EXITFUNC;
00247 return do_ift_cuda();
00248 }
00249 #endif
00250
00251 if (!is_complex()) {
00252 LOGERR("complex image expected. Input image is real image.");
00253 throw ImageFormatException("complex image expected. Input image is real image.");
00254 }
00255
00256 if (!is_ri()) {
00257 LOGWARN("run IFT on AP data, only RI should be used. Converting.");
00258 }
00259
00260 get_data();
00261 EMData* dat = copy_head();
00262 dat->set_size(nx, ny, nz);
00263 ap2ri();
00264
00265 float *d = dat->get_data();
00266 int ndim = get_ndim();
00267
00268
00269
00270
00271 memcpy((char *) d, (char *) rdata, nx * ny * nz * sizeof(float));
00272
00273 int offset = is_fftodd() ? 1 : 2;
00274
00275 if (ndim == 1) {
00276 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00277 } else {
00278 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00279
00280 #ifndef CUDA_FFT
00281 size_t row_size = (nx - offset) * sizeof(float);
00282 for (int i = 1; i < ny * nz; i++) {
00283 memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
00284 }
00285 #endif
00286 }
00287
00288 dat->set_size(nx - offset, ny, nz);
00289 #if defined FFTW2 || defined FFTW3 || defined CUDA_FFT //native fft and ACML already done normalization
00290
00291 float scale = 1.0f / ((nx - offset) * ny * nz);
00292 dat->mult(scale);
00293 #endif //FFTW2 || FFTW3
00294 dat->set_fftodd(false);
00295 dat->set_fftpad(false);
00296 dat->set_complex(false);
00297 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
00298 dat->set_ri(false);
00299 dat->update();
00300
00301
00302 EXITFUNC;
00303 return dat;
00304 }
00305
00306
00307
00308
00309
00310 EMData *EMData::do_ift_inplace()
00311 {
00312 ENTERFUNC;
00313
00314 if (!is_complex()) {
00315 LOGERR("complex image expected. Input image is real image.");
00316 throw ImageFormatException("complex image expected. Input image is real image.");
00317 }
00318
00319 if (!is_ri()) {
00320 LOGWARN("run IFT on AP data, only RI should be used. ");
00321 }
00322 ap2ri();
00323
00324 int offset = is_fftodd() ? 1 : 2;
00325 float* data = get_data();
00326 EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
00327
00328 #if defined FFTW2 || defined FFTW3 || defined CUDA_FFT //native fft and ACML already done normalization
00329
00330 int nxo = nx - offset;
00331 float scale = 1.0f / (nxo * ny * nz);
00332 mult(scale);
00333 #endif //FFTW2 || FFTW3
00334
00335 #ifndef CUDA_FFT
00336 set_fftpad(true);
00337 #else
00338 set_size(nx - offset, ny, nz);
00339 #endif
00340 set_complex(false);
00341 if(ny==1 && nz==1) set_complex_x(false);
00342 set_ri(false);
00343 update();
00344
00345
00346 EXITFUNC;
00347 return this;
00348 }
00349 #undef rdata
00350
00351
00352 std::string EMData::render_ap24(int x0, int y0, int ixsize, int iysize,
00353 int bpl, float scale, int mingray, int maxgray,
00354 float render_min, float render_max,float gamma,int flags)
00355 {
00356 ENTERFUNC;
00357
00358 int asrgb;
00359 int hist=(flags&2)/2;
00360 int invy=(flags&4)?1:0;
00361
00362 if (!is_complex()) throw ImageDimensionException("complex only");
00363
00364 if (get_ndim() != 2) {
00365 throw ImageDimensionException("2D only");
00366 }
00367
00368 if (is_complex()) ri2ap();
00369
00370 if (render_max <= render_min) {
00371 render_max = render_min + 0.01f;
00372 }
00373
00374 if (gamma<=0) gamma=1.0;
00375
00376
00377
00378
00379
00380
00381
00382
00383 static int smg0=0,smg1=0;
00384 static float sgam=0;
00385 static unsigned char gammamap[4096];
00386 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00387 for (int i=0; i<4096; i++) {
00388 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00389 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00390 }
00391 }
00392 smg0=mingray;
00393 smg1=maxgray;
00394 sgam=gamma;
00395
00396 if (flags&8) asrgb=4;
00397 else if (flags&1) asrgb=3;
00398 else throw ImageDimensionException("must set flag 1 or 8");
00399
00400 std::string ret=std::string();
00401
00402 ret.assign(iysize*bpl+hist*1024,char(mingray));
00403 unsigned char *data=(unsigned char *)ret.data();
00404 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00405 if (hist) {
00406 for (int i=0; i<256; i++) histd[i]=0;
00407 }
00408
00409 float rm = render_min;
00410 float inv_scale = 1.0f / scale;
00411 int ysize = iysize;
00412 int xsize = ixsize;
00413
00414 int ymin = 0;
00415 if (iysize * inv_scale > ny) {
00416 ymin = (int) (iysize - ny / inv_scale);
00417 }
00418
00419 float gs = (maxgray - mingray) / (render_max - render_min);
00420 float gs2 = 4095.999f / (render_max - render_min);
00421
00422 if (render_max < render_min) {
00423 gs = 0;
00424 rm = FLT_MAX;
00425 }
00426
00427 int dsx = -1;
00428 int dsy = 0;
00429 int remx = 0;
00430 int remy = 0;
00431 const int scale_n = 100000;
00432
00433 int addi = 0;
00434 int addr = 0;
00435 if (inv_scale == floor(inv_scale)) {
00436 dsx = (int) inv_scale;
00437 dsy = (int) (inv_scale * nx);
00438 }
00439 else {
00440 addi = (int) floor(inv_scale);
00441 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00442 }
00443
00444 int xmin = 0;
00445 if (x0 < 0) {
00446 xmin = (int) (-x0 / inv_scale);
00447 xsize -= (int) floor(x0 / inv_scale);
00448 x0 = 0;
00449 }
00450
00451 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00452 xsize = (int) ((nx - x0) / inv_scale + xmin);
00453 }
00454 int ymax = ysize - 1;
00455 if (y0 < 0) {
00456 ymax = (int) (ysize + y0 / inv_scale - 1);
00457 ymin += (int) floor(y0 / inv_scale);
00458 y0 = 0;
00459 }
00460
00461 if (xmin < 0) xmin = 0;
00462 if (ymin < 0) ymin = 0;
00463 if (xsize > ixsize) xsize = ixsize;
00464 if (ymax > iysize) ymax = iysize;
00465
00466 int lmax = nx * ny - 1;
00467
00468 int mid=nx*ny/2;
00469 float* image_data = get_data();
00470 if (dsx != -1) {
00471 int l = y0 * nx;
00472 for (int j = ymax; j >= ymin; j--) {
00473 int ll = x0;
00474 for (int i = xmin; i < xsize; i++) {
00475 if (l + ll > lmax || ll >= nx - 2) break;
00476
00477 int k = 0;
00478 unsigned char p;
00479 int ph;
00480 if (ll >= nx / 2) {
00481 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00482 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00483 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00484 }
00485 else {
00486 k = nx * ny - (l + 2 * ll) - 2;
00487 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00488 }
00489 if (k>=mid) k-=mid;
00490 else k+=mid;
00491 float t = image_data[k];
00492 if (t <= rm) p = mingray;
00493 else if (t >= render_max) p = maxgray;
00494 else if (gamma!=1.0) {
00495 k=(int)(gs2 * (t-render_min));
00496 p = gammamap[k];
00497 }
00498 else {
00499 p = (unsigned char) (gs * (t - render_min));
00500 p += mingray;
00501 }
00502 if (ph<256) {
00503 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00504 data[i * asrgb + j * bpl+1] = p*ph/256;
00505 data[i * asrgb + j * bpl+2] = 0;
00506 }
00507 else if (ph<512) {
00508 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00509 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00510 data[i * asrgb + j * bpl] = 0;
00511 }
00512 else {
00513 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00514 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00515 data[i * asrgb + j * bpl+1] = 0;
00516 }
00517 if (hist) histd[p]++;
00518 ll += dsx;
00519 }
00520 l += dsy;
00521 }
00522 }
00523 else {
00524 remy = 10;
00525 int l = y0 * nx;
00526 for (int j = ymax; j >= ymin; j--) {
00527 int br = l;
00528 remx = 10;
00529 int ll = x0;
00530 for (int i = xmin; i < xsize - 1; i++) {
00531 if (l + ll > lmax || ll >= nx - 2) {
00532 break;
00533 }
00534 int k = 0;
00535 unsigned char p;
00536 int ph;
00537 if (ll >= nx / 2) {
00538 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00539 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00540 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00541 }
00542 else {
00543 k = nx * ny - (l + 2 * ll) - 2;
00544 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00545 }
00546 if (k>=mid) k-=mid;
00547 else k+=mid;
00548
00549 float t = image_data[k];
00550 if (t <= rm)
00551 p = mingray;
00552 else if (t >= render_max) {
00553 p = maxgray;
00554 }
00555 else if (gamma!=1.0) {
00556 k=(int)(gs2 * (t-render_min));
00557 p = gammamap[k];
00558 }
00559 else {
00560 p = (unsigned char) (gs * (t - render_min));
00561 p += mingray;
00562 }
00563 if (ph<256) {
00564 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00565 data[i * asrgb + j * bpl+1] = p*ph/256;
00566 data[i * asrgb + j * bpl+2] = 0;
00567 }
00568 else if (ph<512) {
00569 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00570 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00571 data[i * asrgb + j * bpl] = 0;
00572 }
00573 else {
00574 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00575 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00576 data[i * asrgb + j * bpl+1] = 0;
00577 }
00578 if (hist) histd[p]++;
00579 ll += addi;
00580 remx += addr;
00581 if (remx > scale_n) {
00582 remx -= scale_n;
00583 ll++;
00584 }
00585 }
00586 l = br + addi * nx;
00587 remy += addr;
00588 if (remy > scale_n) {
00589 remy -= scale_n;
00590 l += nx;
00591 }
00592 }
00593 }
00594
00595
00596 if (asrgb==4) {
00597 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00598 for (int i=xmin; i<xsize*4; i+=4) {
00599 data[i+j+3]=255;
00600 }
00601 }
00602 }
00603
00604 EXITFUNC;
00605
00606
00607 if (invy) {
00608 int x,y;
00609 char swp;
00610 for (y=0; y<iysize/2; y++) {
00611 for (x=0; x<ixsize; x++) {
00612 swp=ret[y*bpl+x];
00613 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00614 ret[(iysize-y-1)*bpl+x]=swp;
00615 }
00616 }
00617 }
00618
00619
00620 return ret;
00621 }
00622
00623
00624 void EMData::render_amp24( int x0, int y0, int ixsize, int iysize,
00625 int bpl, float scale, int mingray, int maxgray,
00626 float render_min, float render_max, void *ref,
00627 void cmap(void *, int coord, unsigned char *tri))
00628 {
00629 ENTERFUNC;
00630
00631 if (get_ndim() != 2) {
00632 throw ImageDimensionException("2D only");
00633 }
00634
00635 if (is_complex()) {
00636 ri2ap();
00637 }
00638
00639 if (render_max <= render_min) {
00640 render_max = render_min + 0.01f;
00641 }
00642
00643 std::string ret=std::string();
00644 ret.resize(iysize*bpl);
00645 unsigned char *data=(unsigned char *)ret.data();
00646
00647 float rm = render_min;
00648 float inv_scale = 1.0f / scale;
00649 int ysize = iysize;
00650 int xsize = ixsize;
00651 const int scale_n = 100000;
00652
00653 int ymin = 0;
00654 if ( iysize * inv_scale > ny) {
00655 ymin = (int) (iysize - ny / inv_scale);
00656 }
00657 float gs = (maxgray - mingray) / (render_max - render_min);
00658 if (render_max < render_min) {
00659 gs = 0;
00660 rm = FLT_MAX;
00661 }
00662 int dsx = -1;
00663 int dsy = 0;
00664 if (inv_scale == floor(inv_scale)) {
00665 dsx = (int) inv_scale;
00666 dsy = (int) (inv_scale * nx);
00667 }
00668 int addi = 0;
00669 int addr = 0;
00670
00671 if (dsx == -1) {
00672 addi = (int) floor(inv_scale);
00673 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00674 }
00675
00676 int remx = 0;
00677 int remy = 0;
00678 int xmin = 0;
00679 if (x0 < 0) {
00680 xmin = (int) (-x0 / inv_scale);
00681 xsize -= (int) floor(x0 / inv_scale);
00682 x0 = 0;
00683 }
00684
00685 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00686 xsize = (int) ((nx - x0) / inv_scale + xmin);
00687 }
00688 int ymax = ysize - 1;
00689 if (y0 < 0) {
00690 ymax = (int) (ysize + y0 / inv_scale - 1);
00691 ymin += (int) floor(y0 / inv_scale);
00692 y0 = 0;
00693 }
00694
00695
00696 if (xmin < 0) {
00697 xmin = 0;
00698 }
00699
00700 if (ymin < 0) {
00701 ymin = 0;
00702 }
00703 if (xsize > ixsize) {
00704 xsize = ixsize;
00705 }
00706 if (ymax > iysize) {
00707 ymax = iysize;
00708 }
00709
00710 int lmax = nx * ny - 1;
00711 unsigned char tri[3];
00712 float* image_data = get_data();
00713 if (is_complex()) {
00714 if (dsx != -1) {
00715 int l = y0 * nx;
00716 for (int j = ymax; j >= ymin; j--) {
00717 int ll = x0;
00718 for (int i = xmin; i < xsize; i++, ll += dsx) {
00719 if (l + ll > lmax || ll >= nx - 2) {
00720 break;
00721 }
00722 int kk = 0;
00723 if (ll >= nx / 2) {
00724 if (l >= (ny - inv_scale) * nx) {
00725 kk = 2 * (ll - nx / 2) + 2;
00726 }
00727 else {
00728 kk = 2 * (ll - nx / 2) + l + 2 + nx;
00729 }
00730 }
00731 else {
00732 kk = nx * ny - (l + 2 * ll) - 2;
00733 }
00734 int k = 0;
00735 float t = image_data[kk];
00736 if (t <= rm) {
00737 k = mingray;
00738 }
00739 else if (t >= render_max) {
00740 k = maxgray;
00741 }
00742 else {
00743 k = (int) (gs * (t - render_min));
00744 k += mingray;
00745 }
00746 tri[0] = static_cast < unsigned char >(k);
00747 cmap(ref, kk, tri);
00748 data[i * 3 + j * bpl] = tri[0];
00749 data[i * 3 + 1 + j * bpl] = tri[1];
00750 data[i * 3 + 2 + j * bpl] = tri[2];
00751 }
00752 l += dsy;
00753 }
00754 }
00755 else {
00756 remy = 10;
00757 for (int j = ymax, l = y0 * nx; j >= ymin; j--) {
00758 int br = l;
00759 remx = 10;
00760 for (int i = xmin, ll = x0; i < xsize - 1; i++) {
00761 if (l + ll > lmax || ll >= nx - 2) {
00762 break;
00763 }
00764 int kk = 0;
00765 if (ll >= nx / 2) {
00766 if (l >= (ny * nx - nx)) {
00767 kk = 2 * (ll - nx / 2) + 2;
00768 }
00769 else {
00770 kk = 2 * (ll - nx / 2) + l + 2 + nx;
00771 }
00772 }
00773 else {
00774 kk = nx * ny - (l + 2 * ll) - 2;
00775 }
00776 int k = 0;
00777 float t = image_data[kk];
00778 if (t <= rm) {
00779 k = mingray;
00780 }
00781 else if (t >= render_max) {
00782 k = maxgray;
00783 }
00784 else {
00785 k = (int) (gs * (t - render_min));
00786 k += mingray;
00787 }
00788 tri[0] = static_cast < unsigned char >(k);
00789 cmap(ref, kk, tri);
00790 data[i * 3 + j * bpl] = tri[0];
00791 data[i * 3 + 1 + j * bpl] = tri[1];
00792 data[i * 3 + 2 + j * bpl] = tri[2];
00793 ll += addi;
00794 remx += addr;
00795 if (remx > scale_n) {
00796 remx -= scale_n;
00797 ll++;
00798 }
00799 }
00800 l = br + addi * nx;
00801 remy += addr;
00802 if (remy > scale_n) {
00803 remy -= scale_n;
00804 l += nx;
00805 }
00806 }
00807 }
00808 }
00809 else {
00810 if (dsx != -1) {
00811 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00812 int br = l;
00813 for (int i = xmin; i < xsize; i++, l += dsx) {
00814 if (l > lmax) {
00815 break;
00816 }
00817 float t = image_data[l];
00818 int k = 0;
00819 if (t <= rm) {
00820 k = mingray;
00821 }
00822 else if (t >= render_max) {
00823 k = maxgray;
00824 }
00825 else {
00826 k = (int) (gs * (t - render_min));
00827 k += mingray;
00828 }
00829 tri[0] = static_cast < unsigned char >(k);
00830 cmap(ref, l, tri);
00831 data[i * 3 + j * bpl] = tri[0];
00832 data[i * 3 + 1 + j * bpl] = tri[1];
00833 data[i * 3 + 2 + j * bpl] = tri[2];
00834 }
00835 l = br + dsy;
00836 }
00837 }
00838 else {
00839 remy = 10;
00840 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00841 int br = l;
00842 remx = 10;
00843 for (int i = xmin; i < xsize; i++) {
00844 if (l > lmax) {
00845 break;
00846 }
00847 float t = image_data[l];
00848 int k = 0;
00849 if (t <= rm) {
00850 k = mingray;
00851 }
00852 else if (t >= render_max) {
00853 k = maxgray;
00854 }
00855 else {
00856 k = (int) (gs * (t - render_min));
00857 k += mingray;
00858 }
00859 tri[0] = static_cast < unsigned char >(k);
00860 cmap(ref, l, tri);
00861 data[i * 3 + j * bpl] = tri[0];
00862 data[i * 3 + 1 + j * bpl] = tri[1];
00863 data[i * 3 + 2 + j * bpl] = tri[2];
00864 l += addi;
00865 remx += addr;
00866 if (remx > scale_n) {
00867 remx -= scale_n;
00868 l++;
00869 }
00870 }
00871 l = br + addi * nx;
00872 remy += addr;
00873 if (remy > scale_n) {
00874 remy -= scale_n;
00875 l += nx;
00876 }
00877 }
00878 }
00879 }
00880
00881 EXITFUNC;
00882 }
00883
00884 void EMData::ap2ri()
00885 {
00886 ENTERFUNC;
00887
00888 if (!is_complex() || is_ri()) {
00889 return;
00890 }
00891
00892 #ifdef EMAN2_USING_CUDA
00893 if (gpu_operation_preferred()) {
00894 EMDataForCuda tmp = get_data_struct_for_cuda();
00895 emdata_ap2ri(&tmp);
00896 set_ri(true);
00897 gpu_update();
00898 EXITFUNC;
00899 return;
00900 }
00901 #endif
00902
00903 Util::ap2ri(get_data(), nx * ny * nz);
00904 set_ri(true);
00905 update();
00906 EXITFUNC;
00907 }
00908
00909 void EMData::ri2inten()
00910 {
00911 ENTERFUNC;
00912
00913 if (!is_complex()) return;
00914 if (!is_ri()) ap2ri();
00915
00916 #ifdef EMAN2_USING_CUDA
00917 if (gpu_operation_preferred()) {
00918 EMDataForCuda tmp = get_data_struct_for_cuda();
00919 emdata_ri2inten(&tmp);
00920 set_attr("is_intensity", int(1));
00921 gpu_update();
00922 EXITFUNC;
00923 return;
00924 }
00925 #endif
00926
00927 float * data = get_data();
00928 size_t size = nx * ny * nz;
00929 for (size_t i = 0; i < size; i += 2) {
00930 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
00931 data[i+1]=0;
00932 }
00933
00934 set_attr("is_intensity", int(1));
00935 update();
00936 EXITFUNC;
00937 }
00938
00939
00940 void EMData::ri2ap()
00941 {
00942 ENTERFUNC;
00943
00944 if (!is_complex() || !is_ri()) {
00945 return;
00946 }
00947 #ifdef EMAN2_USING_CUDA
00948 if (gpu_operation_preferred()) {
00949 EMDataForCuda tmp = get_data_struct_for_cuda();
00950 emdata_ri2ap(&tmp);
00951 set_ri(false);
00952 gpu_update();
00953 EXITFUNC;
00954 return;
00955 }
00956 #endif
00957
00958 float * data = get_data();
00959
00960 size_t size = nx * ny * nz;
00961 for (size_t i = 0; i < size; i += 2) {
00962 #ifdef _WIN32
00963 float f = (float)_hypot(data[i], data[i + 1]);
00964 #else
00965 float f = (float)hypot(data[i], data[i + 1]);
00966 #endif
00967 if (data[i] == 0 && data[i + 1] == 0) {
00968 data[i + 1] = 0;
00969 }
00970 else {
00971 data[i + 1] = atan2(data[i + 1], data[i]);
00972 }
00973 data[i] = f;
00974 }
00975
00976 set_ri(false);
00977 update();
00978 EXITFUNC;
00979 }
00980
00981
00982 float calc_bessel(const int n, const float& x) {
00983 gsl_sf_result result;
00984
00985 gsl_sf_bessel_Jn_e(n,(double)x, &result);
00986 return (float)result.val;
00987 }
00988
00989 EMData* EMData::bispecRotTransInvN(int N, int NK)
00990 {
00991
00992 int EndP = this -> get_xsize();
00993 int Mid = (int) ((1+EndP)/2);
00994 int End = 2*Mid-1;
00995
00996 int CountxyMax = End*End;
00997
00998 int *SortfkInds = new int[CountxyMax];
00999 int *kVecX = new int[CountxyMax];
01000 int *kVecY = new int[CountxyMax];
01001 float *fkVecR = new float[CountxyMax];
01002 float *fkVecI = new float[CountxyMax];
01003 float *absD1fkVec = new float[CountxyMax];
01004 float *absD1fkVecSorted = new float[CountxyMax];
01005
01006 float *jxjyatan2 = new float[End*End];
01007
01008 EMData * ThisCopy = new EMData(End,End);
01009
01010 for (int jx=0; jx <End ; jx++) {
01011 for (int jy=0; jy <End ; jy++) {
01012 float ValNow = this -> get_value_at(jx,jy);
01013 ThisCopy -> set_value_at(jx,jy,ValNow);
01014
01015 }}
01016
01017
01018 EMData* fk = ThisCopy -> do_fft();
01019 fk ->process_inplace("xform.fourierorigin.tocenter");
01020
01021
01022 EMData* fkRCopy = new EMData(End,End);
01023 EMData* fkICopy = new EMData(End,End);
01024 EMData* fkCopy = new EMData(End,End);
01025
01026
01027 for (int jCount= 0; jCount<End*End; jCount++) {
01028
01029 int jx = jCount%End ;
01030 int jy = (jCount-jx)/End ;
01031 jxjyatan2[jCount] = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
01032 }
01033
01034
01035 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01036
01037 int kx = kEx/2;
01038 int kIx = kx+Mid-1;
01039 int kCx = -kx ;
01040 int kCIx = kCx+ Mid-1 ;
01041 for (int kEy= 0 ; kEy<End; kEy++) {
01042 int kIy = kEy ;
01043 int ky = kEy+1-Mid;
01044 float realVal = fk -> get_value_at(kEx ,kEy) ;
01045 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01046 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01047 float fkAng = atan2(imagVal,realVal);
01048
01049 float NewRealVal ;
01050 float NewImagVal ;
01051 float AngMatlab ;
01052
01053 if (kIx==Mid-1) {
01054
01055
01056 }
01057
01058 if (kIx>Mid-1){
01059
01060 }
01061
01062 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01063 NewRealVal = absVal*cos(AngMatlab);
01064 NewImagVal = absVal*sin(AngMatlab);
01065
01066
01067 fkVecR[kIy+kIx *End] = NewRealVal ;
01068 fkVecR[kIy+kCIx*End] = NewRealVal ;
01069 fkVecI[kIy+kIx *End] = NewImagVal ;
01070 fkVecI[kIy+kCIx*End] = -NewImagVal ;
01071 absD1fkVec[kIy + kIx *End] = absVal;
01072 absD1fkVec[kIy + kCIx *End] = absVal;
01073 kVecX[kIy+kIx *End] = kx ;
01074 kVecX[kIy+kCIx *End] = kCx ;
01075 kVecY[kIy+kIx *End] = ky ;
01076 kVecY[kIy+kCIx *End] = ky ;
01077
01078
01079
01080
01081
01082 fkCopy -> set_value_at(kIx ,kIy, absVal);
01083 fkCopy -> set_value_at(kCIx,kIy, absVal);
01084 fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
01085 fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
01086 fkICopy -> set_value_at(kIx, kIy, NewImagVal);
01087 fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
01088
01089 }
01090 }
01091 system("rm -f fkCopy.???");
01092 system("rm -f fk?Copy.???");
01093 fkCopy -> write_image("fkCopy.img");
01094 fkRCopy -> write_image("fkRCopy.img");
01095 fkICopy -> write_image("fkICopy.img");
01096
01097 cout << "Starting the sort "<< endl;
01098
01099 vector< pair<float, int> > absInds;
01100 for(int i = 0; i < CountxyMax; ++i ) {
01101 pair<float,int> p;
01102 p = make_pair(absD1fkVec[i],i);
01103 absInds.push_back( p);
01104 }
01105
01106 std::sort(absInds.begin(),absInds.end());
01107
01108 for(int i = 0; i < CountxyMax; ++i ) {
01109 pair<float,int> p ;
01110 p = absInds[i] ;
01111 absD1fkVecSorted[CountxyMax-1-i] = p.first ;
01112 SortfkInds[CountxyMax-1-i] = p.second ;
01113 }
01114
01115 cout << "Ending the sort "<< endl;
01116
01117
01118
01119
01120 for(int i = 0; i < CountxyMax; ++i ) {
01121 int Si = SortfkInds[i];
01122 int kIx = (int) Si/End; kIx = (int) i/End;
01123 int kIy = Si - kIx*End; kIy = i - kIx*End;
01124 int iC = (End-1-kIx)*End + (End-1-kIy);
01125
01126
01127 }
01128 cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0] << endl;
01129
01130
01131
01132 for(int i = 0; i < NK; ++i ) {
01133 int Si= SortfkInds[i];
01134 int kIx = (int) Si/End;
01135 int kIy = Si - kIx*End;
01136
01137 }
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152 cout << "NK = " << NK << endl;
01153 float frR= 3.0/4.0;
01154 int LradRange= (int) (floor(Mid/frR)) ;
01155
01156 float *radRange = new float[LradRange];
01157 radRange[0]=0;
01158 for (int irad=1; irad < LradRange; irad++){
01159 radRange[irad] = radRange[irad-1] + frR; }
01160
01161
01162
01163
01164 cout << "Starting the calculation of invariants for N= " << N << endl;
01165
01166
01167
01168 EMData* RotTransInv = new EMData();
01169 RotTransInv -> set_size(LradRange,LradRange);
01170
01171
01172
01173
01174
01175
01176
01177 for (int jr1=0; jr1 < LradRange ; jr1++ ) {
01178 float r1= radRange[jr1];
01179
01180 for (int jr2=0; jr2<LradRange; jr2++ ) {
01181 float r2= radRange[jr2];
01182 float RotTransInvTemp=0;
01183 for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
01184 int Countkxy =SortfkInds[jCountkxy] ;
01185 int kx = kVecX[Countkxy] ;
01186 int ky = kVecY[Countkxy] ;
01187 float k2 = (float)(kx*kx+ky*ky);
01188 if (k2==0) { continue;}
01189 float phiK =0;
01190 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; phiK= atan2((float)ky,(float)kx);
01191
01192 float fkR = fkVecR[Countkxy] ;
01193 float fkI = fkVecI[Countkxy] ;
01194
01195
01196 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
01197 int Countqxy =SortfkInds[jCountqxy] ;
01198 int qx = kVecX[Countqxy] ;
01199 int qy = kVecY[Countqxy] ;
01200 int q2 = qx*qx+qy*qy;
01201 if (q2==0) {continue;}
01202 float phiQ =0;
01203 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; phiQ=atan2((float)qy,(float)qx);
01204 float fqR = fkVecR[Countqxy] ;
01205 float fqI = fkVecI[Countqxy] ;
01206 int kCx = (-kx-qx);
01207 int kCy = (-ky-qy);
01208 int kCIx = ((kCx+Mid+2*End)%End);
01209 int kCIy = ((kCy+Mid+2*End)%End);
01210 kCx = kCIx-Mid;
01211 kCy = kCIy-Mid;
01212 int CountCxy = kCIx*End+kCIy;
01213 float fCR = fkVecR[CountCxy];
01214 float fCI = fkVecI[CountCxy];
01215 if (jr1+jr2==-1) {
01216 printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
01217 printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f, kCx=%d, kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01218 }
01219 for (int p=0; p<NK; p++){
01220
01221 if (SortfkInds[p]==CountCxy){
01222 float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
01223 float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
01224
01225
01226 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI))
01227 * cos(N*(phiK-phiQ+M_PI));
01228 bispectemp -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
01229 * sin(N*(phiK-phiQ+M_PI));
01230 float bess1 = calc_bessel(N, Arg1 );
01231 float bess2 = calc_bessel(N, Arg2 );
01232
01233
01234
01235 RotTransInvTemp = RotTransInvTemp + bispectemp * bess1*bess2 ;
01236
01237 }
01238 }
01239 }
01240 }
01241 RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp) ;
01242
01243 }
01244 }
01245
01246
01247 return RotTransInv ;
01248
01249
01250 }
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281 EMData* EMData::bispecRotTransInvDirect(int type)
01282 {
01283
01284 int EndP = this -> get_xsize();
01285 int Mid = (int) ((1+EndP)/2);
01286 int End = 2*Mid-1;
01287
01288 int CountxyMax = End*End;
01289
01290 int *SortfkInds = new int[CountxyMax];
01291 int *kVecX = new int[CountxyMax];
01292 int *kVecY = new int[CountxyMax];
01293 float *fkVecR = new float[CountxyMax];
01294 float *fkVecI = new float[CountxyMax];
01295 float *absD1fkVec = new float[CountxyMax];
01296 float *absD1fkVecSorted = new float[CountxyMax];
01297
01298
01299 float *jxjyatan2 = new float[End*End];
01300
01301
01302 EMData * ThisCopy = new EMData(End,End);
01303
01304 for (int jx=0; jx <End ; jx++) {
01305 for (int jy=0; jy <End ; jy++) {
01306 float ValNow = this -> get_value_at(jx,jy);
01307 ThisCopy -> set_value_at(jx,jy,ValNow);
01308 jxjyatan2[jy*End + jx] = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
01309
01310 }}
01311
01312
01313 EMData* fk = ThisCopy -> do_fft();
01314 fk ->process_inplace("xform.fourierorigin.tocenter");
01315
01316
01317
01318 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01319
01320 int kx = kEx/2;
01321 int kIx = kx+Mid-1;
01322 int kCx = -kx ;
01323 int kCIx = kCx+ Mid-1 ;
01324 for (int kEy= 0 ; kEy<End; kEy++) {
01325 int kIy = kEy ;
01326 int ky = kEy+1-Mid;
01327 float realVal = fk -> get_value_at(kEx ,kEy) ;
01328 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01329 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01330 float fkAng = atan2(imagVal,realVal);
01331
01332 float NewRealVal ;
01333 float NewImagVal ;
01334 float AngMatlab ;
01335
01336 if (kIx==Mid-1) {
01337
01338 }
01339
01340 if (kIx>Mid-1){
01341
01342 }
01343
01344 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01345 NewRealVal = absVal*cos(AngMatlab);
01346 NewImagVal = absVal*sin(AngMatlab);
01347
01348
01349 fkVecR[ kIy +kIx *End] = NewRealVal ;
01350 fkVecR[(End-1-kIy)+kCIx*End] = NewRealVal ;
01351 fkVecI[ kIy +kIx *End] = NewImagVal ;
01352 fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
01353 absD1fkVec[(End-1-kIy) + kIx *End] = absVal;
01354 absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
01355 kVecX[kIy+kIx *End] = kx ;
01356 kVecX[kIy+kCIx *End] = kCx ;
01357 kVecY[kIy+kIx *End] = ky ;
01358 kVecY[kIy+kCIx *End] = ky ;
01359
01360
01361
01362
01363
01364
01365
01366
01367 }
01368 }
01369
01370
01371 for (int TotalInd = 0 ; TotalInd < CountxyMax ; TotalInd++){
01372 int kx = kVecX[TotalInd];
01373 int kIx = kx+Mid-1;
01374 int ky = kVecY[TotalInd];
01375 int kIy = ky+Mid-1;
01376 float fkR = fkVecR[kIy+kIx *End] ;
01377 float fkI = fkVecI[kIy+kIx *End] ;
01378 }
01379
01380 float frR= 3.0/4.0;
01381 frR= 1;
01382 int LradRange= (int) (1+floor(Mid/frR -.1)) ;
01383
01384 float *radRange = new float[LradRange];
01385 for (int irad=0; irad < LradRange; irad++){
01386 radRange[irad] = frR*irad;
01387
01388 }
01389
01390 cout << "Starting the calculation of invariants" << endl;
01391
01392
01393 if (type==0) {
01394 int LthetaRange = 59;
01395 float ftR = (2.0f*M_PI/LthetaRange );
01396 float *thetaRange = new float[LthetaRange];
01397
01398 for (int ith=0; ith < LthetaRange; ith++){
01399 thetaRange[ith] = ftR*ith; }
01400
01401
01402
01403
01404
01405 int TotalVol = LradRange*LradRange*LthetaRange;
01406
01407 float *RotTransInv = new float[TotalVol];
01408 float *WeightInv = new float[TotalVol];
01409
01410 for (int jW=0; jW<TotalVol; jW++) {
01411 RotTransInv[jW] = 0;
01412 WeightInv[jW] = 0;
01413 }
01414
01415
01416 for (int jW=0; jW<TotalVol; jW++) {
01417 RotTransInv[jW] = 0;
01418 WeightInv[jW] = 0;
01419 }
01420
01421
01422
01423 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01424 int kx = kVecX[Countkxy] ;
01425 int ky = kVecY[Countkxy] ;
01426 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01427 float phiK =0;
01428 if (k2>0) phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1];
01429 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01430 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01431
01432
01433 if ((k2==0)|| (k2>Mid) ) { continue;}
01434
01435 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01436 int qx = kVecX[Countqxy] ;
01437 int qy = kVecY[Countqxy] ;
01438 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01439 if ((q2==0)|| (q2>Mid) ) {continue;}
01440 float phiQ =0;
01441 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1];
01442 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01443 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01444 int kCx = (-kx-qx);
01445 int kCy = (-ky-qy);
01446 int kCIx = ((kCx+Mid+2*End)%End);
01447 int kCIy = ((kCy+Mid+2*End)%End);
01448 kCx = ((kCIx+End-1)%End)+1-Mid;
01449 kCy = ((kCIy+End-1)%End)+1-Mid ;
01450
01451 float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
01452 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01453 float fCR = fkVecR[CountCxy];
01454 float fCI = fkVecI[CountCxy];
01455
01456
01457
01458
01459 float phiC = jxjyatan2[ (kCy+Mid-1)*End + kCx+Mid-1];
01460 float phiQK = (4*M_PI+phiQ-phiK);
01461 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
01462
01463
01464
01465 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01466
01467 if ((q2<k2) ) continue;
01468
01469
01470
01471
01472
01473
01474 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01475 int k2IndHi = k2IndLo;
01476 float k2Lo= radRange[k2IndLo];
01477 if (k2IndLo+1< LradRange) {
01478 k2IndHi = k2IndLo+1;
01479 }
01480 float k2Hi= radRange[k2IndHi];
01481
01482 float kCof =k2-k2Lo;
01483
01484 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01485 int q2IndHi=q2IndLo;
01486 float q2Lo= radRange[q2IndLo];
01487 if (q2IndLo+1 < LradRange) {
01488 q2IndHi = q2IndLo+1 ;
01489 }
01490 float qCof = q2-q2Lo;
01491
01492 if ((qCof<0) || (qCof >1) ) {
01493 cout<< "Weird! qCof="<< qCof << " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
01494 int x ;
01495 cin >> x ;
01496 }
01497
01498 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
01499 int thetaIndHi = thetaIndLo;
01500
01501 float thetaLo = thetaRange[thetaIndLo];
01502 float thetaHi = thetaLo;
01503 float thetaCof = 0;
01504
01505 if (thetaIndLo+1< LthetaRange) {
01506 thetaIndHi = thetaIndLo +1;
01507 }else{
01508 thetaIndHi=0;
01509 }
01510
01511 thetaHi = thetaRange[thetaIndHi];
01512
01513 if (thetaHi==thetaLo) {
01514 thetaCof =0 ;
01515 } else {
01516 thetaCof = (phiQK-thetaLo)/(thetaHi-thetaLo);
01517 }
01518
01519 if ((thetaCof>2*M_PI) ) {
01520 cout<< "Weird! thetaCof="<< thetaCof <<endl ;
01521 thetaCof=0;
01522 }
01523
01524
01525
01526
01527
01528 for (int jk =1; jk<=2; jk++){
01529 for (int jq =1; jq<=2; jq++){
01530 for (int jtheta =1; jtheta<=2; jtheta++){
01531
01532 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
01533 * (thetaCof+(1-2*thetaCof)*(jtheta==1));
01534
01535
01536 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01537 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01538 int thetaInd = thetaIndLo*(jtheta==1) + thetaIndHi*(jtheta ==2);
01539 int TotalInd = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
01540
01541
01542
01543
01544
01545
01546 RotTransInv[TotalInd] += Weight*bispectemp;
01547 WeightInv[TotalInd] += Weight;
01548
01549 }}}
01550 }
01551 }
01552
01553 cout << "Finished Main Section " << endl;
01554
01555
01556
01557 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
01558 EMData *RotTransInvF = new EMData(LradRange,LradRange,LthetaRange);
01559 EMData *WeightImage = new EMData(LradRange,LradRange,LthetaRange);
01560
01561
01562
01563
01564
01565
01566
01567 for (int jtheta =0; jtheta < LthetaRange; jtheta++){
01568 for (int jq =0; jq<LradRange; jq++){
01569 for (int jk =0; jk<LradRange ; jk++){
01570
01571 int TotalInd = jtheta*LradRange*LradRange+jq*LradRange+jk;
01572 float Weight = WeightInv[TotalInd];
01573 WeightImage -> set_value_at(jk,jq,jtheta,Weight);
01574 RotTransInvF -> set_value_at(jk,jq,jtheta,0);
01575 if (Weight <= 0) continue;
01576 RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);
01577 int newjtheta = (LthetaRange- jtheta)%LthetaRange;
01578 RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );
01579 }
01580 }
01581 }
01582
01583 cout << " Almost Done " << endl;
01584 system("rm -f WeightImage.???");
01585 WeightImage -> write_image("WeightImage.img");
01586
01587 return RotTransInvF ;
01588 }
01589
01590 if (type==1) {
01591 int TotalVol = LradRange*LradRange;
01592
01593 float *RotTransInv = new float[TotalVol];
01594 float *WeightInv = new float[TotalVol];
01595
01596 for (int jW=0; jW<TotalVol; jW++) {
01597 RotTransInv[jW] = 0;
01598 WeightInv[jW] = 0;
01599 }
01600
01601
01602 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01603 int kx = kVecX[Countkxy] ;
01604 int ky = kVecY[Countkxy] ;
01605 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01606 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01607 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01608
01609
01610 if ((k2==0)|| (k2>Mid) ) { continue;}
01611
01612 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01613
01614
01615 int qx = kVecX[Countqxy] ;
01616 int qy = kVecY[Countqxy] ;
01617 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01618 if ((q2==0)|| (q2>Mid) ) {continue;}
01619 if ((q2<k2) ) continue;
01620
01621 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01622 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01623
01624 int kCx = (-kx-qx);
01625 int kCy = (-ky-qy);
01626 int kCIx = ((kCx+Mid+2*End)%End);
01627 int kCIy = ((kCy+Mid+2*End)%End);
01628 kCx = ((kCIx+End-1)%End)+1-Mid;
01629 kCy = ((kCIy+End-1)%End)+1-Mid ;
01630
01631 float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
01632 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01633 float fCR = fkVecR[CountCxy];
01634 float fCI = fkVecI[CountCxy];
01635
01636
01637 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01638
01639
01640 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01641 int k2IndHi = k2IndLo;
01642 float k2Lo= radRange[k2IndLo];
01643 if (k2IndLo+1< LradRange) {
01644 k2IndHi = k2IndLo+1;
01645 }
01646 float k2Hi= radRange[k2IndHi];
01647
01648 float kCof =k2-k2Lo;
01649
01650 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01651 int q2IndHi=q2IndLo;
01652 float q2Lo= radRange[q2IndLo];
01653 if (q2IndLo+1 < LradRange) {
01654 q2IndHi = q2IndLo+1 ;
01655 }
01656 float qCof = q2-q2Lo;
01657
01658
01659 for (int jk =1; jk<=2; jk++){
01660 for (int jq =1; jq<=2; jq++){
01661
01662 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
01663
01664 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01665 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01666 int TotalInd = q2Ind*LradRange+k2Ind;
01667 RotTransInv[TotalInd] += Weight*bispectemp;
01668 WeightInv[TotalInd] += Weight;
01669
01670 }}
01671 }
01672 }
01673
01674
01675
01676
01677
01678 EMData *RotTransInvF = new EMData(LradRange,LradRange);
01679 EMData *WeightImage = new EMData(LradRange,LradRange);
01680
01681 for (int jk =0; jk<LradRange ; jk++){
01682 for (int jq =jk; jq<LradRange; jq++){
01683 int TotalInd = jq*LradRange+jk;
01684 int TotalIndBar = jq*LradRange+jk;
01685 float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
01686 if (Weight <=0) continue;
01687 WeightImage -> set_value_at(jk,jq,Weight);
01688 WeightImage -> set_value_at(jq,jk,Weight);
01689 #ifdef _WIN32
01690 float ValNow = pow( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight, 1.0f/3.0f ) ;
01691 #else
01692 float ValNow = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight ) ;
01693 #endif //_WIN32
01694 RotTransInvF -> set_value_at(jk,jq,ValNow);
01695 RotTransInvF -> set_value_at(jq,jk,ValNow );
01696 }}
01697
01698 system("rm -f WeightImage.???");
01699 WeightImage -> write_image("WeightImage.img");
01700
01701 return RotTransInvF ;
01702 }
01703 }
01704
01705
01706 void EMData::insert_clip(const EMData * const block, const IntPoint &origin) {
01707 int nx1 = block->get_xsize();
01708 int ny1 = block->get_ysize();
01709 int nz1 = block->get_zsize();
01710
01711 Region area(origin[0], origin[1], origin[2],nx1, ny1, nz1);
01712
01713 int x0 = (int) area.origin[0];
01714 x0 = x0 < 0 ? 0 : x0;
01715
01716 int y0 = (int) area.origin[1];
01717 y0 = y0 < 0 ? 0 : y0;
01718
01719 int z0 = (int) area.origin[2];
01720 z0 = z0 < 0 ? 0 : z0;
01721
01722 int x1 = (int) (area.origin[0] + area.size[0]);
01723 x1 = x1 > nx ? nx : x1;
01724
01725 int y1 = (int) (area.origin[1] + area.size[1]);
01726 y1 = y1 > ny ? ny : y1;
01727
01728 int z1 = (int) (area.origin[2] + area.size[2]);
01729 z1 = z1 > nz ? nz : z1;
01730 if (z1 <= 0) {
01731 z1 = 1;
01732 }
01733
01734 int xd0 = (int) (area.origin[0] < 0 ? -area.origin[0] : 0);
01735 int yd0 = (int) (area.origin[1] < 0 ? -area.origin[1] : 0);
01736 int zd0 = (int) (area.origin[2] < 0 ? -area.origin[2] : 0);
01737
01738 if (x1 < x0 || y1 < y0 || z1 < z0) return;
01739
01740 size_t clipped_row_size = (x1-x0) * sizeof(float);
01741 int src_secsize = nx1 * ny1;
01742 int dst_secsize = nx * ny;
01743
01744 #ifdef EMAN2_USING_CUDA
01745 if (gpu_operation_preferred()) {
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764 get_cuda_data();
01765 CudaDataLock lock(this);
01766 block->get_cuda_data();
01767 CudaDataLock lock2(block);
01768 for (int i = 0; i < (z1-z0); i++) {
01769 float* ldst = get_cuda_data() + (z0+i) * dst_secsize + y0 * nx + x0;
01770 float* lsrc = block->get_cuda_data() + (zd0+i) * src_secsize + yd0 * nx1 + xd0;
01771
01772
01773
01774
01775 cudaError_t error = cudaMemcpy2D( ldst, (nx)*sizeof(float), lsrc,nx1*sizeof(float), clipped_row_size, (y1-y0), cudaMemcpyDeviceToDevice );
01776 if ( error != cudaSuccess) {
01777 string e = cudaGetErrorString(error);
01778 throw UnexpectedBehaviorException("cudaMemcpy2D failed in get_clip with error: " + e);
01779 }
01780 }
01781 gpu_update();
01782 EXITFUNC;
01783 return;
01784 }
01785 #endif
01786 float *src = block->get_data() + zd0 * src_secsize + yd0 * nx1 + xd0;
01787 float *dst = get_data() + z0 * dst_secsize + y0 * nx + x0;
01788
01789
01790
01791
01792
01793 int src_gap = src_secsize - (y1-y0) * nx1;
01794 int dst_gap = dst_secsize - (y1-y0) * nx;
01795 for (int i = z0; i < z1; i++) {
01796 for (int j = y0; j < y1; j++) {
01797 EMUtil::em_memcpy(dst, src, clipped_row_size);
01798 src += nx1;
01799 dst += nx;
01800 }
01801 src += src_gap;
01802 dst += dst_gap;
01803 }
01804
01805 update();
01806 EXITFUNC;
01807 }
01808
01809
01810 void EMData::insert_scaled_sum(EMData *block, const FloatPoint ¢er,
01811 float scale, float)
01812 {
01813 ENTERFUNC;
01814 float * data = get_data();
01815 if (get_ndim()==3) {
01816
01817 int xs=(int)floor(block->get_xsize()*scale/2.0);
01818 int ys=(int)floor(block->get_ysize()*scale/2.0);
01819 int zs=(int)floor(block->get_zsize()*scale/2.0);
01820 int x0=(int)center[0]-xs;
01821 int x1=(int)center[0]+xs;
01822 int y0=(int)center[1]-ys;
01823 int y1=(int)center[1]+ys;
01824 int z0=(int)center[2]-zs;
01825 int z1=(int)center[2]+zs;
01826
01827 if (x1<0||y1<0||z1<0||x0>get_xsize()||y0>get_ysize()||z0>get_zsize()) return;
01828
01829
01830 if (x0<0) x0=0;
01831 if (y0<0) y0=0;
01832 if (z0<0) z0=0;
01833 if (x1>=get_xsize()) x1=get_xsize()-1;
01834 if (y1>=get_ysize()) y1=get_ysize()-1;
01835 if (z1>=get_zsize()) z1=get_zsize()-1;
01836
01837 float bx=block->get_xsize()/2.0f;
01838 float by=block->get_ysize()/2.0f;
01839 float bz=block->get_zsize()/2.0f;
01840
01841 size_t idx;
01842 for (int x=x0; x<=x1; x++) {
01843 for (int y=y0; y<=y1; y++) {
01844 for (int z=z0; z<=z1; z++) {
01845 idx = x + y * nx + z * nx * ny;
01846 data[idx] +=
01847 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by,(z-center[2])/scale+bz);
01848 }
01849 }
01850 }
01851 update();
01852 }
01853 else if (get_ndim()==2) {
01854
01855 int xs=(int)floor(block->get_xsize()*scale/2.0);
01856 int ys=(int)floor(block->get_ysize()*scale/2.0);
01857 int x0=(int)center[0]-xs;
01858 int x1=(int)center[0]+xs;
01859 int y0=(int)center[1]-ys;
01860 int y1=(int)center[1]+ys;
01861
01862 if (x1<0||y1<0||x0>get_xsize()||y0>get_ysize()) return;
01863
01864
01865 if (x0<0) x0=0;
01866 if (y0<0) y0=0;
01867 if (x1>=get_xsize()) x1=get_xsize()-1;
01868 if (y1>=get_ysize()) y1=get_ysize()-1;
01869
01870 float bx=block->get_xsize()/2.0f;
01871 float by=block->get_ysize()/2.0f;
01872
01873 for (int x=x0; x<=x1; x++) {
01874 for (int y=y0; y<=y1; y++) {
01875 data[x + y * nx] +=
01876 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by);
01877 }
01878 }
01879 update();
01880 }
01881 else {
01882 LOGERR("insert_scaled_sum supports only 2D and 3D data");
01883 throw ImageDimensionException("2D/3D only");
01884 }
01885
01886 EXITFUNC;
01887 }
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973