source: Revenant/marslib/include/mars_splines.h

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

Get to a compile state for terrain procedural generation
  • Property mode set to 100755
File size: 25.3 KB
Line 
1#pragma once
2#pragma warning(disable : 4290)
3
4#include <vector>
5#include "mars_calc.h"
6#include "coroutines.h"
7
8#include <assert.h>
9
10namespace mars
11{
12 template <typename V, typename real = float>
13 class SplineBase
14 {
15 public:
16 typedef V VectorType;
17 typedef typename std::vector <V>::const_iterator PointsConstIterator;
18 typedef typename std::vector <V>::iterator PointsIterator;
19 typedef typename std::vector <V>::reverse_iterator PointsRIterator;
20 using Precision = real;
21
22 protected:
23 real _M[4][4];
24
25 struct MatV3x4
26 {
27 V p[4];
28
29 inline const V & operator [] (const unsigned int i) const { return p[i]; }
30 inline V & operator [] (const unsigned int i) { return p[i]; }
31 };
32 std::vector< MatV3x4 > _M4p;
33
34 private:
35 std::vector <V> _points;
36
37 static const unsigned int
38 PAD_FRONT = 1, PAD_BACK = 2,
39 TOTAL_PAD = PAD_FRONT + PAD_BACK;
40
41 void init (const std::vector <V> & points)
42 {
43 // Add padding points:
44 // - one at the beginning same as the first
45 // - two more at the end that toggle between the last two points in the original list
46 // NOTE: Dependency on this algorithm in CRSplineBase::move
47 assert(points.size() >= 2);
48
49 _points.reserve(points.size() + TOTAL_PAD);
50 _points.insert(_points.begin(), 1, points.front());
51 _points.push_back(*(points.rbegin() + 1));
52 _points.push_back(points.back());
53
54 _M4p.insert(_M4p.begin(), size(), MatV3x4());
55 computeMatrices();
56 }
57
58 void computeMatrices ()
59 {
60 const unsigned int nSize = size();
61 for (unsigned int i = 0; i < nSize; ++i)
62 {
63 _M4p[i][0] = _M[0][0]*_points[i] + _M[0][1]*_points[i + 1] + _M[0][2]*_points[i + 2] + _M[0][3]*_points[i + 3];
64 _M4p[i][1] = _M[1][0]*_points[i] + _M[1][1]*_points[i + 1] + _M[1][2]*_points[i + 2] + _M[1][3]*_points[i + 3];
65 _M4p[i][2] = _M[2][0]*_points[i] + _M[2][1]*_points[i + 1] + _M[2][2]*_points[i + 2] + _M[2][3]*_points[i + 3];
66 _M4p[i][3] = _M[3][0]*_points[i] + _M[3][1]*_points[i + 1] + _M[3][2]*_points[i + 2] + _M[3][3]*_points[i + 3];
67 }
68 }
69
70 public:
71
72 SplineBase (const SplineBase & copy) : _points (copy._points), _M4p(copy._M4p)
73 { memcpy (_M, copy._M, sizeof (real) * 4 * 4); }
74
75 SplineBase () : _M4p(NULL)
76 { _points.reserve(TOTAL_PAD + 2); }
77 SplineBase (const std::vector <V> & points) : _points (points)
78 { init(points); }
79
80 const V compute( const unsigned int i, const real t0, const real t1, const real t2, const real t3) const
81 {
82 assert(i < size());
83 return t0 * _M4p[i][0] +
84 t1 * _M4p[i][1] +
85 t2 * _M4p[i][2] +
86 t3 * _M4p[i][3];
87 }
88 template <typename V3D>
89 void move (const real radius, const V3D & unormal)
90 {
91 const V3D & normal = unormal * radius;
92 unsigned int c = 0;
93 PointsIterator i = _points.begin() + 1;
94 PointsRIterator r = _points.rbegin();
95 V t1, t2;
96 const unsigned int nCount = size() - 1;
97
98 t1 = U(compute (c++, 0, 1, 0, 0)) & normal;
99 while (c < nCount)
100 {
101 t2 = t1;
102 t1 = U(compute (c++, 0, 1, 0, 0)) & normal;
103 *i++ += t2;
104 }
105 t2 = t1;
106 t1 = U(compute(--c, 0, 1, 2, 3)) & normal;
107 *i++ += t2;
108 *i += t1;
109
110 // Account for padding points, see class constructor for description
111 // NOTE: Dependency of point-padding algorithm in constructor
112 _points[0] = _points[1];
113 *r++ = *(i - 1);
114 *r++ = *i;
115
116 // Recompute the cached matrices for each point
117 computeMatrices();
118 }
119
120 // For the following points-list manipulation methods, the 2 and 3
121 // constants represent padding at the beginning and at the end
122 // TODO: Extract spline point padding offsets into constants
123 inline unsigned int size () const { return _points.size() - TOTAL_PAD; }
124 inline V & front () { return *begin(); }
125 inline const V & front () const { return *begin(); }
126 inline V & back () { return *(end() - 1); }
127 inline const V & back () const { return *(end() - 1); }
128 inline PointsConstIterator begin() const { return _points.begin() + PAD_FRONT; }
129 inline PointsConstIterator end() const { return _points.end() - PAD_BACK; }
130 inline PointsIterator begin() { return _points.begin() + PAD_FRONT; }
131 inline PointsIterator end() { return _points.end() - TOTAL_PAD; }
132 inline const V & operator [] (const unsigned int i) const
133 {
134 assert(i + PAD_FRONT < _points.size() - PAD_BACK);
135 return _points[i + PAD_FRONT];
136 }
137
138 void operator += (const V & ptOffset)
139 {
140 for (PointsIterator i = _points.begin(); i != _points.end(); ++i)
141 *i += ptOffset;
142
143 // Recompute the cached matrices for each point
144 computeMatrices();
145 }
146 };
147
148 template <typename V, typename real = float>
149 class CubicBSplineBase : public SplineBase <V, real>
150 {
151 private:
152 // TODO: Optimize by initializing this matrix only once rather than each time a spline object is created
153 void init ()
154 {
155 SplineBase<V, real>::_M[0][0] = static_cast <real> (1.0/6.0);
156 SplineBase<V, real>::_M[0][1] = static_cast <real> (4.0/6.0);
157 SplineBase<V, real>::_M[0][2] = static_cast <real> (1/6.0);
158 SplineBase<V, real>::_M[0][3] = 0;
159
160 SplineBase<V, real>::_M[1][0] = static_cast <real> (-3.0/6.0);
161 SplineBase<V, real>::_M[1][1] = static_cast <real> (0.0);
162 SplineBase<V, real>::_M[1][2] = static_cast <real> (3/6.0);
163 SplineBase<V, real>::_M[1][3] = 0;
164
165 SplineBase<V, real>::_M[2][0] = static_cast <real> (3.0/6.0);
166 SplineBase<V, real>::_M[2][1] = static_cast <real> (-6.0/6.0);
167 SplineBase<V, real>::_M[2][2] = static_cast <real> (3/6.0);
168 SplineBase<V, real>::_M[2][3] = 0;
169
170 SplineBase<V, real>::_M[3][0] = static_cast <real> (-1.0/6.0);
171 SplineBase<V, real>::_M[3][1] = static_cast <real> (3.0/6.0);
172 SplineBase<V, real>::_M[3][2] = static_cast <real> (-3.0/6.0);
173 SplineBase<V, real>::_M[3][3] = static_cast <real> (1.0/6.0);
174 }
175 public:
176 CubicBSplineBase () : SplineBase<V, real> ()
177 {
178 init();
179 }
180 CubicBSplineBase (const std::vector <V> & points) : SplineBase<V, real> (points)
181 {
182 init();
183 }
184 };
185
186 template <typename V, typename real = float>
187 class CRSplineBase : public SplineBase <V, real>
188 {
189 private:
190 // TODO: Optimize by initializing this matrix only once rather than each time a spline object is created
191 void init (const real fT)
192 {
193 using S = SplineBase<V, real>;
194
195 S::_M[0][0] = 0; S::_M[0][1] = 1; S::_M[0][2] = 0; S::_M[0][3] = 0;
196 S::_M[1][0] = -fT; S::_M[1][1] = 0; S::_M[1][2] = fT; S::_M[1][3] = 0;
197 S::_M[2][0] = 2*fT; S::_M[2][1] = fT-3; S::_M[2][2] = 3-2*fT; S::_M[2][3] = -fT;
198 S::_M[3][0] = -fT; S::_M[3][1] = 2-fT; S::_M[3][2] = fT-2; S::_M[3][3] = fT;
199 }
200
201 public:
202 CRSplineBase (const real fT = 0.5f) : SplineBase<V, real>()
203 {
204 init (fT);
205 }
206 CRSplineBase (const std::vector <V> & points, const real fT = 0.5f) : SplineBase<V, real>(points)
207 {
208 init(fT);
209 }
210 };
211
212 template <typename SB, typename real>
213 class ArcLengthFnPiece : public GaussLegendreInt <real, ArcLengthFnPiece <SB, real>, 5>
214 {
215 private:
216 const SB & _link;
217 const unsigned int _i;
218 mutable real _cache100;
219
220 // Copying will break the reference with SB
221 ArcLengthFnPiece (const ArcLengthFnPiece & copy) {}
222
223 public:
224 ArcLengthFnPiece (const SB & link, const unsigned int i)
225 : _link(link), _i(i), _cache100(-1) {}
226
227 inline real f (const real t) const
228 {
229 return MAG(_link.compute (_i, 0, 1, 2*t, 3*t*t));
230 }
231 inline real length (const real a, const real b) const
232 {
233#ifdef _DEBUG
234 assert (b >= a);
235 assert (b <= 1.0 && a >= 0.0 );
236#endif
237 return GaussLegendreInt<real, ArcLengthFnPiece <SB, real>, 5>::compute(a, b);
238 }
239
240 inline real length () const
241 { return _cache100 < 0 ? (_cache100 = length(0, 1)) : _cache100; }
242
243 inline void clearCache () const
244 { _cache100 = -1; }
245 };
246
247 template <typename SB>
248 class Spline
249 {
250 public:
251 typedef SB BaseType;
252 typedef Spline< SB > SplineType;
253 typedef typename SB::Precision real;
254 typedef typename SB::VectorType VectorType;
255
256 const real ONE;
257
258 private:
259 SB _sp;
260 ArcLengthFnPiece<SB, real> ** _lenfns;
261
262 // Approximates the "time" along the curve at the specified physical distance
263 // This is the inverse of arclength
264 class ArcPositionFn
265 {
266 private:
267 const ArcLengthFnPiece<SB, real> & _fnLen;
268 real _fTest;
269
270 ArcPositionFn (const ArcPositionFn & copy) {}
271
272 public:
273 ArcPositionFn (const ArcLengthFnPiece<SB, real> & fnLen, const real fPosition)
274 : _fnLen(fnLen), _fTest(fPosition) {}
275
276 inline const real f(const real t) const { return _fTest - _fnLen.length(0, t); }
277 };
278
279 public:
280 typedef typename std::vector <VectorType>::const_iterator PointsConstIterator;
281 typedef typename std::vector <VectorType>::iterator PointsIterator;
282 typedef typename std::vector <VectorType>::reverse_iterator PointsRIterator;
283
284 private:
285 template <typename L>
286 inline static void populateVector (const typename L::const_iterator ia, const typename L::const_iterator ib, std::vector <VectorType> & v)
287 { v.insert(v.begin(), ia, ib); }
288
289 void init ()
290 {
291 const unsigned int nSegCount = segCount();
292
293 _lenfns = new ArcLengthFnPiece< SB, real >* [nSegCount];
294 for (unsigned int i = 0; i < nSegCount; ++i)
295 _lenfns[i] = new ArcLengthFnPiece<SB, real> (_sp, i);
296 }
297 inline const ArcLengthFnPiece<SB, real> & lenfn (const unsigned int i) const
298 {
299#ifdef _DEBUG
300 assert(i < segCount());
301 assert(_lenfns != NULL);
302 assert(_lenfns[i] != NULL);
303#endif
304 return *_lenfns[i];
305 }
306 void clearCache () const
307 {
308 for (unsigned int i = 0; i < segCount(); ++i)
309 lenfn(i).clearCache();
310 }
311
312 public:
313 Spline (const Spline <SB> & copy)
314 : _lenfns(NULL), _sp (copy._sp), ONE(static_cast< real > (1)) { init (); }
315 Spline (const std::vector <VectorType> & points)
316 : _lenfns(NULL), _sp(points), ONE(static_cast< real > (1)) { init(); }
317 Spline (const real fT, const std::vector <VectorType> & points)
318 : _lenfns(NULL), _sp(points, fT), ONE(static_cast< real > (1)) { init(); }
319 Spline ()
320 : _lenfns(NULL), _sp(), ONE(static_cast< real > (1)) { init (); }
321 Spline (const real fT)
322 : _lenfns(NULL), _sp(fT), ONE(static_cast< real > (1)) { init (); }
323
324 template <typename L>
325 static Spline createFrom (const typename L::const_iterator ia, const typename L::const_iterator ib)
326 {
327 std::vector <VectorType> v;
328 populateVector <L> (ia, ib, v);
329 return Spline(v);
330 }
331
332 template <typename L>
333 static Spline createFrom (const real fT, const typename L::const_iterator ia, const typename L::const_iterator ib)
334 {
335 std::vector <VectorType> v;
336 populateVector <L> (ia, ib, v);
337 return Spline(fT, v);
338 }
339
340 // Returns the length of the specified curve segment [i, i+1]
341 inline real len (const unsigned int i) const { return lenfn(i).length(); }
342 // Returns the length of the specified curve segment between times tA and tB [i+tA, i+tB]
343 inline real len (const unsigned int i, const real tA, const real tB) const { return lenfn(i).length(tA, tB); }
344 // Returns the total length of the curve
345 real len () const
346 {
347 const unsigned nSegs = segCount();
348 real nLen = 0;
349
350 for (unsigned i = 0; i < nSegs; ++i)
351 nLen += len(i);
352
353 return nLen;
354 }
355 real distanceFast (const unsigned int i, const real t) const
356 {
357 real fLen = 0;
358
359 for (unsigned int c = 0; c < i; ++c)
360 fLen += len(c);
361
362 return fLen += len(i) * t;
363 }
364 real distance (const unsigned int i, const real t) const
365 {
366 real fLen = 0;
367
368 for (unsigned int c = 0; c < i; ++c)
369 fLen += len(c);
370
371 return fLen += lenfn(i).length(0, t);
372 }
373 inline real timeAt (const unsigned int i, const real & fPos) const
374 {
375 return mars::BisectMethod< real, ArcPositionFn > :: solve ( ArcPositionFn(lenfn(i), fPos), 0, 1 );
376 }
377 inline const VectorType compute (const unsigned int i, const real t) const
378 {
379 const real
380 t2 = t * t,
381 t3 = t2 * t;
382
383 return _sp.compute(i, 1, t, t2, t3);
384 }
385 inline const VectorType gradient (const unsigned int i, const real t) const
386 { return _sp.compute(i, 0, 1, 2*t, 3*t*t); }
387 inline const VectorType laplace (const unsigned int i, const real t) const
388 { return _sp.compute(i, 0, 0, 2, 6*t); }
389
390 template <typename V3D>
391 void move (const real radius, const V3D & unormal)
392 {
393 _sp.move (radius, unormal);
394 clearCache();
395 }
396 inline unsigned int segCount() const { return _sp.size() - 1; }
397 inline unsigned int size () const { return _sp.size(); }
398 inline PointsConstIterator begin() const { return _sp.begin(); }
399 inline PointsConstIterator end() const { return _sp.end(); }
400 inline PointsIterator begin() { return _sp.begin(); }
401 inline PointsIterator end() { return _sp.end(); }
402 inline const VectorType & operator [] (const unsigned int i) const { return _sp[i]; }
403 inline VectorType & operator [] (const unsigned int i) { return _sp[i]; }
404 inline VectorType & front () { return _sp.front(); }
405 inline const VectorType & front () const { return _sp.front(); }
406 inline VectorType & back () { return _sp.back(); }
407 inline const VectorType & back () const { return _sp.back(); }
408
409 inline void operator += (const VectorType & ptOffset)
410 { _sp += ptOffset; }
411
412 ~Spline()
413 {
414 for (unsigned int i = 0; i < segCount(); ++i)
415 delete _lenfns[i];
416
417 delete [] _lenfns;
418 }
419
420 };
421
422 template <typename S, typename Vs, typename integer = short>
423 class SplineRasterIterator : public std::iterator < std::input_iterator_tag, Vs >
424 {
425 public:
426 typedef S SplineType;
427 typedef Vs RasterVecType;
428 typedef integer RasterPrecision;
429 typedef typename SplineType::real real;
430 typedef typename SplineType::VectorType VectorType;
431
432 private:
433 const SplineType & _sp;
434 Vs _zp, // Resulting iterative raster point (also used to compare guess points against)
435 _tzp; // Guess point (integral) (current)
436 VectorType _trp; // Guess point (rational) (current)
437 unsigned int _i, // Segment number (current)
438 _k; // Iteration count
439 real _accel, // Acceleration (current) (added to current time for next guess point)
440 _t, // Time (current)
441 _t0; // Time (last)
442 const unsigned int
443 _segCount;
444
445 CR_CONTEXT;
446
447 int calcZDiff (const real & t)
448 {
449 using namespace mars;
450
451 _trp = _sp.compute(_i, t);
452 _tzp = _trp;
453
454 const Vs dzp = _tzp - _zp;
455
456 return (SQ(dzp.x) + SQ(dzp.y));
457 }
458 inline void assignAccel (const real acc)
459 {
460 assert(acc >= 0);
461 if (acc > 0)
462 _accel = acc;
463 }
464 bool satisfy (real & acc, unsigned int & zdiff)
465 {
466 using namespace mars;
467 const real t = _t + acc;
468
469 // Calculate integral distance from last position to guessed position
470 zdiff = calcZDiff(t);
471 if (zdiff == 1)
472 {
473 _zp = _tzp;
474 _t = t;
475 assignAccel(acc);
476 return true;
477 } else if (zdiff <= 2)
478 {
479 const VectorType drp = _trp - _zp;
480
481 if (drp.x > abs(drp.y) && drp.x > 0) // Slope is bias to the positive x-direction
482 ++_zp.x;
483 else if (drp.y > abs(drp.x) && drp.y > 0) // Slope is bias to the positive y-direction
484 ++_zp.y;
485 else if (drp.y < -abs(drp.x) && drp.y < 0) // Slope is bias to the negative y-direction
486 --_zp.y;
487 else if (drp.x < -abs(drp.y) && drp.x < 0) // Slope is bias to the negative x-direction
488 --_zp.x;
489 else // Slope is perfectly 1:1 or 0:0
490 {
491 if (zdiff == 0)
492 {
493 if (drp.x == 0 || drp.y == 0)
494 return false;
495
496 if (drp.x < 0)
497 --_zp.x;
498 else
499 ++_zp.x;
500
501 if (drp.y < 0)
502 --_zp.y;
503 else
504 ++_zp.y;
505 } else
506 _zp = _tzp;
507 assignAccel(acc);
508 _t = t;
509 return true;
510 }
511
512 if (zdiff == 2)
513 {
514 assignAccel(acc);
515 _t += _accel / 2;
516 } else
517 {
518 assignAccel(acc);
519 _t = t;
520 }
521 return true;
522 } else
523 return false;
524 }
525 template <int N> struct StaticIteration {};
526 template <int N> void drilldown (real & acc, unsigned int & zdiff) {
527 drilldown(acc, zdiff, StaticIteration<N>());
528 }
529 template <int N> void drilldown (real & acc, unsigned int & zdiff, StaticIteration<N>)
530 {
531 using namespace mars;
532
533 // If the guessed position is not satisfied then it is will have an integral non-zero distance of at least +4 from the last position
534 if (!satisfy (acc, zdiff))
535 {
536 if (_trp != _zp)
537 {
538 const real accbak = acc;
539
540 acc /= MAG(_trp - _zp);
541 if (acc >= 1.0)
542 acc = accbak * 2;
543 } else
544 acc *= 2;
545
546 drilldown <N - 1> (acc, zdiff);
547 }
548 }
549 inline void drilldown (real & acc, unsigned int & zdiff, StaticIteration<0>)
550 {
551 if (!satisfy (acc, zdiff))
552 {
553 drilldownBeta(acc);
554 }
555 }
556 void drilldownBeta(real & acc)
557 {
558 using namespace mars;
559 static const real SPLINE_PREC_EPSILON = std::numeric_limits <real>::epsilon();
560
561 real dacc = acc;
562 signed short ddd = 3;
563 unsigned int zdiff;
564
565 do
566 {
567 if (!satisfy(acc, zdiff))
568 {
569 if (zdiff == 0)
570 {
571 if (ddd ^ 1)
572 {
573 dacc /= 2;
574 ddd = 1;
575 }
576 acc += dacc;
577 } else
578 {
579 if (ddd ^ 0)
580 {
581 dacc /= 2;
582 ddd = 0;
583 }
584 acc -= dacc;
585 }
586 } else
587 return;
588
589 // Acc < 1.0 prevents an index-out-of-bounds exception and checking delta-acc here ensures numerical stability
590 } while (acc < 1.0 && acc > 0 && (dacc > SPLINE_PREC_EPSILON || dacc < -SPLINE_PREC_EPSILON));
591
592 // Catch-all for acceleration values that are negative or too small, accounts for numerical instability as well as monotonic exceptions
593 if (acc > SPLINE_PREC_EPSILON)
594 _accel = acc;
595
596 _zp = _tzp;
597 _t += _accel;
598 }
599
600 void findAdjacentCell ()
601 {
602 using namespace mars;
603
604 real
605 acc = _accel;
606 unsigned int
607 zdiff;
608
609 drilldown <6> (acc, zdiff);
610 }
611 void process ()
612 {
613 CR_START();
614
615 _t = 0;
616 _zp = _sp.compute (0, 0);
617
618 for (_i = 0; _i < _segCount; ++_i)
619 {
620 _accel = static_cast <real> (1.0) / _sp.len(_i);
621 _t0 = _t = 0;
622
623 while (_t < 1.0)
624 {
625 findAdjacentCell();
626 CR_RETURN_VOID;
627 _t0 = _t;
628 ++_k;
629 }
630 }
631
632 CR_END();
633 }
634
635 public:
636 SplineRasterIterator (const SplineRasterIterator & copy);
637 SplineRasterIterator (const SplineType & sp)
638 : _sp (sp), _segCount (sp.segCount()), _k(0)
639 {
640 CR_INIT();
641
642 process();
643 }
644
645 // XTODO: Optimize!!!
646 SplineRasterIterator & operator += (const real rDist)
647 {
648 if (rDist > 0)
649 {
650 const unsigned int nSegs = _sp.segCount();
651 const real fDist = distanceFast();
652 const real fTot = static_cast <real> (rDist) + fDist;
653 real fAcc = fDist;
654
655 while (_i < nSegs)
656 {
657 const real fTesterz = _sp.len(_i);
658
659 if (fAcc + fTesterz > fTot)
660 {
661 _t0 = _t = _sp.timeAt (_i, fTot - fAcc);
662 _zp = _sp.compute(_i, _t);
663 _accel = static_cast <real> (1.0) / _sp.len(_i);
664 _k += static_cast <unsigned int> (rDist);
665 return *this;
666 }
667 else
668 fAcc += fTesterz;
669 ++_i;
670 }
671
672 // Terminate the loop if we get here
673 _t0 = _t = 1.0f;
674 }
675
676 return *this;
677 }
678
679 inline SplineRasterIterator & operator++()
680 {
681 process();
682 return *this;
683 }
684 inline SplineRasterIterator operator++(int)
685 {
686 SplineRasterIterator old = *this;
687
688 process();
689 return old;
690 }
691 // Retrieves the number of pixels traversed thus far, or the current pixel in the iteration
692 // NOTE: This is not the same as the distance of spline curve traversed thus far, so it cannot be
693 // compared accurately with the spline length since this deals with number of pixels plotted as
694 // opposed to any physical measure of distance from the origin.
695 inline unsigned int index () const { return _k; }
696 inline unsigned int segment () const { return _i; }
697 inline real distanceFast() const { return _sp.distanceFast(_i, _t0); }
698 inline real distance() const { return _sp.distance(_i, _t0); }
699 inline const VectorType gradient () const { return _sp.gradient(_i, _t0); }
700 inline const VectorType laplace () const { return _sp.laplace(_i, _t0); }
701
702 inline operator const Vs & () const { return _zp; }
703 inline operator bool () const { return _i < _segCount; }
704 inline const SplineType & getSpline() const { return _sp; }
705 };
706
707 template <typename S, typename EOFn, typename EIFn, typename RasterPrec, typename Vs, typename V3>
708 class SplineZoneIterator : public std::iterator < std::input_iterator_tag, Vs >
709 {
710 protected:
711 typedef S SplineType;
712 typedef typename SplineType::real real;
713 typedef SplineRasterIterator <S, Vs, RasterPrec> SplineWalker;
714 typedef typename SplineType::VectorType VectorType;
715
716 private:
717
718 mutable struct TotalRadiusCache
719 {
720 unsigned int result, index, radi;
721
722 inline TotalRadiusCache ()
723 : result(0), index(0), radi(0) {}
724
725 } _crad;
726
727 const unsigned int _netradius; // TODO: Genericize unsigned int _netradius
728 const EOFn _fnEaseOut;
729 const EIFn _fnEaseIn;
730 unsigned int _splnum;
731 signed int _side;
732 SplineWalker * _pItSP;
733 unsigned int _radi;
734 real _fBaseLen, _fLimit, _fOffset, _fSpLen;
735 SplineType * _pSplineA, * _pSplineB;
736 V3 _unormal;
737 typename SplineWalker::RasterVecType _p;
738 bool _done;
739
740 CR_CONTEXT;
741
742 public:
743 SplineZoneIterator (const SplineType & spline, const unsigned int nRadius, const EOFn & fnEaseOut, const EIFn & fnEaseIn, const V3 & unormal)
744 : _netradius (nRadius), _unormal (unormal), _done(false), _pItSP(NULL), _splnum(0), _fnEaseOut(fnEaseOut), _fnEaseIn(fnEaseIn), _radi(0)
745 {
746 using namespace mars;
747
748 CR_INIT();
749 _pSplineB = new S (spline);
750 _pSplineA = new S (spline);
751
752 process();
753 }
754 ~SplineZoneIterator ()
755 {
756 delete _pSplineA;
757 delete _pSplineB;
758 delete _pItSP;
759 }
760
761 inline SplineZoneIterator & operator++()
762 {
763 process();
764 return *this;
765 }
766 inline SplineZoneIterator operator++(int)
767 {
768 SplineZoneIterator old = *this;
769
770 process();
771 return old;
772 }
773
774 inline const Vs & operator * () const { return _p; }
775 inline const Vs & operator -> () const { return _p; }
776
777 inline operator const Vs & () const { return _p; }
778 inline const typename SplineWalker::RasterVecType & pos () const { return _p; }
779
780 // Returns current radius as a ratio to the total radius
781 inline real trad () const { return static_cast <real> (radius()) / static_cast <real> (SplineZoneIterator::crad()); }
782 // Returns the present distance from the spine spline (always positive)
783 inline const unsigned int radius () const { return _radi; }
784 // Returns the present distance from the spine spline (negative if on the negative side and positive if on the positive side)
785 inline const signed int sradius () const { return static_cast< signed int > (radius()) * side(); }
786 inline const signed int side() const { return _side; }
787 // Returns the total radius of the zone at the current time
788 const unsigned int crad () const
789 {
790 if (_crad.result == 0 || _pItSP->index() != _crad.index || _radi != _crad.radi)
791 {
792 const real
793 ONE = _pItSP->getSpline().ONE,
794 frDist = static_cast <real> (_pItSP->distanceFast()) / _fSpLen;
795
796 // FIXME: _radi ends-up being more than the calculated total radius causing the assertion to fail (see below).
797 // This may have been fixed as of the change to use accurate spline length calculation instead of approximate
798 _crad.result =
799 std::max(
800 _radi,
801 static_cast <unsigned int> (
802 _fnEaseOut.invf(frDist) * _fnEaseIn.invf(ONE - frDist) * static_cast <real> (_netradius)
803 ) + 1
804 )
805 ;
806 _crad.index = _pItSP->index();
807 _crad.radi = _radi;
808
809 assert(_crad.radi <= _crad.result);
810 }
811
812 return _crad.result;
813 }
814 inline operator bool () const { return !_done; }
815 inline const unsigned int splineNum () const { return _splnum; }
816 inline const unsigned int index () const { return _pItSP->index(); }
817 inline const real distance () const { return _pItSP->distance(); }
818 inline const real distanceFast () const { return _pItSP->distanceFast(); }
819 inline const VectorType gradient () const { return _pItSP->gradient(); }
820 inline const VectorType laplace () const { return _pItSP->laplace(); }
821 inline const real seglen () const { return _pItSP->getSpline().len(_pItSP->segment()); }
822 inline const VectorType & tailCP () const { return _pItSP->getSpline()[_pItSP->segment()]; }
823 inline const VectorType & leadCP () const { return _pItSP->getSpline()[_pItSP->segment() + 1]; }
824
825 private:
826 void process ()
827 {
828 using namespace mars;
829
830 CR_START();
831
832 _fSpLen =
833 _fBaseLen = _pSplineB->len();
834 for (_pItSP = new SplineWalker (*_pSplineB); *_pItSP; ++*_pItSP)
835 {
836 _p = *_pItSP;
837 CR_RETURN_VOID;
838 }
839 delete _pItSP;
840 _pItSP = NULL;
841
842 // Fill to the positive side of the curve
843 _side = +1;
844 while (_radi++ < _netradius)
845 {
846 prepareWalk(_pSplineA);
847 for (*(_pItSP = new SplineWalker (*_pSplineA)) += _fOffset;
848 *_pItSP && _pItSP->distanceFast() < _fLimit;
849 ++*_pItSP)
850 {
851 _p = *_pItSP;
852 CR_RETURN_VOID;
853 }
854 postWalk(_pSplineA);
855 }
856 _pItSP = NULL;
857
858 // Fill to the negative side of the curve
859 _radi = 0;
860 _side = -1;
861 while (_radi++ < _netradius)
862 {
863 prepareWalk(_pSplineB);
864
865 for (*(_pItSP = new SplineWalker (*_pSplineB)) += _fOffset;
866 *_pItSP && _pItSP->distanceFast() < _fLimit;
867 ++*_pItSP)
868 {
869 _p = *_pItSP;
870 CR_RETURN_VOID;
871 }
872 postWalk(_pSplineB);
873 }
874 _pItSP = NULL;
875
876 _done = true;
877
878 CR_END();
879 }
880 inline void postWalk(const S * pSpline)
881 {
882 delete _pItSP;
883 }
884 void prepareWalk (S * pSpline)
885 {
886 pSpline->move (static_cast< float > (_side), _unormal);
887 _fSpLen = pSpline->len();
888 const real
889 ONE = pSpline->ONE,
890 frRadii = static_cast <real> (_radi) / static_cast <real> (_netradius);
891
892 _fLimit = _fnEaseOut.f( frRadii / _fnEaseIn.invf(ONE - _fnEaseIn.f(frRadii) ) ) * _fSpLen;
893 _fOffset = (ONE - _fnEaseIn.f( frRadii / _fnEaseOut.invf(ONE - _fnEaseOut.f(frRadii ) ) )) * _fSpLen;
894
895 // FIXME: Limit and offset calculations are inaccurate, causes _radi to increment out-of-range.
896 ++_splnum;
897 }
898 };
899}
Note: See TracBrowser for help on using the repository browser.