1 /* convolution_search.cpp: convolution search implementation */
2 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
3 //L
4 //L This library is free software and is distributed under the terms
5 //L and conditions of version 2.1 of the GNU Lesser General Public
6 //L Licence (LGPL) with the following additional clause:
7 //L
8 //L `You may also combine or link a "work that uses the Library" to
9 //L produce a work containing portions of the Library, and distribute
10 //L that work under terms of your choice, provided that you give
11 //L prominent notice with each copy of the work that the specified
12 //L version of the Library is used in it, and that you include or
13 //L provide public access to the complete corresponding
14 //L machine-readable source code for the Library including whatever
15 //L changes were used in the work. (i.e. If you make changes to the
16 //L Library you must distribute those, but you do not need to
17 //L distribute source or object code to those portions of the work
18 //L not covered by this licence.)'
19 //L
20 //L Note that this clause grants an additional right and does not impose
21 //L any additional restriction, and so does not affect compatibility
22 //L with the GNU General Public Licence (GPL). If you wish to negotiate
23 //L other terms, please contact the maintainer.
24 //L
25 //L You can redistribute it and/or modify the library under the terms of
26 //L the GNU Lesser General Public License as published by the Free Software
27 //L Foundation; either version 2.1 of the License, or (at your option) any
28 //L later version.
29 //L
30 //L This library is distributed in the hope that it will be useful, but
31 //L WITHOUT ANY WARRANTY; without even the implied warranty of
32 //L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
33 //L Lesser General Public License for more details.
34 //L
35 //L You should have received a copy of the CCP4 licence and/or GNU
36 //L Lesser General Public License along with this library; if not, write
37 //L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
38 //L The GNU Lesser General Public can also be obtained by writing to the
39 //L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
40 //L MA 02111-1307 USA
41
42 #include "convolution_search.h"
43
44 #include "../core/map_interp.h"
45
46
47 namespace clipper {
48
49
50 // general forms
51 // slow Convolution_search
52
operator ()(Xmap<T> & result,const NXmap<T> & srchval) const53 template<class T> bool Convolution_search_slow<T>::operator() ( Xmap<T>& result, const NXmap<T>& srchval ) const
54 {
55 NX_operator nxop( result, srchval, RTop_orth::identity() );
56 return (*this)( result, srchval, nxop );
57 }
58
operator ()(Xmap<T> & result,const NXmap<T> & srchval,const NX_operator & nxop) const59 template<class T> bool Convolution_search_slow<T>::operator() ( Xmap<T>& result, const NXmap<T>& srchval, const NX_operator& nxop ) const
60 {
61 const Grid_sampling& g = result.grid_sampling();
62
63 // calculate extent of mask function in Xmap space
64 Coord_frac cf;
65 Range<ftype> urange, vrange, wrange;
66 typename NXmap<T>::Map_reference_index inx;
67 for ( inx = srchval.first(); !inx.last(); inx.next() )
68 if ( srchval[inx] != 0.0 ) {
69 cf = nxop.coord_frac( inx.coord().coord_map() );
70 urange.include( cf.u() );
71 vrange.include( cf.v() );
72 wrange.include( cf.w() );
73 }
74 cf = Coord_frac( urange.min(), vrange.min(), wrange.min() );
75 Coord_grid g0 = cf.coord_map( g ).coord_grid() - Coord_grid(1,1,1);
76 cf = Coord_frac( urange.max(), vrange.max(), wrange.max() );
77 Coord_grid g1 = cf.coord_map( g ).coord_grid() + Coord_grid(1,1,1);
78 Grid_range gm( g0, g1 );
79
80 // copy provided NXmap in new NXmap shadowing the xtal grid
81 NXmap<T> target( result.cell(), g, gm );
82 target = 0.0;
83 typename NXmap<T>::Map_reference_index in;
84 Coord_map cm;
85 for ( in = target.first(); !in.last(); in.next() ) {
86 cm = nxop.coord_map( in.coord_orth().coord_frac( result.cell() ) );
87 if ( Interp_linear::can_interp( srchval, cm ) ) {
88 Interp_linear::interp( srchval, cm, target[in] );
89 }
90 }
91
92 // now calculate the search function
93 const Xmap<T>& xmap = *xmp;
94 typename Xmap<T>::Map_reference_index pos;
95 typename Xmap<T>::Map_reference_coord i0, iu, iv, iw;
96 T s;
97 for ( pos = result.first(); !pos.last(); pos.next() ) {
98 g0 = pos.coord() + gm.min();
99 g1 = pos.coord() + gm.max();
100 i0 = Xmap_base::Map_reference_coord( xmap, g0 );
101 s = 0.0;
102 for ( iu = i0; iu.coord().u() <= g1.u(); iu.next_u() )
103 for ( iv = iu; iv.coord().v() <= g1.v(); iv.next_v() )
104 for ( iw = iv; iw.coord().w() <= g1.w(); iw.next_w() ) {
105 in.set_coord( iw.coord() - g0 );
106 s += xmap[iw] * target[in];
107 }
108 result[pos] = s;
109 }
110
111 return true;
112 }
113
114
115 // fast Convolution_search
116
init(const Xmap<T> & xmap)117 template<class T> void Convolution_search_fft<T>::init( const Xmap<T>& xmap )
118 {
119 vol = xmap.cell().volume();
120 rho1.init( xmap.grid_sampling() );
121
122 const Grid_sampling& g = xmap.grid_sampling();
123 T r;
124 typename Xmap<T>::Map_reference_coord i0( xmap, Coord_grid(0,0,0) );
125 typename Xmap<T>::Map_reference_coord iu, iv, iw;
126 for ( iu = i0; iu.coord().u() < g.nu(); iu.next_u() )
127 for ( iv = iu; iv.coord().v() < g.nv(); iv.next_v() )
128 for ( iw = iv; iw.coord().w() < g.nw(); iw.next_w() ) {
129 r = xmap[iw];
130 rho1.real_data( iw.coord() ) = r;
131 }
132 rho1.fft_x_to_h( vol );
133 }
134
operator ()(Xmap<T> & result,const NXmap<T> & srchval) const135 template<class T> bool Convolution_search_fft<T>::operator() ( Xmap<T>& result, const NXmap<T>& srchval ) const
136 {
137 NX_operator nxop( result, srchval, RTop_orth::identity() );
138 return (*this)( result, srchval, nxop );
139 }
140
operator ()(Xmap<T> & result,const NXmap<T> & srchval,const NX_operator & nxop) const141 template<class T> bool Convolution_search_fft<T>::operator() ( Xmap<T>& result, const NXmap<T>& srchval, const NX_operator& nxop ) const
142 {
143 const Grid_sampling& g = rho1.grid_real();
144 FFTmap_p1 map1( g );
145
146 // calculate extent of mask function in Xmap space
147 Coord_frac cf;
148 Range<ftype> urange, vrange, wrange;
149 typename NXmap<T>::Map_reference_index inx;
150 for ( inx = srchval.first(); !inx.last(); inx.next() )
151 if ( srchval[inx] != 0.0 ) {
152 cf = nxop.coord_frac( inx.coord().coord_map() );
153 urange.include( cf.u() );
154 vrange.include( cf.v() );
155 wrange.include( cf.w() );
156 }
157 cf = Coord_frac( urange.min(), vrange.min(), wrange.min() );
158 Coord_grid g0 = cf.coord_map( g ).coord_grid() - Coord_grid(1,1,1);
159 cf = Coord_frac( urange.max(), vrange.max(), wrange.max() );
160 Coord_grid g1 = cf.coord_map( g ).coord_grid() + Coord_grid(1,1,1);
161
162 // copy mask into grid, with offset
163 Coord_grid c, cu;
164 Coord_map cm;
165 T f;
166 for ( c.u() = g0.u(); c.u() <= g1.u(); c.u()++ )
167 for ( c.v() = g0.v(); c.v() <= g1.v(); c.v()++ )
168 for ( c.w() = g0.w(); c.w() <= g1.w(); c.w()++ ) {
169 cm = nxop.coord_map( c.coord_frac( g ) );
170 if ( Interp_linear::can_interp( srchval, cm ) ) {
171 Interp_linear::interp( srchval, cm, f );
172 cu = c.unit(g);
173 map1.real_data( cu ) = f;
174 }
175 }
176
177 // fft
178 map1.fft_x_to_h( vol );
179
180 // calculate
181 const Grid& gh = map1.grid_reci();
182 std::complex<ftype32> r1, f1;
183 for ( c.u() = 0; c.u() < gh.nu(); c.u()++ )
184 for ( c.v() = 0; c.v() < gh.nv(); c.v()++ )
185 for ( c.w() = 0; c.w() < gh.nw(); c.w()++ ) {
186 r1 = rho1.cplx_data( c );
187 f1 = map1.cplx_data( c );
188 map1.cplx_data( c ) = std::conj(f1)*r1;
189 }
190 // invert
191 map1.fft_h_to_x( map1.grid_real().size() / pow( vol, 2 ) );
192
193 // store
194 for ( typename Xmap<T>::Map_reference_index ix = result.first();
195 !ix.last(); ix.next() )
196 result[ix] = map1.real_data( ix.coord() );
197
198 return true;
199 }
200
201
202 // compile templates
203
204 template class Convolution_search_slow<ftype32>;
205 template class Convolution_search_slow<ftype64>;
206 template class Convolution_search_fft<ftype32>;
207 template class Convolution_search_fft<ftype64>;
208
209
210 } // namespace clipper
211