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