Clipper
|
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