00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef eman_reconstructor_h__
00037 #define eman_reconstructor_h__ 1
00038 #include <fstream>
00039 #include <boost/shared_ptr.hpp>
00040 #include "emdata.h"
00041 #include "exception.h"
00042 #include "emobject.h"
00043 #include "interp.h"
00044
00045 using std::vector;
00046 using std::map;
00047 using std::string;
00048 using boost::shared_ptr;
00049
00050 using std::cout;
00051 using std::cerr;
00052 using std::endl;
00053
00054 #include <utility>
00055 using std::pair;
00056
00057 #include "reconstructor_tools.h"
00058
00059 namespace EMAN
00060 {
00061
00062 class Transform3D;
00063 class EMData;
00064
00110 class Reconstructor : public FactoryBase
00111 {
00112 public:
00113 Reconstructor() {}
00114 virtual ~Reconstructor() {}
00117 virtual void setup() = 0;
00118
00126 virtual void setup_seed(EMData* seed,float seed_weight) {throw;}
00127
00136 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; }
00137
00146 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0) {throw;}
00147
00161 virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true ) { throw; }
00162
00167 virtual EMData *finish(bool doift=true) { throw; }
00168
00171 void print_params() const
00172 {
00173 std::cout << "Printing reconstructor params" << std::endl;
00174 for ( Dict::const_iterator it = params.begin(); it != params.end(); ++it )
00175 {
00176 std::cout << (it->first) << " " << (it->second).to_str() << std::endl;
00177 }
00178 std::cout << "Done printing reconstructor params" << std::endl;
00179 }
00180
00181
00182 EMObject& operator[]( const string& key ) { return params[key]; }
00183
00184 private:
00185
00186 Reconstructor(const Reconstructor& that);
00187 Reconstructor& operator=(const Reconstructor& );
00188
00189 };
00190
00202 class ReconstructorVolumeData
00203 {
00204 public:
00208 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) {}
00209
00212 virtual ~ReconstructorVolumeData() { free_memory(); }
00213
00216 const EMData* const get_emdata() { return image; }
00217 protected:
00218
00220 EMData* image;
00222 EMData* tmp_data;
00223
00224
00225 int nx,nx2;
00226 int ny,ny2;
00227 int nz,nz2;
00228
00229 int subnx;
00230 int subny;
00231 int subnz;
00232
00233 int subx0;
00234 int suby0;
00235 int subz0;
00236
00237 protected:
00243 void free_memory()
00244 {
00245 if (image != 0) {delete image; image = 0;}
00246 if ( tmp_data != 0 ) { delete tmp_data; tmp_data = 0; }
00247 }
00248
00254 virtual void normalize_threed(const bool sqrt_damp=false);
00255
00259 virtual void zero_memory()
00260 {
00261 if (tmp_data != 0 ) tmp_data->to_zero();
00262 if (image != 0 ) image->to_zero();
00263 }
00264
00265 private:
00267 ReconstructorVolumeData(const ReconstructorVolumeData& that);
00269 ReconstructorVolumeData& operator=(const ReconstructorVolumeData& );
00270
00271 };
00272
00278 class FourierReconstructorSimple2D : public Reconstructor, public ReconstructorVolumeData
00279 {
00280 public:
00281 FourierReconstructorSimple2D() {}
00282
00283 virtual ~FourierReconstructorSimple2D() { }
00284
00285 virtual void setup();
00286
00287 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00288
00289 virtual EMData *finish(bool doift=true);
00290
00291 virtual string get_name() const { return NAME; }
00292
00293 virtual string get_desc() const { return "performs 2D reconstruction"; }
00294
00295 static Reconstructor *NEW()
00296 {
00297 return new FourierReconstructorSimple2D();
00298 }
00299
00300
00301 virtual TypeDict get_param_types() const
00302 {
00303 TypeDict d;
00304 d.put("nx", EMObject::INT, "Necessary. The x dimension of the input images.");
00305
00306 return d;
00307 }
00308
00309 static const string NAME;
00310 };
00311
00312
00313
00363 class FourierReconstructor : public Reconstructor, public ReconstructorVolumeData
00364 {
00365 public:
00369 FourierReconstructor() { load_default_settings(); }
00370
00374 virtual ~FourierReconstructor() { free_memory(); }
00375
00379 virtual void setup();
00380
00389 virtual void setup_seed(EMData* seed,float seed_weight);
00390
00399 virtual EMData* preprocess_slice( const EMData* const slice, const Transform& t = Transform() );
00400
00410 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00411
00412
00429 virtual int determine_slice_agreement(EMData* slice, const Transform &euler, const float weight=1.0, bool sub=true );
00430
00437 virtual EMData *finish(bool doift=true);
00438
00441 virtual string get_name() const
00442 {
00443 return NAME;
00444 }
00445
00448 virtual string get_desc() const
00449 {
00450 return "Reconstruction via direct Fourier methods using one of a variety of different kernels, most of which are Gaussian based";
00451 }
00452
00456 static Reconstructor *NEW()
00457 {
00458 return new FourierReconstructor();
00459 }
00460
00464 virtual TypeDict get_param_types() const
00465 {
00466 TypeDict d;
00467 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.");
00468 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");
00469 d.put("mode", EMObject::STRING, "Optional. Fourier pixel insertion mode name. gauss_2 is the default.");
00470 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.");
00471 d.put("verbose", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
00472 d.put("subvolume",EMObject::INTARRAY, "Optional. (xorigin,yorigin,zorigin,xsize,ysize,zsize) all in Fourier pixels. Useful for parallelism.");
00473 d.put("savenorm",EMObject::STRING, "Debug. Will cause the normalization volume to be written directly to the specified file when finish() is called.");
00474 return d;
00475 }
00476
00477 static const string NAME;
00478
00479 protected:
00480
00483 void load_default_settings()
00484 {
00485 inserter=0;
00486 image=0;
00487 tmp_data=0;
00488 }
00489
00493 void free_memory();
00494
00497 void load_inserter();
00498
00504 void do_insert_slice_work(const EMData* const input_slice, const Transform & euler,const float weight);
00505
00510 void do_compare_slice_work(EMData* input_slice, const Transform & euler,float weight);
00511
00516 bool pixel_at(const float& xx, const float& yy, const float& zz, float *dt);
00517
00519 FourierPixelInserter3D* inserter;
00520
00521 private:
00524 FourierReconstructor( const FourierReconstructor& that );
00527 FourierReconstructor& operator=( const FourierReconstructor& );
00528
00529 };
00530
00531
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00645 class WienerFourierReconstructor:public Reconstructor, public ReconstructorVolumeData
00646 {
00647 public:
00648 WienerFourierReconstructor() { load_default_settings(); };
00649 virtual ~WienerFourierReconstructor() {};
00650
00651 virtual void setup();
00652
00661 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00662
00663 virtual EMData *finish(bool doift=true);
00664
00665 virtual string get_name() const
00666 {
00667 return NAME;
00668 }
00669
00670 virtual string get_desc() const
00671 {
00672 return "Experimental - Direct Fourier reconstruction taking into account the Wiener filtration of the individual images.";
00673 }
00674
00675 static Reconstructor *NEW()
00676 {
00677 return new WienerFourierReconstructor();
00678 }
00679
00680 virtual TypeDict get_param_types() const
00681 {
00682 TypeDict d;
00683
00684 d.put("size", EMObject::INT);
00685 d.put("mode", EMObject::INT);
00686 d.put("weight", EMObject::FLOAT);
00687 d.put("use_weights", EMObject::BOOL);
00688 d.put("dlog", EMObject::BOOL);
00689 d.put("padratio", EMObject::FLOAT);
00690 d.put("snr", EMObject::FLOATARRAY);
00691 return d;
00692 }
00693
00694 static const string NAME;
00695
00696 private:
00697
00698 WienerFourierReconstructor( const WienerFourierReconstructor& that);
00699
00700 WienerFourierReconstructor& operator=( const WienerFourierReconstructor& );
00701
00702 void load_default_settings()
00703 {
00704 params["size"] = 0;
00705 params["mode"] = 2;
00706 params["weight"] = 1.0;
00707 params["use_weights"] = true;
00708 params["dlog"] = false;
00709 params["padratio"] = 1.0;
00710 }
00711 };
00712
00720 class BackProjectionReconstructor:public Reconstructor, public ReconstructorVolumeData
00721 {
00722 public:
00723 BackProjectionReconstructor() { load_default_settings(); }
00724
00725 virtual ~BackProjectionReconstructor() {}
00726
00727 virtual void setup();
00728
00737 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00738
00739 virtual EMData *finish(bool doift=true);
00740
00741 virtual string get_name() const
00742 {
00743 return NAME;
00744 }
00745
00746 virtual string get_desc() const
00747 {
00748 return "Simple (unfiltered) back-projection reconstruction. Weighting by contributing particles in the class average is optional and default behaviour";
00749 }
00750
00751 static Reconstructor *NEW()
00752 {
00753 return new BackProjectionReconstructor();
00754 }
00755
00756 virtual TypeDict get_param_types() const
00757 {
00758 TypeDict d;
00759 d.put("size", EMObject::INT, "Necessary. The x and y dimensions of the input images.");
00760 d.put("weight", EMObject::FLOAT, "Optional. A temporary value set prior to slice insertion, indicative of the inserted slice's weight. Default sis 1.");
00761 d.put("sym", EMObject::STRING, "Optional. The symmetry to impose on the final reconstruction. Default is c1");
00762 d.put("zsample", EMObject::INT, "Optional. The z dimensions of the reconstructed volume.");
00763 return d;
00764 }
00765
00766 static const string NAME;
00767
00768 private:
00769
00770 BackProjectionReconstructor( const BackProjectionReconstructor& that);
00771
00772 BackProjectionReconstructor& operator=( const BackProjectionReconstructor& );
00773
00774 void load_default_settings()
00775 {
00776 params["weight"] = 1.0;
00777 params["use_weights"] = true;
00778 params["size"] = 0;
00779 params["sym"] = "c1";
00780 params["zsample"] = 0;
00781 }
00782
00783 EMData* preprocess_slice(const EMData* const slice, const Transform& t);
00784 };
00785
00786
00790 EMData* padfft_slice( const EMData* const slice, const Transform& t, int npad );
00791
00792 class nn4Reconstructor:public Reconstructor
00793 {
00794 public:
00795 nn4Reconstructor();
00796
00797 nn4Reconstructor( const string& symmetry, int size, int npad );
00798
00799 virtual ~nn4Reconstructor();
00800
00801 virtual void setup();
00802
00811 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00812
00813 virtual EMData *finish(bool doift=true);
00814
00815 virtual string get_name() const
00816 {
00817 return NAME;
00818 }
00819
00820 virtual string get_desc() const
00821 {
00822 return "Direct Fourier inversion routine";
00823 }
00824
00825 static Reconstructor *NEW()
00826 {
00827 return new nn4Reconstructor();
00828 }
00829
00830 virtual TypeDict get_param_types() const
00831 {
00832 TypeDict d;
00833 d.put("size", EMObject::INT);
00834 d.put("npad", EMObject::INT);
00835 d.put("sign", EMObject::INT);
00836 d.put("ndim", EMObject::INT);
00837 d.put("snr", EMObject::FLOAT);
00838 d.put("symmetry", EMObject::STRING);
00839 d.put("snr", EMObject::FLOAT);
00840 d.put("fftvol", EMObject::EMDATA);
00841 d.put("weight", EMObject::EMDATA);
00842 d.put("weighting", EMObject::INT);
00843 return d;
00844 }
00845
00846 void setup( const string& symmetry, int size, int npad );
00847
00848 int insert_padfft_slice( EMData* padded, const Transform& trans, int mult=1 );
00849
00850 static const string NAME;
00851
00852 private:
00853 EMData* m_volume;
00854 EMData* m_wptr;
00855 EMData* m_result;
00856 bool m_delete_volume;
00857 bool m_delete_weight;
00858 string m_symmetry;
00859 int m_weighting;
00860 int m_vnx, m_vny, m_vnz;
00861 int m_npad;
00862 int m_nsym;
00863 int m_ndim;
00864 int m_vnzp, m_vnyp, m_vnxp;
00865 int m_vnzc, m_vnyc, m_vnxc;
00866 void buildFFTVolume();
00867 void buildNormVolume();
00868 float m_wghta;
00869 float m_wghtb;
00870 float m_osnr;
00871 void load_default_settings()
00872 {
00873
00874 }
00875 };
00876
00877
00878
00879
00880
00881
00882 class nnSSNR_Reconstructor:public Reconstructor
00883 {
00884
00885 public:
00886 nnSSNR_Reconstructor();
00887
00888 nnSSNR_Reconstructor( const string& symmetry, int size, int npad);
00889
00890 ~nnSSNR_Reconstructor();
00891
00892 virtual void setup();
00893
00902 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00903
00904 virtual EMData *finish(bool doift=true);
00905
00906 virtual string get_name() const
00907 {
00908 return NAME;
00909 }
00910
00911 virtual string get_desc() const
00912 {
00913 return "Reconstruction by nearest neighbor with 3D SSNR";
00914 }
00915
00916 static Reconstructor *NEW()
00917 {
00918 return new nnSSNR_Reconstructor();
00919 }
00920
00921 virtual TypeDict get_param_types() const
00922 {
00923 TypeDict d;
00924 d.put("size", EMObject::INT);
00925 d.put("npad", EMObject::INT);
00926 d.put("symmetry", EMObject::STRING);
00927 d.put("fftvol", EMObject::EMDATA);
00928 d.put("weight", EMObject::EMDATA);
00929 d.put("weight2", EMObject::EMDATA);
00930 d.put("SSNR", EMObject::EMDATA);
00931 d.put("w", EMObject::FLOAT);
00932 return d;
00933 }
00934
00935 void setup( const string& symmetry, int size, int npad);
00936
00937 int insert_padfft_slice( EMData* padded, const Transform& trans, int mult=1 );
00938
00939 static const string NAME;
00940
00941 private:
00942 EMData* m_volume;
00943 EMData* m_wptr;
00944 EMData* m_wptr2;
00945 EMData* m_result;
00946 bool m_delete_volume;
00947 bool m_delete_weight;
00948 bool m_delete_weight2;
00949 string m_symmetry;
00950 int m_weighting;
00951 int m_vnx, m_vny, m_vnz;
00952 int m_npad;
00953 int m_nsym;
00954 int m_vnzp, m_vnyp, m_vnxp;
00955 int m_vnzc, m_vnyc, m_vnxc;
00956 void buildFFTVolume();
00957 void buildNormVolume();
00958 void buildNorm2Volume();
00959 float m_wghta;
00960 float m_wghtb;
00961 };
00962
00963
00967 class nn4_ctfReconstructor:public Reconstructor
00968 {
00969 public:
00970 nn4_ctfReconstructor();
00971
00972 nn4_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign );
00973
00974 virtual ~nn4_ctfReconstructor();
00975
00976 virtual void setup();
00977
00987 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
00988
00989 virtual EMData *finish(bool doift=true);
00990
00991 virtual string get_name() const
00992 {
00993 return NAME;
00994 }
00995
00996 virtual string get_desc() const
00997 {
00998 return "Direct Fourier inversion reconstruction routine";
00999 }
01000
01001 static Reconstructor *NEW()
01002 {
01003 return new nn4_ctfReconstructor();
01004 }
01005
01006
01007 TypeDict get_param_types() const
01008 {
01009 TypeDict d;
01010 d.put("size", EMObject::INT);
01011 d.put("npad", EMObject::INT);
01012 d.put("sign", EMObject::INT);
01013 d.put("symmetry", EMObject::STRING);
01014 d.put("snr", EMObject::FLOAT);
01015 d.put("fftvol", EMObject::EMDATA);
01016 d.put("weight", EMObject::EMDATA);
01017 d.put("weighting", EMObject::INT);
01018 d.put("varsnr", EMObject::INT);
01019 return d;
01020 }
01021
01022 void setup( const string& symmetry, int size, int npad, float snr, int sign );
01023
01024 int insert_padfft_slice( EMData* padfft, const Transform& trans, int mult=1);
01025
01026 int insert_buffed_slice( const EMData* buffer, int mult );
01027
01028 static const string NAME;
01029
01030 private:
01031 EMData* m_volume;
01032 EMData* m_result;
01033 EMData* m_wptr;
01034 bool m_delete_volume;
01035 bool m_delete_weight;
01036 int m_vnx, m_vny, m_vnz;
01037 int m_vnzp, m_vnyp, m_vnxp;
01038 int m_vnxc, m_vnyc, m_vnzc;
01039 int m_npad;
01040 int m_sign;
01041 int m_varsnr;
01042 int m_weighting;
01043 float m_wghta, m_wghtb;
01044 float m_snr;
01045 string m_symmetry;
01046 int m_nsym;
01047
01048 void buildFFTVolume();
01049 void buildNormVolume();
01050
01051 };
01052
01053
01054
01055
01056
01057
01058 class nnSSNR_ctfReconstructor:public Reconstructor
01059 {
01060
01061 public:
01062 nnSSNR_ctfReconstructor();
01063
01064 nnSSNR_ctfReconstructor( const string& symmetry, int size, int npad, float snr, int sign);
01065
01066 ~nnSSNR_ctfReconstructor();
01067
01068 virtual void setup();
01069
01079 virtual int insert_slice(const EMData* const slice, const Transform & euler,const float weight=1.0);
01080
01081
01082 virtual EMData *finish(bool doift=true);
01083
01084 virtual string get_name() const
01085 {
01086 return NAME;
01087 }
01088
01089 virtual string get_desc() const
01090 {
01091 return "Reconstruction by nearest neighbor with 3D SSNR with CTF";
01092 }
01093
01094 static Reconstructor *NEW()
01095 {
01096 return new nnSSNR_ctfReconstructor();
01097 }
01098
01099 TypeDict get_param_types() const
01100 {
01101 TypeDict d;
01102 d.put("size", EMObject::INT);
01103 d.put("npad", EMObject::INT);
01104 d.put("symmetry", EMObject::STRING);
01105 d.put("fftvol", EMObject::EMDATA);
01106 d.put("fftwvol", EMObject::EMDATA);
01107 d.put("weight", EMObject::EMDATA);
01108 d.put("weight2", EMObject::EMDATA);
01109 d.put("weight3", EMObject::EMDATA);
01110 d.put("SSNR", EMObject::EMDATA);
01111 d.put("w", EMObject::FLOAT);
01112 d.put("sign", EMObject::INT);
01113 d.put("snr", EMObject::FLOAT);
01114 return d;
01115 }
01116 void setup( const string& symmetry, int size, int npad, float snr, int sign);
01117
01118 int insert_padfft_slice( EMData* padded, const Transform& trans, int mult=1 );
01119
01120 static const string NAME;
01121
01122 private:
01123 EMData* m_volume;
01124 EMData* m_wptr;
01125 EMData* m_wptr2;
01126 EMData* m_wptr3;
01127 EMData* m_result;
01128 bool m_delete_volume;
01129 bool m_delete_weight;
01130 bool m_delete_weight2;
01131 bool m_delete_weight3;
01132 string m_symmetry;
01133 int m_weighting;
01134 int m_vnx, m_vny, m_vnz;
01135 int m_npad;
01136 int m_nsym;
01137 int m_vnzp, m_vnyp, m_vnxp;
01138 int m_vnzc, m_vnyc, m_vnxc;
01139 void buildFFTVolume();
01140 void buildNormVolume();
01141 void buildNorm2Volume();
01142 void buildNorm3Volume();
01143 float m_wghta;
01144 float m_wghtb;
01145 int m_sign;
01146 float m_snr;
01147 int wiener;
01148 };
01149
01150 template <> Factory < Reconstructor >::Factory();
01151
01152 void dump_reconstructors();
01153 map<string, vector<string> > dump_reconstructors_list();
01154
01155
01156 struct point_t
01157 {
01158 int pos2;
01159 float real;
01160 float imag;
01161 float ctf2;
01162 };
01163
01164
01165 class newfile_store
01166 {
01167 public:
01168 newfile_store( const string& prefix, int npad, bool ctf );
01169
01170 virtual ~newfile_store();
01171
01172 void add_image( EMData* data, const Transform& tf );
01173
01174 void add_tovol( EMData* fftvol, EMData* wgtvol, const vector<int>& mults, int pbegin, int pend );
01175
01176 void get_image( int id, EMData* buf );
01177
01178 void read( int nprj );
01179
01180 void restart( );
01181
01182 private:
01183 int m_npad;
01184
01185 bool m_ctf;
01186
01187 string m_bin_file;
01188 string m_txt_file;
01189
01190 shared_ptr<std::ofstream> m_bin_of;
01191 shared_ptr<std::ofstream> m_txt_of;
01192 shared_ptr<std::ifstream> m_bin_if;
01193 vector< std::istream::off_type > m_offsets;
01194
01195 vector< point_t > m_points;
01196 };
01197
01198 class file_store
01199 {
01200 public:
01201 file_store(const string& filename, int npad, int write, bool CTF);
01202
01203 virtual ~file_store();
01204
01205 void add_image(EMData* data, const Transform& tf);
01206
01207 void get_image(int id, EMData* padfft);
01208
01209 void restart();
01210 private:
01211 shared_ptr<std::ifstream> m_ihandle;
01212 shared_ptr<std::ofstream> m_bin_ohandle;
01213 shared_ptr<std::ofstream> m_txt_ohandle;
01214 string m_bin_file;
01215 string m_txt_file;
01216 int m_ctf;
01217 int m_npad;
01218 int m_prev;
01219 int m_x_out;
01220 int m_y_out;
01221 int m_z_out;
01222 int m_write;
01223 std::istream::off_type m_totsize;
01224 float m_Cs;
01225 float m_pixel;
01226 float m_voltage;
01227 float m_ctf_applied;
01228 float m_amp_contrast;
01229 vector< float > m_defocuses;
01230 vector< float > m_phis;
01231 vector< float > m_thetas;
01232 vector< float > m_psis;
01233 };
01234
01235 }
01236
01237 #endif
01238
01239