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