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

fundamentals.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 
00034 using namespace EMAN;
00035 using std::complex;
00036 
00037 namespace EMAN {
00038 
00039 EMData* 
00040 periodogram(EMData* f) {
00041         // These are actual dimensions
00042         int nx  = f->get_xsize();
00043         int ny  = f->get_ysize();
00044         int nz  = f->get_zsize();
00045         
00046         
00047                 // We manifestly assume no zero-padding here, just the 
00048                 // necessary extension along x for the fft
00049 
00050         if (f->is_complex()) nx = (nx - 2 + f->is_fftodd()); // nx is the real-space size of the input image
00051         int lsd2 = (nx + 2 - nx%2) / 2; // Extended x-dimension of the complex image
00052 
00053 //  Process f if real
00054         EMData* fp = NULL;
00055         if(f->is_complex()) fp = f->copy(); // we need to make a full copy so that we don't damage the original
00056         else {
00057                 
00058                 fp = f->norm_pad(false, 1); // Extend and do the FFT if f is real
00059                 fp->do_fft_inplace();
00060 
00061                 
00062                 
00063                 
00064         }
00065         fp->set_array_offsets(1,1,1);
00066 
00067         //  Periodogram: fp:=|fp|**2
00068         for (int iz = 1; iz <= nz; iz++) {
00069                 for (int iy = 1; iy <= ny; iy++) {
00070                         for (int ix = 1; ix <= lsd2; ix++) {
00071                                 float fpr = real(fp->cmplx(ix,iy,iz));
00072                                 float fpi = imag(fp->cmplx(ix,iy,iz));
00073                                 fp->cmplx(ix,iy,iz) = fpr*fpr + fpi*fpi;
00074                         }
00075                 }
00076         }
00077         //  Create power as a 3D array (-n/2:n/2+n%2-1)
00078         int nyt, nzt;
00079         int nx2 = nx/2;
00080         int ny2 = ny/2; if(ny2 == 0) nyt =0; else nyt=ny;
00081         int nz2 = nz/2; if(nz2 == 0) nzt =0; else nzt=nz;
00082         int nx2p = nx2+nx%2;
00083         int ny2p = ny2+ny%2;
00084         int nz2p = nz2+nz%2;
00085         EMData& power = *(new EMData()); // output image
00086         power.set_size(nx, ny, nz);
00087         power.set_array_offsets(-nx2,-ny2,-nz2);
00088 //If instead of preservation of the norm one would prefer to have peak of a PW of a single sine wave equal one
00089 //                             multiply power by the scale below, or the other way around.
00090         float scale = 4.0f/float (nx*nx)/float (ny*ny)/float (nz*nz);
00091         for (int iz = 1; iz <= nz; iz++) {
00092                 int jz=iz-1; 
00093                 if(jz>=nz2p) jz=jz-nzt;
00094                 for (int iy = 1; iy <= ny; iy++) {
00095                         int jy=iy-1; 
00096                         if(jy>=ny2p) jy=jy-nyt;
00097                         for (int ix = 1; ix <= lsd2; ix++) {
00098                                 int jx=ix-1; 
00099                                 if(jx>=nx2p) jx=jx-nx;
00100                                 power(jx,jy,jz) = real(fp->cmplx(ix,iy,iz)) * scale;
00101                         }
00102                 }
00103         }
00104 
00105 //  Create the Friedel related half
00106         int  nzb, nze, nyb, nye, nxb, nxe;
00107         nxb =-nx2+(nx+1)%2;
00108         nxe = nx2-(nx+1)%2;
00109         if(ny2 == 0) {nyb =0; nye = 0;} else {nyb =-ny2+(ny+1)%2; nye = ny2-(ny+1)%2;}
00110         if(nz2 == 0) {nzb =0; nze = 0;} else {nzb =-nz2+(nz+1)%2; nze = nz2-(nz+1)%2;}
00111         for (int iz = nzb; iz <= nze; iz++) {
00112                 for (int iy = nyb; iy <= nye; iy++) {
00113                         for (int ix = 1; ix <= nxe; ix++) { // Note this loop begins with 1 - FFT should create correct Friedel related 0 plane
00114                                 power(-ix,-iy,-iz) = power(ix,iy,iz);
00115                         }
00116                 }
00117         }
00118         if(ny2 != 0)  {
00119                 if(nz2 != 0)  {
00120                         if(nz%2 == 0) {  //if nz even, fix the first slice
00121                                 for (int iy = nyb; iy <= nye; iy++) {
00122                                         for (int ix = nxb; ix <= -1; ix++) { 
00123                                                 power(ix,iy,-nz2) = power(-ix,-iy,-nz2);
00124                                         }
00125                                 }
00126                                 if(ny%2 == 0) {  //if ny even, fix the first line
00127                                         for (int ix = nxb; ix <= -1; ix++) { 
00128                                                 power(ix,-ny2,-nz2) = power(-ix,-ny2,-nz2);
00129                                         }
00130                                 }
00131                         }
00132                 }
00133                 if(ny%2 == 0) {  //if ny even, fix the first column
00134                         for (int iz = nzb; iz <= nze; iz++) {
00135                                 for (int ix = nxb; ix <= -1; ix++) {
00136                                         power(ix,-ny2,-iz) = power(-ix,-ny2,iz);
00137                                 }
00138                         }
00139                 }
00140                 
00141         }
00142                 
00143         if( fp ) {
00144                 delete fp; // avoid a memory leak!
00145                 fp = 0;
00146         }
00147         //power[0][0][0]=power[1][0][0];  //Steve requested the original origin.
00148         
00149         int sz[3];
00150         sz[0] = nx;
00151         sz[1] = ny;
00152         sz[2] = nz;
00153         int max_size = *std::max_element(&sz[0],&sz[3]);
00154         // set the pixel size for the power spectrum, only ration of the frequency pixel size is considered     
00155         power.set_attr("apix_x", float(max_size)/nx);
00156         if(ny2 > 0) power.set_attr("apix_y", float(max_size)/ny);
00157         if(nz2 > 0) power.set_attr("apix_z", float(max_size)/nz);
00158         
00159         power.update();
00160         power.set_array_offsets(0,0,0);
00161         return &power;
00162 //OVER AND OUT
00163 }
00164 
00165 
00166 /*
00167 fourierproduct
00168 Purpose: Calculate various fundamental image processing operations 
00169          based on Fourier products for 1-2-3D images. 
00170 Method: Calculate FFT (if needed),
00171         Real f input -> padded fp workspace
00172         Complex f input -> set fp to f
00173         Real g input -> padded gp workspace
00174         Complex g input -> set gp to g
00175         No g input -> set gp tp fp
00176          -> real output -> if f real delete workspace fp and if g real delete workspace gp.
00177 Input: f real or complex 1-2-3D image, g - real or complex 1-2-3D image.
00178 Output: 1-2-3D real image with the result
00179         (correlation F*conjg(G), convolution F*G, 
00180                 autocorrelation |F|^2, self-correlation |F|)
00181         and with the origin at the image center int[n/2]+1.
00182 */  
00183         EMData* 
00184         fourierproduct(EMData* f, EMData* g, fp_flag flag, fp_type ptype, bool center) {
00185                 size_t normfact;
00186                 //std::complex<float> phase_mult;
00187                 // Not only does the value of "flag" determine how we handle
00188                 // periodicity, but it also determines whether or not we should
00189                 // normalize the results.  Here's some convenience bools:
00190                 bool donorm = (0 == flag%2) ? true : false;
00191                 // the 2x padding is hardcoded for now
00192                 int  npad  = (flag >= 3) ? 2 : 1;  // amount of padding used
00193                 // g may be NULL.  If so, have g point to the same object as f.  In that
00194                 // case we need to be careful later on not to try to delete g's workspace
00195                 // as well as f's workspace, since they will be the same.
00196                 bool  gexists = true;
00197                 if (!g) { g = f; gexists = false; }
00198                 if ( f->is_complex() || g->is_complex() ) {
00199                         // Fourier input only allowed for circulant
00200                         if (CIRCULANT != flag) {
00201                                 LOGERR("Cannot perform normalization or padding on Fourier type.");
00202                                 throw InvalidValueException(flag, "Cannot perform normalization or padding on Fourier type.");
00203                         }
00204                 }
00205                 // These are actual dimensions of f (and real-space sizes for ny and nz)
00206                 int nx  = f->get_xsize();
00207                 int ny  = f->get_ysize();
00208                 int nz  = f->get_zsize();
00209                 // We manifestly assume no zero-padding here, just the 
00210                 // necessary extension along x for the fft
00211                 if (!f->is_real()) nx = (nx - 2 + (f->is_fftodd() ? 1 : 0)); 
00212 
00213                 // these are padded dimensions
00214                 const int nxp = npad*nx;
00215                 const int nyp = (ny > 1) ? npad*ny : 1; // don't pad y for 1-d image
00216                 const int nzp = (nz > 1) ? npad*nz : 1; // don't pad z for 2-d image
00217 
00218                 // now one half of the padded, fft-extended size along x
00219                 const int lsd2 = (nxp + 2 - nxp%2) / 2; 
00220 
00221                 EMData* fp = NULL;
00222                 if (f->is_complex()) { 
00223                         // If f is already a fourier object then fp is a copy of f.
00224                         // (The fp workspace is modified, so we copy f to keep f pristine.)
00225                         fp=f->copy();
00226                 } else {
00227                         //  [normalize] [pad] compute fft
00228                         fp = f->norm_pad(donorm, npad);
00229                         fp->do_fft_inplace();
00230                 }
00231                 // The [padded] fft-extended version of g is gp.
00232                 EMData* gp = NULL;
00233                 if(f==g) {
00234                         // g is an alias for f, so gp should be an alias for fp
00235                         gp=fp;
00236                 } else if (g->is_complex()) {
00237                         // g is already a Fourier object, so gp is just an alias for g
00238                         // (The gp workspace is not modified, so we don't need a copy.)
00239                         gp = g;
00240                 } else {
00241                         // normal case: g is real and different from f, so compute gp
00242                         gp = g->norm_pad(donorm, npad);
00243                         gp->do_fft_inplace();
00244                 }
00245                 // Get complex matrix views of fp and gp; matrices start from 1 (not 0)
00246                 fp->set_array_offsets(1,1,1);
00247                 gp->set_array_offsets(1,1,1);
00248                 
00249                 // If the center flag is true, put the center of the correlation in the middle
00250                 // If it is false, put it in (0,0), this approach saves time, but it is diffcult to manage the result
00251                 if (center) {
00252                         //  Multiply two functions (the real work of this routine)
00253                         int itmp = nx/2;
00254                         //float sx  = float(-twopi*float(itmp)/float(nxp));
00255                         float sxn = 2*float(itmp)/float(nxp);
00256                         float sx = -M_PI*sxn;
00257                         itmp = ny/2;
00258                         //float sy  = float(-twopi*float(itmp)/float(nyp));
00259                         float syn = 2*float(itmp)/float(nyp);
00260                         float sy = -M_PI*syn;
00261                         itmp = nz/2;
00262                         //float sz  = float(-twopi*float(itmp)/float(nzp));
00263                         float szn = 2*float(itmp)/float(nzp);
00264                         float sz = -M_PI*szn;
00265                         if ( nx%2==0 && (ny%2==0 || ny==1 ) && (nz%2==0 || nz==1 ) ) {
00266                                 switch (ptype) {
00267                                         case AUTOCORRELATION:
00268                                         // fpmat := |fpmat|^2
00269                                         // Note nxp are padded dimensions
00270                                                 for (int iz = 1; iz <= nzp; iz++) {
00271                                                         for (int iy = 1; iy <= nyp; iy++) {
00272                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00273                                                                         float fpr = real(fp->cmplx(ix,iy,iz));
00274                                                                         float fpi = imag(fp->cmplx(ix,iy,iz));
00275                                                                         fp->cmplx(ix,iy,iz) = complex<float>(fpr*fpr+fpi*fpi, 0.0f);
00276                                                                 }
00277                                                         }
00278                                                 }
00279                                                 break;
00280                                         case SELF_CORRELATION:
00281                                         // fpmat:=|fpmat|
00282                                         // Note nxp are padded dimensions
00283                                                 for (int iz = 1; iz <= nzp; iz++) {
00284                                                         for (int iy = 1; iy <= nyp; iy++) {
00285                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00286                                                                         fp->cmplx(ix,iy,iz) = complex<float>(abs(fp->cmplx(ix,iy,iz)), 0.0f);
00287                                                                 }
00288                                                         }
00289                                                 }
00290                                                 break;
00291                                         case CORRELATION:
00292                                         // fpmat:=fpmat*conjg(gpmat)
00293                                         // Note nxp are padded dimensions
00294                                                 for (int iz = 1; iz <= nzp; iz++) {
00295                                                         for (int iy = 1; iy <= nyp; iy++) {
00296                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00297                                                                         fp->cmplx(ix,iy,iz) *= conj(gp->cmplx(ix,iy,iz));
00298                                                                 }
00299                                                         }
00300                                                 }
00301                                         break;
00302                                         case CONVOLUTION:
00303                                         // fpmat:=fpmat*gpmat
00304                                         // Note nxp are padded dimensions
00305                                                 for (int iz = 1; iz <= nzp; iz++) {
00306                                                         for (int iy = 1; iy <= nyp; iy++) {
00307                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00308                                                                         fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz);
00309                                                                 }
00310                                                         }
00311                                                 }
00312                                                 break;
00313                                         default:
00314                                                 LOGERR("Illegal option in Fourier Product");
00315                                                 throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00316                                 }                                       
00317                                 for (int iz = 1; iz <= nzp; iz++) {
00318                                         for (int iy = 1; iy <= nyp; iy++) {
00319                                                 for (int ix = (iz+iy+1)%2+1; ix <= lsd2; ix+=2) {
00320                                                         fp->cmplx(ix,iy,iz) = -fp->cmplx(ix,iy,iz);
00321                                                 }
00322                                         }
00323                                 }
00324                         } else {                        
00325                                 switch (ptype) {
00326                                         case AUTOCORRELATION:
00327                                         // fpmat := |fpmat|^2
00328                                         // Note nxp are padded dimensions
00329                                                 for (int iz = 1; iz <= nzp; iz++) {
00330                                                 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00331                                                         for (int iy = 1; iy <= nyp; iy++) {
00332                                                         int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00333                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00334                                                                         int jx=ix-1; float arg=sx*jx+argy;
00335                                                                         float fpr = real(fp->cmplx(ix,iy,iz));
00336                                                                         float fpi = imag(fp->cmplx(ix,iy,iz));
00337                                                                         fp->cmplx(ix,iy,iz)= (fpr*fpr + fpi*fpi) *std::complex<float>(cos(arg),sin(arg));
00338                                                                 }
00339                                                         }
00340                                                 }
00341                                                 break;
00342                                         case SELF_CORRELATION:
00343                                         // fpmat:=|fpmat|
00344                                         // Note nxp are padded dimensions
00345                                                 for (int iz = 1; iz <= nzp; iz++) {
00346                                                 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00347                                                         for (int iy = 1; iy <= nyp; iy++) {
00348                                                         int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00349                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00350                                                                         int jx=ix-1; float arg=sx*jx+argy;
00351                                                                         fp->cmplx(ix,iy,iz) = abs(fp->cmplx(ix,iy,iz)) *std::complex<float>(cos(arg),sin(arg));
00352                                                                 }
00353                                                         }
00354                                                 }
00355                                                 break;
00356                                         case CORRELATION:
00357                                         // fpmat:=fpmat*conjg(gpmat)
00358                                         // Note nxp are padded dimensions
00359                                                 for (int iz = 1; iz <= nzp; iz++) {
00360                                                 int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00361                                                         for (int iy = 1; iy <= nyp; iy++) {
00362                                                         int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00363                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00364                                                                         int jx=ix-1; float arg=sx*jx+argy;
00365                                                                         fp->cmplx(ix,iy,iz) *= conj(gp->cmplx(ix,iy,iz)) *std::complex<float>(cos(arg),sin(arg));
00366                                                                 }
00367                                                         }
00368                                                 }
00369                                         break;
00370                                         case CONVOLUTION:
00371                                         // fpmat:=fpmat*gpmat
00372                                         // Note nxp are padded dimensions
00373                                                 if(npad == 1) {
00374                                                         sx -= 4*(nx%2)/float(nx);
00375                                                         sy -= 4*(ny%2)/float(ny);
00376                                                         sz -= 4*(nz%2)/float(nz);
00377                                                 }
00378                                                 for (int iz = 1; iz <= nzp; iz++) {
00379                                                         int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00380                                                         for (int iy = 1; iy <= nyp; iy++) {
00381                                                                 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00382                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00383                                                                         int jx=ix-1; float arg=sx*jx+argy;
00384                                                                         fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz) *std::complex<float>(cos(arg),sin(arg));
00385                                                                 }
00386                                                         }
00387                                                 }
00388                                                 break;
00389                                         default:
00390                                                 LOGERR("Illegal option in Fourier Product");
00391                                                 throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00392                                 }
00393                         }
00394                 } else {
00395                         // If the center flag is false, then just do basic multiplication
00396                         // Here I aterd the method of complex calculation. This method is much faster than the previous one.
00397                         switch (ptype) {
00398                                 case AUTOCORRELATION:
00399                                         for (int iz = 1; iz <= nzp; iz++) {
00400                                                 for (int iy = 1; iy <= nyp; iy++) {
00401                                                         for (int ix = 1; ix <= lsd2; ix++) {
00402                                                                 float fpr = real(fp->cmplx(ix,iy,iz));
00403                                                                 float fpi = imag(fp->cmplx(ix,iy,iz));
00404                                                                 fp->cmplx(ix,iy,iz) = complex<float>(fpr*fpr+fpi*fpi, 0.0f);
00405                                                         }
00406                                                 }
00407                                         }
00408                                         break;
00409                                 case SELF_CORRELATION:
00410                                         for (int iz = 1; iz <= nzp; iz++) {
00411                                                 for (int iy = 1; iy <= nyp; iy++) {
00412                                                         for (int ix = 1; ix <= lsd2; ix++) {
00413                                                                 fp->cmplx(ix,iy,iz) = complex<float>(abs(fp->cmplx(ix,iy,iz)), 0.0f);
00414                                                         }
00415                                                 }
00416                                         }
00417                                         break;
00418                                 case CORRELATION:
00419                                         //phase_mult = 1;
00420                                         for (int iz = 1; iz <= nzp; iz++) {
00421                                                 for (int iy = 1; iy <= nyp; iy++) {
00422                                                         for (int ix = 1; ix <= lsd2; ix++) {
00423                                                                 fp->cmplx(ix,iy,iz)*= conj(gp->cmplx(ix,iy,iz));
00424                                                         }
00425                                                 }
00426                                         }
00427                                         break;
00428                                 case CONVOLUTION:
00429                                         if(npad == 1) {
00430                                                 float sx = -M_PI*2*(nx%2)/float(nx);
00431                                                 float sy = -M_PI*2*(ny%2)/float(ny);
00432                                                 float sz = -M_PI*2*(nz%2)/float(nz);
00433                                                 for (int iz = 1; iz <= nzp; iz++) {
00434                                                         int jz=iz-1; if(jz>nzp/2) jz=jz-nzp; float argz=sz*jz;
00435                                                         for (int iy = 1; iy <= nyp; iy++) {
00436                                                                 int jy=iy-1; if(jy>nyp/2) jy=jy-nyp; float argy=sy*jy+argz;
00437                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00438                                                                         int jx=ix-1; float arg=sx*jx+argy;
00439                                                                         fp->cmplx(ix,iy,iz) *= gp->cmplx(ix,iy,iz) *std::complex<float>(cos(arg),sin(arg));
00440                                                                 }
00441                                                         }
00442                                                 }
00443                                         } else {
00444                                                 for (int iz = 1; iz <= nzp; iz++) {
00445                                                         for (int iy = 1; iy <= nyp; iy++) {
00446                                                                 for (int ix = 1; ix <= lsd2; ix++) {
00447                                                                         fp->cmplx(ix,iy,iz)*= gp->cmplx(ix,iy,iz);
00448                                                                 }
00449                                                         }
00450                                                 }
00451                                         }
00452                                         break;
00453                                 default:
00454                                         LOGERR("Illegal option in Fourier Product");
00455                                         throw InvalidValueException(ptype, "Illegal option in Fourier Product");
00456                         }
00457                 }
00458                 // Now done w/ gp, so let's get rid of it (if it's not an alias of fp or simply g was complex on input);
00459                 if (gexists && (f != g) && (!g->is_complex())) {
00460                         if( gp ) {
00461                                 delete gp;
00462                                 gp = 0;
00463                         }
00464                 }
00465                 // back transform
00466                 fp->do_ift_inplace();
00467                 if(center && npad ==2)  fp->depad();
00468                 else                    fp->depad_corner();
00469 
00470                 //vector<int> saved_offsets = fp->get_array_offsets();  I do not know what the meaning of it was, did not work anyway PAP
00471                 fp->set_array_offsets(1,1,1);
00472 
00473                 normfact = (size_t)(nxp/nx)*(nyp/ny)*(nzp/nz);  // Normalization factor for the padded operations
00474                 if(normfact>1) {
00475                         for (int iz = 1; iz <= nz; iz++) for (int iy = 1; iy <= ny; iy++) for (int ix = 1; ix <= nx; ix++) (*fp)(ix,iy,iz) *= normfact;
00476                 }
00477                 // Lag normalization
00478                 if(flag>4)  {
00479                         normfact = (size_t)nx*ny*nz;  // Normalization factor
00480                         int nxc=nx/2+1, nyc=ny/2+1, nzc=nz/2+1;
00481                         for (int iz = 1; iz <= nz; iz++) {
00482                                 float lagz=float(normfact/(nz-abs(iz-nzc)));
00483                                 for (int iy = 1; iy <= ny; iy++) {
00484                                         float lagyz=lagz/(ny-abs(iy-nyc));
00485                                         for (int ix = 1; ix <= nx; ix++) {
00486                                                 (*fp)(ix,iy,iz) *= lagyz/(nx-abs(ix-nxc));
00487                                         }
00488                                 }
00489                         }       
00490                 }
00491                 //OVER AND OUT
00492                 //fp->set_array_offsets(saved_offsets);  This was strange and did not work, PAP
00493                 fp->set_array_offsets(0,0,0);
00494                 fp->update();
00495                 return fp;
00496         }
00497 }
00498 

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