reconstructor.h

Go to the documentation of this file.

/*
 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
 * Copyright (c) 2000-2006 Baylor College of Medicine
 *
 * This software is issued under a joint BSD/GNU license. You may use the
 * source code in this file under either license. However, note that the
 * complete EMAN2 and SPARX software packages have some GPL dependencies,
 * so you are responsible for compliance with the licenses of these packages
 * if you opt to use BSD licensing. The warranty disclaimer below holds
 * in either instance.
 *
 * This complete copyright notice must be included in any revised version of the
 * source code. Additional authorship citations may be added, but existing
 * author citations must be preserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 * */
 
//#define RECONDEBUG            // This is used to add some data objects for debugging reconstructors, it should normally be commented out!
 
#ifndef eman_reconstructor_h__
#define eman_reconstructor_h__ 1
#include <fstream>
#include <memory>
#include "emdata.h"
#include "exception.h"
#include "emobject.h"
#include "interp.h"
 
using std::vector;
using std::map;
using std::string;
using std::shared_ptr;
 
using std::cout;
using std::cerr;
using std::endl;
 
#include <utility>
using std::pair;
 
#include "reconstructor_tools.h"
 
namespace EMAN
{
 
        class Transform;
        class EMData;
 
        class Reconstructor : public FactoryBase
        {
          public:
                Reconstructor() {}
                virtual ~Reconstructor() {}
                virtual void setup() = 0;
                
                virtual void setup_seed(EMData* seed,float seed_weight) {throw;}
 
                virtual void setup_seedandweights(EMData* seed,EMData* weight) {throw;}
 
                
                virtual EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() ) { EMData *ret=slice->copy(); ret->set_attr("reconstruct_preproc",(int)1); return ret; }
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight) {throw;}
                int insert_slice(const EMData* const slice, const Transform & euler) { return this->insert_slice(slice, euler, 1.0f); }
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true ) { throw; }
 
                virtual EMData* projection(const Transform &euler, int ret_fourier = 1) { throw; }
 
                
                virtual EMData *finish(bool doift=true) { throw; }
                
                virtual void clear() {throw; }
 
                void print_params() const
                {
                        std::cout << "Printing reconstructor params" << std::endl;
                        for ( Dict::const_iterator it = params.begin(); it != params.end(); ++it )
                        {
                                std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
                        }
                        std::cout << "Done printing reconstructor params" << std::endl;
                }
 
 
                EMObject& operator[]( const string& key ) { return params[key]; }
 
          private:
                // Disallow copy construction
                Reconstructor(const Reconstructor& that);
                Reconstructor& operator=(const Reconstructor& );
 
        };
 
        class ReconstructorVolumeData
        {
                public:
                        inline ReconstructorVolumeData() : image(0), tmp_data(0), nx(0), ny(0), nz(0), subnx(0), subny(0), subnz(0), subx0(0), suby0(0), subz0(0) {}
                        
                        virtual ~ReconstructorVolumeData() { free_memory(); }
 
                        const EMData* get_emdata() { return image; }
                protected:
                        //These EMData pointers will most probably be allocated in setup() and released in finish()
                        EMData* image;
                        EMData* tmp_data;
 
                        // nx,ny,nz generally will store the dimensions of image
                        int nx,nx2;
                        int ny,ny2;
                        int nz,nz2;
                        
                        // Note that subvolume capabilities may not be 100% functional
                        // subvolume size
                        int subnx;
                        int subny;
                        int subnz;
                        
                        // subvolume origin
                        int subx0;
                        int suby0;
                        int subz0;
                        
#ifdef RECONDEBUG
                        double *ddata,*dnorm;
#endif
 
                protected:
                        void free_memory()
                        {
                                if (image != 0)  {delete image; image = 0;}
                                if ( tmp_data != 0 ) { delete tmp_data; tmp_data = 0; }
                        }
 
                        virtual void normalize_threed(const bool sqrt_damp=false,const bool wiener=false);
 
                        virtual void zero_memory()
                        {
                                if (tmp_data != 0 ) tmp_data->to_zero();
                                if (image != 0 ) image->to_zero();
                        }
 
                private:
                ReconstructorVolumeData(const ReconstructorVolumeData& that);
                ReconstructorVolumeData& operator=(const ReconstructorVolumeData& );
 
        };
 
        class FourierReconstructorSimple2D : public Reconstructor, public ReconstructorVolumeData
        {
                public:
                        FourierReconstructorSimple2D() {}
 
                        virtual ~FourierReconstructorSimple2D() { }
 
                        virtual void setup();
 
                        virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                        virtual EMData *finish(bool doift=true);
 
                        virtual string get_name() const { return NAME; }
 
                        virtual string get_desc() const { return "performs 2D reconstruction"; }
 
                        static Reconstructor *NEW()
                        {
                                return new FourierReconstructorSimple2D();
                        }
 
 
                        virtual TypeDict get_param_types() const
                        {
                                TypeDict d;
                                d.put("nx", EMObject::INT, "Necessary. The x dimension of the input images.");
//                              d.put("sym", EMObject::STRING, "Symmetry - assumed to be C1 if not specified");
                                return d;
                        }
                        
                        static const string NAME;
        };
 
 
 
        class FourierReconstructor : public Reconstructor, public ReconstructorVolumeData
        {
          public:
                FourierReconstructor() { load_default_settings(); }
 
                virtual ~FourierReconstructor() { free_memory(); }
 
                virtual void setup();
 
                virtual void setup_seed(EMData* seed,float seed_weight);
 
                virtual void setup_seedandweights(EMData* seed,EMData* weight);
 
                virtual EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() );
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual EMData* projection(const Transform &euler, int ret_fourier);
 
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
 
                virtual EMData *finish(bool doift=true);
                
                virtual void clear();
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based";
                }
 
                static Reconstructor *NEW()
                {
                        return new FourierReconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("sym", EMObject::STRING, "Required. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1, ie - an asymmetric object");
                        d.put("mode", EMObject::STRING, "Optional. Fourier pixel insertion mode name (nearest_neighbor, gauss_2, gauss_3, gauss_5, gauss_5_slow, gypergeom_5, experimental) gauss_2 is the default.");
                        d.put("corners", EMObject::BOOL, "Optional. If not set, then reconstruction will cover a spherical volume with a radius of nx/2. If set, the full Fourier volume will be reconstructed, but will be ~2x slower.");
                        d.put("sqrtnorm", EMObject::BOOL, "Optional. When normalizing, additionally divides by the sqrt of the normalization factor to damp exaggerated features. Is this justifyable ? No idea (yet). Default is false.");
                        d.put("usessnr", EMObject::BOOL, "Optional. Looks for and uses the class_ssnr header parameter from each slice to weight each voxel during insertion to the reconstruction.");
                        d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        d.put("quiet", EMObject::BOOL, "Optional. If false, print verbose information.");
                        d.put("subvolume",EMObject::INTARRAY, "Optional. (xorigin,yorigin,zorigin,xsize,ysize,zsize) all in Fourier pixels. Useful for parallelism.");
                        d.put("savenorm",EMObject::STRING, "Debug. Will cause the normalization volume to be written directly to the specified file when finish() is called.");
                        
                        d.put("normout",EMObject::EMDATA, "Will write the normalization volume to the given EMData object file when finish() is called.");
                        return d;
                }
                
                static const string NAME;
 
          protected:
                virtual void load_default_settings();
 
                virtual void free_memory();
 
                virtual void load_inserter();
 
                virtual void do_insert_slice_work(const EMData* const input_slice, const Transform & euler,const float weight, const bool corners=false);
 
                virtual void do_compare_slice_work(EMData* input_slice, const Transform & euler,float weight);
 
                virtual bool pixel_at(const float& xx, const float& yy, const float& zz, float *dt);
 
                FourierPixelInserter3D* inserter;
 
          private:
                FourierReconstructor( const FourierReconstructor& that );
                FourierReconstructor& operator=( const FourierReconstructor& );
 
        };
 
 
 
        class WienerFourierReconstructor : public FourierReconstructor
        {
          public:
                WienerFourierReconstructor() {};
 
                virtual ~WienerFourierReconstructor() { }
 
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based. This version also incorporates a nonisotropic Wiener filter based on SNR estimates stored in the class-average headers by the ctf.auto averager.";
                }
 
                static Reconstructor *NEW()
                {
                        return new WienerFourierReconstructor();
                }
                
                static const string NAME;
 
          protected:
                
                virtual void do_insert_slice_work(const EMData* const input_slice, const Transform & euler,const float weight);
 
                virtual void do_compare_slice_work(EMData* input_slice, const Transform & euler,float weight);
 
                virtual bool pixel_at(const float& xx, const float& yy, const float& zz, float *dt);
 
//              FourierPixelInserter3D* inserter;
 
          private:
                WienerFourierReconstructor( const WienerFourierReconstructor& that );
                WienerFourierReconstructor& operator=( const WienerFourierReconstructor& );
 
        };
 
        class FourierIterReconstructor : public Reconstructor, public ReconstructorVolumeData
        {
          public:
                FourierIterReconstructor(): ref_vol(0) { }
 
                virtual ~FourierIterReconstructor() { free_memory(); }
 
                virtual void setup();
 
                virtual void setup_seed(EMData* seed,float seed_weight);
                virtual void setup_seedandweights(EMData* seed,EMData* weight);
 
                virtual EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() );
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
 
                virtual EMData *finish(bool doift=true);
                
                virtual void clear();
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Reconstruction via Fourier interpolation using an iterative strategy to define the interpolation kernel";
                }
 
                static Reconstructor *NEW()
                {
                        return new FourierIterReconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1, ie - an asymmetric object");
                        d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        return d;
                }
                
                static const string NAME;
 
          protected:
 
                virtual void free_memory();
                
                EMData *ref_vol;
 
          private:
                FourierIterReconstructor( const FourierIterReconstructor& that );
                FourierIterReconstructor& operator=( const FourierIterReconstructor& );
 
        };
 
        
        class BackProjectionReconstructor:public Reconstructor, public ReconstructorVolumeData
        {
          public:
                BackProjectionReconstructor() { load_default_settings();  }
 
                virtual ~BackProjectionReconstructor() {}
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Simple (unfiltered) back-projection reconstruction. Weighting by contributing particles in the class average is optional and default behaviour";
                }
 
                static Reconstructor *NEW()
                {
                        return new BackProjectionReconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("weight", EMObject::FLOAT, "Optional. A temporary value set prior to slice insertion, indicative of the inserted slice's weight. Default sis 1.");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry to impose on the final reconstruction. Default is c1");
                        d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        return d;
                }
                
                static const string NAME;
                
          private:
                // Disallow copy construction
                BackProjectionReconstructor( const BackProjectionReconstructor& that);
                // Disallow assignment
                BackProjectionReconstructor& operator=( const BackProjectionReconstructor& );
 
                void load_default_settings()
                {
                        params["weight"] = 1.0;
                        params["use_weights"] = true;
                        params["sym"] = "c1";
                        params["verbose"] = false;
                }
 
                EMData* preprocess_slice(const EMData* const slice, const Transform& t);
        };
 
        class RealMedianReconstructor:public Reconstructor, public ReconstructorVolumeData
        {
          public:
                RealMedianReconstructor() {  }
 
                virtual ~RealMedianReconstructor() {}
 
                virtual void setup();
 
                virtual EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() );
                
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "A back projection reconstructor with alternative statistics. mode:\n\
0 - median rather than mean per voxel (default)\n\
1 - similar to a mode. progressive subdivision of numeric axis to find dense value cluster\n\
2 - median with a bias towands zero (postive and negative)\n\
3 - median with a bias towards the 1/2 of the values spanning a narrower range\n\
4 - specifically for tomography, requires slices to be in order by tilt\n\
";
                }
 
                static Reconstructor *NEW()
                {
                        return new RealMedianReconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("mode", EMObject::INT, "Default is 0, median. 1 - compute estimated statistical mode.");
                        d.put("size", EMObject::INTARRAY, "Required. The dimensions of the real-space output volume, including any padding (must be handled by the calling application). Assumed that apix x/y/z identical.");
                        d.put("weight", EMObject::FLOAT, "Optional. A temporary value set prior to slice insertion, indicative of the inserted slice's weight. Default sis 1.");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry to impose on the final reconstruction. Default is c1");
                        d.put("verbose", EMObject::INT, "Optional. Toggles writing diagnostic information to standard out. Default is 0.");
                        return d;
                }
                
                static const string NAME;
                
          private:
                // Disallow copy construction
                RealMedianReconstructor( const RealMedianReconstructor& that);
                // Disallow assignment
                RealMedianReconstructor& operator=( const RealMedianReconstructor& );
                
                vector<EMData *> slices;
                vector<Transform> xforms;
        };
 
        EMData* padfft_slice( const EMData* const slice, const Transform& t, int npad );
 
        class nn4Reconstructor:public Reconstructor
        {
          public:
                nn4Reconstructor();
 
                nn4Reconstructor( const string& symmetry, int size, int npad );
 
                virtual ~nn4Reconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4Reconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("ndim",           EMObject::INT);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad );
 
                int insert_padfft_slice( EMData* padded, const Transform& trans, float mult=1 );
 
                static const string NAME;
 
          private:
                EMData* m_volume;
                EMData* m_wptr;
                string  m_symmetry;
                int m_weighting;
                int m_vnx, m_vny, m_vnz;
                int m_npad;
                int m_nsym;
                int m_ndim;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnzc, m_vnyc, m_vnxc;
                void buildFFTVolume();
                void buildNormVolume();
                float m_wghta;
                float m_wghtb;
                float m_osnr;
                void load_default_settings()
                {
                        //params["use_weights"] = false;
                }
        };
 
 
        class nn4_rectReconstructor:public Reconstructor
        {
          public:
                nn4_rectReconstructor();
 
                nn4_rectReconstructor( const string& symmetry, int size, int npad );
 
                virtual ~nn4_rectReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4_rectReconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("sizeprojection", EMObject::INT);
                        d.put("sizex",          EMObject::INT);
                        d.put("sizey",          EMObject::INT);
                        d.put("sizez",          EMObject::INT);
                        d.put("xratio",         EMObject::FLOAT);
                        d.put("yratio",         EMObject::FLOAT);
                        d.put("zratio",         EMObject::FLOAT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("ndim",           EMObject::INT);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad );
 
                int insert_padfft_slice( EMData* padded, const Transform& trans, float mult=1 );
 
                static const string NAME;
 
          private:
                EMData* m_volume;
                EMData* m_wptr;
                string  m_symmetry;
                int m_weighting;
                int m_vnx, m_vny, m_vnz;
                int m_npad;
                int m_nsym;
                int m_ndim;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnzc, m_vnyc, m_vnxc;
                int m_count;
                float m_xratio,m_yratio,m_zratio;//ratio of x,y,z direction in the 3d volume comparing to the cubic case
                float m_xscale,m_yscale;//ratio of x,y direction of 2D FFT after scaling and roatating operations
                int m_sizeofprojection;
                void buildFFTVolume();
                void buildNormVolume();
                float m_wghta;
                float m_wghtb;
                float m_osnr;
                void load_default_settings()
                {
                        //params["use_weights"] = false;
                }
        };
 
 
     /* Fourier Reconstruction by nearest neighbor with 3D SSNR
        Added by Zhengfan Yang on 03/16/07
     */
 
        class nnSSNR_Reconstructor:public Reconstructor
        {
 
          public:
                nnSSNR_Reconstructor();
 
                nnSSNR_Reconstructor( const string& symmetry, int size, int npad);
 
                ~nnSSNR_Reconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Reconstruction by nearest neighbor with 3D SSNR";
                }
 
                static Reconstructor *NEW()
                {
                        return new nnSSNR_Reconstructor();
                }
 
                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size", EMObject::INT);
                        d.put("npad", EMObject::INT);
                        d.put("symmetry", EMObject::STRING);
                        d.put("fftvol", EMObject::EMDATA);
                        d.put("weight", EMObject::EMDATA);
                        d.put("weight2", EMObject::EMDATA);
                        d.put("SSNR", EMObject::EMDATA);
                        d.put("vol_ssnr", EMObject::EMDATA);
                        d.put("w", EMObject::FLOAT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad);
 
                int insert_padfft_slice( EMData* padded, const Transform& trans, float mult=1 );
 
                static const string NAME;
                
          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_wptr2;
                string  m_symmetry;
                int m_weighting;
                int m_vnx, m_vny, m_vnz;
                int m_npad;
                int m_nsym;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnzc, m_vnyc, m_vnxc;
                void buildFFTVolume();
                void buildNormVolume();
                void buildNorm2Volume();
                float m_wghta;
                float m_wghtb;
        };
 
 
        class nn4_ctfReconstructor:public Reconstructor
        {
          public:
                nn4_ctfReconstructor();
 
                nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign );
 
                virtual ~nn4_ctfReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion reconstruction routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4_ctfReconstructor();
                }
 
 
                TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        d.put("varsnr",     EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad, float snr, int sign );
 
                int insert_padfft_slice( EMData* padfft, EMData* ctf2d2, const Transform& trans, float mult=1);
 
                int insert_buffed_slice( const EMData* buffer, float mult );
 
                static const string NAME;
                
          private:
                EMData* m_volume;
                EMData* m_wptr;
                int m_vnx, m_vny, m_vnz;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnxc, m_vnyc, m_vnzc;
                int m_npad;
                int m_sign;
                int m_varsnr;
                int m_weighting;
                float m_wghta, m_wghtb;
                float m_snr;
                string m_symmetry;
                int m_nsym;
 
                void buildFFTVolume();
                void buildNormVolume();
 
        };
 
 
 
        class nn4_ctfwReconstructor:public Reconstructor
        {
          public:
                nn4_ctfwReconstructor();
 
                nn4_ctfwReconstructor( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf );
 
                virtual ~nn4_ctfwReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler, const float weight);
 
                virtual EMData *finish(bool compensate=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion reconstruction routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4_ctfwReconstructor();
                }
 
 
                TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("refvol",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        d.put("varsnr",     EMObject::INT);
                        d.put("do_ctf",     EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf );
 
                int insert_padfft_slice_weighted( EMData* padfft, EMData* ctf2d2, vector<float> bckgnoise, const Transform& trans, const float weight);
 
                //  It appears here one could add extra functions  int insert_buffed_slice( const EMData* buffer, float mult );
                //int insert_buffed_slice( const EMData* buffer, float mult );
 
                static const string NAME;
 
          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_refvol;
                int    m_vnx, m_vny, m_vnz;
                int    m_vnzp, m_vnyp, m_vnxp;
                int    m_vnxc, m_vnyc, m_vnzc;
                int    m_npad;
                int    m_sign;
                int    m_varsnr;
                int    m_weighting;
                float  m_wghta, m_wghtb;
                float  m_snr;
                string m_symmetry;
                int    m_nsym;
                int    m_do_ctf;
 
                void buildFFTVolume();
                void buildNormVolume();
 
        };
        
        
        class nn4_ctfwsReconstructor:public Reconstructor
        {
          public:
                nn4_ctfwsReconstructor();
 
                nn4_ctfwsReconstructor( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf );
 
                virtual ~nn4_ctfwsReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler, const float weight);
 
                virtual EMData *finish(bool compensate=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion reconstruction routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4_ctfwsReconstructor();
                }
 
 
                TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("refvol",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        d.put("varsnr",     EMObject::INT);
                        d.put("do_ctf",     EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int size, int npad, float snr, int sign, int do_ctf );
 
                int insert_padfft_slice_weighted( EMData* padfft, EMData* ctf2d2, vector<float> bckgnoise, const Transform& trans, const float weight);
 
                //  It appears here one could add extra functions  int insert_buffed_slice( const EMData* buffer, float mult );
                //int insert_buffed_slice( const EMData* buffer, float mult );
 
                static const string NAME;
 
          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_refvol;
                int    m_vnx, m_vny, m_vnz;
                int    m_vnzp, m_vnyp, m_vnxp;
                int    m_vnxc, m_vnyc, m_vnzc;
                int    m_npad;
                int    m_sign;
                int    m_varsnr;
                int    m_weighting;
                float  m_wghta, m_wghtb;
                float  m_snr;
                string m_symmetry;
                int    m_nsym;
                int    m_do_ctf;
 
                void buildFFTVolume();
                void buildNormVolume();
 
        };
        
        class nn4_ctf_rectReconstructor:public Reconstructor
        {
          public:
                nn4_ctf_rectReconstructor();
 
                nn4_ctf_rectReconstructor( const string& symmetry, int size, int npad, float snr, int sign );
 
                virtual ~nn4_ctf_rectReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler, const float weight);
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Direct Fourier inversion reconstruction routine";
                }
 
                static Reconstructor *NEW()
                {
                        return new nn4_ctf_rectReconstructor();
                }
 
 
                TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",           EMObject::INT);
                        d.put("npad",           EMObject::INT);
                        d.put("sign",           EMObject::INT);
                        d.put("symmetry",       EMObject::STRING);
                        d.put("snr",            EMObject::FLOAT);
                        d.put("fftvol",         EMObject::EMDATA);
                        d.put("weight",         EMObject::EMDATA);
                        d.put("refvol",         EMObject::EMDATA);
                        d.put("weighting",  EMObject::INT);
                        d.put("varsnr",     EMObject::INT);
                        return d;
                }
 
                void setup( const string& symmetry, int sizeprojection, int npad, float snr, int sign );
 
                int insert_padfft_slice( EMData* padfft, const Transform& trans, float mult=1);
 
                int insert_buffed_slice( const EMData* buffer, float mult );
                
                static const string NAME;
                
          private:
                EMData* m_volume;
                EMData* m_wptr;
                int m_vnx, m_vny, m_vnz;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnxc, m_vnyc, m_vnzc;
                int m_count;
                float m_xratio, m_yratio, m_zratio;//ratio of x,y,z direction in the 3d volume comparing to the cubic case
                int m_sizeofprojection;
                int m_npad;
                int m_sign;
            int m_varsnr;
                int m_weighting;
                float m_wghta, m_wghtb;
                float m_snr;
                string m_symmetry;
                int m_nsym;
 
                void buildFFTVolume();
                void buildNormVolume();
        };
 
 
 
     /* Fourier Reconstruction by nearest neighbor with 3D SSNR and CTF
        Added by Zhengfan Yang on 04/11/07
     */
 
        class nnSSNR_ctfReconstructor:public Reconstructor
        {
 
          public:
                nnSSNR_ctfReconstructor();
 
                nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign);
 
                ~nnSSNR_ctfReconstructor();
 
                virtual void setup();
 
                virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight);
 
 
                virtual EMData *finish(bool doift=true);
 
                virtual string get_name() const
                {
                        return NAME;
                }
 
                virtual string get_desc() const
                {
                        return "Reconstruction by nearest neighbor with 3D SSNR with CTF";
                }
 
                static Reconstructor *NEW()
                {
                        return new nnSSNR_ctfReconstructor();
                }
 
                TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("size",     EMObject::INT);
                        d.put("npad",     EMObject::INT);
                        d.put("symmetry", EMObject::STRING);
                        d.put("fftvol",   EMObject::EMDATA);
                        d.put("fftwvol",  EMObject::EMDATA);
                        d.put("weight",   EMObject::EMDATA);
                        d.put("weight2",  EMObject::EMDATA);
                        d.put("weight3",  EMObject::EMDATA);
                        d.put("SSNR",     EMObject::EMDATA);
                        d.put("vol_ssnr", EMObject::EMDATA);
                        d.put("w",        EMObject::FLOAT);
                        d.put("sign",     EMObject::INT);
                        d.put("snr",      EMObject::FLOAT);
                        return d;
                }
                void setup( const string& symmetry, int size, int npad, float snr, int sign);
 
                int insert_padfft_slice( EMData* padded, const Transform& trans, float mult=1 );
 
                static const string NAME;     
                
          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_wptr2;
                EMData* m_wptr3;
                string  m_symmetry;
                int m_weighting;
                int m_vnx, m_vny, m_vnz;
                int m_npad;
                int m_nsym;
                int m_vnzp, m_vnyp, m_vnxp;
                int m_vnzc, m_vnyc, m_vnxc;
                void buildFFTVolume();
                void buildNormVolume();
                void buildNorm2Volume();
                void buildNorm3Volume();
                float m_wghta;
                float m_wghtb;
                int   m_sign;
                float m_snr;
                int wiener;
        };
 
        template <> Factory < Reconstructor >::Factory();
 
        void dump_reconstructors();
        map<string, vector<string> > dump_reconstructors_list();
 
 
    struct point_t
    {
        int pos2;
        float real;
        float imag;
        float ctf2;
    };
 
 
        class newfile_store
        {
        public:
                newfile_store( const string& prefix, int npad, bool ctf );
 
                virtual ~newfile_store();
 
                void add_image( EMData* data, const Transform& tf );
 
                void add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend );
 
                void get_image( int id, EMData* buf );
 
        void read( int nprj );
 
                void restart( );
 
        private:
                int m_npad;
 
                bool m_ctf;
 
                string m_bin_file;
                string m_txt_file;
 
                shared_ptr<std::ofstream> m_bin_of;
                shared_ptr<std::ofstream> m_txt_of;
                shared_ptr<std::ifstream> m_bin_if;
                vector< std::istream::off_type > m_offsets;
 
        vector< point_t > m_points;
        };
 
        class file_store
        {
          public:
                file_store(const string& filename, int npad, int write, bool CTF);
 
                virtual ~file_store();
 
                void add_image(EMData* data, const Transform& tf);
 
                void get_image(int id, EMData* padfft);
 
                void restart();
          private:
                shared_ptr<std::ifstream> m_ihandle;
                shared_ptr<std::ofstream> m_bin_ohandle;
                shared_ptr<std::ofstream> m_txt_ohandle;
                string m_bin_file;
                string m_txt_file;
                int m_ctf;
                int m_npad;
                int m_prev;
                int m_x_out;
                int m_y_out;
                int m_z_out;
                int m_write;
                std::istream::off_type m_totsize;
                float m_Cs;
                float m_pixel;
                float m_voltage;
                float m_ctf_applied;
                float m_amp_contrast;
                vector< float > m_defocuses;
                vector< float > m_phis;
                vector< float > m_thetas;
                vector< float > m_psis;
        };
 
}
 
#endif