~mok0/clipper/trunk

« back to all changes in this revision

Viewing changes to clipper/core/fftmap.h

  • Committer: Morten Kjeldgaard
  • Date: 2011-02-09 21:25:56 UTC
  • Revision ID: mok@bioxray.dk-20110209212556-3zbon5ukh4zf3x00
Tags: 2.0.3
versionĀ 2.0.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*! \file lib/fftmap.h
 
2
    Header file for P1 fft map
 
3
*/
 
4
//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
 
5
//L
 
6
//L  This library is free software and is distributed under the terms
 
7
//L  and conditions of version 2.1 of the GNU Lesser General Public
 
8
//L  Licence (LGPL) with the following additional clause:
 
9
//L
 
10
//L     `You may also combine or link a "work that uses the Library" to
 
11
//L     produce a work containing portions of the Library, and distribute
 
12
//L     that work under terms of your choice, provided that you give
 
13
//L     prominent notice with each copy of the work that the specified
 
14
//L     version of the Library is used in it, and that you include or
 
15
//L     provide public access to the complete corresponding
 
16
//L     machine-readable source code for the Library including whatever
 
17
//L     changes were used in the work. (i.e. If you make changes to the
 
18
//L     Library you must distribute those, but you do not need to
 
19
//L     distribute source or object code to those portions of the work
 
20
//L     not covered by this licence.)'
 
21
//L
 
22
//L  Note that this clause grants an additional right and does not impose
 
23
//L  any additional restriction, and so does not affect compatibility
 
24
//L  with the GNU General Public Licence (GPL). If you wish to negotiate
 
25
//L  other terms, please contact the maintainer.
 
26
//L
 
27
//L  You can redistribute it and/or modify the library under the terms of
 
28
//L  the GNU Lesser General Public License as published by the Free Software
 
29
//L  Foundation; either version 2.1 of the License, or (at your option) any
 
30
//L  later version.
 
31
//L
 
32
//L  This library is distributed in the hope that it will be useful, but
 
33
//L  WITHOUT ANY WARRANTY; without even the implied warranty of
 
34
//L  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
35
//L  Lesser General Public License for more details.
 
36
//L
 
37
//L  You should have received a copy of the CCP4 licence and/or GNU
 
38
//L  Lesser General Public License along with this library; if not, write
 
39
//L  to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
 
40
//L  The GNU Lesser General Public can also be obtained by writing to the
 
41
//L  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 
42
//L  MA 02111-1307 USA
 
43
 
 
44
 
 
45
#ifndef CLIPPER_FFTMAP
 
46
#define CLIPPER_FFTMAP
 
47
 
 
48
 
 
49
#include "coords.h"
 
50
 
 
51
#include <complex>
 
52
 
 
53
 
 
54
namespace clipper
 
55
{
 
56
  // fft type
 
57
  typedef float ffttype;
 
58
 
 
59
  // forward definition
 
60
  namespace datatypes {
 
61
    template<class T> class F_phi;
 
62
  }
 
63
 
 
64
 
 
65
  //! FFTmap_p1: low level P1 map used for calculating FFTs
 
66
  /*! This is a pure real P1 map, with an extra section in reciprocal
 
67
    space to allow generation of the full set of resiprocal space
 
68
    magnitudes. Access is by Coord_grid in both spaces, and indices
 
69
    must be non-negative and in range. The first and last sections
 
70
    along the half-length direction only have half the elements
 
71
    stored, the contents of the other half is ignored.
 
72
  */
 
73
  class FFTmap_p1
 
74
  {
 
75
  public:
 
76
    enum FFTtype { Default, Measure, Estimate };  //!< optimisation options
 
77
    //! Null constructor
 
78
    FFTmap_p1();
 
79
    //! Copy constructor
 
80
    FFTmap_p1( const FFTmap_p1& other ) { copy(other); }
 
81
    //! Constructor: takes grid
 
82
    FFTmap_p1( const Grid_sampling& grid_sam, const FFTtype type = Default );
 
83
    //! Assignment operator
 
84
    const FFTmap_p1& operator =(const FFTmap_p1& other) { return copy(other); }
 
85
    //! initialiser: takes grid
 
86
    void init( const Grid_sampling& grid_sam, const FFTtype type = Default );
 
87
    //! Reset
 
88
    void reset();
 
89
 
 
90
    //! Return real space grid.
 
91
    const Grid_sampling& grid_real() const { return grid_sam_; }
 
92
    //! Return reciprocal space grid (i.e. half real grid + 1 section).
 
93
    const Grid& grid_reci() const { return grid_reci_; }
 
94
    //! Test whether a coordinate is in the valid part of the recip. grid.
 
95
    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()) ) ) ); }
 
96
 
 
97
    //! Transform to real space
 
98
    void fft_h_to_x( const ftype& scale );
 
99
    //! Transform to reciprocal space
 
100
    void fft_x_to_h( const ftype& scale );
 
101
 
 
102
    //! get reciprocal space data: slow form with hemisphere check
 
103
    std::complex<ffttype> get_hkl( const HKL& hkl ) const;
 
104
    //! set reciprocal space data: slow form with hemisphere check
 
105
    void set_hkl( const HKL& hkl, const std::complex<ffttype>& f );
 
106
    //! get reciprocal space data
 
107
    const std::complex<ffttype>& cplx_data( const Coord_grid& c ) const
 
108
      { return data_c[ grid_reci_.index( c ) ]; }
 
109
    //! set reciprocal space data
 
110
    std::complex<ffttype>& cplx_data( const Coord_grid& c )
 
111
      { return data_c[ grid_reci_.index( c ) ]; }
 
112
    //! get real space data
 
113
    const ffttype& real_data( const Coord_grid& c ) const
 
114
      { return data_r[ grid_real_.index( c ) ]; }
 
115
    //! set real space data
 
116
    ffttype& real_data( const Coord_grid& c )
 
117
      { return data_r[ grid_real_.index( c ) ]; }
 
118
 
 
119
    //! set/get default optimisation type
 
120
    static FFTtype& default_type() { return default_type_; }
 
121
 
 
122
    void debug() const;
 
123
 
 
124
  protected:
 
125
    const FFTmap_p1& copy( const FFTmap_p1& other );  //!< copy function
 
126
 
 
127
    enum FFTmode { NONE, RECI, REAL, OTHER };  //!< space enumeration
 
128
 
 
129
    FFTmode mode;                       //!< real or reciprocal space?
 
130
    FFTtype type_;                      //!< optimisation options
 
131
    Grid_sampling grid_sam_;            //!< unit cell grid
 
132
    Grid grid_reci_;                    //!< reciprocal space grid
 
133
    Grid grid_real_;                    //!< real space grid
 
134
    Grid grid_half_;                    //!< half grid (for marking unique)
 
135
 
 
136
    Matrix<char> req_kl, req_uv;        //!< reci section lookup
 
137
    std::vector<char> req_l, req_u;     //!< real section lookup
 
138
 
 
139
    std::vector<ffttype> datavec;       //!< vector for the data
 
140
    ffttype* data_r;                    //!< pointer to real data
 
141
    std::complex<ffttype>* data_c;      //!< pointer to complex data
 
142
 
 
143
    static FFTtype default_type_;       //!< default optimisation options
 
144
  };
 
145
 
 
146
 
 
147
  //! FFTmap: P1 map with symmetry used for calculating FFTs
 
148
  /*! The FFTmap is represented in P1 in memory. However, it also has
 
149
    a spacegroup, and the contained data remains consistent with this
 
150
    spacegroup at all times. It has three states - unassigned,
 
151
    real-space, and reciprocal space. In real space it contains real
 
152
    map data. In reciprocal space it holds a hemisphere of complex
 
153
    structure factors, with the Friedels duplicated on the zero
 
154
    section.
 
155
 
 
156
    The user should be able to ignore all the issues of spacegroup
 
157
    symmetry, Friedel opposites, and storage order.
 
158
  */
 
159
  class FFTmap : private FFTmap_p1
 
160
  {
 
161
  public:
 
162
    //! Null constructor
 
163
    FFTmap();
 
164
    //! Constructor: takes spacegroup, cell, grid
 
165
    FFTmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
 
166
    //! initialiser
 
167
    void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
 
168
    //! Reset
 
169
    void reset();
 
170
 
 
171
    //! get the cell
 
172
    const Cell& cell() const { return cell_; }
 
173
    //! get the spacegroup
 
174
    const Spacegroup& spacegroup() const { return spacegroup_; }
 
175
    //! get the cell grid
 
176
    const Grid_sampling& grid_sampling() const { return FFTmap_p1::grid_real(); }
 
177
    //! Transform to real space
 
178
    void fft_h_to_x();
 
179
    //! Transform to reciprocal space
 
180
    void fft_x_to_h();
 
181
 
 
182
    //! get reciprocal space data
 
183
    template<class T> void get_recip_data( const HKL& rfl, datatypes::F_phi<T>& fphi ) const;
 
184
    //! set reciprocal space data
 
185
    template<class T> void set_recip_data( const HKL& rfl, const datatypes::F_phi<T>& fphi );
 
186
 
 
187
    //! get real space data
 
188
    template<class T> void get_real_data( const Coord_grid& c, T& datum ) const;
 
189
    //! set real space data
 
190
    template<class T> void set_real_data( const Coord_grid& c, const T& datum );
 
191
 
 
192
    //! get reciprocal space data (No error checking)
 
193
    datatypes::F_phi<ffttype> get_recip_data( const HKL& rfl ) const;
 
194
    //! get real space data (No error checking)
 
195
    const ffttype& get_real_data( const Coord_grid& c ) const
 
196
      { return real_data(c.unit(grid_real())); }
 
197
 
 
198
    //! calculate map-like object from reflection-like object
 
199
    template<class H, class X> void fft_rfl_to_map( const H& h, X& x );
 
200
    //! calculate reflection-like object from map-like object
 
201
    template<class H, class X> void fft_map_to_rfl( const X& x, H& h );
 
202
 
 
203
    void debug() const;
 
204
 
 
205
  protected:
 
206
    Cell cell_;                         //!< unit cell
 
207
    Spacegroup spacegroup_;             //!< spacegroup
 
208
    std::vector<Isymop> isymop;         //!< Integerised symops
 
209
  };
 
210
 
 
211
 
 
212
  // template implementations
 
213
 
 
214
 
 
215
  /*! Fill this FFTmap object from a reflection object, transform it,
 
216
    and fill the given map object from the FFTmap. This will work for
 
217
    any reflection data object which implements a HKL_reference_index,
 
218
    and every map data object which implements a Map_reference_index.
 
219
 
 
220
    For the results to be sensible, the spacegroup, cell and grids
 
221
    should match. (The map will be zeroed if necessary).
 
222
 
 
223
    \param h The source reflection data object.
 
224
    \param x The target map object.
 
225
  */
 
226
  template<class H, class X> void FFTmap::fft_rfl_to_map( const H& h, X& x )
 
227
  {
 
228
    // zero the map
 
229
    reset();
 
230
 
 
231
    // copy from reflection data
 
232
    typename H::HKL_reference_index ih;
 
233
    for ( ih = h.first_data(); !ih.last(); h.next_data( ih ) )
 
234
      set_recip_data( ih.hkl(), h[ih] );
 
235
 
 
236
    // fft
 
237
    fft_h_to_x();
 
238
 
 
239
    // and copy into the map
 
240
    typename X::Map_reference_index ix;
 
241
    for ( ix = x.first(); !ix.last(); ix.next() )
 
242
      get_real_data( ix.coord(), x[ix] );
 
243
  }
 
244
 
 
245
 
 
246
  /*! Fill this FFTmap object from a map object, transform it, and
 
247
    fill the given reflection object from the FFTmap. This will work
 
248
    for any reflection data object which implements a
 
249
    HKL_reference_index, and every map data object which implements a
 
250
    Map_reference_index.
 
251
 
 
252
    For the results to be sensible, the spacegroup, cell and grids
 
253
    should match. (The map will be zeroed if necessary).
 
254
 
 
255
    \param x The source map object.
 
256
    \param h The target reflection data object.
 
257
  */
 
258
  template<class H, class X> void FFTmap::fft_map_to_rfl( const X& x, H& h )
 
259
  {
 
260
    // zero the map
 
261
    reset();
 
262
 
 
263
    // copy into the map
 
264
    typename X::Map_reference_index ix;
 
265
    for ( ix = x.first(); !ix.last(); ix.next() )
 
266
      set_real_data( ix.coord(), x[ix] );
 
267
 
 
268
    // fft
 
269
    fft_x_to_h();
 
270
 
 
271
    // now fill it
 
272
    typename H::HKL_reference_index ih;
 
273
    for ( ih = h.first(); !ih.last(); ih.next() )
 
274
      get_recip_data( ih.hkl(), h[ih] );
 
275
  }
 
276
 
 
277
 
 
278
} // namespace clipper
 
279
 
 
280
#endif