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 UnexpectedBehaviorException("Bad alloc");
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 dat->update();
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 update();
00225
00226 EXITFUNC;
00227 return this;
00228 }
00229
00230 EMData *EMData::do_ift_cuda()
00231 {
00232 ENTERFUNC;
00233
00234 if (!is_complex()) {
00235 LOGERR("complex image expected. Input image is real image.");
00236 throw ImageFormatException("complex image expected. Input image is real image.");
00237 }
00238
00239 if (!is_ri()) {
00240 throw ImageFormatException("complex ri expected. Got amplitude/phase.");
00241 }
00242
00243 int offset = is_fftodd() ? 1 : 2;
00244 EMData* dat = new EMData(0,0,nx-offset,ny,nz,attr_dict);
00245 if(!dat->rw_alloc()) throw UnexpectedBehaviorException("Bad alloc");
00246
00247 if(cudarwdata == 0){copy_to_cuda();}
00248
00249
00250 int ndim = get_ndim();
00251 if ( ndim == 1 ) {
00252 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, nx-offset,1,1,1);
00253 } else if (ndim == 2) {
00254 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, ny,nx-offset,1,1);
00255 } else if (ndim == 3) {
00256 cuda_dd_fft_complex_to_real_nd(cudarwdata,dat->cudarwdata, nz,ny,nx-offset,1);
00257 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00258
00259
00260 float scale = 1.0f/static_cast<float>((dat->get_size()));
00261 dat->mult(scale);
00262
00263 dat->set_fftpad(false);
00264 dat->set_fftodd(false);
00265 dat->set_complex(false);
00266 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
00267 dat->set_ri(false);
00268
00269 dat->update();
00270
00271 EXITFUNC;
00272 return dat;
00273 }
00274
00275
00276
00277
00278
00279
00280 EMData *EMData::do_ift_inplace_cuda()
00281 {
00282 ENTERFUNC;
00283
00284 if (!is_complex()) {
00285 LOGERR("complex image expected. Input image is real image.");
00286 throw ImageFormatException("complex image expected. Input image is real image.");
00287 }
00288
00289 if (!is_ri()) {
00290 LOGWARN("run IFT on AP data, only RI should be used. ");
00291 }
00292
00293 int offset = is_fftodd() ? 1 : 2;
00294
00295 if(cudarwdata == 0){copy_to_cuda();}
00296
00297 int ndim = get_ndim();
00298 if ( ndim == 1 ) {
00299 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, nx-offset,1,1,1);
00300 } else if (ndim == 2) {
00301 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, ny,nx-offset,1,1);
00302 } else if (ndim == 3) {
00303 cuda_dd_fft_complex_to_real_nd(cudarwdata,cudarwdata, nz,ny,nx-offset,1);
00304 } else throw ImageDimensionException("No cuda FFT support of images with dimensions exceeding 3");
00305 #if defined FFTW2 || defined FFTW3 //native fft and ACML already done normalization
00306
00307 int nxo = nx - offset;
00308 float scale = 1.0f / (nxo * ny * nz);
00309 mult(scale);
00310 #endif //FFTW2 || FFTW3
00311
00312 set_fftpad(true);
00313 set_complex(false);
00314
00315 if(ny==1 && nz==1) set_complex_x(false);
00316 set_ri(false);
00317 update();
00318
00319 EXITFUNC;
00320 return this;
00321 }
00322
00323 #endif //EMAN2_USING_CUDA
00324
00325 EMData *EMData::do_ift()
00326 {
00327 ENTERFUNC;
00328
00329 if (!is_complex()) {
00330 LOGERR("complex image expected. Input image is real image.");
00331 throw ImageFormatException("complex image expected. Input image is real image.");
00332 }
00333
00334 if (!is_ri()) {
00335 LOGWARN("run IFT on AP data, only RI should be used. Converting.");
00336 }
00337
00338 get_data();
00339 EMData* dat = copy_head();
00340 dat->set_size(nx, ny, nz);
00341 ap2ri();
00342
00343 float *d = dat->get_data();
00344 int ndim = get_ndim();
00345
00346
00347
00348
00349 memcpy((char *) d, (char *) rdata, nx * ny * nz * sizeof(float));
00350
00351 int offset = is_fftodd() ? 1 : 2;
00352
00353 if (ndim == 1) {
00354 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00355 } else {
00356 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
00357
00358 size_t row_size = (nx - offset) * sizeof(float);
00359 for (int i = 1; i < ny * nz; i++) {
00360 memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
00361 }
00362
00363 }
00364
00365 dat->set_size(nx - offset, ny, nz);
00366 #if defined FFTW2 || defined FFTW3 //native fft and ACML already done normalization
00367
00368 float scale = 1.0f / ((nx - offset) * ny * nz);
00369 dat->mult(scale);
00370 #endif //FFTW2 || FFTW3
00371 dat->set_fftodd(false);
00372 dat->set_fftpad(false);
00373 dat->set_complex(false);
00374 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
00375 dat->set_ri(false);
00376 dat->update();
00377
00378
00379 EXITFUNC;
00380 return dat;
00381 }
00382
00383
00384
00385
00386
00387 EMData *EMData::do_ift_inplace()
00388 {
00389 ENTERFUNC;
00390
00391 if (!is_complex()) {
00392 LOGERR("complex image expected. Input image is real image.");
00393 throw ImageFormatException("complex image expected. Input image is real image.");
00394 }
00395
00396 if (!is_ri()) {
00397 LOGWARN("run IFT on AP data, only RI should be used. ");
00398 }
00399 ap2ri();
00400
00401 int offset = is_fftodd() ? 1 : 2;
00402 float* data = get_data();
00403 EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
00404
00405 #if defined FFTW2 || defined FFTW3 //native fft and ACML already done normalization
00406
00407 int nxo = nx - offset;
00408 float scale = 1.0f / ((size_t)nxo * ny * nz);
00409 mult(scale);
00410 #endif //FFTW2 || FFTW3
00411
00412 set_fftpad(true);
00413 set_complex(false);
00414 if(ny==1 && nz==1) set_complex_x(false);
00415 set_ri(false);
00416 update();
00417
00418 EXITFUNC;
00419 return this;
00420 }
00421 #undef rdata
00422
00423
00424 std::string EMData::render_ap24(int x0, int y0, int ixsize, int iysize,
00425 int bpl, float scale, int mingray, int maxgray,
00426 float render_min, float render_max,float gamma,int flags)
00427 {
00428 ENTERFUNC;
00429
00430 int asrgb;
00431 int hist=(flags&2)/2;
00432 int invy=(flags&4)?1:0;
00433
00434 if (!is_complex()) throw ImageDimensionException("complex only");
00435
00436 if (get_ndim() != 2) {
00437 throw ImageDimensionException("2D only");
00438 }
00439
00440 if (is_complex()) ri2ap();
00441
00442 if (render_max <= render_min) {
00443 render_max = render_min + 0.01f;
00444 }
00445
00446 if (gamma<=0) gamma=1.0;
00447
00448
00449
00450
00451
00452
00453
00454
00455 static int smg0=0,smg1=0;
00456 static float sgam=0;
00457 static unsigned char gammamap[4096];
00458 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00459 for (int i=0; i<4096; i++) {
00460 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00461 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00462 }
00463 }
00464 smg0=mingray;
00465 smg1=maxgray;
00466 sgam=gamma;
00467
00468 if (flags&8) asrgb=4;
00469 else if (flags&1) asrgb=3;
00470 else throw ImageDimensionException("must set flag 1 or 8");
00471
00472 std::string ret=std::string();
00473
00474 ret.assign(iysize*bpl+hist*1024,char(mingray));
00475 unsigned char *data=(unsigned char *)ret.data();
00476 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00477 if (hist) {
00478 for (int i=0; i<256; i++) histd[i]=0;
00479 }
00480
00481 float rm = render_min;
00482 float inv_scale = 1.0f / scale;
00483 int ysize = iysize;
00484 int xsize = ixsize;
00485
00486 int ymin = 0;
00487 if (iysize * inv_scale > ny) {
00488 ymin = (int) (iysize - ny / inv_scale);
00489 }
00490
00491 float gs = (maxgray - mingray) / (render_max - render_min);
00492 float gs2 = 4095.999f / (render_max - render_min);
00493
00494 if (render_max < render_min) {
00495 gs = 0;
00496 rm = FLT_MAX;
00497 }
00498
00499 int dsx = -1;
00500 int dsy = 0;
00501 int remx = 0;
00502 int remy = 0;
00503 const int scale_n = 100000;
00504
00505 int addi = 0;
00506 int addr = 0;
00507 if (inv_scale == floor(inv_scale)) {
00508 dsx = (int) inv_scale;
00509 dsy = (int) (inv_scale * nx);
00510 }
00511 else {
00512 addi = (int) floor(inv_scale);
00513 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00514 }
00515
00516 int xmin = 0;
00517 if (x0 < 0) {
00518 xmin = (int) (-x0 / inv_scale);
00519 xsize -= (int) floor(x0 / inv_scale);
00520 x0 = 0;
00521 }
00522
00523 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00524 xsize = (int) ((nx - x0) / inv_scale + xmin);
00525 }
00526 int ymax = ysize - 1;
00527 if (y0 < 0) {
00528 ymax = (int) (ysize + y0 / inv_scale - 1);
00529 ymin += (int) floor(y0 / inv_scale);
00530 y0 = 0;
00531 }
00532
00533 if (xmin < 0) xmin = 0;
00534 if (ymin < 0) ymin = 0;
00535 if (xsize > ixsize) xsize = ixsize;
00536 if (ymax > iysize) ymax = iysize;
00537
00538 int lmax = nx * ny - 1;
00539
00540 int mid=nx*ny/2;
00541 float* image_data = get_data();
00542 if (dsx != -1) {
00543 int l = y0 * nx;
00544 for (int j = ymax; j >= ymin; j--) {
00545 int ll = x0;
00546 for (int i = xmin; i < xsize; i++) {
00547 if (l + ll > lmax || ll >= nx - 2) break;
00548
00549 int k = 0;
00550 unsigned char p;
00551 int ph;
00552 if (ll >= nx / 2) {
00553 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00554 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00555 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00556 }
00557 else {
00558 k = nx * ny - (l + 2 * ll) - 2;
00559 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00560 }
00561 if (k>=mid) k-=mid;
00562 else k+=mid;
00563 float t = image_data[k];
00564 if (t <= rm) p = mingray;
00565 else if (t >= render_max) p = maxgray;
00566 else if (gamma!=1.0) {
00567 k=(int)(gs2 * (t-render_min));
00568 p = gammamap[k];
00569 }
00570 else {
00571 p = (unsigned char) (gs * (t - render_min));
00572 p += mingray;
00573 }
00574 if (ph<256) {
00575 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00576 data[i * asrgb + j * bpl+1] = p*ph/256;
00577 data[i * asrgb + j * bpl+2] = 0;
00578 }
00579 else if (ph<512) {
00580 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00581 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00582 data[i * asrgb + j * bpl] = 0;
00583 }
00584 else {
00585 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00586 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00587 data[i * asrgb + j * bpl+1] = 0;
00588 }
00589 if (hist) histd[p]++;
00590 ll += dsx;
00591 }
00592 l += dsy;
00593 }
00594 }
00595 else {
00596 remy = 10;
00597 int l = y0 * nx;
00598 for (int j = ymax; j >= ymin; j--) {
00599 int br = l;
00600 remx = 10;
00601 int ll = x0;
00602 for (int i = xmin; i < xsize - 1; i++) {
00603 if (l + ll > lmax || ll >= nx - 2) {
00604 break;
00605 }
00606 int k = 0;
00607 unsigned char p;
00608 int ph;
00609 if (ll >= nx / 2) {
00610 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
00611 else k = 2 * (ll - nx / 2) + l + 2 + nx;
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 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00617 }
00618 if (k>=mid) k-=mid;
00619 else k+=mid;
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
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975 Util::ap2ri(get_data(), (size_t)nx * ny * nz);
00976 set_ri(true);
00977 update();
00978 EXITFUNC;
00979 }
00980
00981 void EMData::ri2inten()
00982 {
00983 ENTERFUNC;
00984
00985 if (!is_complex()) return;
00986 if (!is_ri()) ap2ri();
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999 float * data = get_data();
01000 size_t size = (size_t)nx * ny * nz;
01001 for (size_t i = 0; i < size; i += 2) {
01002 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
01003 data[i+1]=0;
01004 }
01005
01006 set_attr("is_intensity", int(1));
01007 update();
01008 EXITFUNC;
01009 }
01010
01011
01012 void EMData::ri2ap()
01013 {
01014 ENTERFUNC;
01015
01016 if (!is_complex() || !is_ri()) {
01017 return;
01018 }
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030 float * data = get_data();
01031
01032 size_t size = (size_t)nx * ny * nz;
01033 for (size_t i = 0; i < size; i += 2) {
01034 #ifdef _WIN32
01035 float f = (float)_hypot(data[i], data[i + 1]);
01036 #else
01037 float f = (float)hypot(data[i], data[i + 1]);
01038 #endif
01039 if (data[i] == 0 && data[i + 1] == 0) {
01040 data[i + 1] = 0;
01041 }
01042 else {
01043 data[i + 1] = atan2(data[i + 1], data[i]);
01044 }
01045 data[i] = f;
01046 }
01047
01048 set_ri(false);
01049 update();
01050 EXITFUNC;
01051 }
01052
01053
01054 float calc_bessel(const int n, const float& x) {
01055 gsl_sf_result result;
01056
01057 gsl_sf_bessel_Jn_e(n,(double)x, &result);
01058 return (float)result.val;
01059 }
01060
01061 EMData* EMData::bispecRotTransInvN(int N, int NK)
01062 {
01063
01064 int EndP = this -> get_xsize();
01065 int Mid = (int) ((1+EndP)/2);
01066 int End = 2*Mid-1;
01067
01068 int CountxyMax = End*End;
01069
01070 int *SortfkInds = new int[CountxyMax];
01071 int *kVecX = new int[CountxyMax];
01072 int *kVecY = new int[CountxyMax];
01073 float *fkVecR = new float[CountxyMax];
01074 float *fkVecI = new float[CountxyMax];
01075 float *absD1fkVec = new float[CountxyMax];
01076 float *absD1fkVecSorted = new float[CountxyMax];
01077
01078 float *jxjyatan2 = new float[End*End];
01079
01080 EMData * ThisCopy = new EMData(End,End);
01081
01082 for (int jx=0; jx <End ; jx++) {
01083 for (int jy=0; jy <End ; jy++) {
01084 float ValNow = this -> get_value_at(jx,jy);
01085 ThisCopy -> set_value_at(jx,jy,ValNow);
01086
01087 }}
01088
01089
01090 EMData* fk = ThisCopy -> do_fft();
01091 fk ->process_inplace("xform.fourierorigin.tocenter");
01092
01093
01094 EMData* fkRCopy = new EMData(End,End);
01095 EMData* fkICopy = new EMData(End,End);
01096 EMData* fkCopy = new EMData(End,End);
01097
01098
01099 for (int jCount= 0; jCount<End*End; jCount++) {
01100
01101 int jx = jCount%End ;
01102 int jy = (jCount-jx)/End ;
01103 jxjyatan2[jCount] = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
01104 }
01105
01106
01107 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01108
01109 int kx = kEx/2;
01110 int kIx = kx+Mid-1;
01111 int kCx = -kx ;
01112 int kCIx = kCx+ Mid-1 ;
01113 for (int kEy= 0 ; kEy<End; kEy++) {
01114 int kIy = kEy ;
01115 int ky = kEy+1-Mid;
01116 float realVal = fk -> get_value_at(kEx ,kEy) ;
01117 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01118 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01119 float fkAng = atan2(imagVal,realVal);
01120
01121 float NewRealVal ;
01122 float NewImagVal ;
01123 float AngMatlab ;
01124
01125 if (kIx==Mid-1) {
01126
01127
01128 }
01129
01130 if (kIx>Mid-1){
01131
01132 }
01133
01134 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01135 NewRealVal = absVal*cos(AngMatlab);
01136 NewImagVal = absVal*sin(AngMatlab);
01137
01138
01139 fkVecR[kIy+kIx *End] = NewRealVal ;
01140 fkVecR[kIy+kCIx*End] = NewRealVal ;
01141 fkVecI[kIy+kIx *End] = NewImagVal ;
01142 fkVecI[kIy+kCIx*End] = -NewImagVal ;
01143 absD1fkVec[kIy + kIx *End] = absVal;
01144 absD1fkVec[kIy + kCIx *End] = absVal;
01145 kVecX[kIy+kIx *End] = kx ;
01146 kVecX[kIy+kCIx *End] = kCx ;
01147 kVecY[kIy+kIx *End] = ky ;
01148 kVecY[kIy+kCIx *End] = ky ;
01149
01150
01151
01152
01153
01154 fkCopy -> set_value_at(kIx ,kIy, absVal);
01155 fkCopy -> set_value_at(kCIx,kIy, absVal);
01156 fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
01157 fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
01158 fkICopy -> set_value_at(kIx, kIy, NewImagVal);
01159 fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
01160
01161 }
01162 }
01163 #ifdef _WIN32
01164 _unlink("fkCopy.???");
01165 _unlink("fk?Copy.???");
01166 #else
01167 system("rm -f fkCopy.???");
01168 system("rm -f fk?Copy.???");
01169 #endif //_WIN32
01170 fkCopy -> write_image("fkCopy.img");
01171 fkRCopy -> write_image("fkRCopy.img");
01172 fkICopy -> write_image("fkICopy.img");
01173
01174 cout << "Starting the sort "<< endl;
01175
01176 vector< pair<float, int> > absInds;
01177 for(int i = 0; i < CountxyMax; ++i ) {
01178 pair<float,int> p;
01179 p = make_pair(absD1fkVec[i],i);
01180 absInds.push_back( p);
01181 }
01182
01183 std::sort(absInds.begin(),absInds.end());
01184
01185 for(int i = 0; i < CountxyMax; ++i ) {
01186 pair<float,int> p ;
01187 p = absInds[i] ;
01188 absD1fkVecSorted[CountxyMax-1-i] = p.first ;
01189 SortfkInds[CountxyMax-1-i] = p.second ;
01190 }
01191
01192 cout << "Ending the sort "<< endl;
01193
01194
01195
01196
01197 for(int i = 0; i < CountxyMax; ++i ) {
01198 int Si = SortfkInds[i];
01199 int kIx = (int) Si/End; kIx = (int) i/End;
01200
01201
01202
01203
01204 }
01205 cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0] << endl;
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229 cout << "NK = " << NK << endl;
01230 float frR= 3.0/4.0;
01231 int LradRange= (int) (floor(Mid/frR)) ;
01232
01233 float *radRange = new float[LradRange];
01234 radRange[0]=0;
01235 for (int irad=1; irad < LradRange; irad++){
01236 radRange[irad] = radRange[irad-1] + frR; }
01237
01238
01239
01240
01241 cout << "Starting the calculation of invariants for N= " << N << endl;
01242
01243
01244
01245 EMData* RotTransInv = new EMData();
01246 RotTransInv -> set_size(LradRange,LradRange);
01247
01248
01249
01250
01251
01252
01253
01254 for (int jr1=0; jr1 < LradRange ; jr1++ ) {
01255 float r1= radRange[jr1];
01256
01257 for (int jr2=0; jr2<LradRange; jr2++ ) {
01258 float r2= radRange[jr2];
01259 float RotTransInvTemp=0;
01260 for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
01261 int Countkxy =SortfkInds[jCountkxy] ;
01262 int kx = kVecX[Countkxy] ;
01263 int ky = kVecY[Countkxy] ;
01264 float k2 = (float)(kx*kx+ky*ky);
01265 if (k2==0) { continue;}
01266 float phiK =0;
01267 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; phiK= atan2((float)ky,(float)kx);
01268
01269 float fkR = fkVecR[Countkxy] ;
01270 float fkI = fkVecI[Countkxy] ;
01271
01272
01273 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
01274 int Countqxy =SortfkInds[jCountqxy] ;
01275 int qx = kVecX[Countqxy] ;
01276 int qy = kVecY[Countqxy] ;
01277 int q2 = qx*qx+qy*qy;
01278 if (q2==0) {continue;}
01279 float phiQ =0;
01280 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; phiQ=atan2((float)qy,(float)qx);
01281 float fqR = fkVecR[Countqxy] ;
01282 float fqI = fkVecI[Countqxy] ;
01283 int kCx = (-kx-qx);
01284 int kCy = (-ky-qy);
01285 int kCIx = ((kCx+Mid+2*End)%End);
01286 int kCIy = ((kCy+Mid+2*End)%End);
01287 kCx = kCIx-Mid;
01288 kCy = kCIy-Mid;
01289 int CountCxy = kCIx*End+kCIy;
01290 float fCR = fkVecR[CountCxy];
01291 float fCI = fkVecI[CountCxy];
01292 if (jr1+jr2==-1) {
01293 printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
01294 printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f, kCx=%d, kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01295 }
01296 for (int p=0; p<NK; p++){
01297
01298 if (SortfkInds[p]==CountCxy){
01299 float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
01300 float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
01301
01302
01303 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI))
01304 * cos(N*(phiK-phiQ+M_PI));
01305 bispectemp -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
01306 * sin(N*(phiK-phiQ+M_PI));
01307 float bess1 = calc_bessel(N, Arg1 );
01308 float bess2 = calc_bessel(N, Arg2 );
01309
01310
01311
01312 RotTransInvTemp = RotTransInvTemp + bispectemp * bess1*bess2 ;
01313
01314 }
01315 }
01316 }
01317 }
01318 RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp) ;
01319
01320 }
01321 }
01322
01323
01324 return RotTransInv ;
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
01355
01356
01357
01358 EMData* EMData::bispecRotTransInvDirect(int type)
01359 {
01360
01361 int EndP = this -> get_xsize();
01362 int Mid = (int) ((1+EndP)/2);
01363 int End = 2*Mid-1;
01364
01365 int CountxyMax = End*End;
01366
01367
01368 int *kVecX = new int[CountxyMax];
01369 int *kVecY = new int[CountxyMax];
01370 float *fkVecR = new float[CountxyMax];
01371 float *fkVecI = new float[CountxyMax];
01372 float *absD1fkVec = new float[CountxyMax];
01373
01374
01375
01376 float *jxjyatan2 = new float[End*End];
01377
01378
01379 EMData * ThisCopy = new EMData(End,End);
01380
01381 for (int jx=0; jx <End ; jx++) {
01382 for (int jy=0; jy <End ; jy++) {
01383 float ValNow = this -> get_value_at(jx,jy);
01384 ThisCopy -> set_value_at(jx,jy,ValNow);
01385 jxjyatan2[jy*End + jx] = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
01386
01387 }}
01388
01389
01390 EMData* fk = ThisCopy -> do_fft();
01391 fk ->process_inplace("xform.fourierorigin.tocenter");
01392
01393
01394
01395 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01396
01397 int kx = kEx/2;
01398 int kIx = kx+Mid-1;
01399 int kCx = -kx ;
01400 int kCIx = kCx+ Mid-1 ;
01401 for (int kEy= 0 ; kEy<End; kEy++) {
01402 int kIy = kEy ;
01403 int ky = kEy+1-Mid;
01404 float realVal = fk -> get_value_at(kEx ,kEy) ;
01405 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01406 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01407 float fkAng = atan2(imagVal,realVal);
01408
01409 float NewRealVal ;
01410 float NewImagVal ;
01411 float AngMatlab ;
01412
01413 if (kIx==Mid-1) {
01414
01415 }
01416
01417 if (kIx>Mid-1){
01418
01419 }
01420
01421 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01422 NewRealVal = absVal*cos(AngMatlab);
01423 NewImagVal = absVal*sin(AngMatlab);
01424
01425
01426 fkVecR[ kIy +kIx *End] = NewRealVal ;
01427 fkVecR[(End-1-kIy)+kCIx*End] = NewRealVal ;
01428 fkVecI[ kIy +kIx *End] = NewImagVal ;
01429 fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
01430 absD1fkVec[(End-1-kIy) + kIx *End] = absVal;
01431 absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
01432 kVecX[kIy+kIx *End] = kx ;
01433 kVecX[kIy+kCIx *End] = kCx ;
01434 kVecY[kIy+kIx *End] = ky ;
01435 kVecY[kIy+kCIx *End] = ky ;
01436
01437
01438
01439
01440
01441
01442
01443
01444 }
01445 }
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457 float frR= 3.0/4.0;
01458 frR= 1;
01459 int LradRange= (int) (1+floor(Mid/frR -.1)) ;
01460
01461 float *radRange = new float[LradRange];
01462 for (int irad=0; irad < LradRange; irad++){
01463 radRange[irad] = frR*irad;
01464
01465 }
01466
01467 cout << "Starting the calculation of invariants" << endl;
01468
01469
01470 if (type==0) {
01471 int LthetaRange = 59;
01472 float ftR = (2.0f*M_PI/LthetaRange );
01473 float *thetaRange = new float[LthetaRange];
01474
01475 for (int ith=0; ith < LthetaRange; ith++){
01476 thetaRange[ith] = ftR*ith; }
01477
01478 int TotalVol = LradRange*LradRange*LthetaRange;
01479
01480 float *RotTransInv = new float[TotalVol];
01481 float *WeightInv = new float[TotalVol];
01482
01483 for (int jW=0; jW<TotalVol; jW++) {
01484 RotTransInv[jW] = 0;
01485 WeightInv[jW] = 0;
01486 }
01487
01488 for (int jW=0; jW<TotalVol; jW++) {
01489 RotTransInv[jW] = 0;
01490 WeightInv[jW] = 0;
01491 }
01492
01493
01494
01495 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01496 int kx = kVecX[Countkxy] ;
01497 int ky = kVecY[Countkxy] ;
01498 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01499 float phiK =0;
01500 if (k2>0) phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1];
01501 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01502 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01503
01504
01505 if ((k2==0)|| (k2>Mid) ) { continue;}
01506
01507 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01508 int qx = kVecX[Countqxy] ;
01509 int qy = kVecY[Countqxy] ;
01510 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01511 if ((q2==0)|| (q2>Mid) ) {continue;}
01512 float phiQ =0;
01513 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1];
01514 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01515 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01516 int kCx = (-kx-qx);
01517 int kCy = (-ky-qy);
01518 int kCIx = ((kCx+Mid+2*End)%End);
01519 int kCIy = ((kCy+Mid+2*End)%End);
01520 kCx = ((kCIx+End-1)%End)+1-Mid;
01521 kCy = ((kCIy+End-1)%End)+1-Mid ;
01522
01523
01524 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01525 float fCR = fkVecR[CountCxy];
01526 float fCI = fkVecI[CountCxy];
01527
01528
01529
01530
01531
01532 float phiQK = (4*M_PI+phiQ-phiK);
01533 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
01534
01535
01536
01537 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01538
01539 if ((q2<k2) ) continue;
01540
01541
01542
01543
01544
01545
01546 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01547 int k2IndHi = k2IndLo;
01548 float k2Lo= radRange[k2IndLo];
01549 if (k2IndLo+1< LradRange) {
01550 k2IndHi = k2IndLo+1;
01551 }
01552
01553
01554 float kCof =k2-k2Lo;
01555
01556 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01557 int q2IndHi=q2IndLo;
01558 float q2Lo= radRange[q2IndLo];
01559 if (q2IndLo+1 < LradRange) {
01560 q2IndHi = q2IndLo+1 ;
01561 }
01562 float qCof = q2-q2Lo;
01563
01564 if ((qCof<0) || (qCof >1) ) {
01565 cout<< "Weird! qCof="<< qCof << " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
01566 int x ;
01567 cin >> x ;
01568 }
01569
01570 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
01571 int thetaIndHi = thetaIndLo;
01572
01573 float thetaLo = thetaRange[thetaIndLo];
01574 float thetaHi = thetaLo;
01575 float thetaCof = 0;
01576
01577 if (thetaIndLo+1< LthetaRange) {
01578 thetaIndHi = thetaIndLo +1;
01579 }else{
01580 thetaIndHi=0;
01581 }
01582
01583 thetaHi = thetaRange[thetaIndHi];
01584
01585 if (thetaHi==thetaLo) {
01586 thetaCof =0 ;
01587 } else {
01588 thetaCof = (phiQK-thetaLo)/(thetaHi-thetaLo);
01589 }
01590
01591 if ((thetaCof>2*M_PI) ) {
01592 cout<< "Weird! thetaCof="<< thetaCof <<endl ;
01593 thetaCof=0;
01594 }
01595
01596
01597
01598
01599
01600 for (int jk =1; jk<=2; jk++){
01601 for (int jq =1; jq<=2; jq++){
01602 for (int jtheta =1; jtheta<=2; jtheta++){
01603
01604 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
01605 * (thetaCof+(1-2*thetaCof)*(jtheta==1));
01606
01607
01608 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01609 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01610 int thetaInd = thetaIndLo*(jtheta==1) + thetaIndHi*(jtheta ==2);
01611 int TotalInd = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
01612
01613
01614
01615
01616
01617
01618 RotTransInv[TotalInd] += Weight*bispectemp;
01619 WeightInv[TotalInd] += Weight;
01620
01621 }}}
01622 }
01623 }
01624
01625 cout << "Finished Main Section " << endl;
01626
01627
01628
01629 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
01630 EMData *RotTransInvF = new EMData(LradRange,LradRange,LthetaRange);
01631 EMData *WeightImage = new EMData(LradRange,LradRange,LthetaRange);
01632
01633
01634
01635
01636
01637
01638
01639 for (int jtheta =0; jtheta < LthetaRange; jtheta++){
01640 for (int jq =0; jq<LradRange; jq++){
01641 for (int jk =0; jk<LradRange ; jk++){
01642
01643 int TotalInd = jtheta*LradRange*LradRange+jq*LradRange+jk;
01644 float Weight = WeightInv[TotalInd];
01645 WeightImage -> set_value_at(jk,jq,jtheta,Weight);
01646 RotTransInvF -> set_value_at(jk,jq,jtheta,0);
01647 if (Weight <= 0) continue;
01648 RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);
01649 int newjtheta = (LthetaRange- jtheta)%LthetaRange;
01650 RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );
01651 }
01652 }
01653 }
01654
01655 cout << " Almost Done " << endl;
01656 #ifdef _WIN32
01657 _unlink("WeightImage.???");
01658 #else
01659 system("rm -f WeightImage.???");
01660 #endif //_WIN32
01661 WeightImage -> write_image("WeightImage.img");
01662
01663 return RotTransInvF ;
01664 }
01665
01666 if (type==1) {
01667 int TotalVol = LradRange*LradRange;
01668
01669 float *RotTransInv = new float[TotalVol];
01670 float *WeightInv = new float[TotalVol];
01671
01672 for (int jW=0; jW<TotalVol; jW++) {
01673 RotTransInv[jW] = 0;
01674 WeightInv[jW] = 0;
01675 }
01676
01677
01678 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01679 int kx = kVecX[Countkxy] ;
01680 int ky = kVecY[Countkxy] ;
01681 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01682 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01683 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01684
01685
01686 if ((k2==0)|| (k2>Mid) ) { continue;}
01687
01688 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01689
01690
01691 int qx = kVecX[Countqxy] ;
01692 int qy = kVecY[Countqxy] ;
01693 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01694 if ((q2==0)|| (q2>Mid) ) {continue;}
01695 if ((q2<k2) ) continue;
01696
01697 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01698 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01699
01700 int kCx = (-kx-qx);
01701 int kCy = (-ky-qy);
01702 int kCIx = ((kCx+Mid+2*End)%End);
01703 int kCIy = ((kCy+Mid+2*End)%End);
01704 kCx = ((kCIx+End-1)%End)+1-Mid;
01705 kCy = ((kCIy+End-1)%End)+1-Mid ;
01706
01707
01708 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01709 float fCR = fkVecR[CountCxy];
01710 float fCI = fkVecI[CountCxy];
01711
01712
01713 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01714
01715
01716 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01717 int k2IndHi = k2IndLo;
01718 float k2Lo= radRange[k2IndLo];
01719 if (k2IndLo+1< LradRange) {
01720 k2IndHi = k2IndLo+1;
01721 }
01722
01723
01724 float kCof =k2-k2Lo;
01725
01726 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01727 int q2IndHi=q2IndLo;
01728 float q2Lo= radRange[q2IndLo];
01729 if (q2IndLo+1 < LradRange) {
01730 q2IndHi = q2IndLo+1 ;
01731 }
01732 float qCof = q2-q2Lo;
01733
01734
01735 for (int jk =1; jk<=2; jk++){
01736 for (int jq =1; jq<=2; jq++){
01737
01738 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
01739
01740 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01741 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01742 int TotalInd = q2Ind*LradRange+k2Ind;
01743 RotTransInv[TotalInd] += Weight*bispectemp;
01744 WeightInv[TotalInd] += Weight;
01745
01746 }}
01747 }
01748 }
01749
01750
01751
01752
01753
01754 EMData *RotTransInvF = new EMData(LradRange,LradRange);
01755 EMData *WeightImage = new EMData(LradRange,LradRange);
01756
01757 for (int jk =0; jk<LradRange ; jk++){
01758 for (int jq =jk; jq<LradRange; jq++){
01759 int TotalInd = jq*LradRange+jk;
01760 int TotalIndBar = jq*LradRange+jk;
01761 float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
01762 if (Weight <=0) continue;
01763 WeightImage -> set_value_at(jk,jq,Weight);
01764 WeightImage -> set_value_at(jq,jk,Weight);
01765 #ifdef _WIN32
01766 float ValNow = pow( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight, 1.0f/3.0f ) ;
01767 #else
01768 float ValNow = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight ) ;
01769 #endif //_WIN32
01770 RotTransInvF -> set_value_at(jk,jq,ValNow);
01771 RotTransInvF -> set_value_at(jq,jk,ValNow );
01772 }}
01773
01774 #ifdef _WIN32
01775 _unlink("WeightImage.???");
01776 #else
01777 system("rm -f WeightImage.???");
01778 #endif //_WIN32
01779 WeightImage -> write_image("WeightImage.img");
01780
01781 return RotTransInvF ;
01782 }
01783 return 0;
01784 }
01785
01786
01787 void EMData::insert_clip(const EMData * const block, const IntPoint &origin) {
01788 int nx1 = block->get_xsize();
01789 int ny1 = block->get_ysize();
01790 int nz1 = block->get_zsize();
01791
01792 Region area(origin[0], origin[1], origin[2],nx1, ny1, nz1);
01793
01794
01795 int x0 = (int) area.origin[0];
01796 x0 = x0 < 0 ? 0 : x0;
01797
01798 int y0 = (int) area.origin[1];
01799 y0 = y0 < 0 ? 0 : y0;
01800
01801 int z0 = (int) area.origin[2];
01802 z0 = z0 < 0 ? 0 : z0;
01803
01804 int x1 = (int) (area.origin[0] + area.size[0]);
01805 x1 = x1 > nx ? nx : x1;
01806
01807 int y1 = (int) (area.origin[1] + area.size[1]);
01808 y1 = y1 > ny ? ny : y1;
01809
01810 int z1 = (int) (area.origin[2] + area.size[2]);
01811 z1 = z1 > nz ? nz : z1;
01812 if (z1 <= 0) {
01813 z1 = 1;
01814 }
01815
01816 int xd0 = (int) (area.origin[0] < 0 ? -area.origin[0] : 0);
01817 int yd0 = (int) (area.origin[1] < 0 ? -area.origin[1] : 0);
01818 int zd0 = (int) (area.origin[2] < 0 ? -area.origin[2] : 0);
01819
01820 if (x1 < x0 || y1 < y0 || z1 < z0) return;
01821
01822 size_t clipped_row_size = (x1-x0) * sizeof(float);
01823 int src_secsize = nx1 * ny1;
01824 int dst_secsize = nx * ny;
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847 float *src = block->get_data() + zd0 * src_secsize + yd0 * nx1 + xd0;
01848 float *dst = get_data() + z0 * dst_secsize + y0 * nx + x0;
01849
01850 int src_gap = src_secsize - (y1-y0) * nx1;
01851 int dst_gap = dst_secsize - (y1-y0) * nx;
01852
01853 for (int i = z0; i < z1; i++) {
01854 for (int j = y0; j < y1; j++) {
01855 EMUtil::em_memcpy(dst, src, clipped_row_size);
01856 src += nx1;
01857 dst += nx;
01858 }
01859 src += src_gap;
01860 dst += dst_gap;
01861 }
01862
01863 #ifdef EMAN2_USING_CUDA
01864 if(block->cudarwdata){
01865 copy_to_cuda();
01866 }
01867 #endif
01868
01869 update();
01870 EXITFUNC;
01871 }
01872
01873
01874 void EMData::insert_scaled_sum(EMData *block, const FloatPoint ¢er,
01875 float scale, float)
01876 {
01877 ENTERFUNC;
01878 float * data = get_data();
01879 if (get_ndim()==3) {
01880
01881 int xs=(int)floor(block->get_xsize()*scale/2.0);
01882 int ys=(int)floor(block->get_ysize()*scale/2.0);
01883 int zs=(int)floor(block->get_zsize()*scale/2.0);
01884 int x0=(int)center[0]-xs;
01885 int x1=(int)center[0]+xs;
01886 int y0=(int)center[1]-ys;
01887 int y1=(int)center[1]+ys;
01888 int z0=(int)center[2]-zs;
01889 int z1=(int)center[2]+zs;
01890
01891 if (x1<0||y1<0||z1<0||x0>get_xsize()||y0>get_ysize()||z0>get_zsize()) return;
01892
01893
01894 if (x0<0) x0=0;
01895 if (y0<0) y0=0;
01896 if (z0<0) z0=0;
01897 if (x1>=get_xsize()) x1=get_xsize()-1;
01898 if (y1>=get_ysize()) y1=get_ysize()-1;
01899 if (z1>=get_zsize()) z1=get_zsize()-1;
01900
01901 float bx=block->get_xsize()/2.0f;
01902 float by=block->get_ysize()/2.0f;
01903 float bz=block->get_zsize()/2.0f;
01904
01905 size_t idx;
01906 for (int x=x0; x<=x1; x++) {
01907 for (int y=y0; y<=y1; y++) {
01908 for (int z=z0; z<=z1; z++) {
01909 idx = x + y * nx + (size_t)z * nx * ny;
01910 data[idx] +=
01911 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by,(z-center[2])/scale+bz);
01912 }
01913 }
01914 }
01915 update();
01916 }
01917 else if (get_ndim()==2) {
01918
01919 int xs=(int)floor(block->get_xsize()*scale/2.0);
01920 int ys=(int)floor(block->get_ysize()*scale/2.0);
01921 int x0=(int)center[0]-xs;
01922 int x1=(int)center[0]+xs;
01923 int y0=(int)center[1]-ys;
01924 int y1=(int)center[1]+ys;
01925
01926 if (x1<0||y1<0||x0>get_xsize()||y0>get_ysize()) return;
01927
01928
01929 if (x0<0) x0=0;
01930 if (y0<0) y0=0;
01931 if (x1>=get_xsize()) x1=get_xsize()-1;
01932 if (y1>=get_ysize()) y1=get_ysize()-1;
01933
01934 float bx=block->get_xsize()/2.0f;
01935 float by=block->get_ysize()/2.0f;
01936
01937 for (int x=x0; x<=x1; x++) {
01938 for (int y=y0; y<=y1; y++) {
01939 data[x + y * nx] +=
01940 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by);
01941 }
01942 }
01943 update();
01944 }
01945 else {
01946 LOGERR("insert_scaled_sum supports only 2D and 3D data");
01947 throw ImageDimensionException("2D/3D only");
01948 }
01949
01950 EXITFUNC;
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
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037