reconstructor.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "reconstructor.h"
00037 #include "plugins/reconstructor_template.h"
00038 #include "ctf.h"
00039 #include "emassert.h"
00040 #include "symmetry.h"
00041 #include <cstring>
00042 #include <fstream>
00043 #include <iomanip>
00044 #include <boost/bind.hpp>
00045 #include <boost/format.hpp>
00046 
00047 #include <gsl/gsl_statistics_double.h>
00048 #include <gsl/gsl_fit.h>
00049 
00050 using namespace EMAN;
00051 using std::complex;
00052 
00053 
00054 #include <iostream>
00055 using std::cerr;
00056 using std::endl;
00057 using std::cout; // for debug
00058 
00059 #include <algorithm>
00060 // find, for_each
00061 
00062 #include <iomanip>
00063 using std::setprecision;
00064 
00065 template < typename T > void checked_delete( T*& x )
00066 {
00067     typedef char type_must_be_complete[ sizeof(T)? 1: -1 ];
00068     (void) sizeof(type_must_be_complete);
00069     delete x;
00070     x = NULL;
00071 }
00072 
00073 const string FourierReconstructor::NAME = "fourier";
00074 const string FourierReconstructorSimple2D::NAME = "fouriersimple2D";
00075 const string WienerFourierReconstructor::NAME = "wiener_fourier";
00076 const string BackProjectionReconstructor::NAME = "back_projection";
00077 const string nn4Reconstructor::NAME = "nn4";
00078 const string nnSSNR_Reconstructor::NAME = "nnSSNR";
00079 const string nn4_ctfReconstructor::NAME = "nn4_ctf";
00080 const string nnSSNR_ctfReconstructor::NAME = "nnSSNR_ctf";
00081 
00082 template <> Factory < Reconstructor >::Factory()
00083 {
00084         force_add<FourierReconstructor>();
00085         force_add<FourierReconstructorSimple2D>();
00086 //      force_add(&BaldwinWoolfordReconstructor::NEW);
00087 //      force_add(&WienerFourierReconstructor::NEW);
00088         force_add<BackProjectionReconstructor>();
00089         force_add<nn4Reconstructor>();
00090         force_add<nnSSNR_Reconstructor>();
00091         force_add<nn4_ctfReconstructor>();
00092         force_add<nnSSNR_ctfReconstructor>();
00093 //      force_add<XYZReconstructor>();
00094 }
00095 
00096 void FourierReconstructorSimple2D::setup()
00097 {
00098         nx = params.set_default("nx",0);
00099 
00100         if ( nx < 0 ) throw InvalidValueException(nx, "nx must be positive");
00101 
00102         bool is_fftodd = (nx % 2 == 1);
00103 
00104         ny = nx;
00105         nx += 2-is_fftodd;
00106 
00107         image = new EMData();
00108         image->set_size(nx, ny);
00109         image->set_complex(true);
00110         image->set_fftodd(is_fftodd);
00111         image->set_ri(true);
00112 
00113         tmp_data = new EMData();
00114         tmp_data->set_size(nx/2, nx);
00115 }
00116 
00117 int FourierReconstructorSimple2D::insert_slice(const EMData* const slice, const Transform & euler, const float weight)
00118 {
00119 
00120         // Are these exceptions really necessary? (d.woolford)
00121         if (!slice) throw NullPointerException("EMData pointer (input image) is NULL");
00122 
00123         if ( slice->get_ndim() != 1 ) throw ImageDimensionException("Image dimension must be 1");
00124 
00125         // I could also test to make sure the image is the right dimensions...
00126         if (slice->is_complex()) throw ImageFormatException("The image is complex, expecting real");
00127 
00128         EMData* working_slice = slice->process("xform.phaseorigin.tocorner");
00129 
00130         // Fourier transform the slice
00131         working_slice->do_fft_inplace();
00132 
00133         float* rdata = image->get_data();
00134         float* norm = tmp_data->get_data();
00135         float* dat = working_slice->get_data();
00136 
00137         float g[4];
00138         int offset[4];
00139         float dt[2];
00140         offset[0] = 0; offset[1] = 2; offset[2] = nx; offset[3] = nx+2;
00141 
00142         float alt = -((float)(euler.get_rotation("2d"))["alpha"])*M_PI/180.0f;
00143         for (int x = 0; x < working_slice->get_xsize() / 2; x++) {
00144 
00145                 float rx = (float) x;
00146 
00147                 float xx = rx*cos(alt);
00148                 float yy = rx*sin(alt);
00149                 float cc = 1.0;
00150 
00151                 if (xx < 0) {
00152                         xx = -xx;
00153                         yy = -yy;
00154                         cc = -1.0;
00155                 }
00156 
00157                 yy += ny / 2;
00158 
00159 
00160                 dt[0] = dat[2*x];
00161                 dt[1] = cc * dat[2*x+1];
00162 
00163                 // PHIL IS INTERESTED FROM HERE DOWN
00164                 int x0 = (int) floor(xx);
00165                 int y0 = (int) floor(yy);
00166 
00167                 int i = 2*x0 + y0*nx;
00168 
00169                 float dx = xx - x0;
00170                 float dy = yy - y0;
00171 
00172                 g[0] = Util::agauss(1, dx, dy, 0, EMConsts::I2G);
00173                 g[1] = Util::agauss(1, 1 - dx, dy, 0, EMConsts::I2G);
00174                 g[2] = Util::agauss(1, dx, 1 - dy, 0, EMConsts::I2G);
00175                 g[3] = Util::agauss(1, 1 - dx, 1 - dy, 0, EMConsts::I2G);
00176 
00177                 // At the extreme we can only do some much...
00178                 if ( x0 == nx-2 ) {
00179                         int k = i + offset[0];
00180                         rdata[k] += g[0] * dt[0];
00181                         rdata[k + 1] += g[0] * dt[1];
00182                         norm[k/2] += g[0];
00183 
00184                         k = i + offset[2];
00185                         rdata[k] += g[2] * dt[0];
00186                         rdata[k + 1] += g[2] * dt[1];
00187                         norm[k/2] += g[2];
00188                         continue;
00189 
00190                 }
00191                 // capture and accommodate for periodic boundary conditions in the x direction
00192                 if ( x0 > nx-2 ) {
00193                         int dif = x0 - (nx-2);
00194                         x0 -= dif;
00195                 }
00196                 // At the extreme we can only do some much...
00197                 if ( y0 == ny -1 ) {
00198                         int k = i + offset[0];
00199                         rdata[k] += g[0] * dt[0];
00200                         rdata[k + 1] += g[0] * dt[1];
00201                         norm[k/2] += g[0];
00202 
00203                         k = i + offset[1];
00204                         rdata[k] += g[1] * dt[0];
00205                         rdata[k + 1] += g[1] * dt[1];
00206                         norm[k/2] += g[1];
00207                         continue;
00208                 }
00209                 // capture and accommodate for periodic boundary conditions in the y direction
00210                 if ( y0 > ny-1) {
00211                         int dif = y0 - (ny-1);
00212                         y0 -= dif;
00213                 }
00214 
00215                 if (x0 >= nx - 2 || y0 >= ny - 1) continue;
00216 
00217 
00218 
00219 
00220                 for (int j = 0; j < 4; j++)
00221                 {
00222                         int k = i + offset[j];
00223                         rdata[k] += g[j] * dt[0];
00224                         rdata[k + 1] += g[j] * dt[1];
00225                         norm[k/2] += g[j];
00226 
00227                 }
00228         }
00229 
00230         return 0;
00231 
00232 }
00233 
00234 EMData *FourierReconstructorSimple2D::finish(bool doift)
00235 {
00236         normalize_threed();
00237 
00238         image->process_inplace("xform.fourierorigin.tocorner");
00239         image->do_ift_inplace();
00240         image->depad();
00241         image->process_inplace("xform.phaseorigin.tocenter");
00242 
00243         return          image;
00244 }
00245 
00246 void ReconstructorVolumeData::normalize_threed(const bool sqrt_damp)
00247 // normalizes the 3-D Fourier volume. Also imposes appropriate complex conjugate relationships
00248 {
00249         float* norm = tmp_data->get_data();
00250         float* rdata = image->get_data();
00251 
00252 //      printf("%d %d %d %d %d %d\n",subnx,subny,subnz,image->get_xsize(),image->get_ysize(),image->get_zsize());
00253 
00254         // FIXME should throw a sensible error
00255         if ( 0 == norm ) throw NullPointerException("The normalization volume was null!");
00256         if ( 0 == rdata ) throw NullPointerException("The complex reconstruction volume was null!");
00257 
00258         for (int i = 0; i < subnx * subny * subnz; i += 2) {
00259                 float d = norm[i/2];
00260                 if (sqrt_damp) d*=sqrt(d);
00261                 if (d == 0) {
00262                         rdata[i] = 0;
00263                         rdata[i + 1] = 0;
00264                 }
00265                 else {
00266 //                      rdata[i]=1.0/d;
00267 //                      rdata[i+1]=0.0;         // for debugging only
00268                         rdata[i] /= d;
00269                         rdata[i + 1] /= d;
00270                 }
00271         }
00272 
00273         // enforce complex conj, only works on subvolumes if the complete conjugate plane is in the volume
00274         if (subx0==0 && subnx>1 && subny==ny && subnz==nz) {
00275                 for (int z=0; z<=nz/2; z++) {
00276                         for (int y=1; y<=ny; y++) {
00277                                 if (y==0 && z==0) continue;
00278                                 // x=0
00279                                 int i=(y%ny)*subnx+(z%nz)*subnx*subny;
00280                                 int i2=(ny-y)*subnx+((nz-z)%nz)*subnx*subny;
00281                                 float ar=(rdata[i]+rdata[i2])/2.0f;
00282                                 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00283                                 rdata[i]=ar;
00284                                 rdata[i2]=ar;
00285                                 rdata[i+1]=ai;
00286                                 rdata[i2+1]=-ai;
00287                         }
00288                 }
00289         }
00290 
00291         if (subx0+subnx==nx && subnx>1 && subny==ny && subnz==nz) {
00292                 for (int z=0; z<=nz/2; z++) {
00293                         for (int y=1; y<=ny; y++) {
00294                                 if (y==0 && z==0) continue;
00295                                 // x=0
00296                                 int i=(y%ny)*subnx+(z%nz)*subnx*subny+subnx-2;
00297                                 int i2=(ny-y)*subnx+((nz-z)%nz)*subnx*subny+subnx-2;
00298                                 float ar=(rdata[i]+rdata[i2])/2.0f;
00299                                 float ai=(rdata[i+1]-rdata[i2+1])/2.0f;
00300                                 rdata[i]=ar;
00301                                 rdata[i2]=ar;
00302                                 rdata[i+1]=ai;
00303                                 rdata[i2+1]=-ai;
00304                         }
00305                 }
00306         }
00307 }
00308 
00309 void FourierReconstructor::free_memory()
00310 {
00311         if (image) { delete image; image=0; }
00312         if (tmp_data) { delete tmp_data; tmp_data=0; }
00313         if ( inserter != 0 )
00314         {
00315                 delete inserter;
00316                 inserter = 0;
00317         }
00318 }
00319 
00320 #include <sstream>
00321 
00322 void FourierReconstructor::load_inserter()
00323 {
00324 //      ss
00325 //      string mode = (string)params["mode"];
00326         Dict parms;
00327         parms["data"] = image;
00328         parms["norm"] = tmp_data->get_data();
00329         // These aren't necessary because we deal with them before calling the inserter
00330 //      parms["subnx"] = nx;
00331 //      parms["subny"] = ny;
00332 //      parms["subnx"] = nz;
00333 //      parms["subx0"] = x0;
00334 //      parms["suby0"] = y0;
00335 //      parms["subz0"] = z0;
00336 
00337         if ( inserter != 0 )
00338         {
00339                 delete inserter;
00340         }
00341 
00342         inserter = Factory<FourierPixelInserter3D>::get((string)params["mode"], parms);
00343         inserter->init();
00344 }
00345 
00346 void FourierReconstructor::setup()
00347 {
00348         // default setting behavior - does not override if the parameter is already set
00349         params.set_default("mode","gauss_2");
00350 
00351         vector<int> size=params["size"];
00352 
00353         nx = size[0];
00354         ny = size[1];
00355         nz = size[2];
00356         nx2=nx/2-1;
00357         ny2=ny/2;
00358         nz2=nz/2;
00359 
00360 
00361         // Adjust nx if for Fourier transform even odd issues
00362         bool is_fftodd = (nx % 2 == 1);
00363         // The Fourier transform requires one extra pixel in the x direction,
00364         // which is two spaces in memory, one each for the complex and the
00365         // real components
00366         nx += 2-is_fftodd;
00367 
00368         if (params.has_key("subvolume")) {
00369                 vector<int> sub=params["subvolume"];
00370                 subx0=sub[0];
00371                 suby0=sub[1];
00372                 subz0=sub[2];
00373                 subnx=sub[3];
00374                 subny=sub[4];
00375                 subnz=sub[5];
00376 
00377                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00378                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00379 
00380         }
00381         else {
00382                 subx0=suby0=subz0=0;
00383                 subnx=nx;
00384                 subny=ny;
00385                 subnz=nz;
00386         }
00387 
00388 
00389         // Odd dimension support is here atm, but not above.
00390         if (image) delete image;
00391         image = new EMData();
00392         image->set_size(subnx, subny, subnz);
00393         image->set_complex(true);
00394         image->set_fftodd(is_fftodd);
00395         image->set_ri(true);
00396         image->to_zero();
00397 
00398         if (params.has_key("subvolume")) {
00399                 image->set_attr("subvolume_x0",subx0);
00400                 image->set_attr("subvolume_y0",suby0);
00401                 image->set_attr("subvolume_z0",subz0);
00402                 image->set_attr("subvolume_full_nx",nx);
00403                 image->set_attr("subvolume_full_ny",ny);
00404                 image->set_attr("subvolume_full_nz",nz);
00405         }
00406         
00407         if (tmp_data) delete tmp_data;
00408         tmp_data = new EMData();
00409         tmp_data->set_size(subnx/2, subny, subnz);
00410         tmp_data->to_zero();
00411         tmp_data->update();
00412 
00413         load_inserter();
00414 
00415         if ( (bool) params["quiet"] == false )
00416         {
00417                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00418                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00419                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00420         }
00421 }
00422 
00423 void FourierReconstructor::setup_seed(EMData* seed,float seed_weight) {
00424         // default setting behavior - does not override if the parameter is already set
00425         params.set_default("mode","gauss_2");
00426 
00427         vector<int> size=params["size"];
00428 
00429         nx = size[0];
00430         ny = size[1];
00431         nz = size[2];
00432         nx2=nx/2-1;
00433         ny2=ny/2;
00434         nz2=nz/2;
00435 
00436 
00437         // Adjust nx if for Fourier transform even odd issues
00438         bool is_fftodd = (nx % 2 == 1);
00439         // The Fourier transform requires one extra pixel in the x direction,
00440         // which is two spaces in memory, one each for the complex and the
00441         // real components
00442         nx += 2-is_fftodd;
00443 
00444         if (params.has_key("subvolume")) {
00445                 vector<int> sub=params["subvolume"];
00446                 subx0=sub[0];
00447                 suby0=sub[1];
00448                 subz0=sub[2];
00449                 subnx=sub[3];
00450                 subny=sub[4];
00451                 subnz=sub[5];
00452 
00453                 if (subx0<0 || suby0<0 || subz0<0 || subx0+subnx>nx || suby0+subny>ny || subz0+subnz>nz)
00454                         throw ImageDimensionException("The subvolume cannot extend outside the reconstructed volume");
00455 
00456         }
00457         else {
00458                 subx0=suby0=subz0=0;
00459                 subnx=nx;
00460                 subny=ny;
00461                 subnz=nz;
00462         }
00463 
00464         if (seed->get_xsize()!=subnx || seed->get_ysize()!=subny || seed->get_zsize()!=subnz || !seed->is_complex())
00465                 throw ImageDimensionException("The dimensions of the seed volume do not match the reconstruction size");
00466 
00467         // Odd dimension support is here atm, but not above.
00468         image = seed;
00469         if (params.has_key("subvolume")) {
00470                 image->set_attr("subvolume_x0",subx0);
00471                 image->set_attr("subvolume_y0",suby0);
00472                 image->set_attr("subvolume_z0",subz0);
00473                 image->set_attr("subvolume_full_nx",nx);
00474                 image->set_attr("subvolume_full_ny",ny);
00475                 image->set_attr("subvolume_full_nz",nz);
00476         }
00477 
00478         if (tmp_data) delete tmp_data;
00479         tmp_data = new EMData();
00480         tmp_data->set_size(subnx/2, subny, subnz);
00481         tmp_data->to_value(seed_weight);
00482 
00483         load_inserter();
00484 
00485         if ( (bool) params["quiet"] == false )
00486         {
00487                 cout << "Seeded direct Fourier inversion";
00488                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00489                 cout << "3D Fourier subvolume is " << subnx << " " << subny << " " << subnz << endl;
00490                 cout << "You will require approximately " << setprecision(3) << (subnx*subny*subnz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00491         }
00492 }
00493 
00494 
00495 EMData* FourierReconstructor::preprocess_slice( const EMData* const slice,  const Transform& t )
00496 {
00497         // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
00498         EMData* return_slice = 0;
00499         Transform tmp(t);
00500         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly, this way we only do 2d translation,scaling and mirroring
00501 
00502         if (tmp.is_identity()) return_slice=slice->copy();
00503         else return_slice = slice->process("xform",Dict("transform",&tmp));
00504 
00505         return_slice->process_inplace("xform.phaseorigin.tocorner");
00506 
00507 //      printf("@@@ %d\n",(int)return_slice->get_attr("nx"));
00508         // Fourier transform the slice
00509         return_slice->do_fft_inplace();
00510 
00511 //      printf("%d\n",(int)return_slice->get_attr("nx"));
00512 
00513         return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00514 
00515         // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
00516 //      return_slice->process_inplace("xform.fourierorigin.tocenter");
00517 
00518         return_slice->set_attr("reconstruct_preproc",(int)1);
00519         return return_slice;
00520 }
00521 
00522 
00523 int FourierReconstructor::insert_slice(const EMData* const input_slice, const Transform & arg, const float weight)
00524 {
00525         // Are these exceptions really necessary? (d.woolford)
00526         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00527 
00528         Transform * rotation;
00529 /*      if ( input_slice->has_attr("xform.projection") ) {
00530                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00531         } else {*/
00532         rotation = new Transform(arg); // assignment operator
00533 //      }
00534 
00535         EMData *slice;
00536         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00537         else slice = preprocess_slice( input_slice, *rotation);
00538 
00539 
00540         // We must use only the rotational component of the transform, scaling, translation and mirroring
00541         // are not implemented in Fourier space, but are in preprocess_slice
00542         rotation->set_scale(1.0);
00543         rotation->set_mirror(false);
00544         rotation->set_trans(0,0,0);
00545 
00546         // Finally to the pixel wise slice insertion
00547         do_insert_slice_work(slice, *rotation, weight);
00548 
00549         delete rotation; rotation=0;
00550         delete slice;
00551 
00552 //      image->update();
00553         return 0;
00554 }
00555 
00556 void FourierReconstructor::do_insert_slice_work(const EMData* const input_slice, const Transform & arg,const float weight)
00557 {
00558         // Reload the inserter if the mode has changed
00559 //      string mode = (string) params["mode"];
00560 //      if ( mode != inserter->get_name() )     load_inserter();
00561 
00562 //      int y_in = input_slice->get_ysize();
00563 //      int x_in = input_slice->get_xsize();
00564 //      // Adjust the dimensions to account for odd and even ffts
00565 //      if (input_slice->is_fftodd()) x_in -= 1;
00566 //      else x_in -= 2;
00567 
00568         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00569 
00570         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00571         float iny=(float)(input_slice->get_ysize());
00572 
00573         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00574                 Transform t3d = arg*(*it);
00575                 for (int y = -iny/2; y < iny/2; y++) {
00576                         for (int x = 0; x <=  inx/2; x++) {
00577 
00578                                 float rx = (float) x/(inx-2.0f);        // coords relative to Nyquist=.5
00579                                 float ry = (float) y/iny;
00580 
00581                                 Vec3f coord(rx,ry,0);
00582                                 coord = coord*t3d; // transpose multiplication
00583                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00584                                 float yy = coord[1];
00585                                 float zz = coord[2];
00586 
00587                                 // Map back to real pixel coordinates in output volume
00588                                 xx=xx*(nx-2);
00589                                 yy=yy*ny;
00590                                 zz=zz*nz;
00591 
00592 //                              if (x==10 && y==0) printf("10,0 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f) %1.0f %d\n",
00593 //                                      xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2),inx,nx);
00594 //                              if (x==0 && y==10) printf("0,10 -> %1.2f,%1.2f,%1.2f\t(%5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f   %5.2f %5.2f %5.2f)\n",
00595 //                                      xx,yy,zz,t3d.at(0,0),t3d.at(0,1),t3d.at(0,2),t3d.at(1,0),t3d.at(1,1),t3d.at(1,2),t3d.at(2,0),t3d.at(2,1),t3d.at(2,2));
00596 
00597                                 //printf("%3.1f %3.1f %3.1f\t %1.4f %1.4f\t%1.4f\n",xx,yy,zz,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00598 //                              if (floor(xx)==45 && floor(yy)==45 &&floor(zz)==0) printf("%d. 45 45 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00599 //                              if (floor(xx)==21 && floor(yy)==21 &&floor(zz)==0) printf("%d. 21 21 0\t %d %d\t %1.4f %1.4f\t%1.4f\n",(int)input_slice->get_attr("n"),x,y,input_slice->get_complex_at(x,y).real(),input_slice->get_complex_at(x,y).imag(),weight);
00600                                 inserter->insert_pixel(xx,yy,zz,input_slice->get_complex_at(x,y),weight);
00601                         }
00602                 }
00603         }
00604 }
00605 
00606 int FourierReconstructor::determine_slice_agreement(EMData*  input_slice, const Transform & arg, const float weight,bool sub)
00607 {
00608         // Are these exceptions really necessary? (d.woolford)
00609         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00610 
00611         Transform * rotation;
00612         rotation = new Transform(arg); // assignment operator
00613 
00614         EMData *slice;
00615         if (input_slice->get_attr_default("reconstruct_preproc",(int) 0)) slice=input_slice->copy();
00616         else slice = preprocess_slice( input_slice, *rotation);
00617 
00618 
00619         // We must use only the rotational component of the transform, scaling, translation and mirroring
00620         // are not implemented in Fourier space, but are in preprocess_slice
00621         rotation->set_scale(1.0);
00622         rotation->set_mirror(false);
00623         rotation->set_trans(0,0,0);
00624 
00625         // Remove the current slice first (not threadsafe, but otherwise performance would be awful)
00626         if (sub) do_insert_slice_work(slice, *rotation, -weight);
00627 
00628         // Compare
00629         do_compare_slice_work(slice, *rotation,weight);
00630 
00631         input_slice->set_attr("reconstruct_norm",slice->get_attr("reconstruct_norm"));
00632         input_slice->set_attr("reconstruct_absqual",slice->get_attr("reconstruct_absqual"));
00633 //      input_slice->set_attr("reconstruct_qual",slice->get_attr("reconstruct_qual"));
00634         input_slice->set_attr("reconstruct_weight",slice->get_attr("reconstruct_weight"));
00635 
00636         // Now put the slice back
00637         if (sub) do_insert_slice_work(slice, *rotation, weight);
00638 
00639 
00640         delete rotation;
00641         delete slice;
00642 
00643 //      image->update();
00644         return 0;
00645 
00646 }
00647 
00648 void FourierReconstructor::do_compare_slice_work(EMData* input_slice, const Transform & arg,float weight)
00649 {
00650 
00651         float dt[3];    // This stores the complex and weight from the volume
00652         float dt2[2];   // This stores the local image complex
00653         float *dat = input_slice->get_data();
00654         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00655 
00656         float inx=(float)(input_slice->get_xsize());            // x/y dimensions of the input image
00657         float iny=(float)(input_slice->get_ysize());
00658 
00659         double dot=0;           // summed pixel*weight dot product
00660         double vweight=0;               // sum of weights
00661         double power=0;         // sum of inten*weight from volume
00662         double power2=0;                // sum of inten*weight from image
00663         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00664                 Transform t3d = arg*(*it);
00665                 for (int y = -iny/2; y < iny/2; y++) {
00666                         for (int x = 0; x <=  inx/2; x++) {
00667                                 if (x==0 && y==0) continue;             // We don't want to use the Fourier origin
00668 
00669                                 float rx = (float) x/(inx-2);   // coords relative to Nyquist=.5
00670                                 float ry = (float) y/iny;
00671 
00672 //                              if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
00673 //                                      continue;
00674 
00675                                 Vec3f coord(rx,ry,0);
00676                                 coord = coord*t3d; // transpose multiplication
00677                                 float xx = coord[0]; // transformed coordinates in terms of Nyquist
00678                                 float yy = coord[1];
00679                                 float zz = coord[2];
00680 
00681 
00682                                 if (fabs(xx)>0.5 || fabs(yy)>=0.5 || fabs(zz)>=0.5) continue;
00683 
00684                                 // Map back to actual pixel coordinates in output volume
00685                                 xx=xx*(nx-2);
00686                                 yy=yy*ny;
00687                                 zz=zz*nz;
00688 
00689 
00690                                 int idx = (int)(x * 2 + inx*(y<0?iny+y:y));
00691                                 dt2[0] = dat[idx];
00692                                 dt2[1] = dat[idx+1];
00693 
00694                                 // value returned indirectly in dt
00695                                 if (!pixel_at(xx,yy,zz,dt) || dt[2]==0) continue;
00696 
00697 //                              printf("%f\t%f\t%f\t%f\t%f\n",dt[0],dt[1],dt[2],dt2[0],dt2[1]);
00698                                 dot+=(dt[0]*dt2[0]+dt[1]*dt2[1])*dt[2];
00699                                 vweight+=dt[2];
00700                                 power+=(dt[0]*dt[0]+dt[1]*dt[1])*dt[2];
00701                                 power2+=(dt2[0]*dt2[0]+dt2[1]*dt2[1])*dt[2];
00702                         }
00703                 }
00704         }
00705 
00706         dot/=sqrt(power*power2);                // normalize the dot product
00707         input_slice->set_attr("reconstruct_norm",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)));
00708         input_slice->set_attr("reconstruct_absqual",(float)dot);
00709         float rw=weight<=0?1.0f:1.0f/weight;
00710         input_slice->set_attr("reconstruct_qual",(float)(dot*rw/((rw-1.0)*dot+1.0)));   // here weight is a proxy for SNR
00711         input_slice->set_attr("reconstruct_weight",(float)vweight/(float)(subnx*subny*subnz));
00712 //      printf("** %g\t%g\t%g\t%g ##\n",dot,vweight,power,power2);
00713         //printf("** %f %f %f ##\n",(float)(power2<=0?1.0:sqrt(power/power2)/(inx*iny)),(float)dot,(float)(dot*weight/((weight-1.0)*dot+1.0)));
00714 }
00715 
00716 bool FourierReconstructor::pixel_at(const float& xx, const float& yy, const float& zz, float *dt)
00717 {
00718         int x0 = (int) floor(xx);
00719         int y0 = (int) floor(yy);
00720         int z0 = (int) floor(zz);
00721         
00722         float *rdata=image->get_data();
00723         float *norm=tmp_data->get_data();
00724         float normsum=0;
00725 
00726         dt[0]=dt[1]=dt[2]=0.0;
00727 
00728         if (nx==subnx) {                        // normal full reconstruction
00729                 if (x0<-nx2-1 || y0<-ny2-1 || z0<-nz2-1 || x0>nx2 || y0>ny2 || z0>nz2 ) return false;
00730 
00731                 // no error checking on add_complex_fast, so we need to be careful here
00732                 int x1=x0+1;
00733                 int y1=y0+1;
00734                 int z1=z0+1;
00735                 if (x0<-nx2) x0=-nx2;
00736                 if (x1>nx2) x1=nx2;
00737                 if (y0<-ny2) y0=-ny2;
00738                 if (y1>ny2) y1=ny2;
00739                 if (z0<-nz2) z0=-nz2;
00740                 if (z1>nz2) z1=nz2;
00741                 
00742                 size_t idx;
00743                 float r, gg;
00744                 for (int k = z0 ; k <= z1; k++) {
00745                         for (int j = y0 ; j <= y1; j++) {
00746                                 for (int i = x0; i <= x1; i ++) {
00747                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00748                                         idx=image->get_complex_index_fast(i,j,k);
00749                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00750                                         
00751                                         dt[0]+=gg*rdata[idx];
00752                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00753                                         dt[2]+=norm[idx/2];
00754                                         normsum+=gg;                            
00755                                 }
00756                         }
00757                 }
00758                 dt[0]/=normsum;
00759                 dt[1]/=normsum;
00760 //              printf("%1.2f,%1.2f,%1.2f\t%1.3f\t%1.3f\t%1.3f\n",xx,yy,zz,dt[0],dt[1],dt[2]);
00761                 return true;
00762         } 
00763         else {                                  // for subvolumes, not optimized yet
00764                 size_t idx;
00765                 float r, gg;
00766                 for (int k = z0 ; k <= z0 + 1; k++) {
00767                         for (int j = y0 ; j <= y0 + 1; j++) {
00768                                 for (int i = x0; i <= x0 + 1; i ++) {
00769                                         r = Util::hypot3sq((float) i - xx, j - yy, k - zz);
00770                                         idx=image->get_complex_index(i,j,k,subx0,suby0,subz0,nx,ny,nz);
00771                                         gg = Util::fast_exp(-r / EMConsts::I2G)*norm[idx/2];
00772                                         
00773                                         dt[0]+=gg*rdata[idx];
00774                                         dt[1]+=(i<0?-1.0f:1.0f)*gg*rdata[idx+1];
00775                                         dt[2]+=norm[idx/2];
00776                                         normsum+=gg;                            
00777                                 }
00778                         }
00779                 }
00780                 
00781                 if (normsum==0)  return false;
00782                 return true;
00783         }
00784 }
00785 
00786 
00787 EMData *FourierReconstructor::finish(bool doift)
00788 {
00789 //      float *norm = tmp_data->get_data();
00790 //      float *rdata = image->get_data();
00791 
00792         bool sqrtnorm=params.set_default("sqrtnorm",false);
00793         normalize_threed(sqrtnorm);
00794 
00795 //      tmp_data->write_image("density.mrc");
00796 
00797         // we may as well delete the tmp data now... it saves memory and the calling program might
00798         // need memory after it gets the return volume.
00799         // If this delete didn't happen now, it would happen when the deconstructor was called --david
00800         // no longer a good idea with the new iterative scheme -- steve
00801 //      if ( tmp_data != 0 )
00802 //      {
00803 //              delete tmp_data;
00804 //              tmp_data = 0;
00805 //      }
00806 
00807 /*      image->process_inplace("xform.fourierorigin.tocorner");*/
00808 
00809         if (doift) {
00810                 image->do_ift_inplace();
00811                 image->depad();
00812                 image->process_inplace("xform.phaseorigin.tocenter");
00813         }
00814         // If the image was padded it should be the original size, as the client would expect
00815         //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
00816         // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
00817         // This should now be handled in the calling program --steve 11/03/09
00818 //      bool is_fftodd = (nx % 2 == 1);
00819 //      if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
00820 //      {
00821 //              FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
00822 //              FloatSize region_size( output_x, output_y, output_z);
00823 //              Region clip_region( origin, region_size );
00824 //              image->clip_inplace( clip_region );
00825 //      }
00826 
00827         // Should be an "if (verbose)" here or something
00828         //print_stats(quality_scores);
00829 
00830         image->update();
00831         
00832         if (params.has_key("savenorm") && strlen((const char *)params["savenorm"])>0) {
00833                 if (tmp_data->get_ysize()%2==0 && tmp_data->get_zsize()%2==0) tmp_data->process_inplace("xform.fourierorigin.tocenter");
00834                 tmp_data->write_image((const char *)params["savenorm"]);
00835         }
00836 
00837         delete tmp_data;
00838         tmp_data=0;
00839         EMData *ret=image;
00840         image=0;
00841         
00842         return ret;
00843 }
00844 
00845 /*
00846 void BaldwinWoolfordReconstructor::setup()
00847 {
00848         //This is a bit of a hack - but for now it suffices
00849         params.set_default("mode","nearest_neighbor");
00850         FourierReconstructor::setup();
00851         // Set up the Baldwin Kernel
00852         int P = (int)((1.0+0.25)*max_input_dim+1);
00853         float r = (float)(max_input_dim+1)/(float)P;
00854         dfreq = 0.2f;
00855         if (W != 0) delete [] W;
00856         int maskwidth = params.set_default("maskwidth",2);
00857         W = Util::getBaldwinGridWeights(maskwidth, (float)P, r,dfreq,0.5f,0.2f);
00858 }
00859 
00860 EMData* BaldwinWoolfordReconstructor::finish(bool doift)
00861 {
00862         tmp_data->write_image("density.mrc");
00863         image->process_inplace("xform.fourierorigin.tocorner");
00864         image->do_ift_inplace();
00865         image->depad();
00866         image->process_inplace("xform.phaseorigin.tocenter");
00867 
00868         if ( (bool) params.set_default("postmultiply", false) )
00869         {
00870                 cout << "POST MULTIPLYING" << endl;
00871         // now undo the Fourier convolution with real space division
00872                 float* d = image->get_data();
00873                 float N = (float) image->get_xsize()/2.0f;
00874                 N *= N;
00875                 size_t rnx = image->get_xsize();
00876                 size_t rny = image->get_ysize();
00877                 size_t rnxy = rnx*rny;
00878                 int cx = image->get_xsize()/2;
00879                 int cy = image->get_ysize()/2;
00880                 int cz = image->get_zsize()/2;
00881                 size_t idx;
00882                 for (int k = 0; k < image->get_zsize(); ++k ){
00883                         for (int j = 0; j < image->get_ysize(); ++j ) {
00884                                 for (int i =0; i < image->get_xsize(); ++i ) {
00885                                         float xd = (float)(i-cx); xd *= xd;
00886                                         float yd = (float)(j-cy); yd *= yd;
00887                                         float zd = (float)(k-cz); zd *= zd;
00888                                         float weight = exp((xd+yd+zd)/N);
00889                                         idx = k*rnxy + j*rnx + i;
00890                                         d[idx] *=  weight;
00891                                 }
00892                         }
00893                 }
00894         }
00895         image->update();
00896         return  image;
00897 }
00898 
00899 #include <iomanip>
00900 
00901 // int BaldwinWoolfordReconstructor::insert_slice_weights(const Transform& t3d)
00902 // {
00903 //      bool fftodd = image->is_fftodd();
00904 //      int rnx = nx-2*!fftodd;
00905 //
00906 //      float y_scale = 1.0, x_scale = 1.0;
00907 //
00908 //      if ( ny != rnx  )
00909 //      {
00910 //              if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
00911 //              else x_scale = (float) ny / (float) rnx;
00912 //      }
00913 //
00914 //      int tnx = tmp_data->get_xsize();
00915 //      int tny = tmp_data->get_ysize();
00916 //      int tnz = tmp_data->get_zsize();
00917 //
00918 //      vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00919 //      for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00920 //              Transform n3d = t3d*(*it);
00921 //
00922 //              for (int y = 0; y < tny; y++) {
00923 //                      for (int x = 0; x < tnx; x++) {
00924 //
00925 //                              float rx = (float) x;
00926 //                              float ry = (float) y;
00927 //
00928 //                              if ( ny != rnx )
00929 //                              {
00930 //                                      if ( rnx > ny ) ry *= y_scale;
00931 //                                      else rx *= x_scale;
00932 //                              }
00933 // //                           float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
00934 // //                           float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
00935 // //                           float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
00936 //
00937 //                              Vec3f coord(rx,(ry - tny/2),0);
00938 //                              coord = coord*n3d; // transpose multiplication
00939 //                              float xx = coord[0];
00940 //                              float yy = coord[1];
00941 //                              float zz = coord[2];
00942 //
00943 //                              if (xx < 0 ){
00944 //                                      xx = -xx;
00945 //                                      yy = -yy;
00946 //                                      zz = -zz;
00947 //                              }
00948 //
00949 //                              yy += tny/2;
00950 //                              zz += tnz/2;
00951 //                              insert_density_at(xx,yy,zz);
00952 //                      }
00953 //              }
00954 //      }
00955 //
00956 //      return 0;
00957 // }
00958 
00959 void BaldwinWoolfordReconstructor::insert_density_at(const float& x, const float& y, const float& z)
00960 {
00961         int xl = Util::fast_floor(x);
00962         int yl = Util::fast_floor(y);
00963         int zl = Util::fast_floor(z);
00964 
00965         // w is the windowing width
00966         int w = params.set_default("maskwidth",2);
00967         float wsquared = (float) w*w;
00968         float dw = 1.0f/w;
00969         dw *= dw;
00970 
00971         // w minus one - this control the number of
00972         // pixels/voxels to the left of the main pixel
00973         // that will have density
00974         int wmox = w-1;
00975         int wmoy = w-1;
00976         int wmoz = w-1;
00977 
00978         // If any coordinate is incedental with a vertex, then
00979         // make sure there is symmetry in density accruing.
00980         // i.e. the window width must be equal in both directions
00981         if ( ((float) xl) == x ) wmox = w;
00982         if ( ((float) yl) == y ) wmoy = w;
00983         if ( ((float) zl) == z ) wmoz = w;
00984 
00985         float* d = tmp_data->get_data();
00986         int tnx = tmp_data->get_xsize();
00987         int tny = tmp_data->get_ysize();
00988         int tnz = tmp_data->get_zsize();
00989         size_t tnxy = tnx*tny;
00990 
00991         int mode = params.set_default("mode","nearest_neighbor");
00992 
00993         for(int k = zl-wmoz; k <= zl+w; ++k ) {
00994                 for(int j = yl-wmoy; j <= yl+w; ++j) {
00995                         for( int i = xl-wmox; i <= xl+w; ++i) {
00996                                 float fac = 1.0;
00997                                 int ic = i, jc = j, kc = k;
00998 
00999                                 // Fourier space is periodic, which is enforced
01000                                 // by the next 6 if statements. These if statements
01001                                 // assume that the Fourier DC components is at
01002                                 // (0,ny/2,nz/2).
01003                                 if ( i <= 0 ) {
01004 
01005                                         if ( x != 0 && i == 0 ) fac = 1.0;
01006                                         else if ( x == 0 && i < 0) continue;
01007 //                                      if (i < 0 ) ic = -i;
01008                                         if (i < 0 ) {
01009                                                 continue;
01010                                                 ic = -i;
01011                                                 jc = tny-jc;
01012                                                 kc = tnz-kc;
01013                                         }
01014                                 }
01015                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01016                                 if ( jc < 0 ) jc = tny+jc;
01017                                 if ( jc >= tny ) jc = jc-tny;
01018                                 if ( kc < 0 ) kc = tnz+kc;
01019                                 if ( kc >= tnz ) kc = kc-tnz;
01020 //                              if ( ic >= tnx ) continue;
01021 //                              if ( jc < 0 ) continue;
01022 //                              if ( jc >= tny ) continue;
01023 //                              if ( kc < 0 ) continue;
01024 //                              if ( kc >= tnz ) continue;
01025                                 // This shouldn't happen
01026                                 // Debug remove later
01027                                 if ( ic < 0 ) { cout << "wo 1" << endl; }
01028                                 if ( ic >= tnx  ){ cout << "wo 2" << endl; }
01029                                 if ( jc < 0 ) { cout << "wo 3" << endl; }
01030                                 if ( jc >= tny ) { cout << "wo 4" << endl; }
01031                                 if ( kc < 0 ) { cout << "wo 5" << endl; }
01032                                 if ( kc >= tnz ) { cout << "wo 6" << endl; }
01033 
01034 
01035                                 float zd = (z-(float)k);
01036                                 float yd = (y-(float)j);
01037                                 float xd = (x-(float)i);
01038                                 zd *= zd; yd *= yd; xd *= xd;
01039                                 float distsquared = xd+yd+zd;
01040                                 // We enforce a spherical kernel
01041                                 if ( mode == 1 && distsquared > wsquared ) continue;
01042 
01043 //                              float f = fac*exp(-dw*(distsquared));
01044                                 float f = fac*exp(-2.467f*(distsquared));
01045                                 // Debug - this error should never occur.
01046                                 if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in density insertion" );
01047                                 d[kc*tnxy+jc*tnx+ic] += f;
01048                         }
01049                 }
01050         }
01051 }
01052 
01053 int BaldwinWoolfordReconstructor::insert_slice(const EMData* const input_slice, const Transform & t, const float weight)
01054 {
01055         Transform * rotation;
01056         if ( input_slice->has_attr("xform.projection") ) {
01057                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
01058         } else {
01059                 rotation = new Transform(t); // assignment operator
01060         }
01061         Transform tmp(*rotation);
01062         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
01063 
01064         Vec2f trans = tmp.get_trans_2d();
01065         float scale = tmp.get_scale();
01066         bool mirror = tmp.get_mirror();
01067         EMData* slice = 0;
01068         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01069                 slice = input_slice->process("xform",Dict("transform",&tmp));
01070         } else if ( mirror == true ) {
01071                 slice = input_slice->process("xform.flip",Dict("axis","x"));
01072         }
01073         if ( slice == 0 ) {
01074                 slice = input_slice->process("xform.phaseorigin.tocorner");
01075         } else {
01076                 slice->process_inplace("xform.phaseorigin.tocorner");
01077         }
01078 
01079         slice->do_fft_inplace();
01080         slice->process_inplace("xform.fourierorigin.tocenter");
01081         float *dat = slice->get_data();
01082         float dt[2];
01083 
01084         bool fftodd = image->is_fftodd();
01085         int rnx = nx-2*!fftodd;
01086 
01087         float y_scale = 1.0, x_scale = 1.0;
01088 
01089         if ( ny != rnx  )
01090         {
01091                 if ( rnx > ny ) y_scale = (float) rnx / (float) ny;
01092                 else x_scale = (float) ny / (float) rnx;
01093         }
01094 
01095         int tnx = tmp_data->get_xsize();
01096         int tny = tmp_data->get_ysize();
01097         int tnz = tmp_data->get_zsize();
01098 
01099         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01100 //      float weight = params.set_default("weight",1.0f);
01101 
01102         rotation->set_scale(1.0); rotation->set_mirror(false); rotation->set_trans(0,0,0);
01103         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01104                 Transform t3d = (*rotation)*(*it);
01105 
01106                 for (int y = 0; y < tny; y++) {
01107                         for (int x = 0; x < tnx; x++) {
01108                                 float rx = (float) x;
01109                                 float ry = (float) y;
01110 
01111                                 if ( ny != rnx )
01112                                 {
01113                                         if ( rnx > ny ) ry *= y_scale;
01114                                         else rx *= x_scale;
01115                                 }
01116 
01117 //                              float xx = rx * n3d[0][0] + (ry - tny/2) * n3d[1][0];
01118 //                              float yy = rx * n3d[0][1] + (ry - tny/2) * n3d[1][1];
01119 //                              float zz = rx * n3d[0][2] + (ry - tny/2) * n3d[1][2];
01120 
01121                                 Vec3f coord(rx,(ry - tny/2),0);
01122                                 coord = coord*t3d; // transpose multiplication
01123                                 float xx = coord[0];
01124                                 float yy = coord[1];
01125                                 float zz = coord[2];
01126 
01127 
01128                                 float cc = 1;
01129                                 if (xx < 0 ){
01130                                         xx = -xx;
01131                                         yy = -yy;
01132                                         zz = -zz;
01133                                         cc = -1;
01134                                 }
01135 
01136                                 yy += tny/2;
01137                                 zz += tnz/2;
01138 
01139                                 int idx = x * 2 + y * (slice->get_xsize());
01140                                 dt[0] = dat[idx];
01141                                 dt[1] = cc * dat[idx+1];
01142 
01143                                 insert_pixel(xx,yy,zz,dt);
01144                         }
01145                 }
01146         }
01147 
01148         if(rotation) {delete rotation; rotation=0;}
01149         delete slice;
01150 
01151         return 0;
01152 }
01153 
01154 void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
01155 {
01156         int xl = Util::fast_floor(x);
01157         int yl = Util::fast_floor(y);
01158         int zl = Util::fast_floor(z);
01159 
01160         // w is the windowing width
01161         int w = params.set_default("maskwidth",2);
01162         float wsquared = (float) w*w;
01163         float dw = 1.0f/w;
01164         dw *= dw;
01165 
01166         int wmox = w-1;
01167         int wmoy = w-1;
01168         int wmoz = w-1;
01169 
01170         // If any coordinate is incedental with a vertex, then
01171         // make sure there is symmetry in density accruing.
01172         // i.e. the window width must be equal in both directions
01173         if ( ((float) xl) == x ) wmox = w;
01174         if ( ((float) yl) == y ) wmoy = w;
01175         if ( ((float) zl) == z ) wmoz = w;
01176 
01177         float* we = tmp_data->get_data();
01178         int tnx = tmp_data->get_xsize();
01179         int tny = tmp_data->get_ysize();
01180         int tnz = tmp_data->get_zsize();
01181         int tnxy = tnx*tny;
01182 
01183         int rnx = 2*tnx;
01184         int rnxy = 2*tnxy;
01185 
01186         int mode = params.set_default("mode","nearest_neighbor");
01187 
01188         float* d = image->get_data();
01189         for(int k = zl-wmoz; k <= zl+w; ++k ) {
01190                 for(int j = yl-wmoy; j <= yl+w; ++j) {
01191                         for( int i = xl-wmox; i <= xl+w; ++i) {
01192                                 float fac = 1.0;
01193                                 int ic = i, jc = j, kc = k;
01194 
01195                                 // Fourier space is periodic, which is enforced
01196                                 // by the next 6 if statements. These if statements
01197                                 // assume that the Fourier DC component is at
01198                                 // (0,ny/2,nz/2).
01199                                 float negfac=1.0;
01200                                 if ( i <= 0 ) {
01201                                         if ( x != 0 && i == 0 ) fac = 1.0;
01202                                         else if ( x == 0 && i < 0) continue;
01203                                         if (i < 0 ) {
01204                                                 continue;
01205                                                 ic = -i;
01206                                                 jc = tny-jc;
01207                                                 kc = tnz-kc;
01208                                                 negfac=-1.0;
01209                                         }
01210                                 }
01211                                 if ( ic >= tnx ) ic = 2*tnx-ic-1;
01212                                 if ( jc < 0 ) jc = tny+jc;
01213                                 if ( jc >= tny ) jc = jc-tny;
01214                                 if ( kc < 0 ) kc = tnz+kc;
01215                                 if ( kc >= tnz ) kc = kc-tnz;
01216 //                              if ( ic >= tnx ) continue;
01217 //                              if ( jc < 0 ) continue;
01218 //                              if ( jc >= tny ) continue;
01219 //                              if ( kc < 0 ) continue;
01220 //                              if ( kc >= tnz ) continue;
01221 
01222                                 float zd = (z-(float)k);
01223                                 float yd = (y-(float)j);
01224                                 float xd = (x-(float)i);
01225                                 zd *= zd; yd *= yd; xd *= xd;
01226                                 float distsquared = xd+yd+zd;
01227 //                              float f = fac*exp(-dw*(distsquared));
01228                                 float f = fac*exp(-2.467f*(distsquared));
01229                                 float weight = f/we[kc*tnxy+jc*tnx+ic];
01230                                 // debug - this error should never occur
01231                                 if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
01232                                 size_t k = kc*rnxy+jc*rnx+2*ic;
01233 
01234                                 float factor, dist,residual;
01235                                 int sizeW,sizeWmid,idx;
01236                                 switch (mode) {
01237                                         case 0:
01238                                                 d[k] += weight*f*dt[0];
01239                                                 d[k+1] += negfac*weight*f*dt[1];
01240                                                 cout << "hello" << endl;
01241                                         break;
01242 
01243                                         case 1:
01244                                                 // We enforce a spherical kernel
01245                                                 if ( distsquared > wsquared ) continue;
01246 
01247                                                 sizeW = (int)(1+2*w/dfreq);
01248                                                 sizeWmid = sizeW/2;
01249 
01250                                                 dist = sqrtf(distsquared);
01251                                                 idx = (int)(sizeWmid + dist/dfreq);
01252                                                 if (idx >= sizeW) throw InvalidValueException(idx, "idx was greater than or equal to sizeW");
01253                                                 residual = dist/dfreq - (int)(dist/dfreq);
01254                                                 if ( fabs(residual) > 1) throw InvalidValueException(residual, "Residual was too big");
01255 
01256                                                 factor = (W[idx]*(1.0f-residual)+W[idx+1]*residual)*weight;
01257 
01258                                                 d[k] += dt[0]*factor;
01259                                                 d[k+1] += dt[1]*factor;
01260                                         break;
01261 
01262                                         default:
01263                                                 throw InvalidValueException(mode, "The mode was unsupported in BaldwinWoolfordReconstructor::insert_pixel");
01264                                         break;
01265                                 }
01266                         }
01267                 }
01268         }
01269 }
01270 
01271 // void BaldwinWoolfordReconstructor::insert_pixel(const float& x, const float& y, const float& z, const float dt[2])
01272 // {
01273 //      int xl = Util::fast_floor(x);
01274 //      int yl = Util::fast_floor(y);
01275 //      int zl = Util::fast_floor(z);
01276 //
01277 //      // w is the windowing width
01278 //      int w = params.set_default("maskwidth",2);
01279 //      float dw = 1.0/w;
01280 //      dw *= dw;
01281 // //   dw = 2;
01282 // //   cout << w << endl;
01283 //      //      int w = 3;
01284 //      // w minus one - this control the number of
01285 //      // pixels/voxels to the left of the main pixel
01286 //      // that will have density
01287 //      int wmox = w-1;
01288 //      int wmoy = w-1;
01289 //      int wmoz = w-1;
01290 //
01291 //      // If any coordinate is incedental with a vertex, then
01292 //      // make sure there is symmetry in density accruing.
01293 //      // i.e. the window width must be equal in both directions
01294 //      if ( ((float) xl) == x ) wmox = w;
01295 //      if ( ((float) yl) == y ) wmoy = w;
01296 //      if ( ((float) zl) == z ) wmoz = w;
01297 //
01298 //      float* d = tmp_data->get_data();
01299 //      int tnx = tmp_data->get_xsize();
01300 //      int tny = tmp_data->get_ysize();
01301 //      int tnz = tmp_data->get_zsize();
01302 //      int tnxy = tnx*tny;
01303 //
01304 //      float weight = 1.0;
01305 // //
01306 //      for(int k = zl-wmoz; k <= zl+w; ++k ) {
01307 //              for(int j = yl-wmoy; j <= yl+w; ++j) {
01308 //                      for( int i = xl-wmox; i <= xl+w; ++i) {
01309 //                              float fac = 1.0;
01310 //                              int ic = i, jc = j, kc = k;
01311 //
01312 //                              // Fourier space is periodic, which is enforced
01313 //                              // by the next 6 if statements. These if statements
01314 //                              // assume that the Fourier DC components is at
01315 //                              // (0,ny/2,nz/2).
01316 //                              if ( i <= 0 ) {
01317 //                                      if ( x != 0 && i == 0 ) fac = 1.0;
01318 //                                      else if ( x == 0 && i < 0) continue;
01319 // //                                   if (i < 0 ) ic = -i;
01320 //                                      if (i < 0 ) {
01321 //                                              ic = -i;
01322 //                                              jc = tny-jc;
01323 //                                              kc = tnz-kc;
01324 //                                      }
01325 //                              }
01326 //                              if ( ic >= tnx ) ic = 2*tnx-ic-1;
01327 //                              if ( jc < 0 ) jc = tny+jc;
01328 //                              if ( jc >= tny ) jc = jc-tny;
01329 //                              if ( kc < 0 ) kc = tnz+kc;
01330 //                              if ( kc >= tnz ) kc = kc-tnz;
01331 //                              // This shouldn't happen
01332 //                              // Debug remove later
01333 //                              if ( ic < 0 ) { cout << "wo 1" << endl; }
01334 //                              if ( ic >= tnx  ){ cout << "wo 2" << endl; }
01335 //                              if ( jc < 0 ) { cout << "wo 3" << endl; }
01336 //                              if ( jc >= tny ) { cout << "wo 4" << endl; }
01337 //                              if ( kc < 0 ) { cout << "wo 5" << endl; }
01338 //                              if ( kc >= tnz ) { cout << "wo 6" << endl; }
01339 //
01340 //
01341 //                              float zd = (z-(float)k);
01342 //                              float yd = (y-(float)j);
01343 //                              float xd = (x-(float)i);
01344 //                              zd *= zd; yd *= yd; xd *= xd;
01345 //                              // Debug - this error should never occur.
01346 //                              if ( (kc*tnxy+jc*tnx+ic) >= tnxy*tnz ) throw OutofRangeException(0,tnxy*tnz,kc*tnxy+jc*tnx+ic, "in weight determination insertion" );
01347 // //                           float f = fac*exp(-dw*(xd+yd+zd)*0.5);
01348 //                              float f = exp(-2.467*(xd+yd+zd));
01349 //                              weight += f*(d[kc*tnxy+jc*tnx+ic] - f);
01350 //                      }
01351 //              }
01352 //      }
01353 //      weight = 1.0/weight;
01354 //      int rnx = 2*tnx;
01355 //      int rnxy = 2*tnxy;
01356 //      d = image->get_data();
01357 //      for(int k = zl-wmoz; k <= zl+w; ++k ) {
01358 //              for(int j = yl-wmoy; j <= yl+w; ++j) {
01359 //                      for( int i = xl-wmox; i <= xl+w; ++i) {
01360 //                              float fac = 1.0;
01361 //                              int ic = i, jc = j, kc = k;
01362 //
01363 //                              // Fourier space is periodic, which is enforced
01364 //                              // by the next 6 if statements. These if statements
01365 //                              // assume that the Fourier DC components is at
01366 //                              // (0,ny/2,nz/2).
01367 //                              float negfac=1.0;
01368 //                              if ( i <= 0 ) {
01369 //                                      if ( x != 0 && i == 0 ) fac = 1.0;
01370 //                                      else if ( x == 0 && i < 0) continue;
01371 //                                      if (i < 0 ) {
01372 //                                              continue;
01373 //                                              ic = -i;
01374 //                                              jc = tny-jc;
01375 //                                              kc = tnz-kc;
01376 //                                              negfac=-1.0;
01377 //                                      }
01378 //                              }
01379 //                              if ( ic >= tnx ) ic = 2*tnx-ic-1;
01380 //                              if ( jc < 0 ) jc = tny+jc;
01381 //                              if ( jc >= tny ) jc = jc-tny;
01382 //                              if ( kc < 0 ) kc = tnz+kc;
01383 //                              if ( kc >= tnz ) kc = kc-tnz;
01384 //                              // This shouldn't happen
01385 //                              // Debug remove later
01386 //
01387 //
01388 //                              float zd = (z-(float)k);
01389 //                              float yd = (y-(float)j);
01390 //                              float xd = (x-(float)i);
01391 //                              zd *= zd; yd *= yd; xd *= xd;
01392 // //                           float f = fac*exp(-dw*(xd+yd+zd));
01393 //                              float f = exp(-4.934*(xd+yd+zd));
01394 //                              // Debug - this error should never occur.
01395 //                              if ( (kc*rnxy+jc*rnx+2*ic+1) >= rnxy*tnz ) throw OutofRangeException(0,rnxy*tnz,kc*rnxy+jc*rnx+2*ic+1, "in pixel insertion" );
01396 //
01397 //                              d[kc*rnxy+jc*rnx+2*ic] += weight*f*dt[0];
01398 //                              d[kc*rnxy+jc*rnx+2*ic+1] += negfac*weight*f*dt[1];
01399 //                      }
01400 //              }
01401 //      }
01402 // }
01403 */
01404 
01405 void WienerFourierReconstructor::setup()
01406 {
01407         int size = params["size"];
01408         image = new EMData();
01409         image->set_size(size + 2, size, size);
01410         image->set_complex(true);
01411         image->set_ri(true);
01412 
01413         nx = image->get_xsize();
01414         ny = image->get_ysize();
01415         nz = image->get_zsize();
01416 
01417         int n = nx * ny * nz;
01418         float *rdata = image->get_data();
01419 
01420         for (int i = 0; i < n; i += 2) {
01421                 float f = Util::get_frand(0.0, 2.0 * M_PI);
01422                 rdata[i] = 1.0e-10f * sin(f);
01423                 rdata[i + 1] = 1.0e-10f * cos(f);
01424         }
01425         image->update();
01426 
01427         tmp_data = new EMData();
01428         tmp_data->set_size(size + 1, size, size);
01429 }
01430 
01431 
01432 EMData *WienerFourierReconstructor::finish(bool doift)
01433 {
01434         float *norm = tmp_data->get_data();
01435         float *rdata = image->get_data();
01436 
01437         size_t size = nx * ny * nz;
01438         for (size_t i = 0; i < size; i += 2) {
01439                 float d = norm[i];
01440                 if (d == 0) {
01441                         rdata[i] = 0;
01442                         rdata[i + 1] = 0;
01443                 }
01444                 else {
01445                         float w = 1 + 1 / d;
01446                         rdata[i] /= d * w;
01447                         rdata[i + 1] /= d * w;
01448                 }
01449         }
01450 
01451         if( tmp_data ) {
01452                 delete tmp_data;
01453                 tmp_data = 0;
01454         }
01455         image->update();
01456         return image;
01457 }
01458 
01459 
01460 int WienerFourierReconstructor::insert_slice(const EMData* const slice, const Transform & euler, const float weight)
01461 {
01462         if (!slice) {
01463                 LOGERR("try to insert NULL slice");
01464                 return 1;
01465         }
01466 
01467         int mode = params["mode"];
01468         float padratio = params["padratio"];
01469         vector < float >snr = params["snr"];
01470 
01471         if (!slice->is_complex()) {
01472                 LOGERR("Only complex slice can be inserted.");
01473                 return 1;
01474         }
01475         float *gimx = 0;
01476         if (mode == 5) {
01477                 gimx = Interp::get_gimx();
01478         }
01479 
01480         int nxy = nx * ny;
01481         int off[8] = {0,0,0,0,0,0,0,0};
01482         if (mode == 2) {
01483                 off[0] = 0;
01484                 off[1] = 2;
01485                 off[2] = nx;
01486                 off[3] = nx + 2;
01487                 off[4] = nxy;
01488                 off[5] = nxy + 2;
01489                 off[6] = nxy + nx;
01490                 off[7] = nxy + nx + 2;
01491         }
01492 
01493         float *norm = tmp_data->get_data();
01494         float *dat = slice->get_data();
01495         float *rdata = image->get_data();
01496 
01497         int rl = Util::square(ny / 2 - 1);
01498         float dt[2];
01499         float g[8];
01500 
01501         for (int y = 0; y < ny; y++) {
01502                 for (int x = 0; x < nx / 2; x++) {
01503                         if ((x * x + Util::square(y - ny / 2)) >= rl)
01504                         {
01505                                 continue;
01506                         }
01507 
01508 #ifdef _WIN32
01509                         int r = Util::round((float)_hypot(x, (float) y - ny / 2) * Ctf::CTFOS / padratio);
01510 #else
01511                         int r = Util::round((float)hypot(x, (float) y - ny / 2) * Ctf::CTFOS / padratio);
01512 #endif
01513                         if (r >= Ctf::CTFOS * ny / 2) {
01514                                 r = Ctf::CTFOS * ny / 2 - 1;
01515                         }
01516 
01517                         float weight = snr[r];
01518 
01519                         float xx = (x * euler[0][0] + (y - ny / 2) * euler[0][1]);
01520                         float yy = (x * euler[1][0] + (y - ny / 2) * euler[1][1]);
01521                         float zz = (x * euler[2][0] + (y - ny / 2) * euler[2][1]);
01522                         float cc = 1;
01523 
01524                         if (xx < 0) {
01525                                 xx = -xx;
01526                                 yy = -yy;
01527                                 zz = -zz;
01528                                 cc = -1.0;
01529                         }
01530 
01531                         yy += ny / 2;
01532                         zz += nz / 2;
01533 
01534                         dt[0] = dat[x * 2 + y * nx] * (1 + 1.0f / weight);
01535                         dt[1] = cc * dat[x * 2 + 1 + y * nx] * (1 + 1.0f / weight);
01536 
01537                         int x0 = 0;
01538                         int y0 = 0;
01539                         int z0 = 0;
01540                         int i = 0;
01541                         int l = 0;
01542                         float dx = 0;
01543                         float dy = 0;
01544                         float dz = 0;
01545 
01546                         int mx0 = 0;
01547                         int my0 = 0;
01548                         int mz0 = 0;
01549 
01550                         size_t idx;
01551                         switch (mode) {
01552                         case 1:
01553                                 x0 = 2 * (int) floor(xx + 0.5f);
01554                                 y0 = (int) floor(yy + 0.5f);
01555                                 z0 = (int) floor(zz + 0.5f);
01556 
01557                                 rdata[x0 + y0 * nx + z0 * nxy] += weight * dt[0];
01558                                 rdata[x0 + y0 * nx + z0 * nxy + 1] += weight * dt[1];
01559                                 norm[x0 + y0 * nx + z0 * nxy] += weight;
01560                                 break;
01561 
01562                         case 2:
01563                                 x0 = (int) floor(xx);
01564                                 y0 = (int) floor(yy);
01565                                 z0 = (int) floor(zz);
01566 
01567                                 dx = xx - x0;
01568                                 dy = yy - y0;
01569                                 dz = zz - z0;
01570 
01571                                 weight /= (float)pow((float)(EMConsts::I2G * M_PI), 1.5f);
01572 
01573                                 if (x0 > nx - 2 || y0 > ny - 1 || z0 > nz - 1) {
01574                                         break;
01575                                 }
01576 
01577                                 i = (int) (x0 * 2 + y0 * nx + z0 * nxy);
01578 
01579 
01580                                 g[0] = Util::agauss(1, dx, dy, dz, EMConsts::I2G);
01581                                 g[1] = Util::agauss(1, 1 - dx, dy, dz, EMConsts::I2G);
01582                                 g[2] = Util::agauss(1, dx, 1 - dy, dz, EMConsts::I2G);
01583                                 g[3] = Util::agauss(1, 1 - dx, 1 - dy, dz, EMConsts::I2G);
01584                                 g[4] = Util::agauss(1, dx, dy, 1 - dz, EMConsts::I2G);
01585                                 g[5] = Util::agauss(1, 1 - dx, dy, 1 - dz, EMConsts::I2G);
01586                                 g[6] = Util::agauss(1, dx, 1 - dy, 1 - dz, EMConsts::I2G);
01587                                 g[7] = Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, EMConsts::I2G);
01588 
01589                                 for (int j = 0; j < 8; j++) {
01590                                         int k = i + off[j];
01591                                         rdata[k] += weight * dt[0] * g[j];
01592                                         rdata[k + 1] += weight * dt[1] * g[j];
01593                                         norm[k] += weight * g[j];
01594                                 }
01595 
01596                                 break;
01597                         case 3:
01598                                 x0 = 2 * (int) floor(xx + 0.5f);
01599                                 y0 = (int) floor(yy + 0.5f);
01600                                 z0 = (int) floor(zz + 0.5f);
01601 
01602                                 weight /= (float)pow((float)(EMConsts::I3G * M_PI), 1.5f);
01603 
01604                                 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) {
01605                                         break;
01606                                 }
01607 
01608                                 l = x0 - 2;
01609                                 if (x0 == 0) {
01610                                         l = x0;
01611                                 }
01612 
01613                                 for (int k = z0 - 1; k <= z0 + 1; k++) {
01614                                         for (int j = y0 - 1; j <= y0 + 1; j++) {
01615                                                 for (int i = l; i <= x0 + 2; i += 2) {
01616                                                         float r = Util::hypot3sq((float) i / 2 - xx, j - yy, k - zz);
01617                                                         float gg = exp(-r / EMConsts::I3G);
01618 
01619                                                         idx = i + j * nx + k * nxy;
01620                                                         rdata[idx] += weight * gg * dt[0];
01621                                                         rdata[idx + 1] += weight * gg * dt[1];
01622                                                         norm[idx] += weight * gg;
01623                                                 }
01624                                         }
01625                                 }
01626                                 break;
01627 
01628                         case 4:
01629                                 x0 = 2 * (int) floor(xx);
01630                                 y0 = (int) floor(yy);
01631                                 z0 = (int) floor(zz);
01632 
01633                                 weight /= (float)pow((float)(EMConsts::I4G * M_PI), 1.5f);
01634 
01635                                 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2) {
01636                                         break;
01637                                 }
01638 
01639                                 l = x0 - 2;
01640                                 if (x0 == 0) {
01641                                         l = x0;
01642                                 }
01643 
01644                                 for (int k = z0 - 1; k <= z0 + 2; ++k) {
01645                                         for (int j = y0 - 1; j <= y0 + 2; ++j) {
01646                                                 for (int i = l; i <= x0 + 4; i += 2) {
01647                                                         float r = Util::hypot3sq((float) i / 2 - xx, j - yy, k - zz);
01648                                                         float gg = exp(-r / EMConsts::I4G);
01649 
01650                                                         idx = i + j * nx + k * nxy;
01651                                                         rdata[idx] += weight * gg * dt[0];
01652                                                         rdata[idx + 1] += weight * gg * dt[1];
01653                                                         norm[idx] += weight * gg;
01654                                                 }
01655                                         }
01656                                 }
01657                                 break;
01658 
01659                         case 5:
01660                                 x0 = (int) floor(xx + .5);
01661                                 y0 = (int) floor(yy + .5);
01662                                 z0 = (int) floor(zz + .5);
01663 
01664                                 weight /= (float)pow((float)(EMConsts::I5G * M_PI), 1.5f);
01665 
01666                                 mx0 = -(int) floor((xx - x0) * 39.0f + 0.5) - 78;
01667                                 my0 = -(int) floor((yy - y0) * 39.0f + 0.5) - 78;
01668                                 mz0 = -(int) floor((zz - z0) * 39.0f + 0.5) - 78;
01669                                 x0 *= 2;
01670 
01671                                 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01672                                         break;
01673 
01674                                 if (x0 == 0) {
01675                                         l = 0;
01676                                         mx0 += 78;
01677                                 }
01678                                 else if (x0 == 2) {
01679                                         l = 0;
01680                                         mx0 += 39;
01681                                 }
01682                                 else
01683                                         l = x0 - 4;
01684                                 for (int k = z0 - 2, mmz = mz0; k <= z0 + 2; k++, mmz += 39) {
01685                                         for (int j = y0 - 2, mmy = my0; j <= y0 + 2; j++, mmy += 39) {
01686                                                 for (int i = l, mmx = mx0; i <= x0 + 4; i += 2, mmx += 39) {
01687                                                         size_t ii = i + j * nx + k * nxy;
01688                                                         float gg = weight * gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
01689 
01690                                                         rdata[ii] += gg * dt[0];
01691                                                         rdata[ii + 1] += gg * dt[1];
01692                                                         norm[ii] += gg;
01693                                                 }
01694                                         }
01695                                 }
01696 
01697                                 if (x0 <= 2) {
01698                                         xx = -xx;
01699                                         yy = -(yy - ny / 2) + ny / 2;
01700                                         zz = -(zz - nz / 2) + nz / 2;
01701                                         x0 = (int) floor(xx + 0.5f);
01702                                         y0 = (int) floor(yy + 0.5f);
01703                                         z0 = (int) floor(zz + 0.5f);
01704                                         int mx0 = -(int) floor((xx - x0) * 39.0f + .5);
01705                                         x0 *= 2;
01706 
01707                                         if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01708                                                 break;
01709 
01710                                         for (int k = z0 - 2, mmz = mz0; k <= z0 + 2; k++, mmz += 39) {
01711                                                 for (int j = y0 - 2, mmy = my0; j <= y0 + 2; j++, mmy += 39) {
01712                                                         for (int i = 0, mmx = mx0; i <= x0 + 4; i += 2, mmx += 39) {
01713                                                                 size_t ii = i + j * nx + k * nxy;
01714                                                                 float gg =
01715                                                                         weight * gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
01716 
01717                                                                 rdata[ii] += gg * dt[0];
01718                                                                 rdata[ii + 1] -= gg * dt[1];
01719                                                                 norm[ii] += gg;
01720                                                         }
01721                                                 }
01722                                         }
01723                                 }
01724                                 break;
01725                                 // mode 6 is now mode 5 without the fast interpolation
01726                         case 6:
01727                                 x0 = 2 * (int) floor(xx + .5);
01728                                 y0 = (int) floor(yy + .5);
01729                                 z0 = (int) floor(zz + .5);
01730 
01731                                 weight /= (float)pow((float)(EMConsts::I5G * M_PI), 1.5f);
01732 
01733                                 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01734                                         break;
01735 
01736                                 if (x0 <= 2)
01737                                         l = 0;
01738                                 else
01739                                         l = x0 - 4;
01740                                 for (int k = z0 - 2; k <= z0 + 2; ++k) {
01741                                         for (int j = y0 - 2; j <= y0 + 2; ++j) {
01742                                                 for (int i = l; i <= x0 + 4; i += 2) {
01743                                                         size_t ii = i + j * nx + k * nxy;
01744                                                         float r = Util::hypot3sq((float) i / 2 - xx, (float) j - yy,
01745                                                                                                    (float) k - zz);
01746                                                         float gg = weight * exp(-r / EMConsts::I5G);
01747 
01748                                                         rdata[ii] += gg * dt[0];
01749                                                         rdata[ii + 1] += gg * dt[1];
01750                                                         norm[ii] += gg;
01751                                                 }
01752                                         }
01753                                 }
01754 
01755                                 if (x0 <= 2) {
01756                                         xx = -xx;
01757                                         yy = -(yy - ny / 2) + ny / 2;
01758                                         zz = -(zz - nz / 2) + nz / 2;
01759                                         x0 = 2 * (int) floor(xx + 0.5f);
01760                                         y0 = (int) floor(yy + 0.5f);
01761                                         z0 = (int) floor(zz + 0.5f);
01762 
01763                                         if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01764                                                 break;
01765 
01766                                         for (int k = z0 - 2; k <= z0 + 2; ++k) {
01767                                                 for (int j = y0 - 2; j <= y0 + 2; ++j) {
01768                                                         for (int i = 0; i <= x0 + 4; i += 2) {
01769                                                                 size_t ii = i + j * nx + k * nxy;
01770                                                                 float r = Util::hypot3sq((float) i / 2 - xx, (float) j - yy,
01771                                                                                                            (float) k - zz);
01772                                                                 float gg = weight * exp(-r / EMConsts::I5G);
01773 
01774                                                                 rdata[ii] += gg * dt[0];
01775                                                                 rdata[ii + 1] -= gg * dt[1];    // note the -, complex conj.
01776                                                                 norm[ii] += gg;
01777                                                         }
01778                                                 }
01779                                         }
01780                                 }
01781                                 break;
01782 
01783                         case 7:
01784                                 x0 = 2 * (int) floor(xx + .5);
01785                                 y0 = (int) floor(yy + .5);
01786                                 z0 = (int) floor(zz + .5);
01787 
01788                                 if (x0 >= nx - 4 || y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01789                                         break;
01790 
01791                                 if (x0 <= 2)
01792                                         l = 0;
01793                                 else
01794                                         l = x0 - 4;
01795                                 for (int k = z0 - 2; k <= z0 + 2; k++) {
01796                                         for (int j = y0 - 2; j <= y0 + 2; j++) {
01797                                                 for (int i = l; i <= x0 + 4; i += 2) {
01798                                                         size_t ii = i + j * nx + k * nxy;
01799                                                         float r = (float)sqrt(Util::hypot3sq((float) i / 2 - xx,
01800                                                                                                                            (float) j - yy,
01801                                                                                                                            (float) k - zz));
01802                                                         float gg = weight * Interp::hyperg(r);
01803 
01804                                                         rdata[ii] += gg * dt[0];
01805                                                         rdata[ii + 1] += gg * dt[1];
01806                                                         norm[ii] += gg;
01807                                                 }
01808                                         }
01809                                 }
01810 
01811                                 if (x0 <= 2) {
01812                                         xx = -xx;
01813                                         yy = -(yy - ny / 2) + ny / 2;
01814                                         zz = -(zz - nz / 2) + nz / 2;
01815                                         x0 = 2 * (int) floor(xx + .5);
01816                                         y0 = (int) floor(yy + .5);
01817                                         z0 = (int) floor(zz + .5);
01818 
01819                                         if (y0 > ny - 3 || z0 > nz - 3 || y0 < 2 || z0 < 2)
01820                                                 break;
01821 
01822                                         for (int k = z0 - 2; k <= z0 + 2; ++k) {
01823                                                 for (int j = y0 - 2; j <= y0 + 2; ++j) {
01824                                                         for (int i = 0; i <= x0 + 4; i += 2) {
01825                                                                 size_t ii = i + j * nx + k * nxy;
01826                                                                 float r = sqrt(Util::hypot3sq((float) i / 2 - xx, (float) j - yy,
01827                                                                                                                         (float) k - zz));
01828                                                                 float gg = weight * Interp::hyperg(r);
01829 
01830                                                                 rdata[ii] += gg * dt[0];
01831                                                                 rdata[ii + 1] -= gg * dt[1];
01832                                                                 norm[ii] += gg;
01833                                                         }
01834                                                 }
01835                                         }
01836                                 }
01837                                 break;
01838                         }
01839 
01840                 }
01841         }
01842 
01843         image->update();
01844         tmp_data->update();
01845 //      slice->update();
01846 
01847         return 0;
01848 }
01849 
01850 void BackProjectionReconstructor::setup()
01851 {
01852         int size = params["size"];
01853         image = new EMData();
01854         nx = size;
01855         ny = size;
01856         if ( (int) params["zsample"] != 0 ) nz = params["zsample"];
01857         else nz = size;
01858         image->set_size(nx, ny, nz);
01859 }
01860 
01861 EMData* BackProjectionReconstructor::preprocess_slice(const EMData* const slice, const Transform& t)
01862 {
01863 
01864         EMData* return_slice = slice->process("normalize.edgemean");
01865         return_slice->process_inplace("filter.linearfourier");
01866 
01867         Transform tmp(t);
01868         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
01869         Vec2f trans = tmp.get_trans_2d();
01870         float scale = tmp.get_scale();
01871         bool mirror = tmp.get_mirror();
01872         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
01873                 return_slice->transform(tmp);
01874         } else if ( mirror == true ) {
01875                 return_slice = slice->process("xform.flip",Dict("axis","x"));
01876         }
01877 
01878         return return_slice;
01879 }
01880 
01881 int BackProjectionReconstructor::insert_slice(const EMData* const input, const Transform &t, const float sliceweight)
01882 {
01883         if (!input) {
01884                 LOGERR("try to insert NULL slice");
01885                 return 1;
01886         }
01887 
01888         if (input->get_xsize() != input->get_ysize() || input->get_xsize() != nx) {
01889                 LOGERR("tried to insert image that was not correction dimensions");
01890                 return 1;
01891         }
01892 
01893         Transform * transform;
01894         if ( input->has_attr("xform.projection") ) {
01895                 transform = (Transform*) (input->get_attr("xform.projection")); // assignment operator
01896         } else {
01897                 transform = new Transform(t); // assignment operator
01898         }
01899         EMData* slice = preprocess_slice(input, t);
01900 
01901         float weight = params["weight"];
01902         slice->mult(weight);
01903 
01904         EMData *tmp = new EMData();
01905         tmp->set_size(nx, ny, nz);
01906 
01907         float *slice_data = slice->get_data();
01908         float *tmp_data = tmp->get_data();
01909 
01910         size_t nxy = nx * ny;
01911         size_t nxy_size = nxy * sizeof(float);;
01912         for (int i = 0; i < nz; ++i) {
01913                 memcpy(&tmp_data[nxy * i], slice_data, nxy_size);
01914         }
01915 
01916         transform->set_scale(1.0);
01917         transform->set_mirror(false);
01918         transform->set_trans(0,0,0);
01919         transform->invert();
01920 
01921         tmp->transform(*transform);
01922         image->add(*tmp);
01923 
01924         if(transform) {delete transform; transform=0;}
01925         delete tmp;
01926         delete slice;
01927 
01928         return 0;
01929 }
01930 
01931 EMData *BackProjectionReconstructor::finish(bool doift)
01932 {
01933 
01934         Symmetry3D* sym = Factory<Symmetry3D>::get((string)params["sym"]);
01935         vector<Transform> syms = sym->get_syms();
01936 
01937         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01938 
01939                 EMData tmpcopy(*image);
01940                 tmpcopy.transform(*it);
01941                 image->add(tmpcopy);
01942         }
01943 
01944         image->mult(1.0f/(float)sym->get_nsym());
01945         delete sym;
01946         return image;
01947 }
01948 
01949 EMData* EMAN::padfft_slice( const EMData* const slice, const Transform& t, int npad )
01950 {
01951         int nx = slice->get_xsize();
01952         int ny = slice->get_ysize();
01953         int ndim = (ny==1) ? 1 : 2;
01954 
01955         if( ndim==2 && nx!=ny )
01956         {
01957                 // FIXME: What kind of exception should we throw here?
01958                 throw std::runtime_error("Tried to padfft a 2D slice which is not square.");
01959         }
01960 
01961         // process 2D slice or 1D line -- subtract the average outside of the circle, zero-pad, fft extend, and fft
01962         EMData* temp = slice->average_circ_sub();
01963 
01964         Assert( temp != NULL );
01965         EMData* zeropadded = temp->norm_pad( false, npad );
01966         Assert( zeropadded != NULL );
01967         checked_delete( temp );
01968 
01969         zeropadded->do_fft_inplace();
01970         EMData* padfftslice = zeropadded;
01971 
01972         // shift the projection
01973         Vec2f trans = t.get_trans_2d();
01974         float sx = -trans[0];
01975         float sy = -trans[1];
01976         if(sx != 0.0f || sy != 0.0)
01977                 padfftslice->process_inplace("filter.shift", Dict("x_shift", sx, "y_shift", sy, "z_shift", 0.0f));
01978 
01979         int remove = slice->get_attr_default("remove", 0);
01980         padfftslice->set_attr( "remove", remove );
01981 
01982 
01983 
01984         padfftslice->center_origin_fft();
01985         return padfftslice;
01986 }
01987 
01988 nn4Reconstructor::nn4Reconstructor()
01989 {
01990     m_volume = NULL;
01991     m_wptr   = NULL;
01992     m_result = NULL;
01993 }
01994 
01995 nn4Reconstructor::nn4Reconstructor( const string& symmetry, int size, int npad )
01996 {
01997     m_volume = NULL;
01998     m_wptr   = NULL;
01999     m_result = NULL;
02000         setup( symmetry, size, npad );
02001         load_default_settings();
02002         print_params();
02003 }
02004 
02005 nn4Reconstructor::~nn4Reconstructor()
02006 {
02007     if( m_delete_volume )
02008         checked_delete(m_volume);
02009 
02010     if( m_delete_weight )
02011         checked_delete( m_wptr );
02012 
02013     checked_delete( m_result );
02014 }
02015 
02016 enum weighting_method { NONE, ESTIMATE, VORONOI };
02017 
02018 float max2d( int kc, const vector<float>& pow_a )
02019 {
02020         float max = 0.0;
02021         for( int i=-kc; i <= kc; ++i ) {
02022                 for( int j=-kc; j <= kc; ++j ) {
02023                         if( i==0 && j==0 ) continue;
02024                         {
02025                                 int c = 2*kc+1 - std::abs(i) - std::abs(j);
02026                                 max = max + pow_a[c];
02027                         }
02028                 }
02029         }
02030         return max;
02031 }
02032 
02033 float max3d( int kc, const vector<float>& pow_a )
02034 {
02035         float max = 0.0;
02036         for( int i=-kc; i <= kc; ++i ) {
02037                 for( int j=-kc; j <= kc; ++j ) {
02038                         for( int k=-kc; k <= kc; ++k ) {
02039                                 if( i==0 && j==0 && k==0 ) continue;
02040                                 // if( i!=0 )
02041                                 {
02042                                         int c = 3*kc+1 - std::abs(i) - std::abs(j) - std::abs(k);
02043                                         max = max + pow_a[c];
02044                                         // max = max + c * c;
02045                                         // max = max + c;
02046                                 }
02047                         }
02048                 }
02049         }
02050         return max;
02051 }
02052 
02053 
02054 void nn4Reconstructor::setup()
02055 {
02056         int size = params["size"];
02057         int npad = params["npad"];
02058 
02059 
02060         string symmetry;
02061         if( params.has_key("symmetry") ) {
02062                 symmetry = params["symmetry"].to_str();
02063         } else {
02064                 symmetry = "c1";
02065         }
02066 
02067         if( params.has_key("ndim") ) {
02068             m_ndim = params["ndim"];
02069         } else {
02070                 m_ndim = 3;
02071         }
02072 
02073     if( params.has_key( "snr" ) ) {
02074         m_osnr = 1.0f/float( params["snr"] );
02075     } else {
02076         m_osnr = 0.0;
02077     }
02078 
02079         setup( symmetry, size, npad );
02080 }
02081 
02082 void nn4Reconstructor::setup( const string& symmetry, int size, int npad )
02083 {
02084         m_weighting = ESTIMATE;
02085         m_wghta = 0.2f;
02086 
02087         m_symmetry = symmetry;
02088         m_npad = npad;
02089         m_nsym = Transform::get_nsym(m_symmetry);
02090 
02091         m_vnx = size;
02092         m_vny = size;
02093         m_vnz = (m_ndim==3) ? size : 1;
02094 
02095         m_vnxp = size*npad;
02096         m_vnyp = size*npad;
02097         m_vnzp = (m_ndim==3) ? size*npad : 1;
02098 
02099         m_vnxc = m_vnxp/2;
02100         m_vnyc = m_vnyp/2;
02101         m_vnzc = (m_ndim==3) ? m_vnzp/2 : 1;
02102 
02103         buildFFTVolume();
02104         buildNormVolume();
02105 
02106 }
02107 
02108 
02109 void nn4Reconstructor::buildFFTVolume() {
02110         int offset = 2 - m_vnxp%2;
02111 
02112         if( params.has_key("fftvol") ) {
02113                 m_volume = params["fftvol"];
02114                 m_delete_volume = false;
02115         } else {
02116                 m_volume = new EMData();
02117                 m_delete_volume = true;
02118         }
02119 
02120         if( m_volume->get_xsize() != m_vnxp+offset &&
02121             m_volume->get_ysize() != m_vnyp &&
02122             m_volume->get_zsize() != m_vnzp ) {
02123                 m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02124                 m_volume->to_zero();
02125         }
02126         // ----------------------------------------------------------------
02127         // Added by Zhengfan Yang on 03/15/07
02128         // Original author: please check whether my revision is correct and
02129         // other Reconstructor need similiar revision.
02130         if ( m_vnxp % 2 == 0 )  m_volume->set_fftodd(0);
02131         else                    m_volume->set_fftodd(1);
02132         // ----------------------------------------------------------------
02133 
02134         m_volume->set_nxc(m_vnxp/2);
02135         m_volume->set_complex(true);
02136         m_volume->set_ri(true);
02137         m_volume->set_fftpad(true);
02138         m_volume->set_attr("npad", m_npad);
02139         m_volume->set_array_offsets(0,1,1);
02140 }
02141 
02142 void nn4Reconstructor::buildNormVolume() {
02143 
02144         if( params.has_key("weight") ) {
02145                 m_wptr = params["weight"];
02146                 m_delete_weight = false;
02147         } else {
02148                 m_wptr = new EMData();
02149                 m_delete_weight = true;
02150         }
02151 
02152         if( m_wptr->get_xsize() != m_vnxc+1 &&
02153                 m_wptr->get_ysize() != m_vnyp &&
02154                 m_wptr->get_zsize() != m_vnzp ) {
02155                 m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02156                 m_wptr->to_zero();
02157         }
02158 
02159         m_wptr->set_array_offsets(0,1,1);
02160 }
02161 
02162 void printImage( const EMData* line )
02163 {
02164         Assert( line->get_zsize()==1 );
02165 
02166 
02167         int nx = line->get_xsize();
02168         int ny = line->get_ysize();
02169         for( int j=0; j < ny; ++j ) {
02170                 for( int i=0; i < nx; ++i )  printf( "%10.3f ", line->get_value_at(i,j) );
02171                 printf( "\n" );
02172         }
02173 }
02174 
02175 
02176 
02177 int nn4Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight) {
02178         // sanity checks
02179         if (!slice) {
02180                 LOGERR("try to insert NULL slice");
02181                 return 1;
02182         }
02183 
02184         int padffted= slice->get_attr_default( "padffted", 0 );
02185         if( m_ndim==3 ) {
02186                 if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
02187                         // FIXME: Why doesn't this throw an exception?
02188                         LOGERR("Tried to insert a slice that is the wrong size.");
02189                         return 1;
02190                 }
02191         } else {
02192                 Assert( m_ndim==2 );
02193                 if( slice->get_ysize() !=1 ) {
02194                         LOGERR( "for 2D reconstruction, a line is excepted" );
02195                         return 1;
02196                 }
02197         }
02198 
02199         EMData* padfft = NULL;
02200 
02201         if( padffted != 0 ) padfft = new EMData(*slice);
02202         else                padfft = padfft_slice( slice, t,  m_npad );
02203 
02204         int mult= slice->get_attr_default( "mult", 1 );
02205         Assert( mult > 0 );
02206 
02207         if( m_ndim==3 ) {
02208                 insert_padfft_slice( padfft, t, mult );
02209         } else {
02210                 float alpha = padfft->get_attr( "alpha" );
02211                 alpha = alpha/180.0f*M_PI;
02212                 for(int i=0; i < m_vnxc+1; ++i ) {
02213                         float xnew = i*cos(alpha);
02214                         float ynew = -i*sin(alpha);
02215                         float btqr = padfft->get_value_at( 2*i, 0, 0 );
02216                         float btqi = padfft->get_value_at( 2*i+1, 0, 0 );
02217                         if( xnew < 0.0 ) {
02218                                 xnew *= -1;
02219                                 ynew *= -1;
02220                                 btqi *= -1;
02221                         }
02222 
02223                         int ixn = int(xnew+0.5+m_vnxp) - m_vnxp;
02224                         int iyn = int(ynew+0.5+m_vnyp) - m_vnyp;
02225 
02226                         if(iyn < 0 ) iyn += m_vnyp;
02227 
02228                         (*m_volume)( 2*ixn, iyn+1, 1 ) += btqr *float(mult);
02229                         (*m_volume)( 2*ixn+1, iyn+1, 1 ) += btqi * float(mult);
02230                         (*m_wptr)(ixn,iyn+1, 1) += float(mult);
02231                 }
02232 
02233         }
02234         checked_delete( padfft );
02235         return 0;
02236 }
02237 
02238 int nn4Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02239 {
02240         Assert( padfft != NULL );
02241         // insert slice for all symmetry related positions
02242         for (int isym=0; isym < m_nsym; isym++) {
02243                 Transform tsym = t.get_sym(m_symmetry, isym);
02244                 m_volume->nn( m_wptr, padfft, tsym, mult);
02245         }
02246         return 0;
02247 }
02248 
02249 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02250 
02251 void circumference( EMData* win )
02252 {
02253         float *tw = win->get_data();
02254         //  mask and subtract circumference average
02255         int ix = win->get_xsize();
02256         int iy = win->get_ysize();
02257         int iz = win->get_zsize();
02258         int L2 = (ix/2)*(ix/2);
02259         int L2P = (ix/2-1)*(ix/2-1);
02260 
02261         int IP = ix/2+1;
02262         int JP = iy/2+1;
02263         int KP = iz/2+1;
02264 
02265         float  TNR = 0.0f;
02266         size_t m = 0;
02267         for (int k = 1; k <= iz; ++k) {
02268                 for (int j = 1; j <= iy; ++j) {
02269                         for (int i = 1; i <= ix; ++i) {
02270                                 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02271                                 if (LR<=(size_t)L2) {
02272                                         if(LR >= (size_t)L2P && LR <= (size_t)L2) {
02273                                                 TNR += tw(i,j,k);
02274                                                 ++m;
02275                                         }
02276                                 }
02277                         }
02278                 }
02279         }
02280 
02281         TNR /=float(m);
02282         for (int k = 1; k <= iz; ++k) {
02283                 for (int j = 1; j <= iy; ++j) {
02284                         for (int i = 1; i <= ix; ++i) {
02285                                 size_t LR = (k-KP)*(k-KP)+(j-JP)*(j-JP)+(i-IP)*(i-IP);
02286                                 if (LR<=(size_t)L2) tw(i,j,k) -= TNR; else tw(i,j,k) = 0.0f;
02287                         }
02288                 }
02289         }
02290 
02291 }
02292 
02293 
02294 EMData* nn4Reconstructor::finish(bool doift)
02295 {
02296         if( m_ndim==3 ) {
02297                 m_volume->symplane0(m_wptr);
02298         } else {
02299                 for( int i=1; i <= m_vnyp; ++i ) {
02300 
02301                         if( (*m_wptr)(0, i, 1)==0.0 ) {
02302                                 int j = m_vnyp + 1 - i;
02303                                 (*m_wptr)(0, i, 1) = (*m_wptr)(0, j, 1);
02304                                 (*m_volume)(0, i, 1) = (*m_volume)(0, j, 1);
02305                                 (*m_volume)(1, i, 1) = (*m_volume)(1, j, 1);
02306                         }
02307                 }
02308         }
02309 
02310 
02311         int box = 7;
02312         int kc = (box-1)/2;
02313         vector< float > pow_a( m_ndim*kc+1, 1.0 );
02314         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02315         pow_a.back()=0.0;
02316 
02317         float alpha = 0.0;
02318         if( m_ndim==3) {
02319                 int vol = box*box*box;
02320                 float max = max3d( kc, pow_a );
02321                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02322         } else {
02323                 int ara = box*box;
02324                 float max = max2d( kc, pow_a );
02325                 alpha = ( 1.0f - 1.0f/(float)ara ) / max;
02326         }
02327 
02328         int ix,iy,iz;
02329         for (iz = 1; iz <= m_vnzp; iz++) {
02330                 for (iy = 1; iy <= m_vnyp; iy++) {
02331                         for (ix = 0; ix <= m_vnxc; ix++) {
02332                                 if ( (*m_wptr)(ix,iy,iz) > 0) {//(*v) should be treated as complex!!
02333                                         float tmp;
02334                                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+m_osnr);
02335 
02336                                         if( m_weighting == ESTIMATE ) {
02337                                                 int cx = ix;
02338                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02339                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02340                                                 float sum = 0.0;
02341                                                 for( int ii = -kc; ii <= kc; ++ii ) {
02342                                                         int nbrcx = cx + ii;
02343                                                         if( nbrcx >= m_vnxc ) continue;
02344                                                         for( int jj= -kc; jj <= kc; ++jj ) {
02345                                                                 int nbrcy = cy + jj;
02346                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02347 
02348                                                                 int kcz = (m_ndim==3) ? kc : 0;
02349                                                                 for( int kk = -kcz; kk <= kcz; ++kk ) {
02350                                                                         int nbrcz = cz + kk;
02351                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02352                                                                         if( nbrcx < 0 ) {
02353                                                                                 nbrcx = -nbrcx;
02354                                                                                 nbrcy = -nbrcy;
02355                                                                                 nbrcz = -nbrcz;
02356                                                                         }
02357                                                                         int nbrix = nbrcx;
02358                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02359                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02360                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02361                                                                                 int c = m_ndim*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02362                                                                                 sum = sum + pow_a[c];
02363                                                                         }
02364                                                                 }
02365                                                         }
02366                                                 }
02367                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
02368                                                 tmp = tmp * wght;
02369                                         }
02370                                         (*m_volume)(2*ix,iy,iz) *= tmp;
02371                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
02372                                 }
02373                         }
02374                 }
02375         }
02376 
02377         //if(m_ndim==2) printImage( m_volume );
02378 
02379         // back fft
02380         m_volume->do_ift_inplace();
02381 
02382         // EMData* win = m_volume->window_center(m_vnx);
02383         m_volume->depad();
02384         circumference( m_volume );
02385         m_volume->set_array_offsets( 0, 0, 0 );
02386 
02387     m_result = m_volume->copy();
02388         return m_result;
02389 }
02390 #undef  tw
02391 
02392 // Added By Zhengfan Yang on 03/16/07
02393 // Beginning of the addition
02394 // --------------------------------------------------------------------------------
02395 
02396 nnSSNR_Reconstructor::nnSSNR_Reconstructor()
02397 {
02398         m_volume = NULL;
02399         m_wptr   = NULL;
02400         m_wptr2  = NULL;
02401         m_result = NULL;
02402 }
02403 
02404 nnSSNR_Reconstructor::nnSSNR_Reconstructor( const string& symmetry, int size, int npad)
02405 {
02406         m_volume = NULL;
02407         m_wptr   = NULL;
02408         m_wptr2  = NULL;
02409         m_result = NULL;
02410 
02411         setup( symmetry, size, npad );
02412 }
02413 
02414 nnSSNR_Reconstructor::~nnSSNR_Reconstructor()
02415 {
02416         if( m_delete_volume ) checked_delete(m_volume);
02417 
02418         if( m_delete_weight ) checked_delete( m_wptr );
02419 
02420         if( m_delete_weight2 ) checked_delete( m_wptr2 );
02421 
02422         checked_delete( m_result );
02423 }
02424 
02425 void nnSSNR_Reconstructor::setup()
02426 {
02427         int size = params["size"];
02428         int npad = params["npad"];
02429 
02430         string symmetry;
02431         if( params.has_key("symmetry") ) symmetry = params["symmetry"].to_str();
02432         else                             symmetry = "c1";
02433 
02434         setup( symmetry, size, npad );
02435 }
02436 
02437 void nnSSNR_Reconstructor::setup( const string& symmetry, int size, int npad )
02438 {
02439 
02440         m_weighting = ESTIMATE;
02441         m_wghta = 0.2f;
02442         m_wghtb = 0.004f;
02443 
02444         m_symmetry = symmetry;
02445         m_npad = npad;
02446         m_nsym = Transform::get_nsym(m_symmetry);
02447 
02448         m_vnx = size;
02449         m_vny = size;
02450         m_vnz = size;
02451 
02452         m_vnxp = size*npad;
02453         m_vnyp = size*npad;
02454         m_vnzp = size*npad;
02455 
02456         m_vnxc = m_vnxp/2;
02457         m_vnyc = m_vnyp/2;
02458         m_vnzc = m_vnzp/2;
02459 
02460         buildFFTVolume();
02461         buildNormVolume();
02462         buildNorm2Volume();
02463 
02464 }
02465 
02466 void nnSSNR_Reconstructor::buildFFTVolume() {
02467 
02468         if( params.has_key("fftvol") ) {
02469                 m_volume = params["fftvol"]; /* volume should be defined in python when PMI is turned on*/
02470                 m_delete_volume = false;
02471         } else {
02472                 m_volume = new EMData();
02473                 m_delete_volume = true;
02474         }
02475         m_volume->set_size(m_vnxp+ 2 - m_vnxp%2,m_vnyp,m_vnzp);
02476         m_volume->to_zero();
02477         if ( m_vnxp % 2 == 0 ) m_volume->set_fftodd(0);
02478         else                   m_volume->set_fftodd(1);
02479 
02480         m_volume->set_nxc(m_vnxc);
02481         m_volume->set_complex(true);
02482         m_volume->set_ri(true);
02483         m_volume->set_fftpad(true);
02484         m_volume->set_attr("npad", m_npad);
02485         m_volume->set_array_offsets(0,1,1);
02486 }
02487 
02488 void nnSSNR_Reconstructor::buildNormVolume() {
02489         if( params.has_key("weight") ) {
02490                 m_wptr          = params["weight"];
02491                 m_delete_weight = false;
02492         } else {
02493                 m_wptr = new EMData();
02494                 m_delete_weight = true;
02495         }
02496 
02497         m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02498         m_wptr->to_zero();
02499 
02500         m_wptr->set_array_offsets(0,1,1);
02501 }
02502 
02503 void nnSSNR_Reconstructor::buildNorm2Volume() {
02504 
02505         if( params.has_key("weight2") ) {
02506                 m_wptr2          = params["weight2"];
02507                 m_delete_weight2 = false;
02508         } else {
02509                 m_wptr2 = new EMData();
02510                 m_delete_weight2 = true;
02511         }
02512         m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02513         m_wptr2->to_zero();
02514         m_wptr2->set_array_offsets(0,1,1);
02515 }
02516 
02517 
02518 int nnSSNR_Reconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight) {
02519         // sanity checks
02520         if (!slice) {
02521                 LOGERR("try to insert NULL slice");
02522                 return 1;
02523         }
02524 
02525         int padffted=slice->get_attr_default( "padffted", 0 );
02526 
02527         if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
02528                 // FIXME: Why doesn't this throw an exception?
02529                 LOGERR("Tried to insert a slice that has wrong size.");
02530                 return 1;
02531         }
02532 
02533         EMData* padfft = NULL;
02534 
02535         if( padffted != 0 ) padfft = new EMData(*slice);
02536         else                padfft = padfft_slice( slice, t, m_npad );
02537 
02538         int mult = slice->get_attr_default("mult", 1);
02539 
02540         Assert( mult > 0 );
02541         insert_padfft_slice( padfft, t, mult );
02542 
02543         checked_delete( padfft );
02544         return 0;
02545 }
02546 
02547 int nnSSNR_Reconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02548 {
02549         Assert( padfft != NULL );
02550         // insert slice for all symmetry related positions
02551         for (int isym=0; isym < m_nsym; isym++) {
02552                 Transform tsym = t.get_sym(m_symmetry, isym);
02553                 m_volume->nn_SSNR( m_wptr, m_wptr2, padfft, tsym, mult);
02554         }
02555         return 0;
02556 }
02557 
02558 
02559 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02560 EMData* nnSSNR_Reconstructor::finish(bool doift)
02561 {
02562 /*
02563   I changed the code on 05/15 so it only returns variance.
02564   Lines commented out are marked by //#
02565   The version prior to the currect changes is r1.190
02566   PAP
02567 */
02568         int kz, ky;
02569         //#int iix, iiy, iiz;
02570         int box = 7;
02571         int kc = (box-1)/2;
02572         float alpha = 0.0;
02573         float argx, argy, argz;
02574         vector< float > pow_a( 3*kc+1, 1.0 );
02575         float w = params["w"];
02576         EMData* SSNR = params["SSNR"];
02577         //#EMData* vol_ssnr = new EMData();
02578         //#vol_ssnr->set_size(m_vnxp, m_vnyp, m_vnzp);
02579         //#vol_ssnr->to_zero();
02580         //#  new line follow
02581         EMData* vol_ssnr = new EMData();
02582         vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
02583         vol_ssnr->to_zero();
02584         if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
02585         else                   vol_ssnr->set_fftodd(1);
02586         vol_ssnr->set_nxc(m_vnxc);
02587         vol_ssnr->set_complex(true);
02588         vol_ssnr->set_ri(true);
02589         vol_ssnr->set_fftpad(false);
02590         //#
02591 
02592         float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
02593         float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
02594 #ifdef _WIN32
02595         float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
02596         int  inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
02597 #else
02598         float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
02599         int  inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
02600 #endif  //_WIN32
02601         SSNR->set_size(inc+1,4,1);
02602 
02603         float *nom    = new float[inc+1];
02604         float *denom  = new float[inc+1];
02605         int  *nn     = new int[inc+1];
02606         int  *ka     = new int[inc+1];
02607         float wght = 1.0f;
02608         for (int i = 0; i <= inc; i++) {
02609                 nom[i] = 0.0f;
02610                 denom[i] = 0.0f;
02611                 nn[i] = 0;
02612                 ka[i] = 0;
02613         }
02614 
02615         m_volume->symplane1(m_wptr, m_wptr2);
02616 
02617         if ( m_weighting == ESTIMATE ) {
02618                 int vol = box*box*box;
02619                 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02620                 pow_a[3*kc] = 0.0;
02621                 float max = max3d( kc, pow_a );
02622                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02623         }
02624 
02625         for (int iz = 1; iz <= m_vnzp; iz++) {
02626                 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
02627                 argz = float(kz*kz)*dz2;
02628                 for (int iy = 1; iy <= m_vnyp; iy++) {
02629                         if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
02630                         argy = argz + float(ky*ky)*dy2;
02631                         for (int ix = 0; ix <= m_vnxc; ix++) {
02632                                 float Kn = (*m_wptr)(ix,iy,iz);
02633                                 argx = std::sqrt(argy + float(ix*ix)*dx2);
02634                                 int r = Util::round(float(inc)*argx);
02635                                 if ( r >= 0 && Kn > 4.5f ) {
02636                                         if ( m_weighting == ESTIMATE ) {
02637                                                 int cx = ix;
02638                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02639                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02640 
02641                                                 float sum = 0.0;
02642                                                 for( int ii = -kc; ii <= kc; ++ii ) {
02643                                                         int nbrcx = cx + ii;
02644                                                         if( nbrcx >= m_vnxc ) continue;
02645                                                     for ( int jj= -kc; jj <= kc; ++jj ) {
02646                                                                 int nbrcy = cy + jj;
02647                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02648                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
02649                                                                         int nbrcz = cz + jj;
02650                                                                         if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02651                                                                         if( nbrcx < 0 ) {
02652                                                                                 nbrcx = -nbrcx;
02653                                                                                 nbrcy = -nbrcy;
02654                                                                                 nbrcz = -nbrcz;
02655                                                                         }
02656                                                                         int nbrix = nbrcx;
02657                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02658                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02659                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
02660                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02661                                                                                 sum = sum + pow_a[c];
02662                                                                         }
02663                                                                 }
02664                                                         }
02665                                                 }
02666                                                 wght = 1.0f / ( 1.0f - alpha * sum );
02667                                         } // end of ( m_weighting == ESTIMATE )
02668                                         float nominator = std::norm(m_volume->cmplx(ix,iy,iz))/Kn;
02669                                         float denominator = ((*m_wptr2)(ix,iy,iz)-nominator)/(Kn-1.0f);
02670                                         // Skip Friedel related values
02671                                         if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
02672                                                 if ( r <= inc ) {
02673                                                         nom[r]   += nominator*wght;
02674                                                         denom[r] += denominator/Kn*wght;
02675                                                         nn[r]    += 2;
02676                                                         ka[r]    += int(Kn);
02677                                                 }
02678 /*
02679 #ifdef  _WIN32
02680                                                 //#float  tmp = _cpp_max(nominator/denominator/Kn-1.0f,0.0f);
02681 #else
02682                                                 //#float  tmp = std::max(nominator/denominator/Kn-1.0f,0.0f);
02683 #endif  //_WIN32
02684                                                 //  Create SSNR as a 3D array (-n/2:n/2+n%2-1)
02685                                                 iix = m_vnxc + ix; iiy = m_vnyc + ky; iiz = m_vnzc + kz;
02686                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
02687                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
02688                                                 // Friedel part
02689                                                 iix = m_vnxc - ix; iiy = m_vnyc - ky; iiz = m_vnzc - kz;
02690                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
02691                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
02692 */
02693 
02694                                         }
02695                                         (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
02696                                         //(*vol_ssnr)(2*ix, iy-1, iz-1) =  real(m_volume->cmplx(ix,iy,iz))*wght/Kn;
02697                                         //(*vol_ssnr)(2*ix+1, iy-1, iz-1) = imag(m_volume->cmplx(ix,iy,iz))*wght/Kn;
02698                                 } // end of Kn>4.5
02699                         }
02700                 }
02701         }
02702 
02703         for (int i = 0; i <= inc; i++)  {
02704                 (*SSNR)(i,0,0) = nom[i];  
02705                 (*SSNR)(i,1,0) = denom[i];    // variance
02706                 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
02707                 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
02708         }
02709         vol_ssnr->update();
02710         return vol_ssnr;
02711 }
02712 #undef  tw
02713 
02714 // -----------------------------------------------------------------------------------
02715 // End of this addition
02716 
02717 //####################################################################################
02718 //** nn4 ctf reconstructor
02719 
02720 nn4_ctfReconstructor::nn4_ctfReconstructor()
02721 {
02722         m_volume  = NULL;
02723         m_wptr    = NULL;
02724         m_result  = NULL;
02725 }
02726 
02727 nn4_ctfReconstructor::nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign )
02728 {
02729         setup( symmetry, size, npad, snr, sign );
02730 }
02731 
02732 nn4_ctfReconstructor::~nn4_ctfReconstructor()
02733 {
02734         if( m_delete_volume ) checked_delete(m_volume);
02735 
02736         if( m_delete_weight ) checked_delete( m_wptr );
02737 
02738         checked_delete( m_result );
02739 }
02740 
02741 void nn4_ctfReconstructor::setup()
02742 {
02743         if( ! params.has_key("size") ) throw std::logic_error("Error: image size is not given");
02744 
02745         int size = params["size"];
02746         int npad = params.has_key("npad") ? int(params["npad"]) : 4;
02747         // int sign = params.has_key("sign") ? int(params["sign"]) : 1;
02748         int sign = 1;
02749         string symmetry = params.has_key("symmetry")? params["symmetry"].to_str() : "c1";
02750 
02751         float snr = params["snr"];
02752 
02753     m_varsnr = params.has_key("varsnr") ? int(params["varsnr"]) : 0;
02754         setup( symmetry, size, npad, snr, sign );
02755 
02756 }
02757 
02758 void nn4_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
02759 {
02760         m_weighting = ESTIMATE;
02761         if( params.has_key("weighting") ) {
02762                 int tmp = int( params["weighting"] );
02763                 if( tmp==0 ) m_weighting = NONE;
02764         }
02765 
02766 
02767 
02768         m_wghta = 0.2f;
02769         m_wghtb = 0.004f;
02770 
02771         m_symmetry = symmetry;
02772         m_npad = npad;
02773         m_sign = sign;
02774         m_nsym = Transform::get_nsym(m_symmetry);
02775 
02776         m_snr = snr;
02777 
02778         m_vnx = size;
02779         m_vny = size;
02780         m_vnz = size;
02781 
02782         m_vnxp = size*npad;
02783         m_vnyp = size*npad;
02784         m_vnzp = size*npad;
02785 
02786         m_vnxc = m_vnxp/2;
02787         m_vnyc = m_vnyp/2;
02788         m_vnzc = m_vnzp/2;
02789 
02790         buildFFTVolume();
02791         buildNormVolume();
02792 }
02793 
02794 void nn4_ctfReconstructor::buildFFTVolume() {
02795         int offset = 2 - m_vnxp%2;
02796         if( params.has_key("fftvol") ) {
02797             m_volume = params["fftvol"];
02798             m_delete_volume = false;
02799         } else {
02800             m_volume = new EMData();
02801             m_delete_volume = true;
02802         }
02803 
02804         if( m_volume->get_xsize() != m_vnxp+offset &&
02805             m_volume->get_ysize() != m_vnyp &&
02806             m_volume->get_zsize() != m_vnzp )
02807         {
02808             m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
02809             m_volume->to_zero();
02810         }
02811 
02812         m_volume->set_nxc(m_vnxp/2);
02813         m_volume->set_complex(true);
02814         m_volume->set_ri(true);
02815         m_volume->set_fftpad(true);
02816         m_volume->set_attr("npad", m_npad);
02817         m_volume->set_array_offsets(0,1,1);
02818 }
02819 
02820 void nn4_ctfReconstructor::buildNormVolume()
02821 {
02822         if( params.has_key("weight") )
02823         {
02824             m_wptr = params["weight"];
02825             m_delete_weight = false;
02826         } else {
02827             m_wptr = new EMData();
02828             m_delete_weight = true;
02829         }
02830 
02831         if( m_wptr->get_xsize() != m_vnxc+1 &&
02832             m_wptr->get_ysize() != m_vnyp &&
02833             m_wptr->get_zsize() != m_vnzp )
02834         {
02835             m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
02836             m_wptr->to_zero();
02837         }
02838 
02839         m_wptr->set_array_offsets(0,1,1);
02840 
02841 }
02842 
02843 int nn4_ctfReconstructor::insert_slice(const EMData* const slice, const Transform& t, const float weight)
02844 {
02845         // sanity checks
02846         if (!slice) {
02847                 LOGERR("try to insert NULL slice");
02848                 return 1;
02849         }
02850 
02851         int buffed = slice->get_attr_default( "buffed", 0 );
02852         if( buffed > 0 )
02853         {
02854             int mult = slice->get_attr_default( "mult", 1 );
02855             insert_buffed_slice( slice, mult );
02856             return 0;
02857         }
02858 
02859         int padffted= slice->get_attr_default("padffted", 0);
02860         if( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  )
02861         {
02862                 // FIXME: Why doesn't this throw an exception?
02863                 LOGERR("Tried to insert a slice that is the wrong size.");
02864                 return 1;
02865         }
02866 
02867         EMData* padfft = NULL;
02868 
02869         if( padffted != 0 ) padfft = new EMData(*slice);
02870         else                padfft = padfft_slice( slice, t, m_npad );
02871 
02872         int mult= slice->get_attr_default("mult", 1);
02873 
02874         Assert( mult > 0 );
02875         insert_padfft_slice( padfft, t, mult );
02876 
02877         checked_delete( padfft );
02878 
02879         return 0;
02880 }
02881 
02882 int nn4_ctfReconstructor::insert_buffed_slice( const EMData* buffed, int mult )
02883 {
02884         const float* bufdata = buffed->get_data();
02885         float* cdata = m_volume->get_data();
02886         float* wdata = m_wptr->get_data();
02887 
02888         int npoint = buffed->get_xsize()/4;
02889         for( int i=0; i < npoint; ++i ) {
02890 
02891                int pos2 = int( bufdata[4*i] );
02892                int pos1 = pos2 * 2;
02893                cdata[pos1  ] += bufdata[4*i+1]*mult;
02894                cdata[pos1+1] += bufdata[4*i+2]*mult;
02895                wdata[pos2  ] += bufdata[4*i+3]*mult;
02896 /*
02897         std::cout << "pos1, pos2, ctfv1, ctfv2, ctf2: ";
02898         std::cout << pos1 << " " << bufdata[5*i+1] << " " << bufdata[5*i+2] << " ";
02899         std::cout << pos2 << " " << bufdata[5*i+4] << std::endl;
02900  */
02901         }
02902         return 0;
02903 }
02904 
02905 int nn4_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
02906 {
02907         Assert( padfft != NULL );
02908         float tmp = padfft->get_attr("ctf_applied");
02909         int   ctf_applied = (int) tmp;
02910 
02911         for( int isym=0; isym < m_nsym; isym++) {
02912                 Transform tsym = t.get_sym( m_symmetry, isym );
02913 
02914                 if(ctf_applied) m_volume->nn_ctf_applied(m_wptr, padfft, tsym, mult);
02915                 else            m_volume->nn_ctf(m_wptr, padfft, tsym, mult);
02916         }
02917 
02918         return 0;
02919 
02920 }
02921 
02922 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
02923 EMData* nn4_ctfReconstructor::finish(bool doift)
02924 {
02925         m_volume->set_array_offsets(0, 1, 1);
02926         m_wptr->set_array_offsets(0, 1, 1);
02927         m_volume->symplane0_ctf(m_wptr);
02928 
02929         int box = 7;
02930         int vol = box*box*box;
02931         int kc = (box-1)/2;
02932         vector< float > pow_a( 3*kc+1, 1.0 );
02933         for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
02934         pow_a[3*kc]=0.0;
02935 
02936 
02937         float max = max3d( kc, pow_a );
02938         float alpha = ( 1.0f - 1.0f/(float)vol ) / max;
02939         float osnr = 1.0f/m_snr;
02940 
02941         // normalize
02942         int ix,iy,iz;
02943         for (iz = 1; iz <= m_vnzp; iz++) {
02944                 for (iy = 1; iy <= m_vnyp; iy++) {
02945                         for (ix = 0; ix <= m_vnxc; ix++) {
02946                                 if ( (*m_wptr)(ix,iy,iz) > 0.0f) {//(*v) should be treated as complex!!
02947                     int iyp = (iy<=m_vnyc) ? iy - 1 : iy-m_vnyp-1;
02948                     int izp = (iz<=m_vnzc) ? iz - 1 : iz-m_vnzp-1;
02949                     float tmp=0.0;
02950                     if( m_varsnr )
02951                     {
02952                                             float freq = sqrt( (float)(ix*ix+iyp*iyp+izp*izp) );
02953                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+freq*osnr)*m_sign;
02954                     }
02955                     else
02956                     {
02957                         tmp = (-2*((ix+iy+iz)%2)+1)/((*m_wptr)(ix,iy,iz)+osnr)*m_sign;
02958                     }
02959 
02960                                         if( m_weighting == ESTIMATE ) {
02961                                                 int cx = ix;
02962                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
02963                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
02964                                                 float sum = 0.0;
02965                                                 for( int ii = -kc; ii <= kc; ++ii ) {
02966                                                         int nbrcx = cx + ii;
02967                                                         if( nbrcx >= m_vnxc ) continue;
02968                                                         for( int jj= -kc; jj <= kc; ++jj ) {
02969                                                                 int nbrcy = cy + jj;
02970                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
02971                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
02972                                                                         int nbrcz = cz + jj;
02973                                                                         if( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
02974                                                                         if( nbrcx < 0 ) {
02975                                                                                 nbrcx = -nbrcx;
02976                                                                                 nbrcy = -nbrcy;
02977                                                                                 nbrcz = -nbrcz;
02978                                                                         }
02979 
02980                                                                         int nbrix = nbrcx;
02981                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
02982                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
02983                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0.0 ) {
02984                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
02985                                                                                 sum = sum + pow_a[c];
02986                                                                                   // if(ix%20==0 && iy%20==0 && iz%20==0)
02987                                                                                  //   std::cout << boost::format( "%4d %4d %4d %4d %10.3f" ) % nbrix % nbriy % nbriz % c % sum << std::endl;
02988                                                                         }
02989                                                                 }
02990                                                         }
02991                                                 }
02992                                                 float wght = 1.0f / ( 1.0f - alpha * sum );
02993 /*
02994                         if(ix%10==0 && iy%10==0)
02995                         {
02996                             std::cout << boost::format( "%4d %4d %4d " ) % ix % iy %iz;
02997                             std::cout << boost::format( "%10.3f %10.3f %10.3f " )  % tmp % wght % sum;
02998                             std::  << boost::format( "%10.3f %10.3e " ) % pow_b[r] % alpha;
02999                             std::cout << std::endl;
03000                         }
03001  */
03002                                                 tmp = tmp * wght;
03003                                         }
03004                                         (*m_volume)(2*ix,iy,iz) *= tmp;
03005                                         (*m_volume)(2*ix+1,iy,iz) *= tmp;
03006                                 }
03007                         }
03008                 }
03009         }
03010 
03011     // back fft
03012     m_volume->do_ift_inplace();
03013     m_volume->depad();
03014     circumference( m_volume );
03015     m_volume->set_array_offsets( 0, 0, 0 );
03016     m_result = m_volume->copy();
03017     return m_result;
03018 }
03019 #undef  tw
03020 
03021 
03022 
03023 
03024 // Added By Zhengfan Yang on 04/11/07
03025 // Beginning of the addition
03026 // --------------------------------------------------------------------------------
03027 
03028 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor()
03029 {
03030     m_volume  = NULL;
03031     m_wptr    = NULL;
03032     m_wptr2   = NULL;
03033     m_wptr3   = NULL;
03034     m_result  = NULL;
03035 }
03036 
03037 nnSSNR_ctfReconstructor::nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign)
03038 {
03039     m_volume  = NULL;
03040     m_wptr    = NULL;
03041     m_wptr2   = NULL;
03042     m_wptr3   = NULL;
03043     m_result  = NULL;
03044 
03045     setup( symmetry, size, npad, snr, sign );
03046 }
03047 
03048 nnSSNR_ctfReconstructor::~nnSSNR_ctfReconstructor()
03049 {
03050 
03051    if( m_delete_volume )
03052         checked_delete(m_volume);
03053    if( m_delete_weight )
03054         checked_delete( m_wptr );
03055    if ( m_delete_weight2 )
03056         checked_delete( m_wptr2 );
03057    if ( m_delete_weight3 )
03058         checked_delete( m_wptr3 );
03059    checked_delete( m_result );
03060 }
03061 
03062 void nnSSNR_ctfReconstructor::setup()
03063 {
03064     int  size = params["size"];
03065     int  npad = params["npad"];
03066     int  sign = params["sign"];
03067     float snr = params["snr"];
03068     string symmetry;
03069     if( params.has_key("symmetry") ) {
03070                 symmetry = params["symmetry"].to_str();
03071     } else {
03072                 symmetry = "c1";
03073     }
03074     setup( symmetry, size, npad, snr, sign );
03075 }
03076 void nnSSNR_ctfReconstructor::setup( const string& symmetry, int size, int npad, float snr, int sign )
03077 {
03078 
03079     m_weighting = ESTIMATE;
03080     m_wghta     = 0.2f;
03081     m_wghtb     = 0.004f;
03082     wiener      = 1;
03083 
03084     m_symmetry  = symmetry;
03085     m_npad      = npad;
03086     m_nsym      = Transform::get_nsym(m_symmetry);
03087 
03088     m_sign      = sign;
03089     m_snr       = snr;
03090 
03091     m_vnx       = size;
03092     m_vny       = size;
03093     m_vnz       = size;
03094 
03095     m_vnxp      = size*npad;
03096     m_vnyp      = size*npad;
03097     m_vnzp      = size*npad;
03098 
03099     m_vnxc      = m_vnxp/2;
03100     m_vnyc      = m_vnyp/2;
03101     m_vnzc      = m_vnzp/2;
03102 
03103     buildFFTVolume();
03104     buildNormVolume();
03105     buildNorm2Volume();
03106     buildNorm3Volume();
03107 }
03108 
03109 void nnSSNR_ctfReconstructor::buildFFTVolume() {
03110 
03111         int offset = 2 - m_vnxp%2;
03112         if( params.has_key("fftvol") ) {
03113                 m_volume = params["fftvol"]; /* volume should be defined in python when PMI is turned on*/
03114                 m_delete_volume = false;
03115     } else {
03116                 m_volume = new EMData();
03117                 m_delete_volume = true;
03118         }
03119 
03120         m_volume->set_size(m_vnxp+offset,m_vnyp,m_vnzp);
03121         m_volume->to_zero();
03122 
03123         if ( m_vnxp % 2 == 0 ) {
03124                 m_volume->set_fftodd(0);
03125         } else {
03126                 m_volume->set_fftodd(1);
03127         }
03128 
03129         m_volume->set_nxc(m_vnxc);
03130         m_volume->set_complex(true);
03131         m_volume->set_ri(true); //(real, imaginary) instead of polar coordinate
03132         m_volume->set_fftpad(true);
03133         m_volume->set_attr("npad", m_npad);
03134         m_volume->set_array_offsets(0,1,1);
03135 }
03136 
03137 
03138 void nnSSNR_ctfReconstructor::buildNormVolume()
03139 {
03140         if( params.has_key("weight") ) {
03141                  m_wptr          = params["weight"];
03142                  m_delete_weight = false;
03143         } else {
03144                 m_wptr = new EMData();
03145                 m_delete_weight = true;
03146         }
03147         m_wptr->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03148         m_wptr->to_zero();
03149         m_wptr->set_array_offsets(0,1,1);
03150 }
03151 
03152 void nnSSNR_ctfReconstructor::buildNorm2Volume() {
03153 
03154         if( params.has_key("weight2") ) {
03155                 m_wptr2          = params["weight2"];
03156                 m_delete_weight2 = false;
03157         } else {
03158                 m_wptr2 = new EMData();
03159                 m_delete_weight2 = true;
03160         }
03161         m_wptr2->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03162         m_wptr2->to_zero();
03163         m_wptr2->set_array_offsets(0,1,1);
03164 }
03165 
03166 void nnSSNR_ctfReconstructor::buildNorm3Volume() {
03167 
03168         if( params.has_key("weight3") ) {
03169                 m_wptr3          = params["weight3"];
03170                 m_delete_weight3 = false;
03171         } else {
03172                 m_wptr3 = new EMData();
03173                 m_delete_weight3 = true;
03174         }
03175         m_wptr3->set_size(m_vnxc+1,m_vnyp,m_vnzp);
03176         m_wptr3->to_zero();
03177         m_wptr3->set_array_offsets(0,1,1);
03178 }
03179 
03180 int nnSSNR_ctfReconstructor::insert_slice(const EMData *const  slice, const Transform& t, const float weight) {
03181         // sanity checks
03182         if (!slice) {
03183                 LOGERR("try to insert NULL slice");
03184                 return 1;
03185         }
03186         int padffted= slice->get_attr_default("padffted", 0);
03187         if ( padffted==0 && (slice->get_xsize()!=slice->get_ysize() || slice->get_xsize()!=m_vnx)  ) {
03188                 // FIXME: Why doesn't this throw an exception?
03189                 LOGERR("Tried to insert a slice that is the wrong size.");
03190                 return 1;
03191         }
03192     EMData* padfft = NULL;
03193 
03194     if( padffted != 0 ) {
03195         padfft = new EMData(*slice);
03196     } else {
03197         padfft = padfft_slice( slice, t, m_npad );
03198     }
03199 
03200     int mult= slice->get_attr_default("mult", 1);
03201 
03202         Assert( mult > 0 );
03203         insert_padfft_slice( padfft, t, mult );
03204 
03205         checked_delete( padfft );
03206         return 0;
03207 }
03208 int nnSSNR_ctfReconstructor::insert_padfft_slice( EMData* padfft, const Transform& t, int mult )
03209 {
03210         Assert( padfft != NULL );
03211 
03212         // insert slice for all symmetry related positions
03213         for (int isym=0; isym < m_nsym; isym++) {
03214                 Transform tsym = t.get_sym(m_symmetry, isym);
03215                 m_volume->nn_SSNR_ctf(m_wptr, m_wptr2, m_wptr3, padfft, tsym, mult);
03216         }
03217 
03218         return 0;
03219 }
03220 
03221 #define  tw(i,j,k)      tw[ i-1 + (j-1+(k-1)*iy)*ix ]
03222 EMData* nnSSNR_ctfReconstructor::finish(bool doift)
03223 {
03224 /*
03225   I changed the code on 05/15 so it only returns variance.
03226   Lines commented out are marked by //#
03227   The version prior to the currect changes is r1.190
03228   PAP
03229 */
03230         /***
03231             m_volume ctf*(P^2D->3D(F^3D))
03232             m_wptr   ctf^2
03233             m_wptr2  |P^2D->3D(F^3D)|^2
03234             m_wptr3  Kn
03235             nominator = sum_rot [ wght*signal ]
03236             denominator  = sum_rot[ wght*variance ]
03237                                                       ***/
03238         int kz, ky;
03239         int box = 7;
03240         int kc  = (box-1)/2;
03241         float alpha = 0.0;
03242         float argx, argy, argz;
03243         vector< float > pow_a( 3*kc+1, 1.0 );
03244         float w = params["w"];
03245         float dx2 = 1.0f/float(m_vnxc)/float(m_vnxc);
03246         float dy2 = 1.0f/float(m_vnyc)/float(m_vnyc);
03247 #ifdef  _WIN32
03248         float dz2 = 1.0f/_cpp_max(float(m_vnzc),1.0f)/_cpp_max(float(m_vnzc),1.0f);
03249         int inc = Util::round(float(_cpp_max(_cpp_max(m_vnxc,m_vnyc),m_vnzc))/w);
03250 #else
03251         float dz2 = 1.0f/std::max(float(m_vnzc),1.0f)/std::max(float(m_vnzc),1.0f);
03252         int inc = Util::round(float(std::max(std::max(m_vnxc,m_vnyc),m_vnzc))/w);
03253 #endif  //_WIN32
03254 
03255         EMData* SSNR = params["SSNR"];
03256         SSNR->set_size(inc+1,4,1);
03257         //#EMData* vol_ssnr = new EMData();
03258         //#vol_ssnr->set_size(m_vnxp, m_vnyp, m_vnzp);
03259         //#vol_ssnr->to_zero();
03260         //#  new linea follow
03261         EMData* vol_ssnr = new EMData();
03262         vol_ssnr->set_size(m_vnxp+ 2 - m_vnxp%2, m_vnyp ,m_vnzp);
03263         vol_ssnr->to_zero();
03264         if ( m_vnxp % 2 == 0 ) vol_ssnr->set_fftodd(0);
03265         else                   vol_ssnr->set_fftodd(1);
03266         vol_ssnr->set_nxc(m_vnxc);
03267         vol_ssnr->set_complex(true);
03268         vol_ssnr->set_ri(true);
03269         vol_ssnr->set_fftpad(false);
03270         //#
03271         float *nom    = new float[inc+1];
03272         float *denom  = new float[inc+1];
03273         int  *ka     = new int[inc+1];
03274         int  *nn     = new int[inc+1];
03275         float wght = 1.f;
03276         for (int i = 0; i <= inc; i++) {
03277                 nom[i]   = 0.0f;
03278                 denom[i] = 0.0f;
03279                 nn[i]    = 0;
03280                 ka[i]    = 0;
03281         }
03282         m_volume->symplane2(m_wptr, m_wptr2, m_wptr3);
03283         if ( m_weighting == ESTIMATE ) {
03284                 int vol = box*box*box;
03285                 for( unsigned int i=1; i < pow_a.size(); ++i ) pow_a[i] = pow_a[i-1] * exp(m_wghta);
03286                 pow_a[3*kc] = 0.0;
03287                 float max = max3d( kc, pow_a );
03288                 alpha = ( 1.0f - 1.0f/(float)vol ) / max;
03289         }
03290         for (int iz = 1; iz <= m_vnzp; iz++) {
03291                 if ( iz-1 > m_vnzc ) kz = iz-1-m_vnzp; else kz = iz-1;
03292                 argz = float(kz*kz)*dz2;
03293                 for (int iy = 1; iy <= m_vnyp; iy++) {
03294                         if ( iy-1 > m_vnyc ) ky = iy-1-m_vnyp; else ky = iy-1;
03295                         argy = argz + float(ky*ky)*dy2;
03296                         for (int ix = 0; ix <= m_vnxc; ix++) {
03297                                 float Kn = (*m_wptr3)(ix,iy,iz);
03298                                 argx = std::sqrt(argy + float(ix*ix)*dx2);
03299                                 int r = Util::round(float(inc)*argx);
03300                                 if ( r >= 0 && Kn > 4.5f ) {
03301                                         if ( m_weighting == ESTIMATE ) {
03302                                                 int cx = ix;
03303                                                 int cy = (iy<=m_vnyc) ? iy - 1 : iy - 1 - m_vnyp;
03304                                                 int cz = (iz<=m_vnzc) ? iz - 1 : iz - 1 - m_vnzp;
03305                                                 float sum = 0.0;
03306                                                 for( int ii = -kc; ii <= kc; ++ii ) {
03307                                                         int nbrcx = cx + ii;
03308                                                         if( nbrcx >= m_vnxc ) continue;
03309                                                         for ( int jj= -kc; jj <= kc; ++jj ) {
03310                                                                 int nbrcy = cy + jj;
03311                                                                 if( nbrcy <= -m_vnyc || nbrcy >= m_vnyc ) continue;
03312                                                                 for( int kk = -kc; kk <= kc; ++kk ) {
03313                                                                         int nbrcz = cz + jj;
03314                                                                         if ( nbrcz <= -m_vnyc || nbrcz >= m_vnyc ) continue;
03315                                                                         if( nbrcx < 0 ) {
03316                                                                                 nbrcx = -nbrcx;
03317                                                                                 nbrcy = -nbrcy;
03318                                                                                 nbrcz = -nbrcz;
03319                                                                         }
03320                                                                         int nbrix = nbrcx;
03321                                                                         int nbriy = nbrcy >= 0 ? nbrcy + 1 : nbrcy + 1 + m_vnyp;
03322                                                                         int nbriz = nbrcz >= 0 ? nbrcz + 1 : nbrcz + 1 + m_vnzp;
03323                                                                         if( (*m_wptr)( nbrix, nbriy, nbriz ) == 0 ) {
03324                                                                                 int c = 3*kc+1 - std::abs(ii) - std::abs(jj) - std::abs(kk);
03325                                                                                 sum = sum + pow_a[c];
03326                                                                         }
03327                                                                 }
03328                                                         }
03329                                                 }
03330 //                                              int r = std::abs(cx) + std::abs(cy) + std::abs(cz);
03331                                                 wght = 1.0f / ( 1.0f - alpha * sum );
03332                                         } // end of ( m_weighting == ESTIMATE )
03333                                         float nominator   = std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz);
03334                                         float denominator = ((*m_wptr2)(ix,iy,iz)-std::norm(m_volume->cmplx(ix,iy,iz))/(*m_wptr)(ix,iy,iz))/(Kn-1.0f);
03335                                         // Skip Friedel related values
03336                                         if( (ix>0 || (kz>=0 && (ky>=0 || kz!=0)))) {
03337                                                 if ( r <= inc ) {
03338                                                         nom[r]   += nominator*wght;
03339                                                         denom[r] += denominator/(*m_wptr)(ix,iy,iz)*wght;
03340                                                         nn[r]    += 2;
03341                                                         ka[r]    += int(Kn);
03342                                                 }
03343 /*
03344 #ifdef  _WIN32
03345                                                 float  tmp = _cpp_max(nominator/denominator/(*m_wptr)(ix,iy,iz)-1.0f,0.0f);
03346 #else
03347                                                 float  tmp = std::max(nominator/denominator/(*m_wptr)(ix,iy,iz)-1.0f,0.0f);
03348 #endif  //_WIN32
03349                                                 //  Create SSNR as a 3D array (-n/2:n/2+n%2-1)
03350                                                 int iix = m_vnxc + ix; int iiy = m_vnyc + ky; int iiz = m_vnzc + kz;
03351                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
03352                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
03353                                                 // Friedel part
03354                                                 iix = m_vnxc - ix; iiy = m_vnyc - ky; iiz = m_vnzc - kz;
03355                                                 if( iix >= 0 && iix < m_vnxp && iiy >= 0 && iiy < m_vnyp && iiz >= 0 && iiz < m_vnzp )
03356                                                         (*vol_ssnr)(iix, iiy, iiz) = tmp;
03357 */
03358                                         }
03359                                         (*vol_ssnr)(2*ix, iy-1, iz-1) = denominator*wght;
03360                                 } // end of Kn>4.5 or whatever
03361                         }
03362                 }
03363         }
03364         for (int i = 0; i <= inc; i++) {
03365                 (*SSNR)(i,0,0) = nom[i];
03366                 (*SSNR)(i,1,0) = denom[i];
03367                 (*SSNR)(i,2,0) = static_cast<float>(nn[i]);
03368                 (*SSNR)(i,3,0) = static_cast<float>(ka[i]);
03369         }
03370         vol_ssnr->update();
03371         return vol_ssnr;
03372 //      }
03373 }
03374 #undef  tw
03375 // -----------------------------------------------------------------------------------
03376 // End of this addition
03377 
03378 
03379 void EMAN::dump_reconstructors()
03380 {
03381         dump_factory < Reconstructor > ();
03382 }
03383 
03384 map<string, vector<string> > EMAN::dump_reconstructors_list()
03385 {
03386         return dump_factory_list < Reconstructor > ();
03387 }
03388 
03389 
03390 
03391 
03392 
03393 using std::ofstream;
03394 using std::ifstream;
03395 
03396 
03397 newfile_store::newfile_store( const string& filename, int npad, bool ctf )
03398     : m_bin_file( filename + ".bin" ),
03399       m_txt_file( filename + ".txt" )
03400 {
03401     m_npad = npad;
03402     m_ctf = ctf;
03403 }
03404 
03405 newfile_store::~newfile_store( )
03406 {
03407 }
03408 
03409 void newfile_store::add_image( EMData* emdata, const Transform& tf )
03410 {
03411     if( m_bin_of == NULL )
03412     {
03413         m_bin_of = shared_ptr<ofstream>( new ofstream(m_bin_file.c_str(), std::ios::out|std::ios::binary) );
03414         m_txt_of = shared_ptr<ofstream>( new ofstream(m_txt_file.c_str()) );
03415     }
03416 
03417 
03418     EMData* padfft = padfft_slice( emdata, tf, m_npad );
03419 
03420     int nx = padfft->get_xsize();
03421     int ny = padfft->get_ysize();
03422     int n2 = ny / 2;
03423     int n = ny;
03424 
03425     float voltage=0.0f, pixel=0.0f, Cs=0.0f, ampcont=0.0f, bfactor=0.0f, defocus=0.0f;
03426 
03427     if( m_ctf )
03428     {
03429         Ctf* ctf = emdata->get_attr( "ctf" );
03430         Dict params = ctf->to_dict();
03431         voltage = params["voltage"];
03432         pixel   = params["apix"];
03433         Cs      = params["cs"];
03434         ampcont = params["ampcont"];
03435         bfactor = params["bfactor"];
03436         defocus = params["defocus"];
03437         if(ctf) {delete ctf; ctf=0;}
03438     }
03439 
03440     vector<point_t> points;
03441     for( int j=-ny/2+1; j <= ny/2; j++ )
03442     {
03443         int jp = (j>=0) ? j+1 : ny+j+1;
03444         for( int i=0; i <= n2; ++i )
03445         {
03446             int r2 = i*i + j*j;
03447             if( (r2<ny*ny/4) && !( (i==0) && (j<0) ) )
03448             {
03449                 float ctf;
03450                 if( m_ctf )
03451                 {
03452                     float ak = std::sqrt( r2/float(ny*ny) )/pixel;
03453                     ctf = Util::tf( defocus, ak, voltage, Cs, ampcont, bfactor, 1);
03454                 }
03455                 else
03456                 {
03457                     ctf = 1.0;
03458                 }
03459 
03460                 float xnew = i*tf[0][0] + j*tf[1][0];
03461                 float ynew = i*tf[0][1] + j*tf[1][1];
03462                 float znew = i*tf[0][2] + j*tf[1][2];
03463                 std::complex<float> btq;
03464                 if (xnew < 0.)
03465                 {
03466                     xnew = -xnew;
03467                     ynew = -ynew;
03468                     znew = -znew;
03469                     btq = conj(padfft->cmplx(i,jp-1));
03470                 }
03471                 else
03472                 {
03473                     btq = padfft->cmplx(i,jp-1);
03474                 }
03475 
03476                 int ixn = int(xnew + 0.5 + n) - n;
03477                 int iyn = int(ynew + 0.5 + n) - n;
03478                 int izn = int(znew + 0.5 + n) - n;
03479                 if ((ixn <= n2) && (iyn >= -n2) && (iyn <= n2) && (izn >= -n2) && (izn <= n2))
03480                 {
03481                     int ixf, iyf, izf;
03482                     if (ixn >= 0)
03483                     {
03484                         int iza, iya;
03485                         if (izn >= 0)
03486                             iza = izn + 1;
03487                         else
03488                             iza = n + izn + 1;
03489 
03490                         if (iyn >= 0)
03491                             iya = iyn + 1;
03492                         else
03493                             iya = n + iyn + 1;
03494 
03495                         ixf = ixn;
03496                         iyf = iya;
03497                         izf = iza;
03498                     }
03499                     else
03500                     {
03501                         int izt, iyt;
03502                         if (izn > 0)
03503                             izt = n - izn + 1;
03504                         else
03505                             izt = -izn + 1;
03506 
03507                         if (iyn > 0)
03508                             iyt = n - iyn + 1;
03509                         else
03510                             iyt = -iyn + 1;
03511 
03512                         ixf = -ixn;
03513                         iyf = iyt;
03514                         izf = izt;
03515                     }
03516 
03517 
03518                     int pos2 = ixf + (iyf-1)*nx/2 + (izf-1)*ny*nx/2;
03519                     float ctfv1 = btq.real() * ctf;
03520                     float ctfv2 = btq.imag() * ctf;
03521                     float ctf2 = ctf*ctf;
03522 
03523                     point_t p;
03524                     p.pos2 = pos2;
03525                     p.real = ctfv1;
03526                     p.imag = ctfv2;
03527                     p.ctf2 = ctf2;
03528 
03529                     points.push_back( p );
03530                 }
03531             }
03532         }
03533     }
03534 
03535 
03536     int npoint = points.size();
03537     std::istream::off_type offset = (m_offsets.size()==0) ? 0 : m_offsets.back();
03538     offset += npoint*sizeof(point_t);
03539     m_offsets.push_back( offset );
03540 
03541     *m_txt_of << m_offsets.back() << std::endl;
03542     m_bin_of->write( (char*)(&points[0]), sizeof(point_t)*npoint );
03543     checked_delete( padfft );
03544 }
03545 
03546 void newfile_store::get_image( int id, EMData* buf )
03547 {
03548     if( m_offsets.size()==0 )
03549     {
03550         ifstream is( m_txt_file.c_str() );
03551         std::istream::off_type off;
03552         while( is >> off )
03553         {
03554             m_offsets.push_back( off );
03555         }
03556 
03557         m_bin_if = shared_ptr<std::ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
03558     }
03559 
03560     Assert( m_bin_if != NULL );
03561 
03562     std::istream::off_type offset = (id==0) ? 0 : m_offsets[id-1];
03563     Assert( offset >= 0 );
03564     m_bin_if->seekg( offset, std::ios::beg );
03565 
03566 
03567     if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() )
03568     {
03569         std::cout << "bad or fail or eof while fetching id, offset: " << id << " " << offset << std::endl;
03570         throw std::logic_error( "bad happen" );
03571     }
03572 
03573     int bufsize = (m_offsets[id] - offset)/sizeof(float);
03574     if( buf->get_xsize() != bufsize )
03575     {
03576         buf->set_size( bufsize, 1, 1 );
03577     }
03578 
03579     char* data = (char*)(buf->get_data());
03580     m_bin_if->read( data, sizeof(float)*bufsize );
03581     buf->update();
03582 }
03583 
03584 void newfile_store::read( int nprj )
03585 {
03586     if( m_offsets.size()==0 )
03587     {
03588         ifstream is( m_txt_file.c_str() );
03589         std::istream::off_type off;
03590         while( is >> off )
03591         {
03592             m_offsets.push_back( off );
03593         }
03594     }
03595 
03596     if( m_bin_if==NULL )
03597     {
03598         m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
03599     }
03600 
03601 
03602     int npoint = m_offsets[0]/sizeof(point_t);
03603     std::ios::off_type prjsize = m_offsets[0];
03604 
03605     try
03606     {
03607         m_points.resize(nprj * npoint);
03608     }
03609     catch( std::exception& e )
03610     {
03611         std::cout << "Error: " << e.what() << std::endl;
03612     }
03613 
03614     int ip = 0;
03615     for( int i=0; i < nprj; ++i )
03616     {
03617         m_bin_if->read( (char*)(&m_points[ip]), prjsize );
03618             if( m_bin_if->bad() || m_bin_if->fail() || m_bin_if->eof() )
03619         {
03620             std::cout << "Error: file hander bad or fail or eof" << std::endl;
03621             return;
03622         }
03623         ip += npoint;
03624     }
03625 }
03626 
03627 void newfile_store::add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend )
03628 {
03629     float* vdata = fftvol->get_data();
03630     float* wdata = wgtvol->get_data();
03631 
03632     int npoint = m_offsets[0]/sizeof(point_t);
03633 //    Assert( int(mults.size())==nprj );
03634     Assert( int(m_points.size())== (pend - pbegin)*npoint );
03635 
03636     for( int iprj=pbegin; iprj < pend; ++iprj )
03637     {
03638         int m = mults[iprj];
03639         if( m==0 ) continue;
03640 
03641         int ipt = (iprj-pbegin)*npoint;
03642         for( int i=0; i < npoint; ++i )
03643         {
03644             int pos2 = m_points[ipt].pos2;
03645             int pos1 = pos2*2;
03646 
03647             wdata[pos2] += m_points[ipt].ctf2*m;
03648             vdata[pos1] += m_points[ipt].real*m;
03649             vdata[pos1+1]+= m_points[ipt].imag*m;
03650             ++ipt;
03651         }
03652     }
03653 }
03654 
03655 void newfile_store::restart()
03656 {
03657     m_bin_if = shared_ptr< ifstream>( new ifstream(m_bin_file.c_str(), std::ios::in|std::ios::binary) );
03658 }
03659 
03660 file_store::file_store(const string& filename, int npad, int write, bool ctf)
03661     : m_bin_file(filename + ".bin"),
03662       m_txt_file(filename + ".txt")
03663 {
03664     m_ctf = ctf;
03665     m_prev = -1;
03666     m_npad = npad;
03667     m_write = write;
03668 }
03669 
03670 file_store::~file_store()
03671 {
03672 }
03673 
03674 void file_store::add_image( EMData* emdata, const Transform& tf )
03675 {
03676 
03677     EMData* padfft = padfft_slice( emdata, tf, m_npad );
03678 
03679     float* data = padfft->get_data();
03680 
03681     if( m_write && m_bin_ohandle == NULL )
03682     {
03683         m_bin_ohandle = shared_ptr< ofstream >( new ofstream(m_bin_file.c_str(), std::ios::out | std::ios::binary) );
03684         m_txt_ohandle = shared_ptr< ofstream >( new ofstream(m_txt_file.c_str() ) );
03685         if( m_ctf )
03686                  *m_txt_ohandle << "Cs pixel voltage ctf_applied amp_contrast defocus ";
03687 
03688         *m_txt_ohandle << "phi theta psi" << std::endl;
03689     }
03690 
03691     m_x_out = padfft->get_xsize();
03692     m_y_out = padfft->get_ysize();
03693     m_z_out = padfft->get_zsize();
03694     m_totsize = m_x_out*m_y_out*m_z_out;
03695 
03696     if( m_ctf )
03697     {
03698         Ctf* ctf = padfft->get_attr( "ctf" );
03699         Dict ctf_params = ctf->to_dict();
03700 
03701         m_ctf_applied = padfft->get_attr( "ctf_applied" );
03702 
03703         m_Cs = ctf_params["cs"];
03704         m_pixel = ctf_params["apix"];
03705         m_voltage = ctf_params["voltage"];
03706         m_amp_contrast = ctf_params["ampcont"];
03707         m_defocuses.push_back( ctf_params["defocus"] );
03708         if(ctf) {delete ctf; ctf=0;}
03709     }
03710 
03711     Dict params = tf.get_rotation( "spider" );
03712     float phi = params.get( "phi" );
03713     float tht = params.get( "theta" );
03714     float psi = params.get( "psi" );
03715 
03716 
03717     m_phis.push_back( phi );
03718     m_thetas.push_back( tht );
03719     m_psis.push_back( psi );
03720 
03721     if( m_write )
03722     {
03723         m_bin_ohandle->write( (char*)data, sizeof(float)*m_totsize );
03724 
03725         if( m_ctf )
03726         {
03727             *m_txt_ohandle << m_Cs << " ";
03728             *m_txt_ohandle << m_pixel << " ";
03729             *m_txt_ohandle << m_voltage << " ";
03730             *m_txt_ohandle << m_ctf_applied << " ";
03731             *m_txt_ohandle << m_amp_contrast << " ";
03732             *m_txt_ohandle << m_defocuses.back() << " ";
03733         }
03734         *m_txt_ohandle << m_phis.back() << " ";
03735         *m_txt_ohandle << m_thetas.back() << " ";
03736         *m_txt_ohandle << m_psis.back() << " ";
03737         *m_txt_ohandle << m_x_out << " ";
03738         *m_txt_ohandle << m_y_out << " ";
03739         *m_txt_ohandle << m_z_out << " ";
03740         *m_txt_ohandle << m_totsize << std::endl;
03741     }
03742 
03743     checked_delete(padfft);
03744 
03745 }
03746 
03747 void file_store::get_image( int id, EMData* padfft )
03748 {
03749 
03750     if( m_phis.size() == 0 ) {
03751         ifstream m_txt_ifs( m_txt_file.c_str() );
03752 
03753         if( !m_txt_ifs )
03754         {
03755             std::cerr << "Error: file " << m_txt_file << " does not exist" << std::endl;
03756         }
03757 
03758         string line;
03759         std::getline( m_txt_ifs, line );
03760 
03761         float first, defocus, phi, theta, psi;
03762 
03763 
03764 
03765         while( m_txt_ifs >> first ) {
03766 
03767             if( m_ctf )
03768             {
03769                 m_Cs = first;
03770                 m_txt_ifs >> m_pixel >> m_voltage;
03771                 m_txt_ifs >> m_ctf_applied >> m_amp_contrast;
03772                 m_txt_ifs >> defocus >> phi >> theta >> psi;
03773                 m_defocuses.push_back( defocus );
03774             }
03775             else
03776             {
03777                 phi = first;
03778                 m_txt_ifs >> theta >> psi;
03779             }
03780 
03781             m_txt_ifs >> m_x_out >> m_y_out >> m_z_out >> m_totsize;
03782             m_phis.push_back( phi );
03783             m_thetas.push_back( theta );
03784             m_psis.push_back( psi );
03785         }
03786     }
03787 
03788     Assert( m_ihandle != NULL );
03789 
03790     std::istream::off_type offset = id*sizeof(float)*m_totsize;
03791     Assert( offset >= 0 );
03792 
03793     if( offset > 0 )
03794     {
03795         m_ihandle->seekg(offset, std::ios::beg);
03796     }
03797 
03798     if( m_ihandle->bad() )
03799     {
03800         std::cout << "bad while fetching id, offset: " << id << " " << offset << std::endl;
03801         throw std::logic_error( "bad happen" );
03802     }
03803 
03804     if( m_ihandle->fail() )
03805     {
03806         std::cout << "fail while fetching id, offset, curoff: " << id << " " << offset << std::endl;
03807         throw std::logic_error( "fail happen" );
03808     }
03809 
03810     if( m_ihandle->eof() )
03811     {
03812         std::cout << "eof while fetching id, offset: " << id << " " << offset << std::endl;
03813         throw std::logic_error( "eof happen" );
03814     }
03815 
03816     if( padfft->get_xsize() != m_x_out ||
03817         padfft->get_ysize() != m_y_out ||
03818         padfft->get_zsize() != m_z_out )
03819     {
03820         padfft->set_size(m_x_out, m_y_out, m_z_out);
03821     }
03822 
03823     char* data = (char*)(padfft->get_data());
03824     m_ihandle->read( data, sizeof(float)*m_totsize );
03825     padfft->update();
03826 
03827     if( m_ctf )
03828     {
03829         padfft->set_attr( "Cs", m_Cs );
03830         padfft->set_attr( "Pixel_size", m_pixel );
03831         padfft->set_attr( "voltage", m_voltage );
03832         padfft->set_attr( "ctf_applied", m_ctf_applied );
03833         padfft->set_attr( "amp_contrast", m_amp_contrast );
03834         padfft->set_attr( "defocus", m_defocuses[id] );
03835     }
03836 
03837     padfft->set_attr( "padffted", 1 );
03838     padfft->set_attr( "phi", m_phis[id] );
03839     padfft->set_attr( "theta", m_thetas[id] );
03840     padfft->set_attr( "psi", m_psis[id] );
03841 
03842 }
03843 
03844 void file_store::restart( )
03845 {
03846     if( m_ihandle == NULL )
03847     {
03848         m_ihandle = shared_ptr< ifstream >( new ifstream(m_bin_file.c_str(), std::ios::in | std::ios::binary) );
03849     }
03850 
03851     if( m_ihandle->bad() || m_ihandle->fail() || m_ihandle->eof() )
03852     {
03853         m_ihandle->open( m_bin_file.c_str(), std::ios::binary );
03854     }
03855 
03856     m_ihandle->seekg( 0, std::ios::beg );
03857 }
03858 
03859 /* vim: set ts=4 noet: */

Generated on Tue May 25 17:13:36 2010 for EMAN2 by  doxygen 1.4.7