source: Revenant/geoworld/include/sphpack.h@ 7ef8ec4

port/mars-tycoon
Last change on this file since 7ef8ec4 was 80a6a52, checked in by Jonathan Neufeld <support@…>, 3 years ago

Get to a compile state for terrain procedural generation
  • Property mode set to 100644
File size: 18.4 KB
Line 
1#ifndef __SPHEREPACK_H__
2#define __SPHEREPACK_H__
3
4#include <list>
5#include <assert.h>
6
7#include "mars_lookup.h"
8#include "mars_linalg.h"
9#include "mars_util.h"
10#include "mars_calc.h"
11#include "mars_ptr.h"
12
13namespace geoworld
14{
15 template <typename T>
16 class RandomSpherePack;
17
18 template <typename T>
19 class RandomSpherePack
20 {
21 public:
22 class Sphere;
23 typedef typename mars::TreeTag3D< T >::VECTOR VECTOR;
24 typedef typename mars::TreeTag3D< T >::SCALAR SCALAR;
25 typedef std::list <VECTOR> PointList;
26 typedef std::list< mars::ptr< Sphere > > SphereList;
27 typedef typename SphereList::iterator SphereListIterator;
28 typedef typename SphereList::const_iterator SphereListConstIterator;
29
30 struct IteratorTraits
31 {
32 typedef RandomSpherePack SpherePack;
33 using SphereType = Sphere;
34 typedef typename SphereList::iterator SphereListIterator;
35 };
36 struct ConstIteratorTraits
37 {
38 typedef const RandomSpherePack SpherePack;
39 using SphereType = const Sphere;
40 typedef typename SphereList::const_iterator SphereListIterator;
41 };
42
43 class Sphere
44 {
45 private:
46 SphereList _subnodes;
47
48 public:
49 const RandomSpherePack< T > * owner;
50
51 VECTOR p;
52 SCALAR r;
53
54 inline Sphere (const RandomSpherePack< T > * pOwner)
55 : owner(pOwner), r(0) {}
56
57 inline Sphere (const Sphere & copy)
58 : owner(copy.owner), p(copy.p), r(copy.r), _subnodes(copy._subnodes) {}
59
60 inline Sphere (const RandomSpherePack< T > * pOwner, const VECTOR & p, const T & r)
61 : owner(pOwner), p(p), r(r) {}
62
63 template< typename Tag >
64 inline operator mars::RadialRegion< Tag > () const { return mars::RadialRegion< Tag > (p, r); }
65
66 inline operator mars::CylindricalRegion< T > () const { return mars::CylindricalRegion< T > (p, r); }
67
68 inline static float computeCull (const float azimuth, const float zenith)
69 {
70 using namespace mars;
71
72 static const float PIxPI = SQ(static_cast<float> (PI));
73 // Calculates dot product of the displacement vector with the cull plane normal
74 // Spherical analog of Pythagoras: cos(c) = cos(a)cos(b)
75 return (PIxPI - SQ(acos(cos(azimuth)*cos(zenith)))) / PIxPI;
76 }
77
78 inline SphereListConstIterator begin() const { return _subnodes.begin(); }
79 inline SphereListConstIterator end() const { return _subnodes.end(); }
80 inline SphereListIterator begin() { return _subnodes.begin(); }
81 inline SphereListIterator end() { return _subnodes.end(); }
82
83 inline bool operator == (const Sphere & b) const
84 {
85 return p == b.p && r == b.r && owner == b.owner;
86 }
87
88 inline size_t count() { return _subnodes.size(); }
89
90 // Connect an existing sphere to this sphere as a child
91 void attach (const mars::ptr< Sphere > & pSphere)
92 {
93 assert (mars::MAG(pSphere->p) > 1);
94 _subnodes.push_back(pSphere);
95 }
96 Sphere * attachCulled( const mars::RangeX<float> &rads, const mars::RangeX<float> &rrng, const mars::SphericalCoords<T> &spc0 )
97 {
98 using namespace mars;
99
100 const T r = static_cast <T> (rads.next());
101 SphericalCoords <T> spc1(
102 r,
103 rrng.next(),
104 rrng.next()
105 );
106
107 // The closer a child is to the intersection of the two spheres (parent and current) the smaller it will be
108 if (rrng.delta >= owner->_minimum)
109 spc1.r = static_cast< T > (static_cast< float > (spc1.r) * computeCull(spc1.azimuth, spc1.zenith));
110
111 // Ignore any and all children that are too small, just discard the information
112 if (spc1.r >= owner->_minimum)
113 {
114 // Apply new randomly generated angles to the angles between the two spheres and the corresponding radius
115 return & attach(
116 spc1.zenith + spc0.zenith,
117 spc1.azimuth + spc0.azimuth,
118 spc1.r / 2,
119 spc1.r
120 );
121 } else
122 return NULL;
123 }
124 // Create and connect a new sphere with the specified radius at the specified spherical coordinates relative to this sphere
125 Sphere & attach (const float zenith, const float azimuth, const T dr, const T rr)
126 {
127 using namespace mars;
128
129 // Sphere is attached at the radius of the parent sphere plus the specified amount and at the spherical angles specified
130 _subnodes.push_back(mars::ptr< Sphere > (new Sphere(owner, SphericalCoords<T>(r + dr, zenith, azimuth) + p, rr)));
131 assert (mars::MAG(_subnodes.back()->p) > 1);
132 return *_subnodes.back();
133 }
134 template< typename PTR, typename R, typename G >
135 void store (mars::LookupTree< typename PTR::Type, R, G > & store, const PTR & pElem) const
136 {
137 store.add(pElem, *this);
138 for (typename SphereList::const_iterator it = _subnodes.begin(); it != _subnodes.end(); ++it)
139 (*it)->store(store, pElem);
140 }
141 };
142
143 private:
144 SphereList _toptier;
145
146 const float _divisor;
147 const mars::RangeX <float> _sizeflux;
148 const mars::RangeX <unsigned short> _kidflux;
149 const T _minimum, _maxtier;
150
151 T _limit; // Non-const for subdivision
152
153 inline float intersectAngle (const Sphere & sphere0, const Sphere & sphere) const
154 {
155 using namespace mars;
156
157 assert(mars::MAG(sphere.p - sphere0.p) < sphere.r + sphere0.r);
158 return acos(
159 static_cast< float > (MAGSQ(sphere.p - sphere0.p) + SQ(sphere0.r) - SQ(sphere.r)) / static_cast< float > (2 * MAG(sphere.p - sphere0.p) * sphere0.r)
160 );
161 }
162
163 inline T distance (const T r1, const T r2) const
164 {
165 if (r1 > r2)
166 return r1 + r2 / 2;
167 else
168 return r2 + r1 / 2;
169 }
170
171 mars::ptr< Sphere > spawnBetween (const Sphere & sphere0, const Sphere & sphere1)
172 {
173 using namespace mars;
174 using namespace std;
175
176 // The vector between the two precursor spheres
177 const vector3D< T > dp = sphere1.p - sphere0.p;
178
179 // The scalar and the square of the distance between the two precursor spheres
180 const T dsq = MAGSQ(dp), d = sqrt(dsq);
181
182 // New sphere radius must be large enough to sufficiently overlap both the parent and child spheres
183 // but yet not exceed the parent sphere's radius
184 const RangeX <T> rrng(
185 sphere0.r / 2 + (d - sphere0.r - sphere1.r) / 2,
186 std::max(sphere0.r, sphere1.r)
187 );
188 const T r = rrng.next();
189
190 const T
191 // Scalar and square distances between the sphere centers of the parent sphere and the newly spliced sphere
192 br = distance(sphere0.r, r), brsq = SQ(br),
193
194 // Distance between the sphere centers of the newly spliced sphere and the child sphere
195 ar = distance(r, sphere1.r),
196
197 // Projection of the vector between the sphere centers of the parent sphere and the newly spliced sphere onto the vector between the parent sphere and the child sphere
198 ch = (brsq + dsq - SQ(ar)) / (2*d);
199
200 assert (br + ar >= d);
201
202 // Setup the rotation about the axis defined by the aforementioned vector
203 const mars::Quaternion< T >
204 q = mars::Quaternion< T >::rotation(RANDf(PI * 2), mars::U(dp));
205
206 const vector3D< T >
207 // First-half: Make an arbitrarily perpendicular vector to the aforementioned axis
208 // Second-half: Use the law of cosines to determine the height of the triangle
209 // (this would be the "opposite" shared by two imaginary right-triangles that compose the triangle)
210 vb = mars::U(vector3D< T > (dp.y, -(dp.x + dp.z), dp.y)) * sqrt(brsq - SQ(ch)),
211
212 // Apply the rotation
213 vbr = q & vb,
214
215 // Form the baseline vector
216 vch = mars::U(dp) * ch,
217
218 // Calculate the hypotenuse vector
219 // - hypotenuse: distance between sphere centers of parent sphere and newly spliced sphere
220 // - opposite: rotated vector
221 pp = vch + vbr;
222
223 assert(SQ(MAG(pp) - br) < 0.0001);
224 assert(MAG(sphere1.p - (sphere0.p + pp)) < sphere1.r + r);
225 assert(d - (MAG(pp) + r + sphere1.r) < 0);
226
227 return mars::ptr< Sphere > (new Sphere (*this, sphere0.p + pp, r));
228 }
229
230 void spawn (const Sphere & sphere0, Sphere & sphere, const T nTier)
231 {
232 using namespace mars;
233
234 static const float pi = static_cast <float> (PI);
235
236 // Next number of breadth kids to generate for this sphere
237 const unsigned short amt = _kidflux.next();
238
239 // Allowed radii for children
240 const RangeX <float> rads (
241 _sizeflux.minimum * static_cast <float> (sphere.r),
242 _sizeflux.maximum * static_cast <float> (sphere.r)
243 );
244
245 // Spherical coordinate representation of vector between parent sphere and this current sphere
246 const SphericalCoords<T> spc0 (sphere0.p - sphere.p);
247 RangeX<float> rrng = naturalArcRange(sphere0, sphere, pi);
248
249 for (unsigned int c = sphere.count(); c < amt; ++c)
250 {
251 Sphere * neww = sphere.attachCulled(rads, rrng, spc0);
252
253 // Do not recurse if we are as deep as we are allowed to go
254 // or if the size of the sphere is small enough and doesn't need more descendants in this branch
255 if (neww != NULL && neww->r > _limit && nTier < _maxtier)
256 spawn (sphere, *neww, nTier + 1);
257 }
258 }
259
260 void subdivide (const Sphere & sphere0, Sphere & sphere)
261 {
262 using namespace mars;
263 const T r = sphere.r;
264 const size_t len = sphere.count();
265 static const float pi = static_cast <float> (PI);
266
267 // Do not sub-divide this sphere if it is already too small
268 if (sphere.r / _divisor >= _limit)
269 {
270 sphere.r /= _divisor;
271
272 for (typename SphereList::iterator it = sphere.begin(); it != sphere.end(); ++it)
273 {
274 mars::ptr< Sphere > pSphTail = *it;
275 // First divide the current sphere and all successive spheres before adding any new spheres
276 subdivide(sphere, *pSphTail);
277
278 // Splice-in a randomly-placed sphere between this current sphere and the iterator current sphere
279 mars::ptr< Sphere > pSphere1 = spawnBetween(sphere, *pSphTail);
280
281 spawnImmediates(*pSphere1, sphere, *pSphTail);
282
283 pSphere1->attach(pSphTail);
284 *it = pSphere1;
285
286 // Generate allowed radii for additional children of the parent of the spliced sphere
287 // Must come after sphere.r /= 2 because it depends on the new radius
288 const RangeX <T> rads (
289 _sizeflux.minimum * static_cast <float> (pSphTail->r),
290 _sizeflux.maximum * static_cast <float> (pSphTail->r)
291 );
292
293 // Spherical coordinate representation of vector between parent sphere and this current sphere
294 const SphericalCoords<T> spc0 (pSphTail->p - pSphere1->p);
295 RangeX<float> rrng = naturalArcRange(*pSphTail, *pSphere1, pi);
296
297 // Add only as many new children as the original conditions permit for a new random quantity
298 for (size_t n = _kidflux.next(); n > len; --n)
299 {
300 pSphTail->attachCulled(rads, rrng, spc0);
301 }
302 }
303 }
304 }
305 void spawnImmediates( Sphere & sphere, Sphere & sphere0, Sphere & sphere1)
306 {
307 using namespace mars;
308 const short kfn = _kidflux.next();
309
310 // Add new children if and only if the new random amount is within the bounds and accounts for the spliced one's descendant
311 if (kfn > _kidflux.minimum)
312 {
313 // Accounts for the spliced one's descendant
314 const unsigned short amt = kfn - 1;
315
316 // Allowed radii for children
317 const RangeX <T> rads (
318 _sizeflux.minimum * static_cast <float> (sphere.r),
319 _sizeflux.maximum * static_cast <float> (sphere.r)
320 );
321
322 const SphericalCoords<T> spc (sphere0.p - sphere.p);
323
324 // Allowed range where neither the descendant nor ancestor intersect with this current sphere
325 const RangeX<float> rrng (
326 intersectAngle(sphere, sphere1),
327 PI - intersectAngle(sphere0, sphere),
328 0.02f / 2.0f
329 );
330
331 // Adds children to the spliced sphere
332 for (unsigned int c = 0; c < amt; ++c)
333 {
334 const T r = static_cast <T> (rads.next());
335 SphericalCoords <T> spc1(
336 r,
337 rrng.next() * 2.0f,
338 rrng.next() * 2.0f
339 );
340
341 // Since the two allowed ranges are disjoint, we compensate for the non-linear issue
342 if (spc1.azimuth > rrng.maximum)
343 spc1.azimuth = rrng.maximum - spc1.azimuth;
344 if (spc1.zenith > rrng.maximum)
345 spc1.zenith = rrng.maximum - spc1.zenith;
346
347 // The closer to the intersections the child is, the smaller it becomes
348 const float f = Sphere::computeCull(spc1.azimuth, spc1.zenith);
349
350 spc1.r *= (f + (1.0f - f)) / 2.0f;
351
352 // Attach the child
353 if (spc1.r >= _minimum)
354 {
355 sphere.attach(
356 spc.zenith + spc1.zenith,
357 spc.azimuth + spc1.azimuth,
358 spc1.r / 2,
359 spc1.r
360 );
361 }
362 }
363 }
364 }
365
366 inline mars::RangeX<float> naturalArcRange( const Sphere & sphere0, Sphere & sphere, const float pi ) const
367 {
368 // Depending on radii of the two spheres (parent and current), determine the angle of the arc spanning the intersection of the two spheres
369 const float deg = (sphere0 == sphere ? 0 : intersectAngle(sphere0, sphere));
370
371 // Allowed angular range dispositions for children
372 return mars::RangeX<float> (-pi + deg, +pi - deg, 0.02f);
373 }
374
375 template <typename I>
376 struct StackInfo
377 {
378 I i, end;
379
380 inline StackInfo(const I & begin, const I & end)
381 : i(begin), end(end) {}
382 };
383
384 public:
385 template< typename TRAITS = ConstIteratorTraits >
386 class ConstIterator : public std::iterator< std::input_iterator_tag, Sphere >
387 {
388 protected:
389 typename TRAITS::SphereType * _spCurr;
390
391 private:
392 typedef typename TRAITS::SphereListIterator SphereListIterator;
393
394 typename TRAITS::SpherePack * const _pack;
395 SphereListIterator _it;
396 StackInfo< SphereListIterator > * _stCurr;
397 std::stack< StackInfo< SphereListIterator > > _stack;
398
399 CR_CONTEXT;
400
401 void process()
402 {
403 CR_START();
404
405 _stack.push(StackInfo< SphereListIterator > (_pack->_toptier.begin(), _pack->_toptier.end()));
406 while (!_stack.empty())
407 {
408 _stCurr = & _stack.top();
409
410 while (_stCurr->i != _stCurr->end)
411 {
412 _spCurr = & **_stCurr->i++;
413
414 CR_RETURN_VOID;
415
416 if (_spCurr->begin() != _spCurr->end())
417 {
418 _stack.push(StackInfo< SphereListIterator > (_spCurr->begin(), _spCurr->end()));
419 _stCurr = & _stack.top();
420 continue;
421 }
422 }
423
424 _stack.pop();
425 }
426
427 _spCurr = NULL;
428
429 CR_END();
430 }
431
432 public:
433 ConstIterator (typename TRAITS::SpherePack * const pack)
434 : _pack(pack)
435 {
436 CR_INIT();
437 process();
438 }
439 inline operator bool () { return _spCurr != NULL; }
440 inline ConstIterator & operator++()
441 {
442 process();
443 return *this;
444 }
445 inline ConstIterator operator++(int)
446 {
447 ConstIterator old = *this;
448 process();
449 return old;
450 }
451 inline const Sphere * operator -> () const { return &(*_spCurr); }
452 inline const Sphere & operator * () const { return *_spCurr; }
453 };
454
455 class Iterator : public ConstIterator< IteratorTraits >
456 {
457 public:
458 Iterator (typename IteratorTraits::SpherePack * const pack)
459 : ConstIterator< IteratorTraits > (pack) {}
460
461 inline Sphere * operator -> () { return &(*this->_spCurr); }
462 inline Sphere & operator * () { return *this->_spCurr; }
463 };
464
465 typedef Iterator iterator;
466 typedef ConstIterator< ConstIteratorTraits > const_iterator;
467
468 public:
469 inline RandomSpherePack (
470 const T limit, // Size below-which spheres are no-longer generated, essentially characterizes the degree of refined detail
471 const T minimum, // The minimum size of a sphere (usually just 1.0)
472 const T maxtier, // The maximum number depth-iterations (as opposed to breadth), this also affects the level of detail
473 const mars::RangeX<float> & sizeflux, // Flux of size of spheres at starting tier
474 const mars::RangeX<unsigned short> & kidflux // Flux of number of breadth-iterations (as opposed to depth),
475 )
476 : _sizeflux(sizeflux), _kidflux(kidflux), _limit(limit), _minimum(minimum), _maxtier(maxtier),
477 // Numeric stability will result if (iff?) divisor is not less than 2 because the gap between sub-divided spheres would be too great
478 _divisor(4.0f / 3.0f)
479 {}
480
481 inline Sphere createEmptySphere () const
482 { return Sphere(this); }
483
484 inline const Sphere & spawn (const VECTOR & pt, const T size)
485 {
486 assert (size >= _minimum);
487
488 _toptier.push_back(
489 mars::ptr <Sphere> (
490 new Sphere(
491 this,
492 pt,
493 static_cast <T> (static_cast <float> (size / _sizeflux.minimum) * _sizeflux.next())
494 )
495 )
496 );
497
498 return * _toptier.back();
499 }
500
501 inline const Sphere & spawn (const float zenith, const float azimuth, const T size)
502 {
503 assert (size >= _minimum);
504
505 const T r = static_cast <T> (static_cast <float> (size / _sizeflux.minimum) * _sizeflux.next());
506 const Sphere & neighbor = *_toptier.back();
507
508 _toptier.push_back(
509 mars::ptr <Sphere> (
510 new Sphere(
511 this,
512 neighbor.p + mars::SphericalCoords< T > (distance (neighbor.r, r), zenith, azimuth),
513 r
514 )
515 )
516 );
517
518 return * _toptier.back();
519 }
520
521 inline void generateDescendants ()
522 {
523 for (typename SphereList::iterator it = _toptier.begin(); it != _toptier.end(); ++it)
524 spawn(**it, **it, 1);
525 }
526 inline void subdivide ()
527 {
528 if (_limit / _divisor >= _minimum)
529 _limit /= _divisor;
530 for (typename SphereList::iterator it = _toptier.begin(); it != _toptier.end(); ++it)
531 subdivide(**it, **it);
532
533 typename SphereList::iterator itTail, itHead;
534
535 itHead = _toptier.begin();
536 itTail = itHead++;
537
538 while (itHead != _toptier.end())
539 {
540 mars::ptr< Sphere > pSphereA = spawnBetween(**itTail, **itHead);
541
542 spawnImmediates(*pSphereA, **itTail, **itHead);
543 _toptier.insert (itHead, pSphereA);
544
545 itTail = itHead++;
546 }
547 }
548 inline iterator iterate () { return iterator(this); }
549 inline const_iterator iterate () const { return const_iterator(this); }
550 inline unsigned int count () const
551 {
552 std::stack< StackInfo< SphereListConstIterator > > stack;
553 unsigned int c = 1;
554
555 stack.push(StackInfo< SphereListConstIterator > (_toptier.begin(), _toptier.end()));
556
557 while (!stack.empty())
558 {
559 StackInfo< SphereListConstIterator > & curr = stack.top();
560
561 while (curr.i != curr.end)
562 {
563 const Sphere & sphere = **curr.i;
564
565 ++c;
566 ++curr.i;
567 if (sphere.begin() != sphere.end())
568 {
569 stack.push(curr);
570 curr = StackInfo< SphereListConstIterator >(sphere.begin(), sphere.end());
571 continue;
572 }
573 }
574
575 stack.pop();
576 }
577 return c;
578 }
579
580 template< typename PTR, typename R, typename G >
581 void store (mars::LookupTree< typename PTR::Type, R, G > & store, const PTR & pElem) const
582 {
583 for (typename SphereList::const_iterator it = _toptier.begin(); it != _toptier.end(); ++it)
584 (*it)->store(store, pElem);
585 }
586
587 // TODO: Find a purpose for this method, may not be necessary, was primordially conceived
588 template< typename E, typename R, typename G, typename QR >
589 float probability (mars::LookupTree< E, R, G > & store, const QR & search, E * = NULL) const
590 {
591 using namespace mars;
592 float f = 0.0f;
593
594 for (typename mars::LookupTree< E, R, G >::EntityRadialIterator it = store.contents(search); it.hasMore(); ++it)
595 {
596 const R & region = it.region();
597 const float ff = GaussianFn<float>::f(MAGSQ(region.p - region.p), 1.0, SQ(region.r), 0, 1.0f, 0);
598
599 if (ff > f)
600 f = ff;
601 }
602 return f;
603 }
604 };
605}
606
607#endif
Note: See TracBrowser for help on using the repository browser.