LOOS  v2.3.2
DensityGrid.hpp
1 /*
2  Simple Density Grid Class for LOOS
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 
29 
30 #if !defined(LOOS_DENSITYGRID_HPP)
31 #define LOOS_DENSITYGRID_HPP
32 
33 #include <cmath>
34 #include <stdexcept>
35 #include <utility>
36 
37 #include <stdexcept>
38 #include <boost/iterator/iterator_facade.hpp>
39 
40 #include <loos.hpp>
41 #include <Coord.hpp>
42 
43 #include <SimpleMeta.hpp>
44 
45 namespace loos {
46 
48  namespace DensityTools {
49 
51 
52 
53  template<class T> class DensityGrid;
54 
56 
61  template<class T>
63  public:
64  DensityGridRow(const long i, DensityGrid<T>& g) : idx(i), grid(g) { }
65 
66  T& operator[](const int i) {
67  assert(i >= 0 && i < grid.dims[0]);
68  assert(grid.ptr != 0);
69  return(grid.ptr[idx + i]);
70  }
71 
72  private:
73  long idx;
74  DensityGrid<T>& grid;
75  };
76 
77 
79 
84  template<class T>
86  public:
87  DensityGridPlane(const long i, DensityGrid<T>& g) : idx(i), grid(g) { }
88 
89  DensityGridRow<T> operator[](const int j) {
90  assert(j >= 0 && j < grid.dims[1]);
91  assert(grid.ptr != 0);
92  return(DensityGridRow<T>(idx + j*grid.dims[0], grid));
93  }
94 
95  private:
96  long idx;
97  DensityGrid<T>& grid;
98  };
99 
100 
102 
125  // T is the fundamental type (i.e. the type stored in the grid)
126  // R is the dereference type (this allows us to use a template to
127  // handle both const and non-const iterators).
128  //
129  // Internally, the iterator position is stored as an index into the
130  // block of data managed by the DensityGrid, and is a long, hence the
131  // difference type below...
132  template<typename T, typename R>
133  class DensityGridIterator : public boost::iterator_facade<
134  DensityGridIterator<T, R>,
135  T,
136  boost::random_access_traversal_tag,
137  R&,
138  long
139  >
140  {
141  public:
142  DensityGridIterator() : src(0), offset(0) { }
143  explicit DensityGridIterator(const DensityGrid<T>& g, long l) : src(&g), offset(l) { }
144 
145  // Converts from any dereference type...
146  template<typename S> DensityGridIterator(const DensityGridIterator<T, S>& o) : src(o.src), offset(o.offset) { }
147 
148  loos::GCoord world() const { return(src->gridToWorld(src->indexToGrid(offset))); }
149  loos::GCoord coords() const { return(world()); }
150  DensityGridpoint grid() const { return(src->indexToGrid(offset)); }
151 
152  private:
153 
154  friend class boost::iterator_core_access;
155 
156  // This lets any other DensityGridIterator access our internal bits...
157  template<typename,typename> friend class DensityGridIterator;
158 
159  void increment() { ++offset; }
160  void decrement() { --offset; }
161  void advance(const long n) { offset += n; }
162 
163  template<typename S>
164  bool equal(const DensityGridIterator<T, S>& other) const {
165  return(src == other.src && offset == other.offset);
166  }
167 
168  R& dereference() const {
169  if (offset < 0 || offset >= src->dimabc)
170  throw(std::range_error("Index out of bounds"));
171  return(src->ptr[offset]);
172  }
173 
174  long distance_to(const DensityGridIterator<T, R>& other) const {
175  if (src != other.src)
176  throw(std::logic_error("Grid mismatch"));
177  return(other.offset - offset);
178  }
179 
180  const DensityGrid<T>* src;
181  long offset;
182  };
183 
184 
186 
223  template<class T>
224  class DensityGrid {
225  public:
226 
227  typedef T value_type;
228  typedef DensityGridIterator<T, T> iterator;
229  typedef DensityGridIterator<T, const T> const_iterator;
230 
231  friend class DensityGridRow<T>;
232  friend class DensityGridPlane<T>;
233 
234  friend class DensityGridIterator<T, T>;
235  friend class DensityGridIterator<T, const T>;
236 
237  typedef std::pair<int, int> Range;
238 
240  DensityGrid() : ptr(0), _gridmin(loos::GCoord(0,0,0)), _gridmax(loos::GCoord(0,0,0)),
241  dims(DensityGridpoint(0,0,0)) { }
242 
244 
248  DensityGrid(const loos::GCoord& gmin, const loos::GCoord& gmax, const DensityGridpoint& griddims) :
249  ptr(0), _gridmin(gmin), _gridmax(gmax), dims(griddims) { init(); }
250 
253  DensityGrid(const loos::GCoord& gmin, const loos::GCoord& gmax, const int dim) :
254  ptr(0), _gridmin(gmin), _gridmax(gmax), dims(dim, dim, dim) { init(); }
255 
257  DensityGrid(const DensityGrid<T>& g) : ptr(0), _gridmin(g._gridmin), _gridmax(g._gridmax),
258  dims(g.dims) {
259  init();
260  memcpy(ptr, g.ptr, dimabc * sizeof(T));
261  meta_ = g.meta_;
262  }
263 
264  void resize(const loos::GCoord& gmin, const loos::GCoord& gmax, const DensityGridpoint& griddims) {
265  if (ptr != 0)
266  delete[] ptr;
267  _gridmin = gmin;
268  _gridmax = gmax;
269  dims = griddims;
270  init();
271  }
272 
275  if (this == &g)
276  return(*this);
277 
278  if (ptr) {
279  delete[] ptr;
280  ptr = 0;
281  }
282  _gridmin = g._gridmin;
283  _gridmax = g._gridmax;
284  dims = g.dims;
285  init();
286  if (dimabc != 0)
287  memcpy(ptr, g.ptr, dimabc * sizeof(T));
288 
289  meta_ = g.meta_;
290 
291  return(*this);
292  }
293 
294 
295  ~DensityGrid() { delete[] ptr; }
296 
297 
298 
299  DensityGrid<T> subset(const Range& c, const Range& b, const Range& a) {
300  DensityGridpoint dim;
301 
302  dim.x(a.second - a.first + 1);
303  dim.y(b.second - b.first + 1);
304  dim.z(c.second - c.first + 1);
305 
306  loos::GCoord bottom = gridToWorld(DensityGridpoint(a.first, b.first, c.first));
307  loos::GCoord top = gridToWorld(DensityGridpoint(a.second, b.second, c.second));
308 
309  DensityGrid<T> sub(bottom, top, dim);
310  for (int k=0; k<dim.z(); ++k)
311  for (int j=0; j<dim.y(); ++j)
312  for (int i=0; i<dim.x(); ++i)
313  sub(k, j, i) = operator()(k+c.first, j+b.first, i+a.first);
314 
315  return(sub);
316  }
317 
319  void zero(void) { for (long i=0; i<dimabc; i++) ptr[i] = 0; }
320 
321 
324  long gridToIndex(const DensityGridpoint v) const {
325  return( (v.z() * dims[1] + v.y()) * dims[0] + v.x() );
326  }
327 
328 
330  DensityGridpoint gridpoint(const loos::GCoord& x) const {
331  DensityGridpoint v;
332 
333  for (int i=0; i<3; i++) {
334  long k = static_cast<long>(floor( (x[i] - _gridmin[i]) * delta[i] + 0.5 ));
335  v[i] = k;
336  }
337 
338  return(v);
339  }
340 
342  DensityGridpoint gridpoint(const double z, const double y, const double x) const {
343  loos::GCoord c(x,y,z);
344  return(gridpoint(c));
345  }
346 
347 
349  bool inRange(const DensityGridpoint& g) const {
350  for (int i=0; i<3; i++)
351  if (g[i] < 0 || g[i] >= dims[i])
352  return(false);
353 
354  return(true);
355  }
356 
357  bool inRange(const int k, const int j, const int i) const {
358  return(inRange(DensityGridpoint(i, j, k)));
359  }
360 
362 
363  // Bypassing the Coord<>::operator[] saves us a considerable
364  // amount of time...hence the Coord<>::x(), etc...
365 
366  T& operator()(const int k, const int j, const int i) {
367  long x = ((k * dims.y()) + j) * dims.x() + i;
368  assert(x < dimabc);
369  return(ptr[x]);
370  }
371 
373  T& operator()(const DensityGridpoint& v) {
374  for (int i=0; i<3; i++)
375  assert(v[i] >= 0 && v[i] < dims[i]);
376  return(ptr[ (v[2] * dims.y() + v[1]) * dims.x() + v[0] ]);
377  }
378 
381  T& operator()(const long i) {
382  assert(i >= 0 && i < dimabc);
383  return(ptr[i]);
384  }
385 
387  T& operator()(const loos::GCoord& x) {
388  return(operator()(gridpoint(x)));
389  }
390 
391 
392  // Const versions...
393 
394 
395  const T& operator()(const int k, const int j, const int i) const {
396  long x = ((k * dims.y()) + j) * dims.x() + i;
397  assert(x < dimabc);
398  return(ptr[x]);
399  }
400 
401  const T& operator()(const DensityGridpoint& v) const {
402  for (int i=0; i<3; i++)
403  assert(v[i] >= 0 && v[i] < dims[i]);
404  return(ptr[ (v[2] * dims.y() + v[1]) * dims.x() + v[0] ]);
405  }
406 
407  const T& operator()(const long i) const {
408  assert(i >= 0 && i < dimabc);
409  return(ptr[i]);
410  }
411 
412 
413  const T& operator()(const loos::GCoord& x) const {
414  return(operator()(gridpoint(x)));
415  }
416 
417 
418  // Slicin' und dicin'...
419 
422  assert(k >= 0 && k < dims[2]);
423  return(DensityGridPlane<T>(k * dimab, *this));
424  }
425 
426 
428  loos::GCoord gridToWorld(const DensityGridpoint& v) const {
429  loos::GCoord c;
430 
431  for (int i=0; i<3; i++)
432  c[i] = static_cast<loos::greal>(v[i]) / delta[i] + _gridmin[i];
433 
434  return(c);
435  }
436 
438  DensityGridpoint indexToGrid(const long idx) const {
439  int c = idx / dimab;
440  int r = idx % dimab;
441  int b = r / dims[0];
442  int a = r % dims[0];
443 
444  return(DensityGridpoint(a, b, c));
445  }
446 
447 
449  double gridDist2(const DensityGridpoint& u, const DensityGridpoint& v) {
450  loos::GCoord x = gridToWorld(u);
451  loos::GCoord y = gridToWorld(v);
452 
453  return(x.distance2(y));
454  }
455 
457  double gridDist(const DensityGridpoint& u, const DensityGridpoint& v) {
458 
459  return(sqrt(gridDist2(u, v)));
460  }
461 
464  std::vector<DensityGridpoint> withinBoxRadius(const double r, const loos::GCoord& u, const int pad = 0) const {
465 
466  DensityGridpoint a = gridpoint(u - r);
467  DensityGridpoint b = gridpoint(u + r);
468 
469  std::vector<DensityGridpoint> res;
470  for (int k=a[2]-pad; k <= b[2]+pad; k++) {
471  if (k < 0 || k >= dims[2])
472  continue;
473  for (int j=a[1]-pad; j <= b[1]+pad; j++) {
474  if (j < 0 || j >= dims[1])
475  continue;
476  for (int i=a[0]-pad;i <= b[0]+pad; i++) {
477  if (i < 0 || i >= dims[0])
478  continue;
479  res.push_back(DensityGridpoint(i, j, k));
480  }
481  }
482  }
483 
484  return(res);
485  }
486 
487 
490 
497  std::vector<DensityGridpoint> withinRadius(const double r, const loos::GCoord& u) const {
498  std::vector<DensityGridpoint> pts = withinBoxRadius(r, u);
499  std::vector<DensityGridpoint> res;
500  std::vector<DensityGridpoint>::iterator i;
501 
502  double r2 = r*r;
503 
504  for (i=pts.begin(); i != pts.end(); i++) {
505  loos::GCoord v = gridToWorld(*i);
506  if (u.distance2(v) <= r2)
507  res.push_back(*i);
508  }
509 
510  return(res);
511  }
512 
513 
516  /*
517  * Here, a functor \a f is called for all grid points that lie on
518  * or within a sphere of radius \a r about the real-space coords
519  * \a u. The functor take two parameters: the value at the gid
520  * point and the distance from that grid point to the center of
521  * the sphere (in real-space units)
522  */
523  template<typename Func>
524  void applyWithinRadius(const double r, const loos::GCoord& u, const Func& f) {
525  DensityGridpoint a = gridpoint(u - r);
526  DensityGridpoint b = gridpoint(u + r);
527  double r2 = r * r;
528 
529  std::vector<DensityGridpoint> res;
530  for (int k=a[2]; k <= b[2]; k++) {
531  if (k < 0 || k >= dims[2])
532  continue;
533  for (int j=a[1]; j <= b[1]; j++) {
534  if (j < 0 || j >= dims[1])
535  continue;
536  for (int i=a[0];i <= b[0]; i++) {
537  if (i < 0 || i >= dims[0])
538  continue;
539 
540  DensityGridpoint point(i, j, k);
541  loos::GCoord v = gridToWorld(point);
542  double d = u.distance2(v);
543  if (d <= r2)
544  f(operator()(v), d);
545  }
546  }
547  }
548 
549 
550  }
551 
552 
553  void scale(const T val) {
554  for (long i = 0; i < dimabc; ++i)
555  ptr[i] *= val;
556  }
557 
558  void clear(const T val = 0) {
559  for (long i = 0; i < dimabc; ++i)
560  ptr[i] = val;
561  }
562 
563 
564  DensityGridpoint gridDims(void) const { return(dims); }
565  loos::GCoord minCoord(void) const { return(_gridmin); }
566  loos::GCoord maxCoord(void) const { return(_gridmax); }
567  loos::GCoord gridDelta(void) const { return(delta); }
568 
569  long maxGridIndex(void) const { return(dimabc); }
570  long size() const { return(dimabc); }
571 
573 
579  friend std::ostream& operator<<(std::ostream& os, const DensityGrid<T>& grid) {
580  os << "# DensityGrid-1.1\n";
581  os << grid.meta_;
582  os << grid.dims << std::endl;
583  os << grid._gridmin << std::endl;
584  os << grid._gridmax << std::endl;
585  // os << boost::format("(%.8f,%.8f,%.8f)\n") % grid._gridmin[0] % grid._gridmin[1] % grid._gridmin[2];
586  // os << boost::format("(%.8f,%.8f,%.8f)\n") % grid._gridmax[0] % grid._gridmax[1] % grid._gridmax[2];
587 
588  return(os.write(reinterpret_cast<char*>(grid.ptr), sizeof(T) * grid.dimabc));
589  }
590 
592 
596  friend std::istream& operator>>(std::istream& is, DensityGrid<T>& grid) {
597  std::string s;
598 
599  std::getline(is, s);
600  if (s != "# DensityGrid-1.1")
601  throw(std::runtime_error("Bad input format for DensityGrid - " + s));
602 
603  is >> grid.meta_;
604 
605  if (grid.ptr)
606  delete[] grid.ptr;
607 
608  is >> grid.dims;
609  is >> grid._gridmin;
610  is >> grid._gridmax;
611  if (is.get() != '\n') // Pull trailing newline off of input...
612  throw(std::runtime_error("Grid parse error in header"));
613 
614  grid.init();
615  is.read(reinterpret_cast<char*>(grid.ptr), sizeof(T) * grid.dimabc);
616  if (is.fail() || is.eof())
617  throw(std::runtime_error("Grid read error"));
618 
619  return(is);
620  }
621 
622  iterator begin() { return(iterator(*this, 0)); }
623  iterator end() { return(iterator(*this, dimabc)); }
624 
625  const_iterator begin() const { return(const_iterator(*this, 0)); }
626  const_iterator end() const { return(const_iterator(*this, dimabc)); }
627 
628  bool empty() const { return(dimabc == 0); }
629 
630  void setMetadata(const std::string& s) { meta_.set(s); }
631  void addMetadata(const std::string& s) { meta_.add(s); }
632 
633  SimpleMeta metadata() const { return(meta_); }
634  void metadata(const SimpleMeta& m) { meta_ = m; }
635 
636 
637  private:
638  void init(void) {
639  dimab = dims[0]*dims[1];
640  dimabc = dimab * dims[2];
641 
642  for (int i=0; i<3; i++)
643  delta[i] = (dims[i] - 1)/ (_gridmax[i] - _gridmin[i]);
644 
645  if (dimabc != 0) {
646  ptr = new T[dimabc];
647  zero();
648  } else
649  ptr = 0;
650  }
651 
652  private:
653  T* ptr;
654  loos::GCoord _gridmin, _gridmax, delta;
655  DensityGridpoint dims;
656  long dimabc, dimab;
657 
658  SimpleMeta meta_;
659  };
660 
661  };
662 
663 };
664 
665 #endif
DensityGridpoint gridpoint(const loos::GCoord &x) const
Converts a real-space coordinate into grid coords.
Random access iterator using the BOOST facade.
void applyWithinRadius(const double r, const loos::GCoord &u, const Func &f)
bool inRange(const DensityGridpoint &g) const
Checks to make sure the gridpoint lies within the grid boundaries.
std::vector< DensityGridpoint > withinRadius(const double r, const loos::GCoord &u) const
DensityGrid(const loos::GCoord &gmin, const loos::GCoord &gmax, const DensityGridpoint &griddims)
Create a grid with explicit location in realspace and dimensions.
DensityGridPlane< T > operator[](const int k)
Returns the kth plane from the grid.
T & operator()(const loos::GCoord &x)
Converts x into grid coords, then accesses that element.
double distance2(const Coord< T > &o) const
Distance squared between two coordinates.
Definition: Coord.hpp:429
A simple 3D grid class of arbitrary types.
Definition: DensityGrid.hpp:53
T & operator()(const DensityGridpoint &v)
Access the grid element indexed by the DensityGridPoint.
friend std::istream & operator>>(std::istream &is, DensityGrid< T > &grid)
Read in a grid.
double gridDist(const DensityGridpoint &u, const DensityGridpoint &v)
Linear distance (in real-space) between two grid coords.
DensityGrid(const DensityGrid< T > &g)
Copies an existing grid.
DensityGrid(const loos::GCoord &gmin, const loos::GCoord &gmax, const int dim)
std::vector< DensityGridpoint > withinBoxRadius(const double r, const loos::GCoord &u, const int pad=0) const
long gridToIndex(const DensityGridpoint v) const
void add(const std::string &s)
Append metadata.
Definition: SimpleMeta.hpp:82
loos::GCoord gridToWorld(const DensityGridpoint &v) const
Converts grid coords to real-space (world) coords.
double gridDist2(const DensityGridpoint &u, const DensityGridpoint &v)
Squared-distance (in real-space) between two grid coords.
const DensityGrid< T > & operator=(const DensityGrid< T > &g)
This is a "deep" copy of grid.
DensityGridpoint gridpoint(const double z, const double y, const double x) const
Converts a real-space coordinate into grid coords.
Namespace for most things not already encapsulated within a class.
Encapsulates an i,j-plane from an DensityGrid.
Definition: DensityGrid.hpp:85
Encapsulates a j-row from an DensityGrid.
Definition: DensityGrid.hpp:62
void set(const std::string &s)
Set metadata to string (deletes existing metadata)
Definition: SimpleMeta.hpp:79
void zero(void)
Zero out all elements.
T & operator()(const int k, const int j, const int i)
Access the grid element indexed by k, j, i.
DensityGridpoint indexToGrid(const long idx) const
Calculates the grid coords from a linear index.