Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

fourierfilter.cpp

Go to the documentation of this file.
00001 /*
00002  * Author: Pawel A.Penczek, 09/09/2006 (Pawel.A.Penczek@uth.tmc.edu)
00003  * Copyright (c) 2000-2006 The University of Texas - Houston Medical School
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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  *DISCLAIMER
00043  * 03/29/05 P.A.Penczek
00044  * The University of Texas
00045  * Pawel.A.Penczek@uth.tmc.edu
00046  * Please do not modify the content of this document
00047  * without a written permission of the author.
00048  ******************************************************/
00049 /*
00050    Fourier filter
00051 Purpose: Apply selected Fourier space filters to 1-2-3D image.
00052 Method: Calculate FFT (if needed), apply the filter, if the input was real, do the iFFT, otherwise exit. .
00053 Real input -> padded workspace -> Reallocate workspace to real output.
00054 Complex input -> complex output.
00055 Input: f real or complex 1-2-3D image
00056 Output: 1-2-3D filtered image (real or complex).
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         // Set amount of Fourier padding
00071         // dopad should be a bool, but EMObject Dict's can't store bools.
00072         int dopad = params["dopad"];
00073         int npad;
00074         if (0 == dopad) {
00075                 // no padding
00076                 npad = 1;
00077         } else if (1 == dopad) {
00078                 // 2x padding (hard-wired)
00079                 npad = 2;
00080         } else if (2 == dopad) {
00081                 npad = 4;
00082         } else {
00083                 // invalid entry
00084                 LOGERR("The dopad parameter must be 0 (false) or 1 (true)");
00085                 return NULL; // FIXME: replace w/ exception throw
00086         }
00087 
00088         // If the input image is already a Fourier image, then we want to
00089         // have this routine return a Fourier image
00090         complex_input = fimage->is_complex();
00091         if ( complex_input && 1 == dopad ) {
00092                 // Cannot pad Fourier input image
00093                 LOGERR("Cannot pad Fourier input image");
00094                 return NULL; // FIXME: replace w/ exception throw
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                 // We manifestly assume no zero-padding here, just the
00104                 // necessary extension along x for the fft
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; // Extended x-dimension of the complex image
00112         int lsd3 = lsd2 - 1;
00113 
00114         //  Perform padding (if necessary) and fft, if the image is not already an fft image
00115         EMData* fp = NULL; // workspace image
00116         if (complex_input) {
00117                 if (doInPlace) {
00118                         // it's okay to change the original image
00119                         fp = fimage;
00120                 } else {
00121                         // fimage must remain pristine
00122                         fp = fimage->copy();
00123                 }
00124         } else {
00125                 if (doInPlace) {
00126                         if (npad>1) {
00127                                 LOGERR("Cannot pad with inplace filter");
00128                                 return NULL;    // FIXME, exception
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         //  And the filter type is:
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         // For the given type of filter set up any necessary parameters for the
00160         // filter calculation.  FIXME: Need parameter bounds checking!
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                         //origin_type = params["origin_type"];
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;} //ignore B-factor for the binary CTF
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                         }//without this bracket, compiler on water will complain about crosses initialization
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                         // for 2d, sqrt(2) = 1.414, relax a little bit to 1.6
00276                         // for 3d, sqrt(3) = 1.732, relax a little bit to 1.9
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; // FIXME: replace w/ exception throw
00283         }
00284         // Perform actual calculation
00285         //  Gaussian bandpass is the only one with center for frequencies
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:   // K-B filter
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                                                 //if (!kbptr)
00406                                                 //      throw
00407                                                 //              NullPointerException("kbptr null!");
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:   //  Sinh filter
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                                                 //if (!kbptr)
00431                                                 //      throw
00432                                                 //              NullPointerException("kbptr null!");
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                                                                 //float argu = kbptr->sinhwin((float) jx);
00443                                                                 //cout << jx<<"  "<< nux<<"  "<<argu<<endl;
00444                                                                 break;
00445                                                 }
00446                                         }
00447                                 }
00448                         }
00449                         break;
00450                 case KAISER_I0_INVERSE:   // 1./(K-B filter)
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                                         //if (!kbptr)
00461                                         //      throw
00462                                         //              NullPointerException("kbptr null!");
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:  // 1./sinh
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                                                 //if (!kbptr)
00486                                                 //      throw
00487                                                 //              NullPointerException("kbptr null!");
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                                                                 //float argu = kbptr->sinhwin((float) jx);
00498                                                                 //cout << jx<<"  "<< nux<<"  "<<argu<<endl;
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                         //if (origin_type) {
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                         /*} else {
00554                                 for ( iz = 1; iz <= nzp; iz++) {
00555                                         jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00556                                         if  (iz>nzp2) { kz=iz-nzp2; } else { kz=iz+nzp2; }
00557                                         for ( iy = 1; iy <= nyp; iy++) {
00558                                                 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00559                                                 if  (iy>nyp2) { ky=iy-nyp2; } else { ky=iy+nyp2; }
00560                                                 for ( ix = 1; ix <= lsd2; ix++) {
00561                                                         jx=ix-1;
00562                                                         fp->cmplx(ix,ky,kz) *=  exp(-float(twopi)*iimag*(xshift*jx/nx + yshift*jy/ny+ zshift*jz/nz));
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]); // (1-df)*table[ir]+df*table[ir+1];
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                 // Return a complex (Fourier) filtered image
00671                 // Note: fp and fimage are the _same_ object if doInPlace
00672                 // is true, so in that case fimage has been filtered.
00673                 // We always return an image (pointer), but if the function
00674                 // was called with doInPlace == true then the calling function
00675                 // will probably ignore the return value.
00676 
00677                 // ELSE Return a real-space filtered image
00678                 //
00679                 // On 12/15/2006 Wei Zhang comment:
00680                 // If input is reald and doInPlace == true, we might need delete fp after copy its
00681                 // data back to fimage, since fp is allocated inside this function and is ignored
00682                 // by caller if doInPlace == true. As a reminder, the caller is EMFourierFuncInPlace
00683                 //
00684         fp->set_array_offsets(0,0,0);
00685         fp->update();
00686         if (doInPlace && !complex_input) {
00687                 // copy workspace data into the original image
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 }

Generated on Tue Jul 12 13:48:58 2011 for EMAN2 by  doxygen 1.3.9.1