WFMath  1.0.2
ball_funcs.h
1 // ball_funcs.h (n-dimensional ball implementation)
2 //
3 // The WorldForge Project
4 // Copyright (C) 2000, 2001 The WorldForge Project
5 //
6 // This program is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation; either version 2 of the License, or
9 // (at your option) any later version.
10 //
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with this program; if not, write to the Free Software
18 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
19 //
20 // For information about WorldForge and its authors, please contact
21 // the Worldforge Web Site at http://www.worldforge.org.
22 //
23 
24 // Author: Ron Steinke
25 
26 #ifndef WFMATH_BALL_FUNCS_H
27 #define WFMATH_BALL_FUNCS_H
28 
29 #include <wfmath/ball.h>
30 
31 #include <wfmath/axisbox.h>
32 #include <wfmath/miniball.h>
33 
34 #include <cassert>
35 
36 namespace WFMath {
37 
38 template<int dim>
39 AxisBox<dim> Ball<dim>::boundingBox() const
40 {
41  Point<dim> p_low, p_high;
42 
43  for(int i = 0; i < dim; ++i) {
44  p_low[i] = m_center[i] - m_radius;
45  p_high[i] = m_center[i] + m_radius;
46  }
47 
48  bool valid = m_center.isValid();
49 
50  p_low.setValid(valid);
51  p_high.setValid(valid);
52 
53  return AxisBox<dim>(p_low, p_high, true);
54 }
55 
56 template<int dim, template<class, class> class container>
57 Ball<dim> BoundingSphere(const container<Point<dim>, std::allocator<Point<dim> > >& c)
58 {
59  _miniball::Miniball<dim> m;
60  _miniball::Wrapped_array<dim> w;
61 
62  typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator i, end = c.end();
63  bool valid = true;
64 
65  for(i = c.begin(); i != end; ++i) {
66  valid = valid && i->isValid();
67  for(int j = 0; j < dim; ++j)
68  w[j] = (*i)[j];
69  m.check_in(w);
70  }
71 
72  m.build();
73 
74 #ifndef NDEBUG
75  double dummy;
76 #endif
77  assert("Check that bounding sphere is good to library accuracy" &&
78  m.accuracy(dummy) < numeric_constants<CoordType>::epsilon());
79 
80  w = m.center();
81  Point<dim> center;
82 
83  for(int j = 0; j < dim; ++j)
84  center[j] = w[j];
85 
86  center.setValid(valid);
87 
88  return Ball<dim>(center, std::sqrt(m.squared_radius()));
89 }
90 
91 template<int dim, template<class, class> class container>
92 Ball<dim> BoundingSphereSloppy(const container<Point<dim>, std::allocator<Point<dim> > >& c)
93 {
94  // This is based on the algorithm given by Jack Ritter
95  // in Volume 2, Number 4 of Ray Tracing News
96  // <http://www.acm.org/tog/resources/RTNews/html/rtnews7b.html>
97 
98  typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator i = c.begin(),
99  end = c.end();
100  if (i == end) {
101  return Ball<dim>();
102  }
103 
104  CoordType min[dim], max[dim];
105  typename container<Point<dim>, std::allocator<Point<dim> > >::const_iterator min_p[dim], max_p[dim];
106  bool valid = i->isValid();
107 
108  for(int j = 0; j < dim; ++j) {
109  min[j] = max[j] = (*i)[j];
110  min_p[j] = max_p[j] = i;
111  }
112 
113  while(++i != end) {
114  valid = valid && i->isValid();
115  for(int j = 0; j < dim; ++j) {
116  if(min[j] > (*i)[j]) {
117  min[j] = (*i)[j];
118  min_p[j] = i;
119  }
120  if(max[j] < (*i)[j]) {
121  max[j] = (*i)[j];
122  max_p[j] = i;
123  }
124  }
125  }
126 
127  CoordType span = -1;
128  int direction = -1;
129 
130  for(int j = 0; j < dim; ++j) {
131  CoordType new_span = max[j] - min[j];
132  if(new_span > span) {
133  span = new_span;
134  direction = j;
135  }
136  }
137 
138  assert("Have a direction of maximum size" && direction != -1);
139 
140  Point<dim> center = Midpoint(*(min_p[direction]), *(max_p[direction]));
141  CoordType dist = SloppyDistance(*(min_p[direction]), center);
142 
143  for(i = c.begin(); i != end; ++i) {
144  if(i == min_p[direction] || i == max_p[direction])
145  continue; // We already have these
146 
147  CoordType new_dist = SloppyDistance(*i, center);
148 
149  if(new_dist > dist) {
150  CoordType delta_dist = (new_dist - dist) / 2;
151  // Even though new_dist may be too large, delta_dist / new_dist
152  // always gives enough of a shift to include the new point.
153  center += (*i - center) * delta_dist / new_dist;
154  dist += delta_dist;
155  assert("Shifted ball contains new point" &&
156  SquaredDistance(*i, center) <= dist * dist);
157  }
158  }
159 
160  center.setValid(valid);
161 
162  return Ball<dim>(center, dist);
163 }
164 
165 // These two are here, instead of defined in the class, to
166 // avoid include order problems
167 
168 template<int dim>
169 inline Ball<dim> Point<dim>::boundingSphere() const
170 {
171  return Ball<dim>(*this, 0);
172 }
173 
174 template<int dim>
175 inline Ball<dim> Point<dim>::boundingSphereSloppy() const
176 {
177  return Ball<dim>(*this, 0);
178 }
179 
180 } // namespace WFMath
181 
182 #endif // WFMATH_BALL_FUNCS_H
Generic library namespace.
Definition: atlasconv.h:45
Ball< dim > BoundingSphereSloppy(const container< Point< dim >, std::allocator< Point< dim > > > &c)
get a bounding sphere for a set of points
Definition: ball_funcs.h:92
void setValid(bool valid=true)
make isValid() return true if you&#39;ve initialized the point by hand
Definition: point.h:129
Point< dim > Midpoint(const Point< dim > &p1, const Point< dim > &p2, CoordType dist=0.5)
Definition: point_funcs.h:240
float CoordType
Basic floating point type.
Definition: const.h:140
Ball< dim > BoundingSphere(const container< Point< dim >, std::allocator< Point< dim > > > &c)
get the minimal bounding sphere for a set of points
Definition: ball_funcs.h:57
A dim dimensional point.
Definition: const.h:50
A dim dimensional ball.
Definition: ball.h:34