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