EMAN2
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, dza = 0.0f, azz = 0.0f;
00068         float  tf=0.0f;
00069         if (!fimage)  return NULL;
00070         const int ndim = fimage->get_ndim();
00071         // Set amount of Fourier padding
00072         // dopad should be a bool, but EMObject Dict's can't store bools.
00073         int dopad = params["dopad"];
00074         int npad;
00075         if (0 == dopad) {
00076                 // no padding
00077                 npad = 1;
00078         } else if (1 == dopad) {
00079                 // 2x padding (hard-wired)
00080                 npad = 2;
00081         } else if (2 == dopad) {
00082         npad = 4;
00083         } else {
00084                 // invalid entry
00085                 LOGERR("The dopad parameter must be 0 (false) or 1 (true)");
00086                 return NULL; // FIXME: replace w/ exception throw
00087         }
00088 
00089         // If the input image is already a Fourier image, then we want to
00090         // have this routine return a Fourier image
00091         complex_input = fimage->is_complex();
00092         if ( complex_input && 1 == dopad ) {
00093                 // Cannot pad Fourier input image
00094                 LOGERR("Cannot pad Fourier input image");
00095                 return NULL; // FIXME: replace w/ exception throw
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                 // We manifestly assume no zero-padding here, just the
00105                 // necessary extension along x for the fft
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; // Extended x-dimension of the complex image
00113         int lsd3 = lsd2 - 1;
00114 
00115         //  Perform padding (if necessary) and fft, if the image is not already an fft image
00116         EMData* fp = NULL; // workspace image
00117         if (complex_input) {
00118                 if (doInPlace) {
00119                         // it's okay to change the original image
00120                         fp = fimage;
00121                 } else {
00122                         // fimage must remain pristine
00123                         fp = fimage->copy();
00124                 }
00125         } else {
00126                 if (doInPlace) {
00127                         if (npad>1) {
00128                                 LOGERR("Cannot pad with inplace filter");
00129                                 return NULL;    // FIXME, exception
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         //  And the filter type is:
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         // For the given type of filter set up any necessary parameters for the
00161         // filter calculation.  FIXME: Need parameter bounds checking!
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                         //origin_type = params["origin_type"];
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;} //ignore B-factor for the binary CTF
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                         }//without this bracket, compiler on water will complain about crosses initialization
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                         // for 2d, sqrt(2) = 1.414, relax a little bit to 1.6
00279                         // for 3d, sqrt(3) = 1.732, relax a little bit to 1.9
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; // FIXME: replace w/ exception throw
00286         }
00287         // Perform actual calculation
00288         //  Gaussian bandpass is the only one with center for frequencies
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.0f;
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.0f;
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.0f;
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:   // K-B filter
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                                                 //if (!kbptr)
00409                                                 //      throw
00410                                                 //              NullPointerException("kbptr null!");
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:   //  Sinh filter
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                                                 //if (!kbptr)
00434                                                 //      throw
00435                                                 //              NullPointerException("kbptr null!");
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                                                                 //float argu = kbptr->sinhwin((float) jx);
00446                                                                 //cout << jx<<"  "<< nux<<"  "<<argu<<endl;
00447                                                                 break;
00448                                                 }
00449                                         }
00450                                 }
00451                         }
00452                         break;
00453                 case KAISER_I0_INVERSE:   // 1./(K-B filter)
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                                         //if (!kbptr)
00464                                         //      throw
00465                                         //              NullPointerException("kbptr null!");
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:  // 1./sinh
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                                                 //if (!kbptr)
00489                                                 //      throw
00490                                                 //              NullPointerException("kbptr null!");
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                                                                 //float argu = kbptr->sinhwin((float) jx);
00501                                                                 //cout << jx<<"  "<< nux<<"  "<<argu<<endl;
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                         //if (origin_type) {
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                         /*} else {
00557                                 for ( iz = 1; iz <= nzp; iz++) {
00558                                         jz=iz-1; if (jz>nzp2) jz=jz-nzp;
00559                                         if  (iz>nzp2) { kz=iz-nzp2; } else { kz=iz+nzp2; }
00560                                         for ( iy = 1; iy <= nyp; iy++) {
00561                                                 jy=iy-1; if (jy>nyp2) jy=jy-nyp;
00562                                                 if  (iy>nyp2) { ky=iy-nyp2; } else { ky=iy+nyp2; }
00563                                                 for ( ix = 1; ix <= lsd2; ix++) {
00564                                                         jx=ix-1;
00565                                                         fp->cmplx(ix,ky,kz) *=  exp(-float(twopi)*iimag*(xshift*jx/nx + yshift*jy/ny+ zshift*jz/nz));
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]); // (1-df)*table[ir]+df*table[ir+1];
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                                                         //  astigmatism makes sense only on 2D
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                 // Return a complex (Fourier) filtered image
00686                 // Note: fp and fimage are the _same_ object if doInPlace
00687                 // is true, so in that case fimage has been filtered.
00688                 // We always return an image (pointer), but if the function
00689                 // was called with doInPlace == true then the calling function
00690                 // will probably ignore the return value.
00691 
00692                 // ELSE Return a real-space filtered image
00693                 //
00694                 // On 12/15/2006 Wei Zhang comment:
00695                 // If input is reald and doInPlace == true, we might need delete fp after copy its
00696                 // data back to fimage, since fp is allocated inside this function and is ignored
00697                 // by caller if doInPlace == true. As a reminder, the caller is EMFourierFuncInPlace
00698                 //
00699         fp->set_array_offsets(0,0,0);
00700         fp->update();
00701         if (doInPlace && !complex_input) {
00702                 // copy workspace data into the original image
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 }