Clipper
fftmap.h
00001 
00004 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
00005 //L
00006 //L  This library is free software and is distributed under the terms
00007 //L  and conditions of version 2.1 of the GNU Lesser General Public
00008 //L  Licence (LGPL) with the following additional clause:
00009 //L
00010 //L     `You may also combine or link a "work that uses the Library" to
00011 //L     produce a work containing portions of the Library, and distribute
00012 //L     that work under terms of your choice, provided that you give
00013 //L     prominent notice with each copy of the work that the specified
00014 //L     version of the Library is used in it, and that you include or
00015 //L     provide public access to the complete corresponding
00016 //L     machine-readable source code for the Library including whatever
00017 //L     changes were used in the work. (i.e. If you make changes to the
00018 //L     Library you must distribute those, but you do not need to
00019 //L     distribute source or object code to those portions of the work
00020 //L     not covered by this licence.)'
00021 //L
00022 //L  Note that this clause grants an additional right and does not impose
00023 //L  any additional restriction, and so does not affect compatibility
00024 //L  with the GNU General Public Licence (GPL). If you wish to negotiate
00025 //L  other terms, please contact the maintainer.
00026 //L
00027 //L  You can redistribute it and/or modify the library under the terms of
00028 //L  the GNU Lesser General Public License as published by the Free Software
00029 //L  Foundation; either version 2.1 of the License, or (at your option) any
00030 //L  later version.
00031 //L
00032 //L  This library is distributed in the hope that it will be useful, but
00033 //L  WITHOUT ANY WARRANTY; without even the implied warranty of
00034 //L  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00035 //L  Lesser General Public License for more details.
00036 //L
00037 //L  You should have received a copy of the CCP4 licence and/or GNU
00038 //L  Lesser General Public License along with this library; if not, write
00039 //L  to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
00040 //L  The GNU Lesser General Public can also be obtained by writing to the
00041 //L  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00042 //L  MA 02111-1307 USA
00043 
00044 
00045 #ifndef CLIPPER_FFTMAP
00046 #define CLIPPER_FFTMAP
00047 
00048 
00049 #include "coords.h"
00050 
00051 #include <complex>
00052 
00053 
00054 namespace clipper
00055 {
00056   // fft type
00057   typedef float ffttype;
00058 
00059   // forward definition
00060   namespace datatypes {
00061     template<class T> class F_phi;
00062   }
00063 
00064   // base class for FFT classes
00065   class FFTmap_base {
00066   public:
00067     enum FFTtype { Default, Measure, Estimate };  
00068   protected:
00069     static Mutex mutex;                 
00070   };
00071 
00073 
00080   class FFTmap_p1 : public FFTmap_base
00081   {
00082   public:
00084     FFTmap_p1();
00086     FFTmap_p1( const FFTmap_p1& other ) { copy(other); }
00088     FFTmap_p1( const Grid_sampling& grid_sam, const FFTtype type = Default );
00090     const FFTmap_p1& operator =(const FFTmap_p1& other) { return copy(other); }
00092     void init( const Grid_sampling& grid_sam, const FFTtype type = Default );
00094     void reset();
00095 
00097     const Grid_sampling& grid_real() const { return grid_sam_; }
00099     const Grid& grid_reci() const { return grid_reci_; }
00101     bool uniq_reci( const Coord_grid& c ) const { return ( (c.w()>0 && c.w()<grid_half_.nw()) || (c.w()<=grid_half_.nw() && ( (c.v()>0 && c.v()<grid_half_.nv()) || (c.v()<=grid_half_.nv() && c.u()<=grid_half_.nu()) ) ) ); }
00102 
00104     void fft_h_to_x( const ftype& scale );
00106     void fft_x_to_h( const ftype& scale );
00107 
00109     std::complex<ffttype> get_hkl( const HKL& hkl ) const;
00111     void set_hkl( const HKL& hkl, const std::complex<ffttype>& f );
00113     const std::complex<ffttype>& cplx_data( const Coord_grid& c ) const
00114       { return data_c[ grid_reci_.index( c ) ]; }
00116     std::complex<ffttype>& cplx_data( const Coord_grid& c )
00117       { return data_c[ grid_reci_.index( c ) ]; }
00119     const ffttype& real_data( const Coord_grid& c ) const
00120       { return data_r[ grid_real_.index( c ) ]; }
00122     ffttype& real_data( const Coord_grid& c )
00123       { return data_r[ grid_real_.index( c ) ]; }
00124 
00126     static FFTtype& default_type() { return default_type_; }
00127 
00128     void debug() const;
00129 
00130   protected:
00131     const FFTmap_p1& copy( const FFTmap_p1& other );  
00132 
00133     enum FFTmode { NONE, RECI, REAL, OTHER };  
00134 
00135     FFTmode mode;                       
00136     FFTtype type_;                      
00137     Grid_sampling grid_sam_;            
00138     Grid grid_reci_;                    
00139     Grid grid_real_;                    
00140     Grid grid_half_;                    
00141 
00142     Matrix<char> req_kl, req_uv;        
00143     std::vector<char> req_l, req_u;     
00144 
00145     std::vector<ffttype> datavec;       
00146     ffttype* data_r;                    
00147     std::complex<ffttype>* data_c;      
00148 
00149     static FFTtype default_type_;       
00150   };
00151 
00152 
00154 
00165   class FFTmap : private FFTmap_p1
00166   {
00167   public:
00169     FFTmap();
00171     FFTmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00173     void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00175     void reset();
00176 
00178     const Cell& cell() const { return cell_; }
00180     const Spacegroup& spacegroup() const { return spacegroup_; }
00182     const Grid_sampling& grid_sampling() const { return FFTmap_p1::grid_real(); }
00184     void fft_h_to_x();
00186     void fft_x_to_h();
00187 
00189     template<class T> void get_recip_data( const HKL& rfl, datatypes::F_phi<T>& fphi ) const;
00191     template<class T> void set_recip_data( const HKL& rfl, const datatypes::F_phi<T>& fphi );
00192 
00194     template<class T> void get_real_data( const Coord_grid& c, T& datum ) const;
00196     template<class T> void set_real_data( const Coord_grid& c, const T& datum );
00197 
00199     datatypes::F_phi<ffttype> get_recip_data( const HKL& rfl ) const;
00201     const ffttype& get_real_data( const Coord_grid& c ) const
00202       { return real_data(c.unit(grid_real())); }
00203 
00205     template<class H, class X> void fft_rfl_to_map( const H& h, X& x );
00207     template<class H, class X> void fft_map_to_rfl( const X& x, H& h );
00208 
00209     void debug() const;
00210 
00211   protected:
00212     Cell cell_;                         
00213     Spacegroup spacegroup_;             
00214     std::vector<Isymop> isymop;         
00215   };
00216 
00217 
00218   // template implementations
00219 
00220 
00232   template<class H, class X> void FFTmap::fft_rfl_to_map( const H& h, X& x )
00233   {
00234     // zero the map
00235     reset();
00236 
00237     // copy from reflection data
00238     typename H::HKL_reference_index ih;
00239     for ( ih = h.first_data(); !ih.last(); h.next_data( ih ) )
00240       set_recip_data( ih.hkl(), h[ih] );
00241 
00242     // fft
00243     fft_h_to_x();
00244 
00245     // and copy into the map
00246     typename X::Map_reference_index ix;
00247     for ( ix = x.first(); !ix.last(); ix.next() )
00248       get_real_data( ix.coord(), x[ix] );
00249   }
00250 
00251 
00264   template<class H, class X> void FFTmap::fft_map_to_rfl( const X& x, H& h )
00265   {
00266     // zero the map
00267     reset();
00268 
00269     // copy into the map
00270     typename X::Map_reference_index ix;
00271     for ( ix = x.first(); !ix.last(); ix.next() )
00272       set_real_data( ix.coord(), x[ix] );
00273 
00274     // fft
00275     fft_x_to_h();
00276 
00277     // now fill it
00278     typename H::HKL_reference_index ih;
00279     for ( ih = h.first(); !ih.last(); ih.next() )
00280       get_recip_data( ih.hkl(), h[ih] );
00281   }
00282 
00283 
00284 } // namespace clipper
00285 
00286 #endif