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