#include "emdata.h"Include dependency graph for fundamentals.cpp:

Go to the source code of this file.
Namespaces | |
| namespace | EMAN |
Functions | |
| EMData * | periodogram (EMData *f) |
| EMData * | fourierproduct (EMData *f, EMData *g, fp_flag flag, fp_type ptype, bool center) |
| Fourier product of two images. | |
|
||||||||||||||||||||||||
|
Fourier product of two images.
Definition at line 184 of file fundamentals.cpp. References abs, EMAN::AUTOCORRELATION, EMAN::CIRCULANT, EMAN::EMData::cmplx(), EMAN::CONVOLUTION, EMAN::EMData::copy(), EMAN::CORRELATION, EMAN::EMData::depad(), EMAN::EMData::depad_corner(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::do_ift_inplace(), flag, EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), imag(), InvalidValueException, EMAN::EMData::is_complex(), EMAN::EMData::is_fftodd(), EMAN::EMData::is_real(), LOGERR, EMAN::EMData::norm_pad(), nx, ny, real(), EMAN::SELF_CORRELATION, EMAN::EMData::set_array_offsets(), and EMAN::EMData::update(). Referenced by EMAN::autocorrelation(), EMAN::convolution(), EMAN::correlation(), EMAN::ConvolutionProcessor::process_inplace(), and EMAN::self_correlation(). 00184 {
00185 size_t normfact;
00186 //std::complex<float> phase_mult;
00187 // Not only does the value of "flag" determine how we handle
00188 // periodicity, but it also determines whether or not we should
00189 // normalize the results. Here's some convenience bools:
00190 bool donorm = (0 == flag%2) ? true : false;
00191 // the 2x padding is hardcoded for now
00192 int npad = (flag >= 3) ? 2 : 1; // amount of padding used
00193 // g may be NULL. If so, have g point to the same object as f. In that
00194 // case we need to be careful later on not to try to delete g's workspace
00195 // as well as f's workspace, since they will be the same.
00196 bool gexists = true;
00197 if (!g) { g = f; gexists = false; }
00198 if ( f->is_complex() || g->is_complex() ) {
00199 // Fourier input only allowed for circulant
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 // These are actual dimensions of f (and real-space sizes for ny and nz)
00206 int nx = f->get_xsize();
00207 int ny = f->get_ysize();
00208 int nz = f->get_zsize();
00209 // We manifestly assume no zero-padding here, just the
00210 // necessary extension along x for the fft
00211 if (!f->is_real()) nx = (nx - 2 + (f->is_fftodd() ? 1 : 0));
00212
00213 // these are padded dimensions
00214 const int nxp = npad*nx;
00215 const int nyp = (ny > 1) ? npad*ny : 1; // don't pad y for 1-d image
00216 const int nzp = (nz > 1) ? npad*nz : 1; // don't pad z for 2-d image
00217
00218 // now one half of the padded, fft-extended size along x
00219 const int lsd2 = (nxp + 2 - nxp%2) / 2;
00220
00221 EMData* fp = NULL;
00222 if (f->is_complex()) {
00223 // If f is already a fourier object then fp is a copy of f.
00224 // (The fp workspace is modified, so we copy f to keep f pristine.)
00225 fp=f->copy();
00226 } else {
00227 // [normalize] [pad] compute fft
00228 fp = f->norm_pad(donorm, npad);
00229 fp->do_fft_inplace();
00230 }
00231 // The [padded] fft-extended version of g is gp.
00232 EMData* gp = NULL;
00233 if(f==g) {
00234 // g is an alias for f, so gp should be an alias for fp
00235 gp=fp;
00236 } else if (g->is_complex()) {
00237 // g is already a Fourier object, so gp is just an alias for g
00238 // (The gp workspace is not modified, so we don't need a copy.)
00239 gp = g;
00240 } else {
00241 // normal case: g is real and different from f, so compute gp
00242 gp = g->norm_pad(donorm, npad);
00243 gp->do_fft_inplace();
00244 }
00245 // Get complex matrix views of fp and gp; matrices start from 1 (not 0)
00246 fp->set_array_offsets(1,1,1);
00247 gp->set_array_offsets(1,1,1);
00248
00249 // If the center flag is true, put the center of the correlation in the middle
00250 // If it is false, put it in (0,0), this approach saves time, but it is diffcult to manage the result
00251 if (center) {
00252 // Multiply two functions (the real work of this routine)
00253 int itmp = nx/2;
00254 //float sx = float(-twopi*float(itmp)/float(nxp));
00255 float sxn = 2*float(itmp)/float(nxp);
00256 float sx = -M_PI*sxn;
00257 itmp = ny/2;
00258 //float sy = float(-twopi*float(itmp)/float(nyp));
00259 float syn = 2*float(itmp)/float(nyp);
00260 float sy = -M_PI*syn;
00261 itmp = nz/2;
00262 //float sz = float(-twopi*float(itmp)/float(nzp));
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 // fpmat := |fpmat|^2
00269 // Note nxp are padded dimensions
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 // fpmat:=|fpmat|
00282 // Note nxp are padded dimensions
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 // fpmat:=fpmat*conjg(gpmat)
00293 // Note nxp are padded dimensions
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 // fpmat:=fpmat*gpmat
00304 // Note nxp are padded dimensions
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 // fpmat := |fpmat|^2
00328 // Note nxp are padded dimensions
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 // fpmat:=|fpmat|
00344 // Note nxp are padded dimensions
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 // fpmat:=fpmat*conjg(gpmat)
00358 // Note nxp are padded dimensions
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 // fpmat:=fpmat*gpmat
00372 // Note nxp are padded dimensions
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 // If the center flag is false, then just do basic multiplication
00396 // Here I aterd the method of complex calculation. This method is much faster than the previous one.
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 //phase_mult = 1;
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 // Now done w/ gp, so let's get rid of it (if it's not an alias of fp or simply g was complex on input);
00459 if (gexists && (f != g) && (!g->is_complex())) {
00460 if( gp ) {
00461 delete gp;
00462 gp = 0;
00463 }
00464 }
00465 // back transform
00466 fp->do_ift_inplace();
00467 if(center && npad ==2) fp->depad();
00468 else fp->depad_corner();
00469
00470 //vector<int> saved_offsets = fp->get_array_offsets(); I do not know what the meaning of it was, did not work anyway PAP
00471 fp->set_array_offsets(1,1,1);
00472
00473 normfact = (size_t)(nxp/nx)*(nyp/ny)*(nzp/nz); // Normalization factor for the padded operations
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 // Lag normalization
00478 if(flag>4) {
00479 normfact = (size_t)nx*ny*nz; // Normalization factor
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 //OVER AND OUT
00492 //fp->set_array_offsets(saved_offsets); This was strange and did not work, PAP
00493 fp->set_array_offsets(0,0,0);
00494 fp->update();
00495 return fp;
00496 }
|
|
|
Definition at line 40 of file fundamentals.cpp. References EMAN::EMData::cmplx(), EMAN::EMData::copy(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), imag(), EMAN::EMData::is_complex(), EMAN::EMData::is_fftodd(), EMAN::EMData::norm_pad(), nx, ny, power(), real(), EMAN::EMData::set_array_offsets(), EMAN::EMData::set_attr(), EMAN::EMData::set_size(), and EMAN::EMData::update(). 00040 {
00041 // These are actual dimensions
00042 int nx = f->get_xsize();
00043 int ny = f->get_ysize();
00044 int nz = f->get_zsize();
00045
00046
00047 // We manifestly assume no zero-padding here, just the
00048 // necessary extension along x for the fft
00049
00050 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd()); // nx is the real-space size of the input image
00051 int lsd2 = (nx + 2 - nx%2) / 2; // Extended x-dimension of the complex image
00052
00053 // Process f if real
00054 EMData* fp = NULL;
00055 if(f->is_complex()) fp = f->copy(); // we need to make a full copy so that we don't damage the original
00056 else {
00057
00058 fp = f->norm_pad(false, 1); // Extend and do the FFT if f is real
00059 fp->do_fft_inplace();
00060
00061
00062
00063
00064 }
00065 fp->set_array_offsets(1,1,1);
00066
00067 // Periodogram: fp:=|fp|**2
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 // Create power as a 3D array (-n/2:n/2+n%2-1)
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()); // output image
00086 power.set_size(nx, ny, nz);
00087 power.set_array_offsets(-nx2,-ny2,-nz2);
00088 //If instead of preservation of the norm one would prefer to have peak of a PW of a single sine wave equal one
00089 // multiply power by the scale below, or the other way around.
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 // Create the Friedel related half
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++) { // Note this loop begins with 1 - FFT should create correct Friedel related 0 plane
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) { //if nz even, fix the first slice
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) { //if ny even, fix the first line
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) { //if ny even, fix the first column
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; // avoid a memory leak!
00145 fp = 0;
00146 }
00147 //power[0][0][0]=power[1][0][0]; //Steve requested the original origin.
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 // set the pixel size for the power spectrum, only ration of the frequency pixel size is considered
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 //OVER AND OUT
00163 }
|
1.3.9.1