#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 167 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(). 00167 { 00168 int normfact; 00169 //std::complex<float> phase_mult; 00170 // Not only does the value of "flag" determine how we handle 00171 // periodicity, but it also determines whether or not we should 00172 // normalize the results. Here's some convenience bools: 00173 bool donorm = (0 == flag%2) ? true : false; 00174 // the 2x padding is hardcoded for now 00175 int npad = (flag >= 3) ? 2 : 1; // amount of padding used 00176 // g may be NULL. If so, have g point to the same object as f. In that 00177 // case we need to be careful later on not to try to delete g's workspace 00178 // as well as f's workspace, since they will be the same. 00179 bool gexists = true; 00180 if (!g) { g = f; gexists = false; } 00181 if ( f->is_complex() || g->is_complex() ) { 00182 // Fourier input only allowed for circulant 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 // These are actual dimensions of f (and real-space sizes for ny and nz) 00189 int nx = f->get_xsize(); 00190 int ny = f->get_ysize(); 00191 int nz = f->get_zsize(); 00192 // We manifestly assume no zero-padding here, just the 00193 // necessary extension along x for the fft 00194 if (!f->is_real()) nx = (nx - 2 + (f->is_fftodd() ? 1 : 0)); 00195 00196 // these are padded dimensions 00197 const int nxp = npad*nx; 00198 const int nyp = (ny > 1) ? npad*ny : 1; // don't pad y for 1-d image 00199 const int nzp = (nz > 1) ? npad*nz : 1; // don't pad z for 2-d image 00200 00201 // now one half of the padded, fft-extended size along x 00202 const int lsd2 = (nxp + 2 - nxp%2) / 2; 00203 // The [padded] fft-extended fourier version of f is fp. 00204 00205 EMData* fp = NULL; 00206 if (f->is_complex()) { 00207 // If f is already a fourier object then fp is a copy of f. 00208 // (The fp workspace is modified, so we copy f to keep f pristine.) 00209 fp=f->copy(); 00210 } else { 00211 // [normalize] [pad] compute fft 00212 fp = f->norm_pad(donorm, npad); 00213 fp->do_fft_inplace(); 00214 } 00215 // The [padded] fft-extended version of g is gp. 00216 EMData* gp = NULL; 00217 if(f==g) { 00218 // g is an alias for f, so gp should be an alias for fp 00219 gp=fp; 00220 } else if (g->is_complex()) { 00221 // g is already a Fourier object, so gp is just an alias for g 00222 // (The gp workspace is not modified, so we don't need a copy.) 00223 gp = g; 00224 } else { 00225 // normal case: g is real and different from f, so compute gp 00226 gp = g->norm_pad(donorm, npad); 00227 gp->do_fft_inplace(); 00228 } 00229 // Get complex matrix views of fp and gp; matrices start from 1 (not 0) 00230 fp->set_array_offsets(1,1,1); 00231 gp->set_array_offsets(1,1,1); 00232 00233 // If the center flag is true, put the center of the correlation in the middle 00234 // If it is false, put it in (0,0), this approach saves time, but it is diffcult to manage the result 00235 if (center) { 00236 // Multiply two functions (the real work of this routine) 00237 int itmp = nx/2; 00238 //float sx = float(-twopi*float(itmp)/float(nxp)); 00239 float sxn = 2*float(itmp)/float(nxp); 00240 float sx = -M_PI*sxn; 00241 itmp = ny/2; 00242 //float sy = float(-twopi*float(itmp)/float(nyp)); 00243 float syn = 2*float(itmp)/float(nyp); 00244 float sy = -M_PI*syn; 00245 itmp = nz/2; 00246 //float sz = float(-twopi*float(itmp)/float(nzp)); 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 // fpmat := |fpmat|^2 00253 // Note nxp are padded dimensions 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 // fpmat:=|fpmat| 00266 // Note nxp are padded dimensions 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 // fpmat:=fpmat*conjg(gpmat) 00277 // Note nxp are padded dimensions 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 // fpmat:=fpmat*gpmat 00288 // Note nxp are padded dimensions 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 // fpmat := |fpmat|^2 00312 // Note nxp are padded dimensions 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 // fpmat:=|fpmat| 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 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 // fpmat:=fpmat*conjg(gpmat) 00342 // Note nxp are padded dimensions 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 // fpmat:=fpmat*gpmat 00356 // Note nxp are padded dimensions 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 // If the center flag is false, then just do basic multiplication 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 //phase_mult = 1; 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 // 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); 00442 if (gexists && (f != g) && (!g->is_complex())) { 00443 if( gp ) { 00444 delete gp; 00445 gp = 0; 00446 } 00447 } 00448 // back transform 00449 fp->do_ift_inplace(); 00450 if(center && npad ==2) fp->depad(); 00451 else fp->depad_corner(); 00452 00453 //vector<int> saved_offsets = fp->get_array_offsets(); I do not know what the meaning of it was, did not work anyway PAP 00454 fp->set_array_offsets(1,1,1); 00455 00456 normfact = (nxp/nx)*(nyp/ny)*(nzp/nz); // Normalization factor for the padded operations 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 // Lag normalization 00461 if(flag>4) { 00462 normfact = nx*ny*nz; // Normalization factor 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 //OVER AND OUT 00475 //fp->set_array_offsets(saved_offsets); This was strange and did not work, PAP 00476 fp->set_array_offsets(0,0,0); 00477 fp->update(); 00478 return fp; 00479 }
|
|
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_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 // We manifestly assume no zero-padding here, just the 00046 // necessary extension along x for the fft 00047 00048 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd()); // nx is the real-space size of the input image 00049 int lsd2 = (nx + 2 - nx%2) / 2; // Extended x-dimension of the complex image 00050 00051 // Process f if real 00052 EMData* fp = NULL; 00053 if(f->is_complex()) fp = f->copy(); // we need to make a full copy so that we don't damage the original 00054 else { 00055 fp = f->norm_pad(false, 1); // Extend and do the FFT if f is real 00056 fp->do_fft_inplace(); 00057 } 00058 fp->set_array_offsets(1,1,1); 00059 00060 // Periodogram: fp:=|fp|**2 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 // Create power as a 3D array (-n/2:n/2+n%2-1) 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()); // output image 00079 power.set_size(nx, ny, nz); 00080 power.set_array_offsets(-nx2,-ny2,-nz2); 00081 //If instead of preservation of the norm one would prefer to have peak of a PW of a single sine wave equal one 00082 // multiply power by the scale below, or the other way around. 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 // Create the Friedel related half 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++) { // Note this loop begins with 1 - FFT should create correct Friedel related 0 plane 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) { //if nz even, fix the first slice 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) { //if ny even, fix the first line 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) { //if ny even, fix the first column 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; // avoid a memory leak! 00138 fp = 0; 00139 } 00140 //power[0][0][0]=power[1][0][0]; //Steve requested the original origin. 00141 00142 power.update(); 00143 power.set_array_offsets(0,0,0); 00144 return &power; 00145 //OVER AND OUT 00146 }
|