#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 181 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::div(), 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(). 00181 { 00182 //std::complex<float> phase_mult; 00183 // Not only does the value of "flag" determine how we handle 00184 // periodicity, but it also determines whether or not we should 00185 // normalize the results. Here's some convenience bools: 00186 bool donorm = (0 == flag%2) ? true : false; 00187 // the 2x padding is hardcoded for now 00188 int npad = (flag >= 3) ? 2 : 1; // amount of padding used 00189 // g may be NULL. If so, have g point to the same object as f. In that 00190 // case we need to be careful later on not to try to delete g's workspace 00191 // as well as f's workspace, since they will be the same. 00192 bool gexists = true; 00193 if (!g) { g = f; gexists = false; } 00194 if ( f->is_complex() || g->is_complex() ) { 00195 // Fourier input only allowed for circulant 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 // These are actual dimensions of f (and real-space sizes for ny and nz) 00202 int nx = f->get_xsize(); 00203 int ny = f->get_ysize(); 00204 int nz = f->get_zsize(); 00205 // We manifestly assume no zero-padding here, just the 00206 // necessary extension along x for the fft 00207 if (!f->is_real()) nx = (nx - 2 + (f->is_fftodd() ? 1 : 0)); 00208 00209 // these are padded dimensions 00210 const int nxp = npad*nx; 00211 const int nyp = (ny > 1) ? npad*ny : 1; // don't pad y for 1-d image 00212 const int nzp = (nz > 1) ? npad*nz : 1; // don't pad z for 2-d image 00213 00214 // now one half of the padded, fft-extended size along x 00215 const int lsd2 = (nxp + 2 - nxp%2) / 2; 00216 00217 EMData* fp = NULL; 00218 if (f->is_complex()) { 00219 // If f is already a fourier object then fp is a copy of f. 00220 // (The fp workspace is modified, so we copy f to keep f pristine.) 00221 fp=f->copy(); 00222 } else { 00223 // [normalize] [pad] compute fft 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 // The [padded] fft-extended version of g is gp. 00231 EMData* gp = NULL; 00232 if(f==g) { 00233 // g is an alias for f, so gp should be an alias for fp 00234 gp=fp; 00235 } else if (g->is_complex()) { 00236 // g is already a Fourier object, so gp is just an alias for g 00237 // (The gp workspace is not modified, so we don't need a copy.) 00238 gp = g; 00239 } else { 00240 // normal case: g is real and different from f, so compute gp 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 // Get complex matrix views of fp and gp; matrices start from 1 (not 0) 00248 fp->set_array_offsets(1,1,1); 00249 gp->set_array_offsets(1,1,1); 00250 00251 // If the center flag is true, put the center of the correlation in the middle 00252 // If it is false, put it in (0,0), this approach saves time, but it is diffcult to manage the result 00253 if (center) { 00254 // Multiply two functions (the real work of this routine) 00255 int itmp = nx/2; 00256 //float sx = float(-twopi*float(itmp)/float(nxp)); 00257 float sxn = 2*float(itmp)/float(nxp); 00258 float sx = -M_PI*sxn; 00259 itmp = ny/2; 00260 //float sy = float(-twopi*float(itmp)/float(nyp)); 00261 float syn = 2*float(itmp)/float(nyp); 00262 float sy = -M_PI*syn; 00263 itmp = nz/2; 00264 //float sz = float(-twopi*float(itmp)/float(nzp)); 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)) ) { // padded images have always even size: if ( (padded) || (even) ) ... 00268 switch (ptype) { 00269 case AUTOCORRELATION: 00270 // fpmat := |fpmat|^2 00271 // Note nxp are padded dimensions 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 // fpmat:=|fpmat| 00284 // Note nxp are padded dimensions 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 // fpmat:=fpmat*conjg(gpmat) 00295 // Note nxp are padded dimensions 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 // fpmat:=fpmat*gpmat 00306 // Note nxp are padded dimensions 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 // fpmat := |fpmat|^2 00330 // Note nxp are padded dimensions 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 // fpmat:=|fpmat| 00346 // Note nxp are padded dimensions 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 // fpmat:=fpmat*conjg(gpmat) 00360 // Note nxp are padded dimensions 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 // fpmat:=fpmat*gpmat 00374 // Note nxp are padded dimensions 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 // If the center flag is false, then just do basic multiplication 00398 // Here I aterd the method of complex calculation. This method is much faster than the previous one. 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 //phase_mult = 1; 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 // 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); 00461 if (gexists && (f != g) && (!g->is_complex())) { 00462 if( gp ) { 00463 delete gp; 00464 gp = 0; 00465 } 00466 } 00467 // back transform 00468 fp->do_ift_inplace(); 00469 if(center && npad ==2) fp->depad(); 00470 else fp->depad_corner(); 00471 00472 //vector<int> saved_offsets = fp->get_array_offsets(); I do not know what the meaning of it was, did not work anyway PAP 00473 fp->set_array_offsets(1,1,1); 00474 00475 // Lag normalization 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 //OVER AND OUT 00489 //fp->set_array_offsets(saved_offsets); This was strange and did not work, PAP 00490 fp->set_array_offsets(0,0,0); 00491 fp->update(); 00492 return fp; 00493 }
|
|
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 // We manifestly assume no zero-padding here, just the 00047 // necessary extension along x for the fft 00048 00049 if (f->is_complex()) nx = (nx - 2 + f->is_fftodd()); // nx is the real-space size of the input image 00050 int lsd2 = (nx + 2 - nx%2) / 2; // Extended x-dimension of the complex image 00051 00052 // Process f if real 00053 EMData* fp = NULL; 00054 if(f->is_complex()) fp = f->copy(); // we need to make a full copy so that we don't damage the original 00055 else { 00056 00057 fp = f->norm_pad(false, 1); // Extend and do the FFT if f is real 00058 fp->do_fft_inplace(); 00059 00060 00061 } 00062 fp->set_array_offsets(1,1,1); 00063 00064 // Periodogram: fp:=|fp|**2 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 // Create power as a 3D array (-n/2:n/2+n%2-1) 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()); // output image 00083 power.set_size(nx, ny, nz); 00084 power.set_array_offsets(-nx2,-ny2,-nz2); 00085 //If instead of preservation of the norm one would prefer to have peak of a PW of a single sine wave equal one 00086 // multiply power by the scale below, or the other way around. 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 // Create the Friedel related half 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++) { // Note this loop begins with 1 - FFT should create correct Friedel related 0 plane 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) { //if nz even, fix the first slice 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) { //if ny even, fix the first line 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) { //if ny even, fix the first column 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; // avoid a memory leak! 00142 fp = 0; 00143 } 00144 //power[0][0][0]=power[1][0][0]; //Steve requested the original origin. 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 // set the pixel size for the power spectrum, only ration of the frequency pixel size is considered 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 //OVER AND OUT 00160 }
|