#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 }
|