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