00001
00002
00003
00004
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 #include "emdata.h"
00033
00034 using namespace EMAN;
00035 using std::complex;
00036
00037 namespace EMAN {
00038
00039 EMData*
00040 periodogram(EMData* f) {
00041
00042 int nx = f->get_xsize();
00043 int ny = f->get_ysize();
00044 int nz = f->get_zsize();
00045
00046
00047
00048
00049
00050 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd());
00051 int lsd2 = (nx + 2 - nx%2) / 2;
00052
00053
00054 EMData* fp = NULL;
00055 if(f->is_complex()) fp = f->copy();
00056 else {
00057
00058 fp = f->norm_pad(false, 1);
00059 fp->do_fft_inplace();
00060
00061
00062
00063
00064 }
00065 fp->set_array_offsets(1,1,1);
00066
00067
00068 for (int iz = 1; iz <= nz; iz++) {
00069 for (int iy = 1; iy <= ny; iy++) {
00070 for (int ix = 1; ix <= lsd2; ix++) {
00071 float fpr = real(fp->cmplx(ix,iy,iz));
00072 float fpi = imag(fp->cmplx(ix,iy,iz));
00073 fp->cmplx(ix,iy,iz) = fpr*fpr + fpi*fpi;
00074 }
00075 }
00076 }
00077
00078 int nyt, nzt;
00079 int nx2 = nx/2;
00080 int ny2 = ny/2; if(ny2 == 0) nyt =0; else nyt=ny;
00081 int nz2 = nz/2; if(nz2 == 0) nzt =0; else nzt=nz;
00082 int nx2p = nx2+nx%2;
00083 int ny2p = ny2+ny%2;
00084 int nz2p = nz2+nz%2;
00085 EMData& power = *(new EMData());
00086 power.set_size(nx, ny, nz);
00087 power.set_array_offsets(-nx2,-ny2,-nz2);
00088
00089
00090 float scale = 4.0f/float (nx*nx)/float (ny*ny)/float (nz*nz);
00091 for (int iz = 1; iz <= nz; iz++) {
00092 int jz=iz-1;
00093 if(jz>=nz2p) jz=jz-nzt;
00094 for (int iy = 1; iy <= ny; iy++) {
00095 int jy=iy-1;
00096 if(jy>=ny2p) jy=jy-nyt;
00097 for (int ix = 1; ix <= lsd2; ix++) {
00098 int jx=ix-1;
00099 if(jx>=nx2p) jx=jx-nx;
00100 power(jx,jy,jz) = real(fp->cmplx(ix,iy,iz)) * scale;
00101 }
00102 }
00103 }
00104
00105
00106 int nzb, nze, nyb, nye, nxb, nxe;
00107 nxb =-nx2+(nx+1)%2;
00108 nxe = nx2-(nx+1)%2;
00109 if(ny2 == 0) {nyb =0; nye = 0;} else {nyb =-ny2+(ny+1)%2; nye = ny2-(ny+1)%2;}
00110 if(nz2 == 0) {nzb =0; nze = 0;} else {nzb =-nz2+(nz+1)%2; nze = nz2-(nz+1)%2;}
00111 for (int iz = nzb; iz <= nze; iz++) {
00112 for (int iy = nyb; iy <= nye; iy++) {
00113 for (int ix = 1; ix <= nxe; ix++) {
00114 power(-ix,-iy,-iz) = power(ix,iy,iz);
00115 }
00116 }
00117 }
00118 if(ny2 != 0) {
00119 if(nz2 != 0) {
00120 if(nz%2 == 0) {
00121 for (int iy = nyb; iy <= nye; iy++) {
00122 for (int ix = nxb; ix <= -1; ix++) {
00123 power(ix,iy,-nz2) = power(-ix,-iy,-nz2);
00124 }
00125 }
00126 if(ny%2 == 0) {
00127 for (int ix = nxb; ix <= -1; ix++) {
00128 power(ix,-ny2,-nz2) = power(-ix,-ny2,-nz2);
00129 }
00130 }
00131 }
00132 }
00133 if(ny%2 == 0) {
00134 for (int iz = nzb; iz <= nze; iz++) {
00135 for (int ix = nxb; ix <= -1; ix++) {
00136 power(ix,-ny2,-iz) = power(-ix,-ny2,iz);
00137 }
00138 }
00139 }
00140
00141 }
00142
00143 if( fp ) {
00144 delete fp;
00145 fp = 0;
00146 }
00147
00148
00149 int sz[3];
00150 sz[0] = nx;
00151 sz[1] = ny;
00152 sz[2] = nz;
00153 int max_size = *std::max_element(&sz[0],&sz[3]);
00154
00155 power.set_attr("apix_x", float(max_size)/nx);
00156 if(ny2 > 0) power.set_attr("apix_y", float(max_size)/ny);
00157 if(nz2 > 0) power.set_attr("apix_z", float(max_size)/nz);
00158
00159 power.update();
00160 power.set_array_offsets(0,0,0);
00161 return &power;
00162
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 EMData*
00184 fourierproduct(EMData* f, EMData* g, fp_flag flag, fp_type ptype, bool center) {
00185 size_t normfact;
00186
00187
00188
00189
00190 bool donorm = (0 == flag%2) ? true : false;
00191
00192 int npad = (flag >= 3) ? 2 : 1;
00193
00194
00195
00196 bool gexists = true;
00197 if (!g) { g = f; gexists = false; }
00198 if ( f->is_complex() || g->is_complex() ) {
00199
00200 if (CIRCULANT != flag) {
00201 LOGERR("Cannot perform normalization or padding on Fourier type.");
00202 throw InvalidValueException(flag, "Cannot perform normalization or padding on Fourier type.");
00203 }
00204 }
00205
00206 int nx = f->get_xsize();
00207 int ny = f->get_ysize();
00208 int nz = f->get_zsize();
00209
00210
00211 if (!f->is_real()) nx = (nx - 2 + (f->is_fftodd() ? 1 : 0));
00212
00213
00214 const int nxp = npad*nx;
00215 const int nyp = (ny > 1) ? npad*ny : 1;
00216 const int nzp = (nz > 1) ? npad*nz : 1;
00217
00218
00219 const int lsd2 = (nxp + 2 - nxp%2) / 2;
00220
00221 EMData* fp = NULL;
00222 if (f->is_complex()) {
00223
00224
00225 fp=f->copy();
00226 } else {
00227
00228 fp = f->norm_pad(donorm, npad);
00229 fp->do_fft_inplace();
00230 }
00231
00232 EMData* gp = NULL;
00233 if(f==g) {
00234
00235 gp=fp;
00236 } else if (g->is_complex()) {
00237
00238
00239 gp = g;
00240 } else {
00241
00242 gp = g->norm_pad(donorm, npad);
00243 gp->do_fft_inplace();
00244 }
00245
00246 fp->set_array_offsets(1,1,1);
00247 gp->set_array_offsets(1,1,1);
00248
00249
00250
00251 if (center) {
00252
00253 int itmp = nx/2;
00254
00255 float sxn = 2*float(itmp)/float(nxp);
00256 float sx = -M_PI*sxn;
00257 itmp = ny/2;
00258
00259 float syn = 2*float(itmp)/float(nyp);
00260 float sy = -M_PI*syn;
00261 itmp = nz/2;
00262
00263 float szn = 2*float(itmp)/float(nzp);
00264 float sz = -M_PI*szn;
00265 if ( nx%2==0 && (ny%2==0 || ny==1 ) && (nz%2==0 || nz==1 ) ) {
00266 switch (ptype) {
00267 case AUTOCORRELATION:
00268
00269
00270 for (int iz = 1; iz <= nzp; iz++) {
00271 for (int iy = 1; iy <= nyp; iy++) {
00272 for (int ix = 1; ix <= lsd2; ix++) {
00273 float fpr = real(fp->cmplx(ix,iy,iz));
00274 float fpi = imag(fp->cmplx(ix,iy,iz));
00275 fp->cmplx(ix,iy,iz) = complex<float>(fpr*fpr+fpi*fpi, 0.0f);
00276 }
00277 }
00278 }
00279 break;
00280 case SELF_CORRELATION:
00281
00282
00283 for (int iz = 1; iz <= nzp; iz++) {
00284 for (int iy = 1; iy <= nyp; iy++) {
00285 for (int ix = 1; ix <= lsd2; ix++) {
00286 fp->cmplx(ix,iy,iz) = complex<float>(abs(fp->cmplx(ix,iy,iz)), 0.0f);
00287 }
00288 }
00289 }
00290 break;
00291 case CORRELATION:
00292
00293
00294 for (int iz = 1; iz <= nzp; iz++) {
00295 for (int iy = 1; iy <= nyp; iy++) {
00296 for (int ix = 1; ix <= lsd2; ix++) {
00297 fp->cmplx(ix,iy,iz) *= conj(gp->cmplx(ix,iy,iz));
00298 }
00299 }
00300 }
00301 break;
00302 case CONVOLUTION:
00303
00304
00305 for (int iz = 1; iz <= nzp; iz++) {
00306 for (int iy = 1; iy <= nyp; iy++) {
00307 for (int ix = 1; ix <= lsd2; ix++) {
00308 fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz);
00309 }
00310 }
00311 }
00312 break;
00313 default:
00314 LOGERR("Illegal option in Fourier Product");
00315 throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00316 }
00317 for (int iz = 1; iz <= nzp; iz++) {
00318 for (int iy = 1; iy <= nyp; iy++) {
00319 for (int ix = (iz+iy+1)%2+1; ix <= lsd2; ix+=2) {
00320 fp->cmplx(ix,iy,iz) = -fp->cmplx(ix,iy,iz);
00321 }
00322 }
00323 }
00324 } else {
00325 switch (ptype) {
00326 case AUTOCORRELATION:
00327
00328
00329 for (int iz = 1; iz <= nzp; iz++) {
00330 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00331 for (int iy = 1; iy <= nyp; iy++) {
00332 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00333 for (int ix = 1; ix <= lsd2; ix++) {
00334 int jx=ix-1; float arg=sx*jx+argy;
00335 float fpr = real(fp->cmplx(ix,iy,iz));
00336 float fpi = imag(fp->cmplx(ix,iy,iz));
00337 fp->cmplx(ix,iy,iz)= (fpr*fpr + fpi*fpi) *std::complex<float>(cos(arg),sin(arg));
00338 }
00339 }
00340 }
00341 break;
00342 case SELF_CORRELATION:
00343
00344
00345 for (int iz = 1; iz <= nzp; iz++) {
00346 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00347 for (int iy = 1; iy <= nyp; iy++) {
00348 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00349 for (int ix = 1; ix <= lsd2; ix++) {
00350 int jx=ix-1; float arg=sx*jx+argy;
00351 fp->cmplx(ix,iy,iz) = abs(fp->cmplx(ix,iy,iz)) *std::complex<float>(cos(arg),sin(arg));
00352 }
00353 }
00354 }
00355 break;
00356 case CORRELATION:
00357
00358
00359 for (int iz = 1; iz <= nzp; iz++) {
00360 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00361 for (int iy = 1; iy <= nyp; iy++) {
00362 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00363 for (int ix = 1; ix <= lsd2; ix++) {
00364 int jx=ix-1; float arg=sx*jx+argy;
00365 fp->cmplx(ix,iy,iz) *= conj(gp->cmplx(ix,iy,iz)) *std::complex<float>(cos(arg),sin(arg));
00366 }
00367 }
00368 }
00369 break;
00370 case CONVOLUTION:
00371
00372
00373 if(npad == 1) {
00374 sx -= 4*(nx%2)/float(nx);
00375 sy -= 4*(ny%2)/float(ny);
00376 sz -= 4*(nz%2)/float(nz);
00377 }
00378 for (int iz = 1; iz <= nzp; iz++) {
00379 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00380 for (int iy = 1; iy <= nyp; iy++) {
00381 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00382 for (int ix = 1; ix <= lsd2; ix++) {
00383 int jx=ix-1; float arg=sx*jx+argy;
00384 fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz) *std::complex<float>(cos(arg),sin(arg));
00385 }
00386 }
00387 }
00388 break;
00389 default:
00390 LOGERR("Illegal option in Fourier Product");
00391 throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00392 }
00393 }
00394 } else {
00395
00396
00397 switch (ptype) {
00398 case AUTOCORRELATION:
00399 for (int iz = 1; iz <= nzp; iz++) {
00400 for (int iy = 1; iy <= nyp; iy++) {
00401 for (int ix = 1; ix <= lsd2; ix++) {
00402 float fpr = real(fp->cmplx(ix,iy,iz));
00403 float fpi = imag(fp->cmplx(ix,iy,iz));
00404 fp->cmplx(ix,iy,iz) = complex<float>(fpr*fpr+fpi*fpi, 0.0f);
00405 }
00406 }
00407 }
00408 break;
00409 case SELF_CORRELATION:
00410 for (int iz = 1; iz <= nzp; iz++) {
00411 for (int iy = 1; iy <= nyp; iy++) {
00412 for (int ix = 1; ix <= lsd2; ix++) {
00413 fp->cmplx(ix,iy,iz) = complex<float>(abs(fp->cmplx(ix,iy,iz)), 0.0f);
00414 }
00415 }
00416 }
00417 break;
00418 case CORRELATION:
00419
00420 for (int iz = 1; iz <= nzp; iz++) {
00421 for (int iy = 1; iy <= nyp; iy++) {
00422 for (int ix = 1; ix <= lsd2; ix++) {
00423 fp->cmplx(ix,iy,iz)*= conj(gp->cmplx(ix,iy,iz));
00424 }
00425 }
00426 }
00427 break;
00428 case CONVOLUTION:
00429 if(npad == 1) {
00430 float sx = -M_PI*2*(nx%2)/float(nx);
00431 float sy = -M_PI*2*(ny%2)/float(ny);
00432 float sz = -M_PI*2*(nz%2)/float(nz);
00433 for (int iz = 1; iz <= nzp; iz++) {
00434 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00435 for (int iy = 1; iy <= nyp; iy++) {
00436 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00437 for (int ix = 1; ix <= lsd2; ix++) {
00438 int jx=ix-1; float arg=sx*jx+argy;
00439 fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz) *std::complex<float>(cos(arg),sin(arg));
00440 }
00441 }
00442 }
00443 } else {
00444 for (int iz = 1; iz <= nzp; iz++) {
00445 for (int iy = 1; iy <= nyp; iy++) {
00446 for (int ix = 1; ix <= lsd2; ix++) {
00447 fp->cmplx(ix,iy,iz)*= gp->cmplx(ix,iy,iz);
00448 }
00449 }
00450 }
00451 }
00452 break;
00453 default:
00454 LOGERR("Illegal option in Fourier Product");
00455 throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00456 }
00457 }
00458
00459 if (gexists && (f != g) && (!g->is_complex())) {
00460 if( gp ) {
00461 delete gp;
00462 gp = 0;
00463 }
00464 }
00465
00466 fp->do_ift_inplace();
00467 if(center && npad ==2) fp->depad();
00468 else fp->depad_corner();
00469
00470
00471 fp->set_array_offsets(1,1,1);
00472
00473 normfact = (size_t)(nxp/nx)*(nyp/ny)*(nzp/nz);
00474 if(normfact>1) {
00475 for (int iz = 1; iz <= nz; iz++) for (int iy = 1; iy <= ny; iy++) for (int ix = 1; ix <= nx; ix++) (*fp)(ix,iy,iz) *= normfact;
00476 }
00477
00478 if(flag>4) {
00479 normfact = (size_t)nx*ny*nz;
00480 int nxc=nx/2+1, nyc=ny/2+1, nzc=nz/2+1;
00481 for (int iz = 1; iz <= nz; iz++) {
00482 float lagz=float(normfact/(nz-abs(iz-nzc)));
00483 for (int iy = 1; iy <= ny; iy++) {
00484 float lagyz=lagz/(ny-abs(iy-nyc));
00485 for (int ix = 1; ix <= nx; ix++) {
00486 (*fp)(ix,iy,iz) *= lagyz/(nx-abs(ix-nxc));
00487 }
00488 }
00489 }
00490 }
00491
00492
00493 fp->set_array_offsets(0,0,0);
00494 fp->update();
00495 return fp;
00496 }
00497 }
00498