LOOS  v2.3.2
internal-water-filter.cpp
1 /*
2  Internal water filter library
3 */
4 
5 
6 /*
7  This file is part of LOOS.
8 
9  LOOS (Lightweight Object-Oriented Structure library)
10  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
11  Department of Biochemistry and Biophysics
12  School of Medicine & Dentistry, University of Rochester
13 
14  This package (LOOS) is free software: you can redistribute it and/or modify
15  it under the terms of the GNU General Public License as published by
16  the Free Software Foundation under version 3 of the License.
17 
18  This package is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  GNU General Public License for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with this program. If not, see <http://www.gnu.org/licenses/>.
25 */
26 
27 
28 #include <internal-water-filter.hpp>
29 
30 
31 using namespace std;
32 
33 namespace loos {
34  namespace DensityTools {
35 
36 
37  string WaterFilterBox::name(void) const {
38  stringstream s;
39  s << boost::format("WaterFilterBox(pad=%f)") % pad_;
40  return(s.str());
41  }
42 
43 
44  vector<int> WaterFilterBox::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
45  bdd_ = boundingBox(prot);
46  vector<int> result(solv.size());
47 
48  AtomicGroup::const_iterator ai;
49 
50  uint j = 0;
51  for (ai = solv.begin(); ai != solv.end(); ++ai) {
52  bool flag = true;
53  GCoord c = (*ai)->coords();
54  for (int i=0; i<3; ++i)
55  if (c[i] < bdd_[0][i] || c[i] > bdd_[1][i]) {
56  flag = false;
57  break;
58  }
59 
60  result[j++] = flag;
61  }
62  return(result);
63  }
64 
65 
66  double WaterFilterBox::volume(void) {
67  GCoord v = bdd_[1] - bdd_[0];
68  return(v[0] * v[1] * v[2]);
69  }
70 
71 
72  vector<GCoord> WaterFilterBox::boundingBox(const AtomicGroup& grp) {
73  vector<GCoord> bdd = grp.boundingBox();
74  bdd[0] = bdd[0] - pad_;
75  bdd[1] = bdd[1] + pad_;
76 
77  return(bdd);
78  }
79 
80 
81  // --------------------------------------------------------------------------------
82  string WaterFilterRadius::name(void) const {
83  stringstream s;
84  s << boost::format("WaterFilterRadius(radius=%f)") % radius_;
85  return(s.str());
86  }
87 
88 
89  vector<int> WaterFilterRadius::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
90  bdd_ = boundingBox(prot);
91  vector<int> result(solv.size());
92 
93  double r2 = radius_ * radius_;
94  for (uint j=0; j<prot.size(); ++j)
95  for (uint i=0; i<solv.size(); ++i)
96  if (prot[j]->coords().distance2(solv[i]->coords()) <= r2) {
97  result[i] = 1;
98  break;
99  }
100 
101  return(result);
102  }
103 
104 
105  double WaterFilterRadius::volume(void) {
106  GCoord v = bdd_[1] - bdd_[0];
107  return(v[0] * v[1] * v[2]);
108  }
109 
110 
111  vector<GCoord> WaterFilterRadius::boundingBox(const AtomicGroup& grp) {
112  vector<GCoord> bdd = grp.boundingBox();
113  bdd[0] = bdd[0] - radius_;
114  bdd[1] = bdd[1] + radius_;
115 
116  return(bdd);
117  }
118 
119 
120  // --------------------------------------------------------------------------------
121 
122 
123  string WaterFilterAxis::name(void) const {
124  stringstream s;
125  s << boost::format("WaterFilterAxis(radius=%f)") % sqrt(radius_);
126  return(s.str());
127  }
128 
129 
130  vector<int> WaterFilterAxis::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
131  bdd_ = boundingBox(prot);
132 
133  vector<int> result(solv.size());
134  AtomicGroup::const_iterator ai;
135  uint j = 0;
136 
137  for (ai = solv.begin(); ai != solv.end(); ++ai) {
138  GCoord a = (*ai)->coords();
139  if (a.z() < bdd_[0][2] || a.z() > bdd_[1][2]) {
140  result[j++] = 0;
141  continue;
142  }
143 
144  // Find the nearest point on the axis to the atom...
145  a -= orig_;
146  double k = (axis_ * a) / axis_.length2();
147  GCoord ah = orig_ + k * axis_;
148  GCoord v = (*ai)->coords() - ah;
149 
150  // Are we within the radius cutoff?
151  double d = v.length2();
152  if (d <= radius_)
153  result[j++] = true;
154  else
155  result[j++] = false;
156  }
157  return(result);
158  }
159 
160  double WaterFilterAxis::volume(void) {
161  return((bdd_[1][2] - bdd_[0][2]) * PI * radius_);
162  }
163 
164 
165  vector<GCoord> WaterFilterAxis::boundingBox(const AtomicGroup& grp) {
166 
167  // Set the principal axis...
168  orig_ = grp.centroid();
169  vector<GCoord> axes = grp.principalAxes();
170  axis_ = axes[0];
171  vector<GCoord> bdd = grp.boundingBox();
172 
173  // Calculate the extents of the box containing the principal axis cylinder...
174  double r = sqrt(radius_);
175  GCoord lbd = orig_ - axis_ - GCoord(r,r,0.0);
176  GCoord ubd = orig_ + axis_ + GCoord(r,r,0.0);
177 
178  // Set the z-bounds to the protein bounding box...
179  lbd[2] = bdd[0][2];
180  ubd[2] = bdd[1][2];
181 
182  // Replace...
183  bdd[0] = lbd;
184  bdd[1] = ubd;
185 
186  return(bdd);
187  }
188 
189  // --------------------------------------------------------------------------------
190 
191 
192  string WaterFilterCore::name(void) const {
193  stringstream s;
194  s << boost::format("WaterFilterCore(radius=%f)") % sqrt(radius_);
195  return(s.str());
196  }
197 
198 
199  GCoord WaterFilterCore::calculateAxis(const AtomicGroup& bundle)
200  {
201  if (!bundle.hasBonds())
202  throw(runtime_error("WaterFilterCore requires model connectivity (bonds)"));
203 
204  vector<AtomicGroup> segments = bundle.splitByMolecule();
205  GCoord axis(0,0,0);
206 
207  for (vector<AtomicGroup>::iterator i = segments.begin(); i != segments.end(); ++i) {
208  vector<GCoord> axes = (*i).principalAxes();
209  if (axes[0].z() < 0.0)
210  axis -= axes[0];
211  else
212  axis += axes[0];
213  }
214 
215  axis /= axis.length();
216  return(axis);
217  }
218 
219 
220 
221  vector<int> WaterFilterCore::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
222  bdd_ = boundingBox(prot);
223 
224  vector<int> result(solv.size());
225  AtomicGroup::const_iterator ai;
226  uint j = 0;
227 
228  for (ai = solv.begin(); ai != solv.end(); ++ai) {
229  GCoord a = (*ai)->coords();
230  if (a.z() < bdd_[0][2] || a.z() > bdd_[1][2]) {
231  result[j++] = 0;
232  continue;
233  }
234 
235  // Find the nearest point on the axis to the atom...
236  a -= orig_;
237  double k = (axis_ * a) / axis_.length2();
238  GCoord ah = orig_ + k * axis_;
239  GCoord v = (*ai)->coords() - ah;
240 
241  // Are we within the radius cutoff?
242  double d = v.length2();
243  if (d <= radius_)
244  result[j++] = true;
245  else
246  result[j++] = false;
247  }
248  return(result);
249  }
250 
251  // TODO: Fix!
252  double WaterFilterCore::volume(void) {
253  return((bdd_[1][2] - bdd_[0][2]) * PI * radius_);
254  }
255 
256 
257  vector<GCoord> WaterFilterCore::boundingBox(const AtomicGroup& grp) {
258 
259  // Set the principal axis...
260  orig_ = grp.centroid();
261  axis_ = calculateAxis(grp);
262  vector<GCoord> bdd = grp.boundingBox();
263 
264  // Calculate the extents of the box containing the principal axis cylinder...
265  double r = sqrt(radius_);
266  GCoord lbd = orig_ - axis_ - GCoord(r,r,0.0);
267  GCoord ubd = orig_ + axis_ + GCoord(r,r,0.0);
268 
269  // Set the z-bounds to the protein bounding box...
270  lbd[2] = bdd[0][2];
271  ubd[2] = bdd[1][2];
272 
273  // Replace...
274  bdd[0] = lbd;
275  bdd[1] = ubd;
276 
277  return(bdd);
278  }
279 
280 
281  // --------------------------------------------------------------------------------
282 
283 
284  string WaterFilterBlob::name(void) const {
285  stringstream s;
286  GCoord min = blob_.minCoord();
287  GCoord max = blob_.maxCoord();
288  DensityGridpoint dim = blob_.gridDims();
289  s << "WaterFilterBlob(" << dim << ":" << min << "x" << max << ")";
290  return(s.str());
291  }
292 
293 
294  double WaterFilterBlob::volume(void) {
295  if (vol >= 0.0)
296  return(vol);
297 
298  GCoord d = blob_.gridDelta();
299  double delta = d[0] * d[1] * d[2];
300  long n = blob_.maxGridIndex();
301  long c = 0;
302  for (long i=0; i<n; ++i)
303  if (blob_(i))
304  ++c;
305 
306  vol = c * delta;
307  return(vol);
308  }
309 
310  vector<int> WaterFilterBlob::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
311  vector<int> result(solv.size());
312  AtomicGroup::const_iterator ci;
313  uint j = 0;
314  for (ci = solv.begin(); ci != solv.end(); ++ci) {
315  GCoord c = (*ci)->coords();
316  DensityGridpoint probe = blob_.gridpoint(c);
317  if (blob_.inRange(probe))
318  result[j++] = (blob_(c) != 0);
319  else
320  result[j++] = 0;
321  }
322 
323  return(result);
324  }
325 
326 
327  // This ignores the protein bounding box...
328  vector<GCoord> WaterFilterBlob::boundingBox(const AtomicGroup& prot) {
329  if (bdd_set)
330  return(bdd_);
331 
332  DensityGridpoint dim = blob_.gridDims();
333  DensityGridpoint min = dim;
334  DensityGridpoint max(0,0,0);
335 
336  for (int k=0; k<dim[2]; ++k)
337  for (int j=0; j<dim[1]; ++j)
338  for (int i=0; i<dim[0]; ++i) {
339  DensityGridpoint probe(i,j,k);
340  if (blob_(probe) != 0)
341  for (int x=0; x<3; ++x) {
342  if (probe[x] < min[x])
343  min[x] = probe[x];
344  if (probe[x] > max[x])
345  max[x] = probe[x];
346  }
347  }
348 
349  vector<GCoord> bdd(2);
350  bdd[0] = blob_.gridToWorld(min);
351  bdd[1] = blob_.gridToWorld(max);
352 
353  return(bdd);
354  }
355 
356 
357 
358  // --------------------------------------------------------------------------------
359 
360 
361  string ZClippedWaterFilter::name(void) const {
362  stringstream s;
363 
364  s << boost::format("ZClippedWaterFilter(%s, %f, %f)") % WaterFilterDecorator::name() % zmin_ % zmax_;
365  return(s.str());
366  }
367 
368 
369  vector<int> ZClippedWaterFilter::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
370  vector<int> result = WaterFilterDecorator::filter(solv, prot);
371 
372  for (uint i=0; i<result.size(); ++i)
373  if (result[i]) {
374  GCoord c = solv[i]->coords();
375  if (c[2] < zmin_ || c[2] > zmax_)
376  result[i] = 0;
377  }
378 
379  return(result);
380  }
381 
382 
383  vector<GCoord> ZClippedWaterFilter::boundingBox(const AtomicGroup& grp) {
384  vector<GCoord> bdd = WaterFilterDecorator::boundingBox(grp);
385  bdd[0][2] = zmin_;
386  bdd[1][2] = zmax_;
387 
388  return(bdd);
389  }
390 
391  // --------------------------------------------------------------------------------
392 
393  string BulkedWaterFilter::name(void) const {
394  stringstream s;
395 
396  s << boost::format("BulkedWaterFilter(%s, %f, %f, %f)") % WaterFilterDecorator::name() % pad_ % zmin_ % zmax_;
397  return(s.str());
398  }
399 
400 
401  vector<int> BulkedWaterFilter::filter(const AtomicGroup& solv, const AtomicGroup& prot) {
402  vector<int> result = WaterFilterDecorator::filter(solv, prot);
403  vector<GCoord> bdd = boundingBox(prot);
404 
405  for (uint i=0; i<result.size(); ++i)
406  if (!result[i]) {
407  GCoord c = solv[i]->coords();
408  if ( ((c[0] >= bdd[0][0] && c[0] <= bdd[1][0]) &&
409  (c[1] >= bdd[0][1] && c[1] <= bdd[1][1]) &&
410  (c[2] >= bdd[0][2] && c[2] <= zmin_))
411  ||
412  ((c[0] >= bdd[0][0] && c[0] <= bdd[1][0]) &&
413  (c[1] >= bdd[0][1] && c[1] <= bdd[1][1]) &&
414  (c[2] <= bdd[1][2] && c[2] >= zmax_)) )
415  result[i] = true;
416  }
417 
418  return(result);
419  }
420 
421 
422  vector<GCoord> BulkedWaterFilter::boundingBox(const AtomicGroup& grp) {
423  vector<GCoord> bdd = grp.boundingBox();
424 
425  bdd[0] -= pad_;
426  bdd[1] += pad_;
427  return(bdd);
428  }
429 
430  };
431 
432 
433 };
std::vector< GCoord > boundingBox(void) const
Bounding box for the group...
bool hasBonds(void) const
Does any atom in the group have bond information???
STL namespace.
std::vector< AtomicGroup > splitByMolecule(void) const
Returns a vector of AtomicGroups split based on bond connectivity.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
double length2(void) const
Length of the Coord (as a vector) squared.
Definition: Coord.hpp:413
std::vector< GCoord > principalAxes(void) const
Compute the principal axes of a group.
Definition: AG_linalg.cpp:104
GCoord centroid(void) const
Centroid of atoms (ignores mass, operates in group coordinates)