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