LOOS  v2.3.2
GridUtils.hpp
1 /*
2  Utility functions/classes for LOOS grids
3 */
4 
5 
6 /*
7  This file is part of LOOS.
8 
9  LOOS (Lightweight Object-Oriented Structure library)
10  Copyright (c) 2009, 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 
29 
30 #if !defined(LOOS_GRID_UTILS_HPP)
31 #define LOOS_GRID_UTILS_HPP
32 
33 #include <DensityGrid.hpp>
34 
35 namespace loos {
36 
37  namespace DensityTools {
38 
39 
40  // Functors for grid operations
41 
43  template<typename T>
44  class Threshold {
45  public:
46  Threshold(const T& t) : thresh(t) { }
47  bool operator()(const T& t) const { return(t >= thresh); }
48  private:
49  T thresh;
50  };
51 
53  template<typename T>
55  public:
56  ThresholdRange(const T& l, const T& h) : lo(l), hi(h) { }
57  bool operator()(const T& t) const { return( t >= lo && t <= hi ); }
58  private:
59  T lo, hi;
60  };
61 
63  template<typename T>
65  public:
66  bool operator()(const T& t) const { return(t > 0.0); }
67  };
68 
69 
71 
80  template<typename T, class Functor>
81  std::vector<DensityGridpoint> floodFill(const DensityGridpoint seed, const DensityGrid<T>& data_grid,
82  const int id, DensityGrid<int>& blob_grid, const Functor& op)
83  {
84  std::vector<DensityGridpoint> stack;
85  stack.push_back(seed);
86  std::vector<DensityGridpoint> list;
87  list.push_back(seed);
88  blob_grid(seed) = id;
89 
90  while(!stack.empty()) {
91  DensityGridpoint point = stack.back();
92  stack.pop_back();
93  for (int k=-1; k<=1; ++k)
94  for (int j=-1; j<=1; ++j)
95  for (int i=-1; i<=1; ++i) {
96  if (i == j && j == k)
97  continue;
98  DensityGridpoint probe = point + DensityGridpoint(i, j, k);
99  if (!data_grid.inRange(probe))
100  continue;
101  if (blob_grid(probe) == 0 && op(data_grid(probe))) {
102  blob_grid(probe) = id;
103  stack.push_back(probe);
104  list.push_back(probe);
105  }
106  }
107 
108  }
109 
110  return(list);
111  }
112 
114  template<typename T, class Functor>
115  int floodFill(const DensityGridpoint seed, const DensityGrid<T>& data_grid, const Functor& op) {
116  DensityGrid<int> blob_grid(data_grid.minCoord(), data_grid.maxCoord(), data_grid.gridDims());
117 
118  std::vector<DensityGridpoint> result = floodFill(seed, data_grid, 1, blob_grid, op);
119  return(result.size());
120  }
121 
122 
124 
135  template<typename T, class Functor>
136  std::vector<loos::GCoord> findPeaks(const DensityGrid<T>& grid, DensityGrid<int>& blobs, const Functor& op) {
137  std::vector<loos::GCoord> peaks;
138 
139  DensityGridpoint dims = grid.gridDims();
140 
141  int id = 0;
142  for (int k=0; k<dims.z(); ++k)
143  for (int j=0; j<dims.y(); ++j)
144  for (int i=0; i<dims.x(); ++i) {
145  DensityGridpoint p (i, j, k);
146  if (!blobs(p) && op(grid(p))) {
147  std::vector<DensityGridpoint> points = floodFill(p, grid, ++id, blobs, op);
148  if (!points.empty()) {
149  loos::GCoord center(0,0,0);
150  double mass = 0.0;
151  for (std::vector<DensityGridpoint>::iterator i = points.begin(); i != points.end(); ++i) {
152  double m = grid(*i);
153  center += m * grid.gridToWorld(*i);
154  mass += m;
155  }
156  center /= mass;
157  peaks.push_back(center);
158  }
159  }
160  }
161 
162  return(peaks);
163  }
164 
165 
167  template<typename T, class Functor>
168  std::vector<loos::GCoord> findPeaks(const DensityGrid<T>& grid, const Functor& op) {
169 
170  DensityGridpoint dims = grid.gridDims();
171  DensityGrid<int> blobs(grid.minCoord(), grid.maxCoord(), dims);
172  return(findPeaks(grid, blobs, op));
173  }
174 
175 
177  template<class T, class Functor>
178  loos::AtomicGroup gridToAtomicGroup(const DensityGrid<T>& grid, const Functor& op) {
179  loos::AtomicGroup group;
180  DensityGridpoint dims = grid.gridDims();
181 
182  int id = 0;
183  for (int k=0; k<dims.z(); ++k)
184  for (int j=0; j<dims.y(); ++j)
185  for (int i=0; i<dims.x(); ++i) {
186  DensityGridpoint p(i, j, k);
187  if (op(grid(p))) {
188  loos::pAtom atom(new loos::Atom(++id, "UNK", grid.gridToWorld(p)));
189  atom->resid(id);
190  atom->resname("GRD");
191  atom->mass(grid(p));
192  group.append(atom);
193  }
194  }
195  return(group);
196  }
197 
198 
199 
201  template<class T>
203  DensityGrid<T> tmp(grid);
204  DensityGridpoint gdim = grid.gridDims();
205  DensityGridpoint kdim = kernel.gridDims();
206 
207  int kkc = kdim.z() / 2;
208  int kjc = kdim.y() / 2;
209  int kic = kdim.x() / 2;
210 
211  for (int k=0; k<gdim.z(); ++k)
212  for (int j=0; j<gdim.y(); ++j)
213  for (int i=0; i<gdim.x(); ++i) {
214  T sum = 0;
215 
216  for (int kk=0; kk<kdim.z(); ++kk)
217  for (int jj=0; jj<kdim.y(); ++jj)
218  for (int ii=0; ii<kdim.x(); ++ii) {
219  int gk = k + (kk - kkc);
220  if (gk < 0 || gk > gdim.z())
221  continue;
222 
223  int gj = j + (jj - kjc);
224  if (gj < 0 || gj > gdim.y())
225  continue;
226 
227  int gi = i + (ii - kic);
228  if (gi < 0 || gi > gdim.x())
229  continue;
230 
231  sum += grid(k,j,i) * kernel(kk, jj, ii);
232  }
233 
234  tmp(k, j, i) = sum;
235  }
236 
237  grid = tmp;
238  }
239 
240 
242  template<class T>
243  void gridConvolve(DensityGrid<T>& grid, std::vector<T>& kernel) {
244  DensityGridpoint gdim = grid.gridDims();
245  int kn = kernel.size();
246  int kc = kn / 2;
247 
248  DensityGrid<T> tmp(grid.minCoord(), grid.maxCoord(), grid.gridDims());
249  tmp.metadata(grid.metadata());
250 
251  // First convolve along the k axis...
252  for (int j=0; j<gdim.y(); ++j)
253  for (int i=0; i<gdim.x(); ++i)
254  for (int k = 0; k<gdim.z(); ++k) {
255  T sum = 0;
256  for (int ii = 0; ii<kn; ++ii) {
257  int idx = k + ii - kc;
258  if (idx < 0 || idx >= gdim.z())
259  continue;
260  sum += grid(idx, j, i) * kernel[ii];
261  }
262  tmp(k, j, i) = sum;
263  }
264 
265  // Convolve along the j axis
266  DensityGrid<T> tmp2(grid.minCoord(), grid.maxCoord(), grid.gridDims());
267  tmp2.metadata(grid.metadata());
268  for (int k=0; k<gdim.z(); ++k)
269  for (int i=0; i<gdim.x(); ++i)
270  for (int j=0; j<gdim.y(); ++j) {
271  T sum = 0;
272  for (int ii=0; ii<kn; ++ii) {
273  int idx = j + ii - kc;
274  if (idx < 0 || idx >= gdim.y())
275  continue;
276  sum += tmp(k, idx, i) * kernel[ii];
277  }
278  tmp2(k, j, i) = sum;
279  }
280 
281  // Finally, convolve along the i axis
282  for (int k=0; k<gdim.z(); ++k)
283  for (int j=0; j<gdim.y(); ++j)
284  for (int i=0; i<gdim.x(); ++i) {
285  T sum = 0;
286  for (int ii=0; ii<kn; ++ii) {
287  int idx = i + ii - kc;
288  if (idx < 0 || idx >= gdim.x())
289  continue;
290  sum += tmp2(k, j, idx) * kernel[ii];
291  }
292  tmp(k, j, i) = sum;
293  }
294 
295  grid = tmp;
296  }
297 
299  std::vector<double> gaussian1d(const int, const double);
300 
301  };
302 
303 };
304 
305 
306 
307 #endif
bool inRange(const DensityGridpoint &g) const
Checks to make sure the gridpoint lies within the grid boundaries.
AtomicGroup & append(pAtom pa)
Append the atom onto the group.
loos::AtomicGroup gridToAtomicGroup(const DensityGrid< T > &grid, const Functor &op)
Converts grid points (determined by functor) into an AtomicGroup of pseudo-atoms. ...
Definition: GridUtils.hpp:178
A simple 3D grid class of arbitrary types.
Definition: DensityGrid.hpp:53
Basic Atom class for handling atom properties.
Definition: Atom.hpp:50
std::vector< double > gaussian1d(const int w, const double sigma)
Construct a 1D gaussian.
Definition: GridUtils.cpp:38
Functor that is true if value is greater than or equal to threshold.
Definition: GridUtils.hpp:44
Functor that is true if value is between low and hi, inclusive.
Definition: GridUtils.hpp:54
void gridConvolve(DensityGrid< T > &grid, DensityGrid< T > &kernel)
Convolve a grid with another grid (kernel)
Definition: GridUtils.hpp:202
std::vector< loos::GCoord > findPeaks(const DensityGrid< T > &grid, DensityGrid< int > &blobs, const Functor &op)
Find peaks in a grid given the criteria defined by the passed functor.
Definition: GridUtils.hpp:136
loos::GCoord gridToWorld(const DensityGridpoint &v) const
Converts grid coords to real-space (world) coords.
Functor that is true for any non-zero density.
Definition: GridUtils.hpp:64
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
std::vector< DensityGridpoint > floodFill(const DensityGridpoint seed, const DensityGrid< T > &data_grid, const int id, DensityGrid< int > &blob_grid, const Functor &op)
Flood-fill a grid.
Definition: GridUtils.hpp:81