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
 *
 * */

#ifndef eman_reconstructor_h__
#define eman_reconstructor_h__ 1
#include <fstream>
#include <boost/shared_ptr.hpp>
#include "emdata.h"
#include "exception.h"
#include "emobject.h"
#include "interp.h"

using std::vector;
using std::map;
using std::string;
using boost::shared_ptr;

using std::cout;
using std::cerr;
using std::endl;

#include <utility>
using std::pair;

#include "reconstructor_tools.h"

namespace EMAN
{

        class Transform3D;
        class EMData;

        class Reconstructor : public FactoryBase
        {
          public:
                Reconstructor() {}
                virtual ~Reconstructor() {}
                virtual void setup() = 0;

                virtual int insert_slice(const EMData* const slice, const Transform3D & euler) {throw;};

                // This is pure virtual yet, but it should be... changed by Wei
                virtual int insert_slice(const EMData* const slice, const Transform & euler) {throw;};

                virtual int determine_slice_agreement(const EMData* const, const Transform3D &, const unsigned int)
                {
                        cout << "You called determine slice agreement but nothing happened - there is no functionality for determing slice agreement using this " << get_name() << " reconstructor" << endl;
                        return 0;
                }

                virtual int determine_slice_agreement(const EMData* const, const Transform &, const unsigned int)
                {
                        throw;
                }

                virtual EMData *finish() = 0;

                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]; }

                // A function used only by the Fourier reconstructor for testing and writing to std out purposes in e2make3d.py
                virtual float get_score(const unsigned int) { return 0; }
                // A function used only by the Fourier reconstructor for testing and writing to std out purposes in e2make3d.py
                virtual float get_norm(const unsigned int) { return 0; }


                virtual int insert_slice_weights(const Transform&) {return 0;}

          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) {}
                        
                        virtual ~ReconstructorVolumeData() { free_memory(); }

                        const EMData* const 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;
                        int ny;
                        int nz;

                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();

                        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 Transform3D & euler) {
                                throw UnexpectedBehaviorException("Transform3D interface is redundant");
                        }

                        virtual int insert_slice(const EMData* const slice, const Transform & euler);

                        virtual EMData *finish();

                        virtual string get_name() const { return "fouriersimple2D"; }

                        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;
                        }
        };



        class FourierReconstructor : public Reconstructor, public ReconstructorVolumeData
        {
          public:
                FourierReconstructor() : image_idx(0), inserter(0), interp_FRC_calculator(0), slice_insertion_flag(true),
                slice_agreement_flag(false), x_scale_factor(0.0), y_scale_factor(0.0), z_scale_factor(0.0),
                max_input_dim(0), output_x(0), output_y(0), output_z(0) { load_default_settings(); }

                virtual ~FourierReconstructor() { free_memory(); }

                virtual void setup();

                virtual int insert_slice(const EMData* const slice, const Transform & t3d = Transform());


                virtual int determine_slice_agreement(const EMData* const input_slice, const Transform & t3d, const unsigned int  num_particles_in_slice = 1);


                virtual int insert_slice(const EMData* const slice, const Transform3D & t3d) {
                        throw UnexpectedBehaviorException("No support for Transform3D in insert_slice, use Transform instead");
                };

                virtual int determine_slice_agreement(const EMData* const input_slice, const Transform3D & t3d, const unsigned int  num_particles_in_slice = 1) {
                        throw UnexpectedBehaviorException("No support for Transform3D in determine_slice_agreement, use Transform instead");
                };

                virtual EMData *finish();

                virtual string get_name() const
                {
                        return "fourier";
                }

                virtual string get_desc() const
                {
                        return "Reconstruction via direct Fourier methods using a combination of nearest neighbour and Gaussian kernels";
                }

                static Reconstructor *NEW()
                {
                        return new FourierReconstructor();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("x_in", EMObject::INT, "Necessary. The x dimension of the input images.");
                        d.put("y_in", EMObject::INT, "Necessary. The y dimension of the input images.");
                        d.put("zsample", EMObject::INT, "Optional. The z dimension (Fourier sampling) of the reconstructed volume, very useful for tomographic reconstruction. Works for general volumes.");
                        d.put("ysample", EMObject::INT, "Optional. The y dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
                        d.put("xsample", EMObject::INT, "Optional. The x dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
                        d.put("mode", EMObject::INT, "Optional. Fourier pixel insertion mode [1-8] - mode 2 is default.");
                        d.put("hard", EMObject::FLOAT, "Optional. The quality metric threshold. Default is 0 (off).");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1");
                        d.put("quiet", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        d.put("3damp", EMObject::BOOL, "Optional. Toggles writing the 3D FFT amplitude image. Default is false.");
                        d.put("weight", EMObject::FLOAT, "Weight of the slice that is being inserted. Default is 1.0.");
                        d.put("start_model", EMObject::EMDATA, "Start model");
                        d.put("start_model_weight", EMObject::FLOAT, "start_model_weight");
                        return d;
                }

                virtual float get_score(const unsigned int idx) { if ( quality_scores.size() > idx ) return quality_scores[idx].get_snr_normed_frc_integral(); else throw UnexpectedBehaviorException("The requested index was beyond the length of the quality scores vector."); }

                virtual float get_norm(const unsigned int idx) { if ( quality_scores.size() > idx ) return quality_scores[idx].get_norm();  else throw UnexpectedBehaviorException("The requested index was beyond the length of the quality scores vector."); }

                virtual void zero_memory();

          protected:
                EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() );
                
                void load_default_settings()
                {
//                      params.put("x_in", = 0;
//                      params.put("y_in", = 0;
//                      params["zsample"] = 0;
//                      params["ysample"] = 0;
//                      params["xsample"] = 0;
//                      params["mode"] = 2;
//                      params["hard"] = 0.05;
//                      params["sym"] = "c1";
//                      params["quiet"] = true;
//                      params["3damp"] = false;
                }

                void free_memory();

                void load_inserter();

                void load_interp_FRC_calculator();

                void do_insert_slice_work(const EMData* const input_slice, const Transform & euler);

                void print_stats(const vector<InterpolatedFRC::QualityScores>& scores);

                vector<InterpolatedFRC::QualityScores> quality_scores, prev_quality_scores;

                unsigned int image_idx;

                FourierPixelInserter3D* inserter;

                InterpolatedFRC* interp_FRC_calculator;

                bool slice_insertion_flag;
                bool slice_agreement_flag;

                float x_scale_factor, y_scale_factor, z_scale_factor;

                int max_input_dim;

                int output_x, output_y, output_z;
          private:
                FourierReconstructor( const FourierReconstructor& that );
                FourierReconstructor& operator=( const FourierReconstructor& );

        };

        class BaldwinWoolfordReconstructor : public FourierReconstructor
        {
                public:
                BaldwinWoolfordReconstructor() : W(0) {}

                virtual ~BaldwinWoolfordReconstructor() {
                        if ( W != 0 )
                                delete [] W;
                }

                virtual string get_name() const
                {
                        return "baldwinwoolford";
                }

                virtual string get_desc() const
                {
                        return "Reconstruction via direct Fourier inversion using gridding and delta function based weights";
                }

                static Reconstructor *NEW()
                {
                        return new BaldwinWoolfordReconstructor();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        d.put("mode", EMObject::INT, "Optional. Fourier pixel insertion mode [1-8] - mode 2 is default.");
                        d.put("x_in", EMObject::INT, "Necessary. The x dimension of the input images.");
                        d.put("y_in", EMObject::INT, "Necessary. The y dimension of the input images.");
                        d.put("zsample", EMObject::INT, "Optional. The z dimension (Fourier sampling) of the reconstructed volume, very useful for tomographic reconstruction. Works for general volumes.");
                        d.put("ysample", EMObject::INT, "Optional. The y dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
                        d.put("xsample", EMObject::INT, "Optional. The x dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
                        d.put("sym", EMObject::STRING, "Optional. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1");
                        d.put("maskwidth", EMObject::INT, "The width of the Fourier space kernel used to interpolate data to grid points" );
                        d.put("postmultiply", EMObject::BOOL, "A flag that controls whether or not the reconstructed volume is post multiplied in real space by the IFT of the weighting function. Default is on");
                        // Currently redundant
                        d.put("3damp", EMObject::BOOL, "this doesn't work, fixme dsaw");
                        d.put("hard", EMObject::FLOAT, "Optional. The quality metric threshold. Default is 0 (off).");
                        d.put("quiet", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
                        return d;
                }
                virtual EMData *finish();

                virtual int insert_slice_weights(const Transform& t3d);

                virtual int insert_slice(const EMData* const image, const Transform3D& t3d) {
                        throw UnexpectedBehaviorException("Transform3D interface is redundant");
                }
                virtual int insert_slice(const EMData* const image, const Transform& t3d);
                void insert_pixel(const float& x, const float& y, const float& z, const float dt[2]);

                void insert_density_at(const float& x, const float& y, const float& z);

                virtual void setup();

                protected:
                void load_default_settings()
                {
                        params["mode"] = 1;
                        params["x_in"] = 0;
                        params["y_in"] = 0;
                        params["zsample"] = 0;
                        params["ysample"] = 0;
                        params["xsample"] = 0;
                        params["sym"] = "c1";
                        params["maskwidth"] = 3;

                        // Currently redundant
                        params["3damp"] = false;
                        params["hard"] = 0.05;
                        params["quiet"] = false;
                }

                private:
                BaldwinWoolfordReconstructor( const BaldwinWoolfordReconstructor& that );
                BaldwinWoolfordReconstructor& operator=( const BaldwinWoolfordReconstructor& );

                float* W;
                float dfreq;

        };

        class WienerFourierReconstructor:public Reconstructor, public ReconstructorVolumeData
        {
          public:
                WienerFourierReconstructor() { load_default_settings(); };
                virtual ~WienerFourierReconstructor() {};

                virtual void setup();

                virtual int insert_slice(const EMData* const slice, const Transform3D & euler);

                virtual EMData *finish();

                virtual string get_name() const
                {
                        return "wiener_fourier";
                }

                virtual string get_desc() const
                {
                        return "Experimental - Direct Fourier reconstruction taking into account the Wiener filtration of the individual images.";
                }

                static Reconstructor *NEW()
                {
                        return new WienerFourierReconstructor();
                }

                virtual TypeDict get_param_types() const
                {
                        TypeDict d;
                        // FIXME:: double check all of thes are need, expecially dlog and weight
                        d.put("size", EMObject::INT);
                        d.put("mode", EMObject::INT);
                        d.put("weight", EMObject::FLOAT);
                        d.put("use_weights", EMObject::BOOL);
                        d.put("dlog", EMObject::BOOL);
                        d.put("padratio", EMObject::FLOAT);
                        d.put("snr", EMObject::FLOATARRAY);
                        return d;
                }
          private:
                // Disallow copy construction
                WienerFourierReconstructor( const WienerFourierReconstructor& that);
                // Disallow assignment
                WienerFourierReconstructor& operator=( const WienerFourierReconstructor& );

                void load_default_settings()
                {
                        params["size"] = 0;
                        params["mode"] = 2;
                        params["weight"] = 1.0;
                        params["use_weights"] = true;
                        params["dlog"] = false;
                        params["padratio"] = 1.0;
                }
        };

        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 Transform3D & euler) {
                        throw UnexpectedBehaviorException("Use of Tranform3D with BackProjection is deprecated"); }

                virtual int insert_slice(const EMData* const slice, const Transform & euler);

                virtual EMData *finish();

                virtual string get_name() const
                {
                        return "back_projection";
                }

                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::INT, "Necessary. The x and y dimensions of the input images.");
                        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("zsample", EMObject::INT, "Optional. The z dimensions of the reconstructed volume.");
                        return d;
                }
          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["size"] = 0;
                        params["sym"] = "c1";
                        params["zsample"] = 0;
                }

                EMData* preprocess_slice(const EMData* const slice, const Transform& t);
        };


        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);

               virtual EMData *finish();

                virtual string get_name() const
                {
                        return "nn4";
                }

                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, int mult=1 );


          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_result;
                bool m_delete_volume;
                bool m_delete_weight;
                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;
                }
        };


     /* 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);

                virtual EMData* finish();

                virtual string get_name() const
                {
                        return "nnSSNR";
                }

                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("w", EMObject::FLOAT);
                        return d;
                }

                void setup( const string& symmetry, int size, int npad);

                int insert_padfft_slice( EMData* padded, const Transform& trans, int mult=1 );

          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_wptr2;
                EMData* m_result;
                bool m_delete_volume;
                bool m_delete_weight;
                bool m_delete_weight2;
                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);

                virtual EMData *finish();

                virtual string get_name() const
                {
                        return "nn4_ctf";
                }

                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, const Transform& trans, int mult=1);

                int insert_buffed_slice( const EMData* buffer, int mult );
          private:
                EMData* m_volume;
                EMData* m_result;
                EMData* m_wptr;
                bool m_delete_volume;
                bool m_delete_weight;
                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();

        };


     /* 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);


                virtual EMData* finish();

                virtual string get_name() const
                {
                        return "nnSSNR_ctf";
                }

                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("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, int mult=1 );

          private:
                EMData* m_volume;
                EMData* m_wptr;
                EMData* m_wptr2;
                EMData* m_wptr3;
                EMData* m_result;
                bool m_delete_volume;
                bool m_delete_weight;
                bool m_delete_weight2;
                bool m_delete_weight3;
                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

/* vim: set ts=4 noet: */