#pragma once
#pragma warning(disable : 4290)

#include <vector>
#include "mars_calc.h"
#include "coroutines.h"

#include <assert.h>

namespace mars
{
	template <typename V, typename real = float>
	class SplineBase
	{
	public:
		typedef V VectorType;
		typedef typename std::vector <V>::const_iterator PointsConstIterator;
		typedef typename std::vector <V>::iterator PointsIterator;
		typedef typename std::vector <V>::reverse_iterator PointsRIterator;
        using Precision = real;

	protected:
		real _M[4][4];

		struct MatV3x4
		{
			V p[4];

			inline const V & operator [] (const unsigned int i) const { return p[i]; }
			inline V & operator [] (const unsigned int i) { return p[i]; }
		};
		std::vector< MatV3x4 > _M4p;

	private:
		std::vector <V> _points;

		static const unsigned int 
			PAD_FRONT = 1, PAD_BACK = 2, 
			TOTAL_PAD = PAD_FRONT + PAD_BACK;

		void init (const std::vector <V> & points)
		{
			// Add padding points:
			// - one at the beginning same as the first
			// - two more at the end that toggle between the last two points in the original list
			// NOTE: Dependency on this algorithm in CRSplineBase::move
			assert(points.size() >= 2);

			_points.reserve(points.size() + TOTAL_PAD);
			_points.insert(_points.begin(), 1, points.front());
			_points.push_back(*(points.rbegin() + 1));
			_points.push_back(points.back());

			_M4p.insert(_M4p.begin(), size(), MatV3x4());
			computeMatrices();
		}

		void computeMatrices ()
		{
			const unsigned int nSize = size();
			for (unsigned int i = 0; i < nSize; ++i)
			{
				_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];
				_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];
				_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];
				_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];
			}
		}

	public:

		SplineBase (const SplineBase & copy) : _points (copy._points), _M4p(copy._M4p)
			{ memcpy (_M, copy._M, sizeof (real) * 4 * 4); }

		SplineBase () : _M4p(NULL)
			{ _points.reserve(TOTAL_PAD + 2); }
		SplineBase (const std::vector <V> & points) : _points (points)
			{ init(points); }

		const V compute( const unsigned int i, const real t0, const real t1, const real t2, const real t3) const
		{
			assert(i < size());
			return  t0 * _M4p[i][0] +
				t1 * _M4p[i][1] + 
				t2 * _M4p[i][2] + 
				t3 * _M4p[i][3];
		}
		template <typename V3D>
		void move (const real radius, const V3D & unormal)
		{
			const V3D & normal = unormal * radius;
			unsigned int c = 0;
			PointsIterator i = _points.begin() + 1;
			PointsRIterator r = _points.rbegin();
			V t1, t2;
			const unsigned int nCount = size() - 1;

			t1 = U(compute (c++, 0, 1, 0, 0)) & normal;
			while (c < nCount)
			{
				t2 = t1;
				t1 = U(compute (c++, 0, 1, 0, 0)) & normal;
				*i++ += t2;
			}
			t2 = t1;
			t1 = U(compute(--c, 0, 1, 2, 3)) & normal;
			*i++ += t2;
			*i += t1;

			// Account for padding points, see class constructor for description
			// NOTE: Dependency of point-padding algorithm in constructor
			_points[0] = _points[1];
			*r++ = *(i - 1);
			*r++ = *i;

			// Recompute the cached matrices for each point
			computeMatrices();
		}

		// For the following points-list manipulation methods, the 2 and 3 
		// constants represent padding at the beginning and at the end
		// TODO: Extract spline point padding offsets into constants
		inline unsigned int size () const { return _points.size() - TOTAL_PAD; }
		inline V & front () { return *begin(); }
		inline const V & front () const { return *begin(); }
		inline V & back () { return *(end() - 1); }
		inline const V & back () const { return *(end() - 1); }
		inline PointsConstIterator begin() const { return _points.begin() + PAD_FRONT; }
		inline PointsConstIterator end() const { return _points.end() - PAD_BACK; }
		inline PointsIterator begin() { return _points.begin() + PAD_FRONT; }
		inline PointsIterator end() { return _points.end() - TOTAL_PAD; }
		inline const V & operator [] (const unsigned int i) const 
		{
			assert(i + PAD_FRONT < _points.size() - PAD_BACK);
			return _points[i + PAD_FRONT];
		}

		void operator += (const V & ptOffset)
		{
			for (PointsIterator i = _points.begin(); i != _points.end(); ++i)
				*i += ptOffset;

			// Recompute the cached matrices for each point
			computeMatrices();
		}
	};

	template <typename V, typename real = float>
	class CubicBSplineBase : public SplineBase <V, real>
	{
	private:
		// TODO: Optimize by initializing this matrix only once rather than each time a spline object is created
		void init ()
		{
			SplineBase<V, real>::_M[0][0] = static_cast <real> (1.0/6.0);
            SplineBase<V, real>::_M[0][1] = static_cast <real> (4.0/6.0);	
            SplineBase<V, real>::_M[0][2] = static_cast <real> (1/6.0);			
            SplineBase<V, real>::_M[0][3] = 0;
            
			SplineBase<V, real>::_M[1][0] = static_cast <real> (-3.0/6.0);	
            SplineBase<V, real>::_M[1][1] = static_cast <real> (0.0);		
            SplineBase<V, real>::_M[1][2] = static_cast <real> (3/6.0);			
            SplineBase<V, real>::_M[1][3] = 0;
            
            SplineBase<V, real>::_M[2][0] = static_cast <real> (3.0/6.0);	
            SplineBase<V, real>::_M[2][1] = static_cast <real> (-6.0/6.0);	
            SplineBase<V, real>::_M[2][2] = static_cast <real> (3/6.0);			
            SplineBase<V, real>::_M[2][3] = 0;
            
			SplineBase<V, real>::_M[3][0] = static_cast <real> (-1.0/6.0);	
            SplineBase<V, real>::_M[3][1] = static_cast <real> (3.0/6.0);	
            SplineBase<V, real>::_M[3][2] = static_cast <real> (-3.0/6.0);		
            SplineBase<V, real>::_M[3][3] = static_cast <real> (1.0/6.0);
		}
	public:
		CubicBSplineBase () : SplineBase<V, real> ()
		{
			init();
		}
		CubicBSplineBase (const std::vector <V> & points) : SplineBase<V, real> (points)
		{
			init();
		}
	};

	template <typename V, typename real = float>
	class CRSplineBase : public SplineBase <V, real>
	{
	private:
		// TODO: Optimize by initializing this matrix only once rather than each time a spline object is created
		void init (const real fT)
		{
            using S = SplineBase<V, real>;
            
			S::_M[0][0] = 0;     S::_M[0][1] = 1;    S::_M[0][2] = 0;		S::_M[0][3] = 0;
			S::_M[1][0] = -fT;   S::_M[1][1] = 0;    S::_M[1][2] = fT;		S::_M[1][3] = 0;
			S::_M[2][0] = 2*fT;  S::_M[2][1] = fT-3; S::_M[2][2] = 3-2*fT;	S::_M[2][3] = -fT;
			S::_M[3][0] = -fT;   S::_M[3][1] = 2-fT; S::_M[3][2] = fT-2;	S::_M[3][3] = fT;
		}

	public:
		CRSplineBase (const real fT = 0.5f) : SplineBase<V, real>()
		{
			init (fT);
		}
		CRSplineBase (const std::vector <V> & points, const real fT = 0.5f) : SplineBase<V, real>(points)
		{
			init(fT);
		}
	};

	template <typename SB, typename real>
	class ArcLengthFnPiece : public GaussLegendreInt <real, ArcLengthFnPiece <SB, real>, 5>
	{
	private:
		const SB & _link;
		const unsigned int _i;
		mutable real _cache100;

		// Copying will break the reference with SB
		ArcLengthFnPiece (const ArcLengthFnPiece & copy) {}

	public:
		ArcLengthFnPiece (const SB & link, const unsigned int i)
			: _link(link), _i(i), _cache100(-1) {}

		inline real f (const real t) const 
		{
			return MAG(_link.compute (_i, 0, 1, 2*t, 3*t*t));
		}
		inline real length (const real a, const real b) const
		{ 
#ifdef _DEBUG
			assert (b >= a);
			assert (b <= 1.0 && a >= 0.0 );
#endif
			return GaussLegendreInt<real, ArcLengthFnPiece <SB, real>, 5>::compute(a, b); 
		}

		inline real length () const
			{ return _cache100 < 0 ? (_cache100 = length(0, 1)) : _cache100; }

		inline void clearCache () const
			{ _cache100 = -1; }
	};

	template <typename SB>
	class Spline
	{
	public:
		typedef SB BaseType;
		typedef Spline< SB > SplineType;
		typedef typename SB::Precision real;
		typedef typename SB::VectorType VectorType;

		const real ONE;

	private:
		SB _sp;
		ArcLengthFnPiece<SB, real> ** _lenfns;

		// Approximates the "time" along the curve at the specified physical distance
		// This is the inverse of arclength
		class ArcPositionFn
		{
		private:
			const ArcLengthFnPiece<SB, real>	&	_fnLen;
			real									_fTest;

			ArcPositionFn (const ArcPositionFn & copy) {}

		public:
			ArcPositionFn (const ArcLengthFnPiece<SB, real> & fnLen, const real fPosition)
				: _fnLen(fnLen), _fTest(fPosition) {}

			inline const real f(const real t) const { return _fTest - _fnLen.length(0, t); }
		};

	public:
		typedef typename std::vector <VectorType>::const_iterator PointsConstIterator;
		typedef typename std::vector <VectorType>::iterator PointsIterator;
		typedef typename std::vector <VectorType>::reverse_iterator PointsRIterator;

	private:
		template <typename L>
		inline static void populateVector (const typename L::const_iterator ia, const typename L::const_iterator ib, std::vector <VectorType> & v)
			{ v.insert(v.begin(), ia, ib); }

		void init ()
		{
			const unsigned int nSegCount = segCount();

			_lenfns = new ArcLengthFnPiece< SB, real >* [nSegCount];
			for (unsigned int i = 0; i < nSegCount; ++i)
				_lenfns[i] = new ArcLengthFnPiece<SB, real> (_sp, i);
		}
		inline const ArcLengthFnPiece<SB, real> & lenfn (const unsigned int i) const
		{
#ifdef _DEBUG
			assert(i < segCount());
			assert(_lenfns != NULL);
			assert(_lenfns[i] != NULL);
#endif
			return *_lenfns[i];
		}
		void clearCache () const
		{
			for (unsigned int i = 0; i < segCount(); ++i)
				lenfn(i).clearCache();
		}

	public:
		Spline (const Spline <SB> & copy)
			: _lenfns(NULL), _sp (copy._sp), ONE(static_cast< real > (1)) { init (); }
		Spline (const std::vector <VectorType> & points)
			: _lenfns(NULL), _sp(points), ONE(static_cast< real > (1)) { init(); }
		Spline (const real fT, const std::vector <VectorType> & points)
			: _lenfns(NULL), _sp(points, fT), ONE(static_cast< real > (1)) { init(); }
		Spline ()
			: _lenfns(NULL), _sp(), ONE(static_cast< real > (1)) { init (); }
		Spline (const real fT)
			: _lenfns(NULL), _sp(fT), ONE(static_cast< real > (1)) { init (); }

		template <typename L>
		static Spline createFrom (const typename L::const_iterator ia, const typename L::const_iterator ib)
		{
			 std::vector <VectorType> v;
			 populateVector <L> (ia, ib, v);
			 return Spline(v);
		} 
		
		template <typename L>
		static Spline createFrom (const real fT, const typename L::const_iterator ia, const typename L::const_iterator ib)
		{
			std::vector <VectorType> v;
			populateVector <L> (ia, ib, v);
			return Spline(fT, v);
		} 

		// Returns the length of the specified curve segment [i, i+1]
		inline real len (const unsigned int i) const { return lenfn(i).length(); }
		// Returns the length of the specified curve segment between times tA and tB [i+tA, i+tB]
		inline real len (const unsigned int i, const real tA, const real tB) const { return lenfn(i).length(tA, tB); }
		// Returns the total length of the curve
		real len () const
		{
			const unsigned nSegs = segCount();
			real nLen = 0;

			for (unsigned i = 0; i < nSegs; ++i)
				nLen += len(i);

			return nLen;
		}
		real distanceFast (const unsigned int i, const real t) const
		{
			real fLen = 0;

			for (unsigned int c = 0; c < i; ++c)
				fLen += len(c);

			return fLen += len(i) * t;
		}
		real distance (const unsigned int i, const real t) const
		{
			real fLen = 0;

			for (unsigned int c = 0; c < i; ++c)
				fLen += len(c);

			return fLen += lenfn(i).length(0, t);
		}
		inline real timeAt (const unsigned int i, const real & fPos) const
		{
			return mars::BisectMethod< real, ArcPositionFn > :: solve ( ArcPositionFn(lenfn(i), fPos), 0, 1 );
		}
		inline const VectorType compute (const unsigned int i, const real t) const
		{
			const real 
				t2 = t * t,
				t3 = t2 * t;

			return _sp.compute(i, 1, t, t2, t3);
		}
		inline const VectorType gradient (const unsigned int i, const real t) const
			{ return _sp.compute(i, 0, 1, 2*t, 3*t*t); }
		inline const VectorType laplace (const unsigned int i, const real t) const
			{ return _sp.compute(i, 0, 0, 2, 6*t); }

		template <typename V3D>
		void move (const real radius, const V3D & unormal)
		{
			_sp.move (radius, unormal);
			clearCache();
		}
		inline unsigned int segCount() const { return _sp.size() - 1; }
		inline unsigned int size () const { return _sp.size(); }
		inline PointsConstIterator begin() const { return _sp.begin(); }
		inline PointsConstIterator end() const { return _sp.end(); }
		inline PointsIterator begin() { return _sp.begin(); }
		inline PointsIterator end() { return _sp.end(); }
		inline const VectorType & operator [] (const unsigned int i) const { return _sp[i]; }
		inline VectorType & operator [] (const unsigned int i) { return _sp[i]; }
		inline VectorType & front () { return _sp.front(); }
		inline const VectorType & front () const { return _sp.front(); }
		inline VectorType & back () { return _sp.back(); }
		inline const VectorType & back () const { return _sp.back(); }

		inline void operator += (const VectorType & ptOffset)
			{ _sp += ptOffset; }

		~Spline()
		{
			for (unsigned int i = 0; i < segCount(); ++i)
				delete _lenfns[i];

			delete [] _lenfns;
		}

	};

	template <typename S, typename Vs, typename integer = short>
	class SplineRasterIterator : public std::iterator < std::input_iterator_tag, Vs >
	{
	public:
		typedef S SplineType;
		typedef Vs RasterVecType;
		typedef integer RasterPrecision;
		typedef typename SplineType::real real;
		typedef typename SplineType::VectorType VectorType;

	private:
		const SplineType & _sp;
		Vs	_zp,				// Resulting iterative raster point (also used to compare guess points against)
			_tzp;				// Guess point (integral) (current)
		VectorType	_trp;		// Guess point (rational) (current)
		unsigned int	_i,		// Segment number (current)
			_k;		// Iteration count
		real	_accel,			// Acceleration (current) (added to current time for next guess point)
			_t,				// Time (current)
			_t0;			// Time (last)
		const unsigned int 
			_segCount;

		CR_CONTEXT;

		int calcZDiff (const real & t)
		{
			using namespace mars;

			_trp = _sp.compute(_i, t);
			_tzp = _trp;

			const Vs dzp = _tzp - _zp;

			return (SQ(dzp.x) + SQ(dzp.y));
		}
		inline void assignAccel (const real acc)
		{
			assert(acc >= 0);
			if (acc > 0)
				_accel = acc;
		}
		bool satisfy (real & acc, unsigned int & zdiff)
		{
			using namespace mars;
			const real t = _t + acc;

			// Calculate integral distance from last position to guessed position
			zdiff = calcZDiff(t);
			if (zdiff == 1)
			{
				_zp = _tzp;
				_t = t;
				assignAccel(acc);
				return true;
			} else if (zdiff <= 2)
			{
				const VectorType drp = _trp - _zp;

				if (drp.x > abs(drp.y) && drp.x > 0)	// Slope is bias to the positive x-direction
					++_zp.x;
				else if (drp.y > abs(drp.x) && drp.y > 0)	// Slope is bias to the positive y-direction
					++_zp.y;
				else if (drp.y < -abs(drp.x) && drp.y < 0)	// Slope is bias to the negative y-direction
					--_zp.y;
				else if (drp.x < -abs(drp.y) && drp.x < 0)	// Slope is bias to the negative x-direction
					--_zp.x;
				else	// Slope is perfectly 1:1 or 0:0
				{
					if (zdiff == 0)
					{
						if (drp.x == 0 || drp.y == 0)
							return false;

						if (drp.x < 0)
							--_zp.x;
						else
							++_zp.x;

						if (drp.y < 0)
							--_zp.y;
						else
							++_zp.y;
					} else
						_zp = _tzp;
					assignAccel(acc);
					_t = t;
					return true;
				}

				if (zdiff == 2)
				{
					assignAccel(acc);
					_t += _accel / 2;
				} else
				{
					assignAccel(acc);
					_t = t;
				}
				return true;
			} else
				return false;
		}
		template <int N> struct StaticIteration {};
        template <int N> void drilldown (real & acc, unsigned int & zdiff) {
            drilldown(acc, zdiff, StaticIteration<N>());
        }
		template <int N> void drilldown (real & acc, unsigned int & zdiff, StaticIteration<N>)
		{
			using namespace mars;

			// 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
			if (!satisfy (acc, zdiff))
			{
				if (_trp != _zp)
				{
					const real accbak = acc;

					acc /= MAG(_trp - _zp);
					if (acc >= 1.0)
						acc = accbak * 2;
				} else
					acc *= 2;

				drilldown <N - 1> (acc, zdiff);
			}
		}
		inline void drilldown (real & acc, unsigned int & zdiff, StaticIteration<0>) 
		{
			if (!satisfy (acc, zdiff))
			{
				drilldownBeta(acc);
			}
		}
		void drilldownBeta(real & acc)
		{
			using namespace mars;
			static const real SPLINE_PREC_EPSILON = std::numeric_limits <real>::epsilon();

			real dacc = acc;
			signed short ddd = 3;
			unsigned int zdiff;

			do
			{
				if (!satisfy(acc, zdiff))
				{
					if (zdiff == 0)
					{
						if (ddd ^ 1)
						{
							dacc /= 2;
							ddd = 1;
						}
						acc += dacc;
					} else
					{
						if (ddd ^ 0)
						{
							dacc /= 2;
							ddd = 0;
						}
						acc -= dacc;
					}
				} else
					return;

				// Acc < 1.0 prevents an index-out-of-bounds exception and checking delta-acc here ensures numerical stability
			} while (acc < 1.0 && acc > 0 && (dacc > SPLINE_PREC_EPSILON || dacc < -SPLINE_PREC_EPSILON));

			// Catch-all for acceleration values that are negative or too small, accounts for numerical instability as well as monotonic exceptions
			if (acc > SPLINE_PREC_EPSILON)
				_accel = acc;

			_zp = _tzp;
			_t += _accel;
		}

		void findAdjacentCell ()
		{
			using namespace mars;

			real 
				acc = _accel;
			unsigned int
				zdiff;

			drilldown <6> (acc, zdiff);
		}
		void process ()
		{
			CR_START();

			_t = 0;
			_zp = _sp.compute (0, 0);

			for (_i = 0; _i < _segCount; ++_i)
			{
				_accel = static_cast <real> (1.0) / _sp.len(_i);
				_t0 = _t = 0;

				while (_t < 1.0)
				{
					findAdjacentCell();
					CR_RETURN_VOID;
					_t0 = _t;
					++_k;
				}
			}

			CR_END();
		}

	public:
		SplineRasterIterator (const SplineRasterIterator & copy);
		SplineRasterIterator (const SplineType & sp)
			: _sp (sp), _segCount (sp.segCount()), _k(0)
		{
			CR_INIT();

			process();
		}

		// XTODO: Optimize!!!
		SplineRasterIterator & operator += (const real rDist)
		{
			if (rDist > 0)
			{
				const unsigned int nSegs = _sp.segCount();
				const real fDist = distanceFast();
				const real fTot = static_cast <real> (rDist) + fDist;
				real fAcc = fDist;

				while (_i < nSegs)
				{
					const real fTesterz = _sp.len(_i);

					if (fAcc + fTesterz > fTot)
					{
						_t0 = _t = _sp.timeAt (_i, fTot - fAcc);
						_zp = _sp.compute(_i, _t);
						_accel = static_cast <real> (1.0) / _sp.len(_i);
						_k += static_cast <unsigned int> (rDist);
						return *this;
					}
					else
						fAcc += fTesterz;
					++_i;
				}

				// Terminate the loop if we get here
				_t0 = _t = 1.0f;
			}

			return *this;
		}

		inline SplineRasterIterator & operator++()
		{
			process();
			return *this;
		}
		inline SplineRasterIterator operator++(int)
		{
			SplineRasterIterator old = *this;

			process();
			return old;
		}
		// Retrieves the number of pixels traversed thus far, or the current pixel in the iteration
		// NOTE: This is not the same as the distance of spline curve traversed thus far, so it cannot be
		// compared accurately with the spline length since this deals with number of pixels plotted as 
		// opposed to any physical measure of distance from the origin.
		inline unsigned int index () const { return _k; }
		inline unsigned int segment () const { return _i; }
		inline real distanceFast() const { return _sp.distanceFast(_i, _t0); }
		inline real distance() const { return _sp.distance(_i, _t0); }
		inline const VectorType gradient () const { return _sp.gradient(_i, _t0); }
		inline const VectorType laplace () const { return _sp.laplace(_i, _t0); }

		inline operator const Vs & () const { return _zp; }
		inline operator bool () const { return _i < _segCount; }
		inline const SplineType & getSpline() const { return _sp; }
	};

	template <typename S, typename EOFn, typename EIFn, typename RasterPrec, typename Vs, typename V3>
	class SplineZoneIterator : public std::iterator < std::input_iterator_tag, Vs >
	{
	protected:
		typedef S SplineType;
		typedef typename SplineType::real real;
		typedef SplineRasterIterator <S, Vs, RasterPrec> SplineWalker;
		typedef typename SplineType::VectorType VectorType;

	private:

		mutable struct TotalRadiusCache 
		{
			unsigned int result, index, radi;

			inline TotalRadiusCache ()
				: result(0), index(0), radi(0) {}

		} _crad;

		const unsigned int _netradius;	// TODO: Genericize unsigned int _netradius
		const EOFn _fnEaseOut;
		const EIFn _fnEaseIn;
		unsigned int _splnum;
		signed int _side;
		SplineWalker * _pItSP;
		unsigned int _radi;
		real _fBaseLen, _fLimit, _fOffset, _fSpLen;
		SplineType * _pSplineA, * _pSplineB;
		V3 _unormal;
		typename SplineWalker::RasterVecType _p;
		bool _done;

		CR_CONTEXT;

	public:
		SplineZoneIterator (const SplineType & spline, const unsigned int nRadius, const EOFn & fnEaseOut, const EIFn & fnEaseIn, const V3 & unormal)
			: _netradius (nRadius), _unormal (unormal), _done(false), _pItSP(NULL), _splnum(0), _fnEaseOut(fnEaseOut), _fnEaseIn(fnEaseIn), _radi(0)
		{
			using namespace mars;

			CR_INIT();
			_pSplineB = new S (spline);
			_pSplineA = new S (spline);

			process();
		}
		~SplineZoneIterator ()
		{
			delete _pSplineA;
			delete _pSplineB;
			delete _pItSP;
		}

		inline SplineZoneIterator & operator++()
		{
			process();
			return *this;
		}
		inline SplineZoneIterator operator++(int)
		{
			SplineZoneIterator old = *this;

			process();
			return old;
		}

		inline const Vs & operator * () const { return _p; }
		inline const Vs & operator -> () const { return _p; }

		inline operator const Vs & () const { return _p; }
		inline const typename SplineWalker::RasterVecType & pos () const { return _p; }

		// Returns current radius as a ratio to the total radius
		inline real trad () const { return static_cast <real> (radius()) / static_cast <real> (SplineZoneIterator::crad()); }
		// Returns the present distance from the spine spline (always positive)
		inline const unsigned int radius () const { return _radi; }
		// Returns the present distance from the spine spline (negative if on the negative side and positive if on the positive side)
		inline const signed int sradius () const { return static_cast< signed int > (radius()) * side(); }
		inline const signed int side() const { return _side; }
		// Returns the total radius of the zone at the current time
		const unsigned int crad () const 
		{ 
			if (_crad.result == 0 || _pItSP->index() != _crad.index || _radi != _crad.radi)
			{
				const real 
					ONE = _pItSP->getSpline().ONE,
					frDist = static_cast <real> (_pItSP->distanceFast()) / _fSpLen;

				// FIXME: _radi ends-up being more than the calculated total radius causing the assertion to fail (see below).  
				// This may have been fixed as of the change to use accurate spline length calculation instead of approximate
				_crad.result = 
					std::max(
						_radi,
						static_cast <unsigned int> (
							_fnEaseOut.invf(frDist) * _fnEaseIn.invf(ONE - frDist) * static_cast <real> (_netradius)
						) + 1
					)
				;
				_crad.index = _pItSP->index();
				_crad.radi = _radi;

				assert(_crad.radi <= _crad.result);
			}

			return _crad.result;
		}
		inline operator bool () const { return !_done; }
		inline const unsigned int splineNum () const { return _splnum; }
		inline const unsigned int index () const { return _pItSP->index(); }
		inline const real distance () const { return _pItSP->distance(); }
		inline const real distanceFast () const { return _pItSP->distanceFast(); }
		inline const VectorType gradient () const { return _pItSP->gradient(); }
		inline const VectorType laplace () const { return _pItSP->laplace(); }
		inline const real seglen () const { return _pItSP->getSpline().len(_pItSP->segment()); }
		inline const VectorType & tailCP () const { return _pItSP->getSpline()[_pItSP->segment()]; }
		inline const VectorType & leadCP () const { return _pItSP->getSpline()[_pItSP->segment() + 1]; }

	private:
		void process ()
		{
			using namespace mars;

			CR_START();

			_fSpLen = 
				_fBaseLen = _pSplineB->len();
			for (_pItSP = new SplineWalker (*_pSplineB); *_pItSP; ++*_pItSP)
			{
				_p = *_pItSP;
				CR_RETURN_VOID;
			}
			delete _pItSP;
			_pItSP = NULL;

			// Fill to the positive side of the curve
			_side = +1;
			while (_radi++ < _netradius)
			{
				prepareWalk(_pSplineA);
				for (*(_pItSP = new SplineWalker (*_pSplineA)) += _fOffset; 
					*_pItSP && _pItSP->distanceFast() < _fLimit;
					++*_pItSP)
				{
					_p = *_pItSP;
					CR_RETURN_VOID;
				}
				postWalk(_pSplineA);
			}
			_pItSP = NULL;

			// Fill to the negative side of the curve
			_radi = 0;
			_side = -1;
			while (_radi++ < _netradius)
			{
				prepareWalk(_pSplineB);

				for (*(_pItSP = new SplineWalker (*_pSplineB)) += _fOffset; 
					*_pItSP && _pItSP->distanceFast() < _fLimit; 
					++*_pItSP)
				{
					_p = *_pItSP;
					CR_RETURN_VOID;
				}
				postWalk(_pSplineB);
			}
			_pItSP = NULL;

			_done = true;

			CR_END();
		}
		inline void postWalk(const S * pSpline)
		{
			delete _pItSP;
		}
		void prepareWalk (S * pSpline)
		{
			pSpline->move (static_cast< float > (_side), _unormal);
			_fSpLen = pSpline->len();
			const real 
				ONE = pSpline->ONE,
				frRadii = static_cast <real> (_radi) / static_cast <real> (_netradius);

			_fLimit = _fnEaseOut.f( frRadii / _fnEaseIn.invf(ONE - _fnEaseIn.f(frRadii) ) ) * _fSpLen;
			_fOffset = (ONE - _fnEaseIn.f( frRadii / _fnEaseOut.invf(ONE - _fnEaseOut.f(frRadii ) ) )) * _fSpLen;

			// FIXME: Limit and offset calculations are inaccurate, causes _radi to increment out-of-range.
			++_splnum;
		}
	};
}