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 #include "processor.h"
00034 #include <algorithm>
00035 #include <cstdlib>
00036
00037 using namespace EMAN;
00038 using namespace std;
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 EMData* Processor::EMFourierFilterFunc(EMData * fimage, Dict params, bool doInPlace)
00059 {
00060 int nx, ny, nz, nyp2, nzp2, ix, iy, iz, jx, jy, jz;
00061 float dx, dy, dz, omega=0, omegaL=0, omegaH=0;
00062 float center=0, gamma=0, argx, argy, argz;
00063 float aa, eps, ord=0, cnst=0, aL, aH, cnstL=0, cnstH=0;
00064 bool complex_input;
00065 vector<float> table;
00066 int undoctf=0;
00067 float voltage=100.0f, ak=0.0f, cs=2.0f, ps=1.0f, b_factor=0.0f, wgh=0.1f, sign=-1.0f;
00068 if (!fimage) return NULL;
00069 const int ndim = fimage->get_ndim();
00070
00071
00072 int dopad = params["dopad"];
00073 int npad;
00074 if (0 == dopad) {
00075
00076 npad = 1;
00077 } else if (1 == dopad) {
00078
00079 npad = 2;
00080 } else if (2 == dopad) {
00081 npad = 4;
00082 } else {
00083
00084 LOGERR("The dopad parameter must be 0 (false) or 1 (true)");
00085 return NULL;
00086 }
00087
00088
00089
00090 complex_input = fimage->is_complex();
00091 if ( complex_input && 1 == dopad ) {
00092
00093 LOGERR("Cannot pad Fourier input image");
00094 return NULL;
00095 }
00096
00097 Util::KaiserBessel* kbptr = 0;
00098
00099
00100 nx = fimage->get_xsize();
00101 ny = fimage->get_ysize();
00102 nz = fimage->get_zsize();
00103
00104
00105 if (complex_input) nx = (nx - 2 + fimage->is_fftodd());
00106
00107 const int nxp = npad*nx;
00108 const int nyp = (ny > 1) ? npad*ny : 1;
00109 const int nzp = (nz > 1) ? npad*nz : 1;
00110
00111 int lsd2 = (nxp + 2 - nxp%2) / 2;
00112 int lsd3 = lsd2 - 1;
00113
00114
00115 EMData* fp = NULL;
00116 if (complex_input) {
00117 if (doInPlace) {
00118
00119 fp = fimage;
00120 } else {
00121
00122 fp = fimage->copy();
00123 }
00124 } else {
00125 if (doInPlace) {
00126 if (npad>1) {
00127 LOGERR("Cannot pad with inplace filter");
00128 return NULL;
00129 }
00130 fp=fimage;
00131 fp->do_fft_inplace();
00132 } else {
00133 fp = fimage->norm_pad( false, npad, 1);
00134 fp->do_fft_inplace();
00135 }
00136 }
00137 fp->set_array_offsets(1,1,1);
00138
00139
00140 int filter_type = params["filter_type"];
00141
00142 nyp2 = nyp/2; nzp2 = nzp/2;
00143 dx = 1.0f/float(nxp);
00144 #ifdef _WIN32
00145 dy = 1.0f/_cpp_max(float(nyp),1.0f);
00146 dz = 1.0f/_cpp_max(float(nzp),1.0f);
00147 #else
00148 dy = 1.0f/std::max(float(nyp),1.0f);
00149 dz = 1.0f/std::max(float(nzp),1.0f);
00150 #endif //_WIN32
00151 float dx2 = dx*dx, dy2 = dy*dy, dz2 = dz*dz;
00152
00153 vector<float>::size_type tsize;
00154 float sz[3];
00155 float szmax;
00156 vector<float>::size_type maxsize;
00157 float xshift=0.0, yshift=0.0, zshift=0.0;
00158
00159
00160
00161 switch (filter_type) {
00162 case TOP_HAT_LOW_PASS:
00163 case TOP_HAT_HIGH_PASS:
00164 omega = params["cutoff_abs"];
00165 omega = 1.0f/omega/omega;
00166 break;
00167 case TOP_HAT_BAND_PASS:
00168 omegaL = params["low_cutoff_frequency"];
00169 omegaH = params["high_cutoff_frequency"];
00170 omegaL = 1.0f/omegaL/omegaL;
00171 omegaH = 1.0f/omegaH/omegaH;
00172 break;
00173 case TOP_HOMOMORPHIC:
00174 omegaL = params["low_cutoff_frequency"];
00175 omegaH = params["high_cutoff_frequency"];
00176 gamma = params["value_at_zero_frequency"];
00177 omegaL = 1.0f/omegaL/omegaL;
00178 omegaH = 1.0f/omegaH/omegaH;
00179 break;
00180 case GAUSS_LOW_PASS:
00181 case GAUSS_HIGH_PASS:
00182 case GAUSS_INVERSE:
00183 omega = params["cutoff_abs"];
00184 omega = 0.5f/omega/omega;
00185 break;
00186 case GAUSS_BAND_PASS:
00187 omega = params["cutoff_abs"];
00188 center = params["center"];
00189 omega = 0.5f/omega/omega;
00190 break;
00191 case GAUSS_HOMOMORPHIC:
00192 omega = params["cutoff_abs"];
00193 gamma = params["value_at_zero_frequency"];
00194 omega = 0.5f/omega/omega;
00195 gamma = 1.0f-gamma;
00196 break;
00197 case BUTTERWORTH_LOW_PASS:
00198 case BUTTERWORTH_HIGH_PASS:
00199 omegaL = params["low_cutoff_frequency"];
00200 omegaH = params["high_cutoff_frequency"];
00201 eps = 0.882f;
00202 aa = 10.624f;
00203 ord = 2.0f*log10(eps/sqrt(aa*aa-1.0f))/log10(omegaL/omegaH);
00204 omegaL = omegaL/pow(eps,2.0f/ord);
00205 break;
00206 case BUTTERWORTH_HOMOMORPHIC:
00207 omegaL = params["low_cutoff_frequency"];
00208 omegaH = params["high_cutoff_frequency"];
00209 gamma = params["value_at_zero_frequency"];
00210 eps = 0.882f;
00211 aa = 10.624f;
00212 ord = 2.0f*log10(eps/sqrt(pow(aa,2)-1.0f))/log10(omegaL/omegaH);
00213 omegaL = omegaL/pow(eps,2.0f/ord);
00214 gamma = 1.0f-gamma;
00215 break;
00216 case SHIFT:
00217 xshift = params["x_shift"];
00218 yshift = params["y_shift"];
00219 zshift = params["z_shift"];
00220
00221 break;
00222 case TANH_LOW_PASS:
00223 case TANH_HIGH_PASS:
00224 omega = params["cutoff_abs"];
00225 aa = params["fall_off"];
00226 cnst = float(pihalf/aa/omega);
00227 break;
00228 case TANH_HOMOMORPHIC:
00229 omega = params["cutoff_abs"];
00230 aa = params["fall_off"];
00231 gamma = params["value_at_zero_frequency"];
00232 cnst = float(pihalf/aa/omega);
00233 gamma=1.0f-gamma;
00234 break;
00235 case TANH_BAND_PASS:
00236 omegaL = params["low_cutoff_frequency"];
00237 aL = params["Low_fall_off"];
00238 omegaH = params["high_cutoff_frequency"];
00239 aH = params["high_fall_off"];
00240 cnstL = float(pihalf/aL/(omegaH-omegaL));
00241 cnstH = float(pihalf/aH/(omegaH-omegaL));
00242 break;
00243 case CTF_:
00244 dz = params["defocus"];
00245 cs = params["Cs"];
00246 voltage = params["voltage"];
00247 ps = params["Pixel_size"];
00248 b_factor = params["B_factor"];
00249 wgh = params["amp_contrast"];
00250 sign = params["sign"];
00251 undoctf = params["undo"];
00252 ix = params["binary"];
00253 if(ix == 1) {undoctf = 2; b_factor=0.0;}
00254 break;
00255 case KAISER_I0:
00256 case KAISER_SINH:
00257 case KAISER_I0_INVERSE:
00258 case KAISER_SINH_INVERSE:
00259 {
00260 float alpha = params["alpha"];
00261 int K = params["K"];
00262 float r = params["r"];
00263 float v = params["v"];
00264 int N = params["N"];
00265 kbptr = new Util::KaiserBessel(alpha, K, r, v, N);
00266 break;
00267 }
00268 case RADIAL_TABLE:
00269 table = params["table"];
00270 tsize = table.size();
00271 sz[0] = static_cast<float>(lsd2);
00272 sz[1] = static_cast<float>(nyp2);
00273 sz[2] = static_cast<float>(nzp2);
00274 szmax = *max_element(&sz[0],&sz[3]);
00275
00276
00277 if (nzp > 1) {maxsize = vector<float>::size_type(1.9*szmax);} else {maxsize = vector<float>::size_type(1.6*szmax);}
00278 for (vector<float>::size_type i = tsize+1; i < maxsize; i++) table.push_back(0.f);
00279 break;
00280 default:
00281 LOGERR("Unknown Fourier Filter type");
00282 return NULL;
00283 }
00284
00285
00286 switch (filter_type) {
00287 case GAUSS_BAND_PASS:
00288 for ( iz = 1; iz <= nzp; iz++) {
00289 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00290 for ( iy = 1; iy <= nyp; iy++) {
00291 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00292 for ( ix = 1; ix <= lsd2; ix++) {
00293 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00294 fp->cmplx(ix,iy,iz) *= exp(-pow(sqrt(argx)-center,2)*omega);
00295 }
00296 }
00297 }
00298 break;
00299 case TOP_HAT_LOW_PASS:
00300 for ( iz = 1; iz <= nzp; iz++) {
00301 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00302 for ( iy = 1; iy <= nyp; iy++) {
00303 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00304 for ( ix = 1; ix <= lsd2; ix++) {
00305 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00306 if (argx*omega>1.0f) fp->cmplx(ix,iy,iz) = 0;
00307 }
00308 }
00309 }
00310 break;
00311 case TOP_HAT_HIGH_PASS:
00312 for ( iz = 1; iz <= nzp; iz++) {
00313 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00314 for ( iy = 1; iy <= nyp; iy++) {
00315 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00316 for ( ix = 1; ix <= lsd2; ix++) {
00317 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00318 if (argx*omega<=1.0f) fp->cmplx(ix,iy,iz) = 0;
00319 }
00320 }
00321 } break;
00322 case TOP_HAT_BAND_PASS:
00323 for ( iz = 1; iz <= nzp; iz++) {
00324 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00325 for ( iy = 1; iy <= nyp; iy++) {
00326 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00327 for ( ix = 1; ix <= lsd2; ix++) {
00328 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00329 if (argx*omegaL<1.0f || argx*omegaH>=1.0f) fp->cmplx(ix,iy,iz) = 0;
00330 }
00331 }
00332 }
00333 break;
00334 case TOP_HOMOMORPHIC:
00335 for ( iz = 1; iz <= nzp; iz++) {
00336 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00337 for ( iy = 1; iy <= nyp; iy++) {
00338 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00339 for ( ix = 1; ix <= lsd2; ix++) {
00340 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00341 if (argx*omegaH>1.0f) fp->cmplx(ix,iy,iz) = 0.0f;
00342 else if (argx*omegaL<=1.0f) fp->cmplx(ix,iy,iz) *= gamma;
00343 }
00344 }
00345 }
00346 break;
00347 case GAUSS_LOW_PASS :
00348 for ( iz = 1; iz <= nzp; iz++) {
00349 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00350 for ( iy = 1; iy <= nyp; iy++) {
00351 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00352 for ( ix = 1; ix <= lsd2; ix++) {
00353 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00354 fp->cmplx(ix,iy,iz) *= exp(-argx*omega);
00355 }
00356 }
00357 }
00358 break;
00359 case GAUSS_HIGH_PASS:
00360 for ( iz = 1; iz <= nzp; iz++) {
00361 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00362 for ( iy = 1; iy <= nyp; iy++) {
00363 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00364 for ( ix = 1; ix <= lsd2; ix++) {
00365 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00366 fp->cmplx(ix,iy,iz) *= 1.0f-exp(-argx*omega);
00367 }
00368 }
00369 }
00370 break;
00371 case GAUSS_HOMOMORPHIC:
00372 for ( iz = 1; iz <= nzp; iz++) {
00373 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00374 for ( iy = 1; iy <= nyp; iy++) {
00375 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00376 for ( ix = 1; ix <= lsd2; ix++) {
00377 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00378 fp->cmplx(ix,iy,iz) *= 1.0f-gamma*exp(-argx*omega);
00379 }
00380 }
00381 }
00382 break;
00383 case GAUSS_INVERSE :
00384 for ( iz = 1; iz <= nzp; iz++) {
00385 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00386 for ( iy = 1; iy <= nyp; iy++) {
00387 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00388 for ( ix = 1; ix <= lsd2; ix++) {
00389 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00390 fp->cmplx(ix,iy,iz) *= exp(argx*omega);
00391 }
00392 }
00393 }
00394 break;
00395 case KAISER_I0:
00396 for ( iz = 1; iz <= nzp; iz++) {
00397 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00398 float nuz = jz*dz;
00399 for ( iy = 1; iy <= nyp; iy++) {
00400 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00401 float nuy = jy*dy;
00402 for ( ix = 1; ix <= lsd2; ix++) {
00403 jx=ix-1;
00404 float nux = jx*dx;
00405
00406
00407
00408 switch (ndim) {
00409 case 3:
00410 fp->cmplx(ix,iy,iz) *= kbptr->i0win(nux)*kbptr->i0win(nuy)*kbptr->i0win(nuz);
00411 break;
00412 case 2:
00413 fp->cmplx(ix,iy,iz) *= kbptr->i0win(nux)*kbptr->i0win(nuy);
00414 break;
00415 case 1:
00416 fp->cmplx(ix,iy,iz)*= kbptr->i0win(nux);
00417 break;
00418 }
00419 }
00420 }
00421 }
00422 break;
00423 case KAISER_SINH:
00424 for ( iz = 1; iz <= nzp; iz++) {
00425 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00426 for ( iy = 1; iy <= nyp; iy++) {
00427 jy=iy-1; if(jy>nyp2) jy=jy-nyp;
00428 for ( ix = 1; ix <= lsd2; ix++) {
00429 jx=ix-1;
00430
00431
00432
00433 switch (ndim) {
00434 case 3:
00435 fp->cmplx(ix,iy,iz)*= kbptr->sinhwin((float)jx)*kbptr->sinhwin((float)jy)*kbptr->sinhwin((float)jz);
00436 break;
00437 case 2:
00438 fp->cmplx(ix,iy,iz)*= kbptr->sinhwin((float)jx)*kbptr->sinhwin((float)jy);
00439 break;
00440 case 1:
00441 fp->cmplx(ix,iy,iz)*= kbptr->sinhwin((float)jx);
00442
00443
00444 break;
00445 }
00446 }
00447 }
00448 }
00449 break;
00450 case KAISER_I0_INVERSE:
00451 for ( iz = 1; iz <= nzp; iz++) {
00452 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00453 float nuz = jz*dz;
00454 for ( iy = 1; iy <= nyp; iy++) {
00455 jy=iy-1; if(jy>nyp2) jy=jy-nyp;
00456 float nuy = jy*dy;
00457 for ( ix = 1; ix <= lsd2; ix++) {
00458 jx=ix-1;
00459 float nux = jx*dx;
00460
00461
00462
00463 switch (ndim) {
00464 case 3:
00465 fp->cmplx(ix,iy,iz) /= (kbptr->i0win(nux)*kbptr->i0win(nuy)*kbptr->i0win(nuz));
00466 break;
00467 case 2:
00468 fp->cmplx(ix,iy,iz) /= (kbptr->i0win(nux)*kbptr->i0win(nuy));
00469 break;
00470 case 1:
00471 fp->cmplx(ix,iy,iz) /= kbptr->i0win(nux);
00472 break;
00473 }
00474 }
00475 }
00476 }
00477 break;
00478 case KAISER_SINH_INVERSE:
00479 for ( iz = 1; iz <= nzp; iz++) {
00480 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00481 for ( iy = 1; iy <= nyp; iy++) {
00482 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00483 for ( ix = 1; ix <= lsd2; ix++) {
00484 jx=ix-1;
00485
00486
00487
00488 switch (ndim) {
00489 case 3:
00490 fp->cmplx(ix,iy,iz) /= (kbptr->sinhwin((float)jx)*kbptr->sinhwin((float)jy)*kbptr->sinhwin((float)jz));
00491 break;
00492 case 2:
00493 fp->cmplx(ix,iy,iz) /= (kbptr->sinhwin((float)jx)*kbptr->sinhwin((float)jy));
00494 break;
00495 case 1:
00496 fp->cmplx(ix,iy,iz) /= kbptr->sinhwin((float)jx);
00497
00498
00499 break;
00500 }
00501 }
00502 }
00503 }
00504 break;
00505 case BUTTERWORTH_LOW_PASS:
00506 for ( iz = 1; iz <= nzp; iz++) {
00507 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00508 for ( iy = 1; iy <= nyp; iy++) {
00509 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00510 for ( ix = 1; ix <= lsd2; ix++) {
00511 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00512 fp->cmplx(ix,iy,iz) *= sqrt(1.0f/(1.0f+pow(sqrt(argx)/omegaL,ord)));
00513 }
00514 }
00515 }
00516 break;
00517 case BUTTERWORTH_HIGH_PASS:
00518 for ( iz = 1; iz <= nzp; iz++) {
00519 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00520 for ( iy = 1; iy <= nyp; iy++) {
00521 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00522 for ( ix = 1; ix <= lsd2; ix++) {
00523 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00524 fp->cmplx(ix,iy,iz) *= 1.0f-sqrt(1.0f/(1.0f+pow(sqrt(argx)/omegaL,ord)));
00525 }
00526 }
00527 }
00528 break;
00529 case BUTTERWORTH_HOMOMORPHIC:
00530 for ( iz = 1; iz <= nzp; iz++) {
00531 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00532 for ( iy = 1; iy <= nyp; iy++) {
00533 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00534 for ( ix = 1; ix <= lsd2; ix++) {
00535 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00536 fp->cmplx(ix,iy,iz) *= 1.0f-gamma*sqrt(1.0f/(1.0f+pow(sqrt(argx)/omegaL,ord)));
00537 }
00538 }
00539 }
00540 break;
00541 case SHIFT:
00542
00543 for ( iz = 1; iz <= nzp; iz++) {
00544 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00545 for ( iy = 1; iy <= nyp; iy++) {
00546 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00547 for ( ix = 1; ix <= lsd2; ix++) {
00548 jx=ix-1;
00549 fp->cmplx(ix,iy,iz) *= exp(-float(twopi)*iimag*(xshift*jx/nx + yshift*jy/ny+ zshift*jz/nz));
00550 }
00551 }
00552 }
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567 break;
00568 case TANH_LOW_PASS:
00569 for ( iz = 1; iz <= nzp; iz++) {
00570 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00571 for ( iy = 1; iy <= nyp; iy++) {
00572 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00573 for ( ix = 1; ix <= lsd2; ix++) {
00574 jx=ix-1; argx = sqrt(argy + float(jx*jx)*dx2);
00575 fp->cmplx(ix,iy,iz) *= 0.5f*(tanh(cnst*(argx+omega))-tanh(cnst*(argx-omega)));
00576 }
00577 }
00578 }
00579 break;
00580 case TANH_HIGH_PASS:
00581 for ( iz = 1; iz <= nzp; iz++) {
00582 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00583 for ( iy = 1; iy <= nyp; iy++) {
00584 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00585 for ( ix = 1; ix <= lsd2; ix++) {
00586 jx=ix-1; sqrt(argx = argy + float(jx*jx)*dx2);
00587 fp->cmplx(ix,iy,iz) *= 1.0f-0.5f*(tanh(cnst*(argx+omega))-tanh(cnst*(argx-omega)));
00588 }
00589 }
00590 }
00591 break;
00592 case TANH_HOMOMORPHIC:
00593 for ( iz = 1; iz <= nzp; iz++) {
00594 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00595 for ( iy = 1; iy <= nyp; iy++) {
00596 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00597 for ( ix = 1; ix <= lsd2; ix++) {
00598 jx=ix-1; argx = sqrt(argy + float(jx*jx)*dx2);
00599 fp->cmplx(ix,iy,iz) *= 1.0f-gamma*0.5f*(tanh(cnst*(argx+omega))-tanh(cnst*(argx-omega)));
00600 }
00601 }
00602 }
00603 break;
00604 case TANH_BAND_PASS:
00605 for ( iz = 1; iz <= nzp; iz++) {
00606 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00607 for ( iy = 1; iy <= nyp; iy++) {
00608 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00609 for ( ix = 1; ix <= lsd2; ix++) {
00610 jx=ix-1; argx = sqrt(argy + float(jx*jx)*dx2);
00611 fp->cmplx(ix,iy,iz) *= 0.5f*(tanh(cnstH*(argx+omegaH))-tanh(cnstH*(argx-omegaH))-tanh(cnstL*(argx+omegaL))+tanh(cnstL*(argx-omegaL)));
00612 }
00613 }
00614 }
00615 break;
00616 case RADIAL_TABLE:
00617 for ( iz = 1; iz <= nzp; iz++) {
00618 jz=iz-1; if (jz>nzp2) jz=jz-nzp; argz = float(jz*jz)*dz2;
00619 for ( iy = 1; iy <= nyp; iy++) {
00620 jy=iy-1; if (jy>nyp2) jy=jy-nyp; argy = argz + float(jy*jy)*dy2;
00621 for ( ix = 1; ix <= lsd2; ix++) {
00622 jx=ix-1; argx = argy + float(jx*jx)*dx2;
00623 float rf = sqrt( argx )*nxp;
00624 int ir = int(rf);
00625 float df = rf - float(ir);
00626 float f = table[ir] + df * (table[ir+1] - table[ir]);
00627 fp->cmplx(ix,iy,iz) *= f;
00628 }
00629 }
00630 }
00631 break;
00632 case CTF_:
00633 for ( iz = 1; iz <= nzp; iz++) {
00634 jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00635 for ( iy = 1; iy <= nyp; iy++) {
00636 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00637 for ( ix = 1; ix <= lsd2; ix++) {
00638 jx=ix-1;
00639 if(ny>1 && nz<=1 ) ak=sqrt(static_cast<float>(jx)/lsd3*static_cast<float>(jx)/lsd3 +
00640 static_cast<float>(jy)/nyp2*static_cast<float>(jy)/nyp2)/ps/2.0f;
00641 else if(ny<=1) ak=sqrt(static_cast<float>(jx)/lsd3*static_cast<float>(jx)/lsd3)/ps/2.0f;
00642 else if(nz>1) ak=sqrt(static_cast<float>(jx)/lsd3*static_cast<float>(jx)/lsd3 +
00643 static_cast<float>(jy)/nyp2*static_cast<float>(jy)/nyp2 +
00644 static_cast<float>(jz)/nzp2*static_cast<float>(jz)/nzp2)/ps/2.0f;
00645 float tf=Util::tf(dz, ak, voltage, cs, wgh, b_factor, sign);
00646 switch (undoctf) {
00647 case 0:
00648 fp->cmplx(ix,iy,iz) *= tf;
00649 break;
00650 case 1:
00651 if( tf>0 && tf < 1e-5 ) tf = 1e-5f;
00652 if( tf<0 && tf > -1e-5 ) tf = -1e-5f;
00653 fp->cmplx(ix,iy,iz) /= tf;
00654 break;
00655 case 2:
00656 if(tf < 0.0f) fp->cmplx(ix,iy,iz) *= -1.0f;
00657 break;
00658 }
00659 }
00660 }
00661 }
00662 break;
00663 }
00664 delete kbptr; kbptr = 0;
00665 if (!complex_input) {
00666 fp->do_ift_inplace();
00667 fp->depad();
00668 }
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 fp->set_array_offsets(0,0,0);
00685 fp->update();
00686 if (doInPlace && !complex_input) {
00687
00688 float* orig = fimage->get_data();
00689 float* work = fp->get_data();
00690 for (size_t i = 0; i < (size_t)nx*ny*nz; i++) orig[i] = work[i];
00691 fimage->update();
00692 }
00693 EXITFUNC;
00694 return fp;
00695 }