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, (size_t)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 (size_t i = 1; i < (size_t)ny * nz; i++) {
00360 memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
00361 }
00362 }
00363
00364 dat->set_size(nx - offset, ny, nz);
00365 #if defined FFTW2 || defined FFTW3 //native fft and ACML already done normalization
00366
00367 float scale = 1.0f / ((nx - offset) * ny * nz);
00368 dat->mult(scale);
00369 #endif //FFTW2 || FFTW3
00370 dat->set_fftodd(false);
00371 dat->set_fftpad(false);
00372 dat->set_complex(false);
00373 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
00374 dat->set_ri(false);
00375 dat->update();
00376
00377
00378 EXITFUNC;
00379 return dat;
00380 }
00381
00382
00383
00384
00385
00386 EMData *EMData::do_ift_inplace()
00387 {
00388 ENTERFUNC;
00389
00390 if (!is_complex()) {
00391 LOGERR("complex image expected. Input image is real image.");
00392 throw ImageFormatException("complex image expected. Input image is real image.");
00393 }
00394
00395 if (!is_ri()) {
00396 LOGWARN("run IFT on AP data, only RI should be used. ");
00397 }
00398 ap2ri();
00399
00400 int offset = is_fftodd() ? 1 : 2;
00401 float* data = get_data();
00402 EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
00403
00404 #if defined FFTW2 || defined FFTW3 //native fft and ACML already done normalization
00405
00406 int nxo = nx - offset;
00407 float scale = 1.0f / ((size_t)nxo * ny * nz);
00408 mult(scale);
00409 #endif //FFTW2 || FFTW3
00410
00411 set_fftpad(true);
00412 set_complex(false);
00413 if(ny==1 && nz==1) set_complex_x(false);
00414 set_ri(false);
00415 update();
00416
00417 EXITFUNC;
00418 return this;
00419 }
00420 #undef rdata
00421
00422
00423 std::string EMData::render_ap24(int x0, int y0, int ixsize, int iysize,
00424 int bpl, float scale, int mingray, int maxgray,
00425 float render_min, float render_max,float gamma,int flags)
00426 {
00427 ENTERFUNC;
00428
00429 int asrgb;
00430 int hist=(flags&2)/2;
00431 int invy=(flags&4)?1:0;
00432
00433 if (!is_complex()) throw ImageDimensionException("complex only");
00434
00435 if (get_ndim() != 2) {
00436 throw ImageDimensionException("2D only");
00437 }
00438
00439 if (is_complex()) ri2ap();
00440
00441 if (render_max <= render_min) {
00442 render_max = render_min + 0.01f;
00443 }
00444
00445 if (gamma<=0) gamma=1.0;
00446
00447
00448
00449
00450
00451
00452
00453
00454 static int smg0=0,smg1=0;
00455 static float sgam=0;
00456 static unsigned char gammamap[4096];
00457 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
00458 for (int i=0; i<4096; i++) {
00459 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00460 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
00461 }
00462 }
00463 smg0=mingray;
00464 smg1=maxgray;
00465 sgam=gamma;
00466
00467 if (flags&8) asrgb=4;
00468 else if (flags&1) asrgb=3;
00469 else throw ImageDimensionException("must set flag 1 or 8");
00470
00471 std::string ret=std::string();
00472
00473 ret.assign(iysize*bpl+hist*1024,char(mingray));
00474 unsigned char *data=(unsigned char *)ret.data();
00475 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
00476 if (hist) {
00477 for (int i=0; i<256; i++) histd[i]=0;
00478 }
00479
00480 float rm = render_min;
00481 float inv_scale = 1.0f / scale;
00482 int ysize = iysize;
00483 int xsize = ixsize;
00484
00485 int ymin = 0;
00486 if (iysize * inv_scale > ny) {
00487 ymin = (int) (iysize - ny / inv_scale);
00488 }
00489
00490 float gs = (maxgray - mingray) / (render_max - render_min);
00491 float gs2 = 4095.999f / (render_max - render_min);
00492
00493 if (render_max < render_min) {
00494 gs = 0;
00495 rm = FLT_MAX;
00496 }
00497
00498 int dsx = -1;
00499 int dsy = 0;
00500 int remx = 0;
00501 int remy = 0;
00502 const int scale_n = 100000;
00503
00504 int addi = 0;
00505 int addr = 0;
00506 if (inv_scale == floor(inv_scale)) {
00507 dsx = (int) inv_scale;
00508 dsy = (int) (inv_scale * nx);
00509 }
00510 else {
00511 addi = (int) floor(inv_scale);
00512 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00513 }
00514
00515 int xmin = 0;
00516 if (x0 < 0) {
00517 xmin = (int) (-x0 / inv_scale);
00518 xsize -= (int) floor(x0 / inv_scale);
00519 x0 = 0;
00520 }
00521
00522 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00523 xsize = (int) ((nx - x0) / inv_scale + xmin);
00524 }
00525 int ymax = ysize - 1;
00526 if (y0 < 0) {
00527 ymax = (int) (ysize + y0 / inv_scale - 1);
00528 ymin += (int) floor(y0 / inv_scale);
00529 y0 = 0;
00530 }
00531
00532 if (xmin < 0) xmin = 0;
00533 if (ymin < 0) ymin = 0;
00534 if (xsize > ixsize) xsize = ixsize;
00535 if (ymax > iysize) ymax = iysize;
00536
00537 int lmax = nx * ny - 1;
00538
00539 int mid=nx*ny/2;
00540 float* image_data = get_data();
00541 if (dsx != -1) {
00542 int l = y0 * nx;
00543 for (int j = ymax; j >= ymin; j--) {
00544 int ll = x0;
00545 for (int i = xmin; i < xsize; i++) {
00546 if (l + ll > lmax || ll >= nx - 2) break;
00547
00548 int k = 0;
00549 unsigned char p;
00550 int ph;
00551 if (ll >= nx / 2) {
00552 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
00553 else k = 2 * (ll - nx / 2) + l + 2 + nx;
00554 if (k>=mid) k-=mid;
00555 else k+=mid;
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 if (k>=mid) k-=mid;
00562 else k+=mid;
00563 }
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 if (k>=mid) k-=mid;
00614 else k+=mid;
00615 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384;
00616 }
00617 else {
00618 k = nx * ny - (l + 2 * ll) - 2;
00619 if (k>=mid) k-=mid;
00620 else k+=mid;
00621 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384;
00622 }
00623
00624 float t = image_data[k];
00625 if (t <= rm)
00626 p = mingray;
00627 else if (t >= render_max) {
00628 p = maxgray;
00629 }
00630 else if (gamma!=1.0) {
00631 k=(int)(gs2 * (t-render_min));
00632 p = gammamap[k];
00633 }
00634 else {
00635 p = (unsigned char) (gs * (t - render_min));
00636 p += mingray;
00637 }
00638 if (ph<256) {
00639 data[i * asrgb + j * bpl] = p*(255-ph)/256;
00640 data[i * asrgb + j * bpl+1] = p*ph/256;
00641 data[i * asrgb + j * bpl+2] = 0;
00642 }
00643 else if (ph<512) {
00644 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
00645 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
00646 data[i * asrgb + j * bpl] = 0;
00647 }
00648 else {
00649 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
00650 data[i * asrgb + j * bpl] = p*(ph-512)/256;
00651 data[i * asrgb + j * bpl+1] = 0;
00652 }
00653 if (hist) histd[p]++;
00654 ll += addi;
00655 remx += addr;
00656 if (remx > scale_n) {
00657 remx -= scale_n;
00658 ll++;
00659 }
00660 }
00661 l = br + addi * nx;
00662 remy += addr;
00663 if (remy > scale_n) {
00664 remy -= scale_n;
00665 l += nx;
00666 }
00667 }
00668 }
00669
00670
00671 if (asrgb==4) {
00672 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
00673 for (int i=xmin; i<xsize*4; i+=4) {
00674 data[i+j+3]=255;
00675 }
00676 }
00677 }
00678
00679 EXITFUNC;
00680
00681
00682 if (invy) {
00683 int x,y;
00684 char swp;
00685 for (y=0; y<iysize/2; y++) {
00686 for (x=0; x<ixsize; x++) {
00687 swp=ret[y*bpl+x];
00688 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
00689 ret[(iysize-y-1)*bpl+x]=swp;
00690 }
00691 }
00692 }
00693
00694
00695 return ret;
00696 }
00697
00698
00699 void EMData::render_amp24( int x0, int y0, int ixsize, int iysize,
00700 int bpl, float scale, int mingray, int maxgray,
00701 float render_min, float render_max, void *ref,
00702 void cmap(void *, int coord, unsigned char *tri))
00703 {
00704 ENTERFUNC;
00705
00706 if (get_ndim() != 2) {
00707 throw ImageDimensionException("2D only");
00708 }
00709
00710 if (is_complex()) {
00711 ri2ap();
00712 }
00713
00714 if (render_max <= render_min) {
00715 render_max = render_min + 0.01f;
00716 }
00717
00718 std::string ret=std::string();
00719 ret.resize(iysize*bpl);
00720 unsigned char *data=(unsigned char *)ret.data();
00721
00722 float rm = render_min;
00723 float inv_scale = 1.0f / scale;
00724 int ysize = iysize;
00725 int xsize = ixsize;
00726 const int scale_n = 100000;
00727
00728 int ymin = 0;
00729 if ( iysize * inv_scale > ny) {
00730 ymin = (int) (iysize - ny / inv_scale);
00731 }
00732 float gs = (maxgray - mingray) / (render_max - render_min);
00733 if (render_max < render_min) {
00734 gs = 0;
00735 rm = FLT_MAX;
00736 }
00737 int dsx = -1;
00738 int dsy = 0;
00739 if (inv_scale == floor(inv_scale)) {
00740 dsx = (int) inv_scale;
00741 dsy = (int) (inv_scale * nx);
00742 }
00743 int addi = 0;
00744 int addr = 0;
00745
00746 if (dsx == -1) {
00747 addi = (int) floor(inv_scale);
00748 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
00749 }
00750
00751 int remx = 0;
00752 int remy = 0;
00753 int xmin = 0;
00754 if (x0 < 0) {
00755 xmin = (int) (-x0 / inv_scale);
00756 xsize -= (int) floor(x0 / inv_scale);
00757 x0 = 0;
00758 }
00759
00760 if ((xsize - xmin) * inv_scale > (nx - x0)) {
00761 xsize = (int) ((nx - x0) / inv_scale + xmin);
00762 }
00763 int ymax = ysize - 1;
00764 if (y0 < 0) {
00765 ymax = (int) (ysize + y0 / inv_scale - 1);
00766 ymin += (int) floor(y0 / inv_scale);
00767 y0 = 0;
00768 }
00769
00770
00771 if (xmin < 0) {
00772 xmin = 0;
00773 }
00774
00775 if (ymin < 0) {
00776 ymin = 0;
00777 }
00778 if (xsize > ixsize) {
00779 xsize = ixsize;
00780 }
00781 if (ymax > iysize) {
00782 ymax = iysize;
00783 }
00784
00785 int lmax = nx * ny - 1;
00786 unsigned char tri[3];
00787 float* image_data = get_data();
00788 if (is_complex()) {
00789 if (dsx != -1) {
00790 int l = y0 * nx;
00791 for (int j = ymax; j >= ymin; j--) {
00792 int ll = x0;
00793 for (int i = xmin; i < xsize; i++, ll += dsx) {
00794 if (l + ll > lmax || ll >= nx - 2) {
00795 break;
00796 }
00797 int kk = 0;
00798 if (ll >= nx / 2) {
00799 if (l >= (ny - inv_scale) * nx) {
00800 kk = 2 * (ll - nx / 2) + 2;
00801 }
00802 else {
00803 kk = 2 * (ll - nx / 2) + l + 2 + nx;
00804 }
00805 }
00806 else {
00807 kk = nx * ny - (l + 2 * ll) - 2;
00808 }
00809 int k = 0;
00810 float t = image_data[kk];
00811 if (t <= rm) {
00812 k = mingray;
00813 }
00814 else if (t >= render_max) {
00815 k = maxgray;
00816 }
00817 else {
00818 k = (int) (gs * (t - render_min));
00819 k += mingray;
00820 }
00821 tri[0] = static_cast < unsigned char >(k);
00822 cmap(ref, kk, tri);
00823 data[i * 3 + j * bpl] = tri[0];
00824 data[i * 3 + 1 + j * bpl] = tri[1];
00825 data[i * 3 + 2 + j * bpl] = tri[2];
00826 }
00827 l += dsy;
00828 }
00829 }
00830 else {
00831 remy = 10;
00832 for (int j = ymax, l = y0 * nx; j >= ymin; j--) {
00833 int br = l;
00834 remx = 10;
00835 for (int i = xmin, ll = x0; i < xsize - 1; i++) {
00836 if (l + ll > lmax || ll >= nx - 2) {
00837 break;
00838 }
00839 int kk = 0;
00840 if (ll >= nx / 2) {
00841 if (l >= (ny * nx - nx)) {
00842 kk = 2 * (ll - nx / 2) + 2;
00843 }
00844 else {
00845 kk = 2 * (ll - nx / 2) + l + 2 + nx;
00846 }
00847 }
00848 else {
00849 kk = nx * ny - (l + 2 * ll) - 2;
00850 }
00851 int k = 0;
00852 float t = image_data[kk];
00853 if (t <= rm) {
00854 k = mingray;
00855 }
00856 else if (t >= render_max) {
00857 k = maxgray;
00858 }
00859 else {
00860 k = (int) (gs * (t - render_min));
00861 k += mingray;
00862 }
00863 tri[0] = static_cast < unsigned char >(k);
00864 cmap(ref, kk, tri);
00865 data[i * 3 + j * bpl] = tri[0];
00866 data[i * 3 + 1 + j * bpl] = tri[1];
00867 data[i * 3 + 2 + j * bpl] = tri[2];
00868 ll += addi;
00869 remx += addr;
00870 if (remx > scale_n) {
00871 remx -= scale_n;
00872 ll++;
00873 }
00874 }
00875 l = br + addi * nx;
00876 remy += addr;
00877 if (remy > scale_n) {
00878 remy -= scale_n;
00879 l += nx;
00880 }
00881 }
00882 }
00883 }
00884 else {
00885 if (dsx != -1) {
00886 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00887 int br = l;
00888 for (int i = xmin; i < xsize; i++, l += dsx) {
00889 if (l > lmax) {
00890 break;
00891 }
00892 float t = image_data[l];
00893 int k = 0;
00894 if (t <= rm) {
00895 k = mingray;
00896 }
00897 else if (t >= render_max) {
00898 k = maxgray;
00899 }
00900 else {
00901 k = (int) (gs * (t - render_min));
00902 k += mingray;
00903 }
00904 tri[0] = static_cast < unsigned char >(k);
00905 cmap(ref, l, tri);
00906 data[i * 3 + j * bpl] = tri[0];
00907 data[i * 3 + 1 + j * bpl] = tri[1];
00908 data[i * 3 + 2 + j * bpl] = tri[2];
00909 }
00910 l = br + dsy;
00911 }
00912 }
00913 else {
00914 remy = 10;
00915 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
00916 int br = l;
00917 remx = 10;
00918 for (int i = xmin; i < xsize; i++) {
00919 if (l > lmax) {
00920 break;
00921 }
00922 float t = image_data[l];
00923 int k = 0;
00924 if (t <= rm) {
00925 k = mingray;
00926 }
00927 else if (t >= render_max) {
00928 k = maxgray;
00929 }
00930 else {
00931 k = (int) (gs * (t - render_min));
00932 k += mingray;
00933 }
00934 tri[0] = static_cast < unsigned char >(k);
00935 cmap(ref, l, tri);
00936 data[i * 3 + j * bpl] = tri[0];
00937 data[i * 3 + 1 + j * bpl] = tri[1];
00938 data[i * 3 + 2 + j * bpl] = tri[2];
00939 l += addi;
00940 remx += addr;
00941 if (remx > scale_n) {
00942 remx -= scale_n;
00943 l++;
00944 }
00945 }
00946 l = br + addi * nx;
00947 remy += addr;
00948 if (remy > scale_n) {
00949 remy -= scale_n;
00950 l += nx;
00951 }
00952 }
00953 }
00954 }
00955
00956 EXITFUNC;
00957 }
00958
00959 void EMData::ap2ri()
00960 {
00961 ENTERFUNC;
00962
00963 if (!is_complex() || is_ri()) {
00964 return;
00965 }
00966
00967 Util::ap2ri(get_data(), (size_t)nx * ny * nz);
00968 set_ri(true);
00969 update();
00970 EXITFUNC;
00971 }
00972
00973 void EMData::ri2inten()
00974 {
00975 ENTERFUNC;
00976
00977 if (!is_complex()) return;
00978 if (!is_ri()) ap2ri();
00979
00980 float * data = get_data();
00981 size_t size = (size_t)nx * ny * nz;
00982 for (size_t i = 0; i < size; i += 2) {
00983 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
00984 data[i+1]=0;
00985 }
00986
00987 set_attr("is_intensity", int(1));
00988 update();
00989 EXITFUNC;
00990 }
00991
00992
00993 void EMData::ri2ap()
00994 {
00995 ENTERFUNC;
00996
00997 if (!is_complex() || !is_ri()) {
00998 return;
00999 }
01000
01001 float * data = get_data();
01002
01003 size_t size = (size_t)nx * ny * nz;
01004 for (size_t i = 0; i < size; i += 2) {
01005 #ifdef _WIN32
01006 float f = (float)_hypot(data[i], data[i + 1]);
01007 #else
01008 float f = (float)hypot(data[i], data[i + 1]);
01009 #endif
01010 if (data[i] == 0 && data[i + 1] == 0) {
01011 data[i + 1] = 0;
01012 }
01013 else {
01014 data[i + 1] = atan2(data[i + 1], data[i]);
01015 }
01016 data[i] = f;
01017 }
01018
01019 set_ri(false);
01020 update();
01021 EXITFUNC;
01022 }
01023
01024
01025 float calc_bessel(const int n, const float& x) {
01026 gsl_sf_result result;
01027
01028 gsl_sf_bessel_Jn_e(n,(double)x, &result);
01029 return (float)result.val;
01030 }
01031
01032 EMData* EMData::bispecRotTransInvN(int N, int NK)
01033 {
01034
01035 int EndP = this -> get_xsize();
01036 int Mid = (int) ((1+EndP)/2);
01037 int End = 2*Mid-1;
01038
01039 int CountxyMax = End*End;
01040
01041 int *SortfkInds = new int[CountxyMax];
01042 int *kVecX = new int[CountxyMax];
01043 int *kVecY = new int[CountxyMax];
01044 float *fkVecR = new float[CountxyMax];
01045 float *fkVecI = new float[CountxyMax];
01046 float *absD1fkVec = new float[CountxyMax];
01047 float *absD1fkVecSorted = new float[CountxyMax];
01048
01049 float *jxjyatan2 = new float[End*End];
01050
01051 EMData * ThisCopy = new EMData(End,End);
01052
01053 for (int jx=0; jx <End ; jx++) {
01054 for (int jy=0; jy <End ; jy++) {
01055 float ValNow = this -> get_value_at(jx,jy);
01056 ThisCopy -> set_value_at(jx,jy,ValNow);
01057
01058 }}
01059
01060
01061 EMData* fk = ThisCopy -> do_fft();
01062 fk ->process_inplace("xform.fourierorigin.tocenter");
01063
01064
01065 EMData* fkRCopy = new EMData(End,End);
01066 EMData* fkICopy = new EMData(End,End);
01067 EMData* fkCopy = new EMData(End,End);
01068
01069
01070 for (int jCount= 0; jCount<End*End; jCount++) {
01071
01072 int jx = jCount%End ;
01073 int jy = (jCount-jx)/End ;
01074 jxjyatan2[jCount] = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
01075 }
01076
01077
01078 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01079
01080 int kx = kEx/2;
01081 int kIx = kx+Mid-1;
01082 int kCx = -kx ;
01083 int kCIx = kCx+ Mid-1 ;
01084 for (int kEy= 0 ; kEy<End; kEy++) {
01085 int kIy = kEy ;
01086 int ky = kEy+1-Mid;
01087 float realVal = fk -> get_value_at(kEx ,kEy) ;
01088 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01089 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01090 float fkAng = atan2(imagVal,realVal);
01091
01092 float NewRealVal ;
01093 float NewImagVal ;
01094 float AngMatlab ;
01095
01096 if (kIx==Mid-1) {
01097
01098
01099 }
01100
01101 if (kIx>Mid-1){
01102
01103 }
01104
01105 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01106 NewRealVal = absVal*cos(AngMatlab);
01107 NewImagVal = absVal*sin(AngMatlab);
01108
01109
01110 fkVecR[kIy+kIx *End] = NewRealVal ;
01111 fkVecR[kIy+kCIx*End] = NewRealVal ;
01112 fkVecI[kIy+kIx *End] = NewImagVal ;
01113 fkVecI[kIy+kCIx*End] = -NewImagVal ;
01114 absD1fkVec[kIy + kIx *End] = absVal;
01115 absD1fkVec[kIy + kCIx *End] = absVal;
01116 kVecX[kIy+kIx *End] = kx ;
01117 kVecX[kIy+kCIx *End] = kCx ;
01118 kVecY[kIy+kIx *End] = ky ;
01119 kVecY[kIy+kCIx *End] = ky ;
01120
01121
01122
01123
01124
01125 fkCopy -> set_value_at(kIx ,kIy, absVal);
01126 fkCopy -> set_value_at(kCIx,kIy, absVal);
01127 fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
01128 fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
01129 fkICopy -> set_value_at(kIx, kIy, NewImagVal);
01130 fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
01131
01132 }
01133 }
01134 #ifdef _WIN32
01135 _unlink("fkCopy.???");
01136 _unlink("fk?Copy.???");
01137 #else
01138 system("rm -f fkCopy.???");
01139 system("rm -f fk?Copy.???");
01140 #endif //_WIN32
01141 fkCopy -> write_image("fkCopy.img");
01142 fkRCopy -> write_image("fkRCopy.img");
01143 fkICopy -> write_image("fkICopy.img");
01144
01145 cout << "Starting the sort "<< endl;
01146
01147 vector< pair<float, int> > absInds;
01148 for(int i = 0; i < CountxyMax; ++i ) {
01149 pair<float,int> p;
01150 p = make_pair(absD1fkVec[i],i);
01151 absInds.push_back( p);
01152 }
01153
01154 std::sort(absInds.begin(),absInds.end());
01155
01156 for(int i = 0; i < CountxyMax; ++i ) {
01157 pair<float,int> p ;
01158 p = absInds[i] ;
01159 absD1fkVecSorted[CountxyMax-1-i] = p.first ;
01160 SortfkInds[CountxyMax-1-i] = p.second ;
01161 }
01162
01163 cout << "Ending the sort "<< endl;
01164
01165
01166
01167
01168 for(int i = 0; i < CountxyMax; ++i ) {
01169 int Si = SortfkInds[i];
01170 int kIx = (int) Si/End; kIx = (int) i/End;
01171
01172
01173
01174
01175 }
01176 cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0] << endl;
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200 cout << "NK = " << NK << endl;
01201 float frR= 3.0/4.0;
01202 int LradRange= (int) (floor(Mid/frR)) ;
01203
01204 float *radRange = new float[LradRange];
01205 radRange[0]=0;
01206 for (int irad=1; irad < LradRange; irad++){
01207 radRange[irad] = radRange[irad-1] + frR; }
01208
01209
01210
01211
01212 cout << "Starting the calculation of invariants for N= " << N << endl;
01213
01214
01215
01216 EMData* RotTransInv = new EMData();
01217 RotTransInv -> set_size(LradRange,LradRange);
01218
01219
01220
01221
01222
01223
01224
01225 for (int jr1=0; jr1 < LradRange ; jr1++ ) {
01226 float r1= radRange[jr1];
01227
01228 for (int jr2=0; jr2<LradRange; jr2++ ) {
01229 float r2= radRange[jr2];
01230 float RotTransInvTemp=0;
01231 for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
01232 int Countkxy =SortfkInds[jCountkxy] ;
01233 int kx = kVecX[Countkxy] ;
01234 int ky = kVecY[Countkxy] ;
01235 float k2 = (float)(kx*kx+ky*ky);
01236 if (k2==0) { continue;}
01237 float phiK =0;
01238 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; phiK= atan2((float)ky,(float)kx);
01239
01240 float fkR = fkVecR[Countkxy] ;
01241 float fkI = fkVecI[Countkxy] ;
01242
01243
01244 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
01245 int Countqxy =SortfkInds[jCountqxy] ;
01246 int qx = kVecX[Countqxy] ;
01247 int qy = kVecY[Countqxy] ;
01248 int q2 = qx*qx+qy*qy;
01249 if (q2==0) {continue;}
01250 float phiQ =0;
01251 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; phiQ=atan2((float)qy,(float)qx);
01252 float fqR = fkVecR[Countqxy] ;
01253 float fqI = fkVecI[Countqxy] ;
01254 int kCx = (-kx-qx);
01255 int kCy = (-ky-qy);
01256 int kCIx = ((kCx+Mid+2*End)%End);
01257 int kCIy = ((kCy+Mid+2*End)%End);
01258 kCx = kCIx-Mid;
01259 kCy = kCIy-Mid;
01260 int CountCxy = kCIx*End+kCIy;
01261 float fCR = fkVecR[CountCxy];
01262 float fCI = fkVecI[CountCxy];
01263 if (jr1+jr2==-1) {
01264 printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
01265 printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f, kCx=%d, kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
01266 }
01267 for (int p=0; p<NK; p++){
01268
01269 if (SortfkInds[p]==CountCxy){
01270 float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
01271 float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
01272
01273
01274 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI))
01275 * cos(N*(phiK-phiQ+M_PI));
01276 bispectemp -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
01277 * sin(N*(phiK-phiQ+M_PI));
01278 float bess1 = calc_bessel(N, Arg1 );
01279 float bess2 = calc_bessel(N, Arg2 );
01280
01281
01282
01283 RotTransInvTemp = RotTransInvTemp + bispectemp * bess1*bess2 ;
01284
01285 }
01286 }
01287 }
01288 }
01289 RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp) ;
01290
01291 }
01292 }
01293
01294
01295 return RotTransInv ;
01296
01297
01298 }
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329 EMData* EMData::bispecRotTransInvDirect(int type)
01330 {
01331
01332 int EndP = this -> get_xsize();
01333 int Mid = (int) ((1+EndP)/2);
01334 int End = 2*Mid-1;
01335
01336 int CountxyMax = End*End;
01337
01338
01339 int *kVecX = new int[CountxyMax];
01340 int *kVecY = new int[CountxyMax];
01341 float *fkVecR = new float[CountxyMax];
01342 float *fkVecI = new float[CountxyMax];
01343 float *absD1fkVec = new float[CountxyMax];
01344
01345
01346
01347 float *jxjyatan2 = new float[End*End];
01348
01349
01350 EMData * ThisCopy = new EMData(End,End);
01351
01352 for (int jx=0; jx <End ; jx++) {
01353 for (int jy=0; jy <End ; jy++) {
01354 float ValNow = this -> get_value_at(jx,jy);
01355 ThisCopy -> set_value_at(jx,jy,ValNow);
01356 jxjyatan2[jy*End + jx] = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
01357
01358 }}
01359
01360
01361 EMData* fk = ThisCopy -> do_fft();
01362 fk ->process_inplace("xform.fourierorigin.tocenter");
01363
01364
01365
01366 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) {
01367
01368 int kx = kEx/2;
01369 int kIx = kx+Mid-1;
01370 int kCx = -kx ;
01371 int kCIx = kCx+ Mid-1 ;
01372 for (int kEy= 0 ; kEy<End; kEy++) {
01373 int kIy = kEy ;
01374 int ky = kEy+1-Mid;
01375 float realVal = fk -> get_value_at(kEx ,kEy) ;
01376 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
01377 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
01378 float fkAng = atan2(imagVal,realVal);
01379
01380 float NewRealVal ;
01381 float NewImagVal ;
01382 float AngMatlab ;
01383
01384 if (kIx==Mid-1) {
01385
01386 }
01387
01388 if (kIx>Mid-1){
01389
01390 }
01391
01392 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
01393 NewRealVal = absVal*cos(AngMatlab);
01394 NewImagVal = absVal*sin(AngMatlab);
01395
01396
01397 fkVecR[ kIy +kIx *End] = NewRealVal ;
01398 fkVecR[(End-1-kIy)+kCIx*End] = NewRealVal ;
01399 fkVecI[ kIy +kIx *End] = NewImagVal ;
01400 fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
01401 absD1fkVec[(End-1-kIy) + kIx *End] = absVal;
01402 absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
01403 kVecX[kIy+kIx *End] = kx ;
01404 kVecX[kIy+kCIx *End] = kCx ;
01405 kVecY[kIy+kIx *End] = ky ;
01406 kVecY[kIy+kCIx *End] = ky ;
01407
01408
01409
01410
01411
01412
01413
01414
01415 }
01416 }
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428 float frR= 3.0/4.0;
01429 frR= 1;
01430 int LradRange= (int) (1+floor(Mid/frR -.1)) ;
01431
01432 float *radRange = new float[LradRange];
01433 for (int irad=0; irad < LradRange; irad++){
01434 radRange[irad] = frR*irad;
01435
01436 }
01437
01438 cout << "Starting the calculation of invariants" << endl;
01439
01440
01441 if (type==0) {
01442 int LthetaRange = 59;
01443 float ftR = (2.0f*M_PI/LthetaRange );
01444 float *thetaRange = new float[LthetaRange];
01445
01446 for (int ith=0; ith < LthetaRange; ith++){
01447 thetaRange[ith] = ftR*ith; }
01448
01449 int TotalVol = LradRange*LradRange*LthetaRange;
01450
01451 float *RotTransInv = new float[TotalVol];
01452 float *WeightInv = new float[TotalVol];
01453
01454 for (int jW=0; jW<TotalVol; jW++) {
01455 RotTransInv[jW] = 0;
01456 WeightInv[jW] = 0;
01457 }
01458
01459 for (int jW=0; jW<TotalVol; jW++) {
01460 RotTransInv[jW] = 0;
01461 WeightInv[jW] = 0;
01462 }
01463
01464
01465
01466 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01467 int kx = kVecX[Countkxy] ;
01468 int ky = kVecY[Countkxy] ;
01469 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01470 float phiK =0;
01471 if (k2>0) phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1];
01472 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01473 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01474
01475
01476 if ((k2==0)|| (k2>Mid) ) { continue;}
01477
01478 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01479 int qx = kVecX[Countqxy] ;
01480 int qy = kVecY[Countqxy] ;
01481 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01482 if ((q2==0)|| (q2>Mid) ) {continue;}
01483 float phiQ =0;
01484 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1];
01485 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01486 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01487 int kCx = (-kx-qx);
01488 int kCy = (-ky-qy);
01489 int kCIx = ((kCx+Mid+2*End)%End);
01490 int kCIy = ((kCy+Mid+2*End)%End);
01491 kCx = ((kCIx+End-1)%End)+1-Mid;
01492 kCy = ((kCIy+End-1)%End)+1-Mid ;
01493
01494
01495 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01496 float fCR = fkVecR[CountCxy];
01497 float fCI = fkVecI[CountCxy];
01498
01499
01500
01501
01502
01503 float phiQK = (4*M_PI+phiQ-phiK);
01504 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
01505
01506
01507
01508 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01509
01510 if ((q2<k2) ) continue;
01511
01512
01513
01514
01515
01516
01517 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01518 int k2IndHi = k2IndLo;
01519 float k2Lo= radRange[k2IndLo];
01520 if (k2IndLo+1< LradRange) {
01521 k2IndHi = k2IndLo+1;
01522 }
01523
01524
01525 float kCof =k2-k2Lo;
01526
01527 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01528 int q2IndHi=q2IndLo;
01529 float q2Lo= radRange[q2IndLo];
01530 if (q2IndLo+1 < LradRange) {
01531 q2IndHi = q2IndLo+1 ;
01532 }
01533 float qCof = q2-q2Lo;
01534
01535 if ((qCof<0) || (qCof >1) ) {
01536 cout<< "Weird! qCof="<< qCof << " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
01537 int x ;
01538 cin >> x ;
01539 }
01540
01541 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
01542 int thetaIndHi = thetaIndLo;
01543
01544 float thetaLo = thetaRange[thetaIndLo];
01545 float thetaHi = thetaLo;
01546 float thetaCof = 0;
01547
01548 if (thetaIndLo+1< LthetaRange) {
01549 thetaIndHi = thetaIndLo +1;
01550 }else{
01551 thetaIndHi=0;
01552 }
01553
01554 thetaHi = thetaRange[thetaIndHi];
01555
01556 if (thetaHi==thetaLo) {
01557 thetaCof =0 ;
01558 } else {
01559 thetaCof = (phiQK-thetaLo)/(thetaHi-thetaLo);
01560 }
01561
01562 if ((thetaCof>2*M_PI) ) {
01563 cout<< "Weird! thetaCof="<< thetaCof <<endl ;
01564 thetaCof=0;
01565 }
01566
01567
01568
01569
01570
01571 for (int jk =1; jk<=2; jk++){
01572 for (int jq =1; jq<=2; jq++){
01573 for (int jtheta =1; jtheta<=2; jtheta++){
01574
01575 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
01576 * (thetaCof+(1-2*thetaCof)*(jtheta==1));
01577
01578
01579 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01580 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01581 int thetaInd = thetaIndLo*(jtheta==1) + thetaIndHi*(jtheta ==2);
01582 int TotalInd = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
01583
01584
01585
01586
01587
01588
01589 RotTransInv[TotalInd] += Weight*bispectemp;
01590 WeightInv[TotalInd] += Weight;
01591
01592 }}}
01593 }
01594 }
01595
01596 cout << "Finished Main Section " << endl;
01597
01598
01599
01600 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
01601 EMData *RotTransInvF = new EMData(LradRange,LradRange,LthetaRange);
01602 EMData *WeightImage = new EMData(LradRange,LradRange,LthetaRange);
01603
01604
01605
01606
01607
01608
01609
01610 for (int jtheta =0; jtheta < LthetaRange; jtheta++){
01611 for (int jq =0; jq<LradRange; jq++){
01612 for (int jk =0; jk<LradRange ; jk++){
01613
01614 int TotalInd = jtheta*LradRange*LradRange+jq*LradRange+jk;
01615 float Weight = WeightInv[TotalInd];
01616 WeightImage -> set_value_at(jk,jq,jtheta,Weight);
01617 RotTransInvF -> set_value_at(jk,jq,jtheta,0);
01618 if (Weight <= 0) continue;
01619 RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);
01620 int newjtheta = (LthetaRange- jtheta)%LthetaRange;
01621 RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );
01622 }
01623 }
01624 }
01625
01626 cout << " Almost Done " << endl;
01627 #ifdef _WIN32
01628 _unlink("WeightImage.???");
01629 #else
01630 system("rm -f WeightImage.???");
01631 #endif //_WIN32
01632 WeightImage -> write_image("WeightImage.img");
01633
01634 return RotTransInvF ;
01635 }
01636
01637 if (type==1) {
01638 int TotalVol = LradRange*LradRange;
01639
01640 float *RotTransInv = new float[TotalVol];
01641 float *WeightInv = new float[TotalVol];
01642
01643 for (int jW=0; jW<TotalVol; jW++) {
01644 RotTransInv[jW] = 0;
01645 WeightInv[jW] = 0;
01646 }
01647
01648
01649 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
01650 int kx = kVecX[Countkxy] ;
01651 int ky = kVecY[Countkxy] ;
01652 float k2 = ::sqrt((float)(kx*kx+ky*ky));
01653 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
01654 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
01655
01656
01657 if ((k2==0)|| (k2>Mid) ) { continue;}
01658
01659 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){
01660
01661
01662 int qx = kVecX[Countqxy] ;
01663 int qy = kVecY[Countqxy] ;
01664 float q2 = ::sqrt((float)(qx*qx+qy*qy));
01665 if ((q2==0)|| (q2>Mid) ) {continue;}
01666 if ((q2<k2) ) continue;
01667
01668 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
01669 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
01670
01671 int kCx = (-kx-qx);
01672 int kCy = (-ky-qy);
01673 int kCIx = ((kCx+Mid+2*End)%End);
01674 int kCIy = ((kCy+Mid+2*End)%End);
01675 kCx = ((kCIx+End-1)%End)+1-Mid;
01676 kCy = ((kCIy+End-1)%End)+1-Mid ;
01677
01678
01679 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
01680 float fCR = fkVecR[CountCxy];
01681 float fCI = fkVecI[CountCxy];
01682
01683
01684 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
01685
01686
01687 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
01688 int k2IndHi = k2IndLo;
01689 float k2Lo= radRange[k2IndLo];
01690 if (k2IndLo+1< LradRange) {
01691 k2IndHi = k2IndLo+1;
01692 }
01693
01694
01695 float kCof =k2-k2Lo;
01696
01697 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
01698 int q2IndHi=q2IndLo;
01699 float q2Lo= radRange[q2IndLo];
01700 if (q2IndLo+1 < LradRange) {
01701 q2IndHi = q2IndLo+1 ;
01702 }
01703 float qCof = q2-q2Lo;
01704
01705
01706 for (int jk =1; jk<=2; jk++){
01707 for (int jq =1; jq<=2; jq++){
01708
01709 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
01710
01711 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
01712 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
01713 int TotalInd = q2Ind*LradRange+k2Ind;
01714 RotTransInv[TotalInd] += Weight*bispectemp;
01715 WeightInv[TotalInd] += Weight;
01716
01717 }}
01718 }
01719 }
01720
01721
01722
01723
01724
01725 EMData *RotTransInvF = new EMData(LradRange,LradRange);
01726 EMData *WeightImage = new EMData(LradRange,LradRange);
01727
01728 for (int jk =0; jk<LradRange ; jk++){
01729 for (int jq =jk; jq<LradRange; jq++){
01730 int TotalInd = jq*LradRange+jk;
01731 int TotalIndBar = jq*LradRange+jk;
01732 float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
01733 if (Weight <=0) continue;
01734 WeightImage -> set_value_at(jk,jq,Weight);
01735 WeightImage -> set_value_at(jq,jk,Weight);
01736 #ifdef _WIN32
01737 float ValNow = pow( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight, 1.0f/3.0f ) ;
01738 #else
01739 float ValNow = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight ) ;
01740 #endif //_WIN32
01741 RotTransInvF -> set_value_at(jk,jq,ValNow);
01742 RotTransInvF -> set_value_at(jq,jk,ValNow );
01743 }}
01744
01745 #ifdef _WIN32
01746 _unlink("WeightImage.???");
01747 #else
01748 system("rm -f WeightImage.???");
01749 #endif //_WIN32
01750 WeightImage -> write_image("WeightImage.img");
01751
01752 return RotTransInvF ;
01753 }
01754 return 0;
01755 }
01756
01757
01758 void EMData::insert_clip(const EMData * const block, const IntPoint &origin) {
01759 int nx1 = block->get_xsize();
01760 int ny1 = block->get_ysize();
01761 int nz1 = block->get_zsize();
01762
01763 Region area(origin[0], origin[1], origin[2],nx1, ny1, nz1);
01764
01765
01766 int x0 = (int) area.origin[0];
01767 x0 = x0 < 0 ? 0 : x0;
01768
01769 int y0 = (int) area.origin[1];
01770 y0 = y0 < 0 ? 0 : y0;
01771
01772 int z0 = (int) area.origin[2];
01773 z0 = z0 < 0 ? 0 : z0;
01774
01775 int x1 = (int) (area.origin[0] + area.size[0]);
01776 x1 = x1 > nx ? nx : x1;
01777
01778 int y1 = (int) (area.origin[1] + area.size[1]);
01779 y1 = y1 > ny ? ny : y1;
01780
01781 int z1 = (int) (area.origin[2] + area.size[2]);
01782 z1 = z1 > nz ? nz : z1;
01783 if (z1 <= 0) {
01784 z1 = 1;
01785 }
01786
01787 int xd0 = (int) (area.origin[0] < 0 ? -area.origin[0] : 0);
01788 int yd0 = (int) (area.origin[1] < 0 ? -area.origin[1] : 0);
01789 int zd0 = (int) (area.origin[2] < 0 ? -area.origin[2] : 0);
01790
01791 if (x1 < x0 || y1 < y0 || z1 < z0) return;
01792
01793 size_t clipped_row_size = (x1-x0) * sizeof(float);
01794 size_t src_secsize = (size_t)(nx1 * ny1);
01795 size_t dst_secsize = (size_t)(nx * ny);
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818 float *src = block->get_data() + zd0 * src_secsize + yd0 * nx1 + xd0;
01819 float *dst = get_data() + z0 * dst_secsize + y0 * nx + x0;
01820
01821 size_t src_gap = src_secsize - (y1-y0) * nx1;
01822 size_t dst_gap = dst_secsize - (y1-y0) * nx;
01823
01824 for (int i = z0; i < z1; i++) {
01825 for (int j = y0; j < y1; j++) {
01826 EMUtil::em_memcpy(dst, src, clipped_row_size);
01827 src += nx1;
01828 dst += nx;
01829 }
01830 src += src_gap;
01831 dst += dst_gap;
01832 }
01833
01834 #ifdef EMAN2_USING_CUDA
01835 if(block->cudarwdata){
01836 copy_to_cuda();
01837 }
01838 #endif
01839
01840 update();
01841 EXITFUNC;
01842 }
01843
01844
01845 void EMData::insert_scaled_sum(EMData *block, const FloatPoint ¢er,
01846 float scale, float)
01847 {
01848 ENTERFUNC;
01849 float * data = get_data();
01850 if (get_ndim()==3) {
01851
01852 int xs=(int)floor(block->get_xsize()*scale/2.0);
01853 int ys=(int)floor(block->get_ysize()*scale/2.0);
01854 int zs=(int)floor(block->get_zsize()*scale/2.0);
01855 int x0=(int)center[0]-xs;
01856 int x1=(int)center[0]+xs;
01857 int y0=(int)center[1]-ys;
01858 int y1=(int)center[1]+ys;
01859 int z0=(int)center[2]-zs;
01860 int z1=(int)center[2]+zs;
01861
01862 if (x1<0||y1<0||z1<0||x0>get_xsize()||y0>get_ysize()||z0>get_zsize()) return;
01863
01864
01865 if (x0<0) x0=0;
01866 if (y0<0) y0=0;
01867 if (z0<0) z0=0;
01868 if (x1>=get_xsize()) x1=get_xsize()-1;
01869 if (y1>=get_ysize()) y1=get_ysize()-1;
01870 if (z1>=get_zsize()) z1=get_zsize()-1;
01871
01872 float bx=block->get_xsize()/2.0f;
01873 float by=block->get_ysize()/2.0f;
01874 float bz=block->get_zsize()/2.0f;
01875
01876 size_t idx;
01877 for (int x=x0; x<=x1; x++) {
01878 for (int y=y0; y<=y1; y++) {
01879 for (int z=z0; z<=z1; z++) {
01880 idx = x + y * nx + (size_t)z * nx * ny;
01881 data[idx] +=
01882 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by,(z-center[2])/scale+bz);
01883 }
01884 }
01885 }
01886 update();
01887 }
01888 else if (get_ndim()==2) {
01889
01890 int xs=(int)floor(block->get_xsize()*scale/2.0);
01891 int ys=(int)floor(block->get_ysize()*scale/2.0);
01892 int x0=(int)center[0]-xs;
01893 int x1=(int)center[0]+xs;
01894 int y0=(int)center[1]-ys;
01895 int y1=(int)center[1]+ys;
01896
01897 if (x1<0||y1<0||x0>get_xsize()||y0>get_ysize()) return;
01898
01899
01900 if (x0<0) x0=0;
01901 if (y0<0) y0=0;
01902 if (x1>=get_xsize()) x1=get_xsize()-1;
01903 if (y1>=get_ysize()) y1=get_ysize()-1;
01904
01905 float bx=block->get_xsize()/2.0f;
01906 float by=block->get_ysize()/2.0f;
01907
01908 for (int x=x0; x<=x1; x++) {
01909 for (int y=y0; y<=y1; y++) {
01910 data[x + y * nx] +=
01911 block->sget_value_at_interp((x-center[0])/scale+bx,(y-center[1])/scale+by);
01912 }
01913 }
01914 update();
01915 }
01916 else {
01917 LOGERR("insert_scaled_sum supports only 2D and 3D data");
01918 throw ImageDimensionException("2D/3D only");
01919 }
01920
01921 EXITFUNC;
01922 }
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008