// K-3D
// Copyright (c) 1995-2004, Timothy M. Shead
//
// Contact: tshead@k-3d.com
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public
// License along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

/** \file
		\author Romain Behar (romainbehar@yahoo.com)
*/

#include <k3dsdk/classes.h>
#include <k3dsdk/color.h>
#include <k3dsdk/imaterial.h>
#include <k3dsdk/irenderman.h>
#include <k3dsdk/material.h>
#include <k3dsdk/material_collection.h>
#include <k3dsdk/measurement.h>
#include <k3dsdk/mesh_source.h>
#include <k3dsdk/module.h>
#include <k3dsdk/object.h>
#include <k3dsdk/persistence.h>
#include <k3dsdk/plugins.h>
#include <k3dsdk/property.h>
#include <k3dsdk/transform.h>

namespace k3d
{

inline bool operator < (const k3d::vector3& a, const k3d::vector3& b)
{
	if(a[0] < b[0])
		return true;
	if(a[0] > b[0])
		return false;
	if(a[1] < b[1])
		return true;
	if(a[1] > b[1])
		return false;
	if(a[2] < b[2])
		return true;

	return false;
}

} // namespace k3d

namespace libk3dmesh
{

namespace detail
{

// This file is part of a program called sphereEversion:
// http://www.dgp.utoronto.ca/~mjmcguff/eversion/
//
// sphereEversion was created by Michael McGuffin.
// The code in this file was almost entirely taken
// (with slight adaptations) from the source code of
// a program called "evert", which was written by
// Nathaniel Thurston.
// evert's source code can be down loaded from
// http://www.geom.umn.edu/docs/outreach/oi/software.html
//
// Michael McGuffin extends many thanks to Nathaniel Thurston,
// Silvio Levy, and the Geometry Center (University of Minnesota)
// for making evert's source code freely available to the public.

// Axis convention: right handed system
class eversion_point
{
public:
	// Location (x, y, z) and normal vector(nx, ny, nz)
	double vertex[3];
	double normal[3];
};

typedef std::vector<eversion_point> geometry_vector_t;
typedef std::vector<geometry_vector_t> geometry_matrix_t;


class TwoJet D(const class ThreeJet x, int index);
class TwoJet
{
public:
	double f;
	double fu, fv;
	double fuv;

	TwoJet() {}
	TwoJet(double d, double du, double dv)
	{ f = d; fu = du; fv = dv; fuv = 0; }
	TwoJet(double d, double du, double dv, double duv)
	{ f = d; fu = du; fv = dv; fuv = duv; }

	bool operator < (double d) { return f < d; }
	bool operator > (double d) { return f > d; }
	bool operator <= (double d) { return f <= d; }
	bool operator >= (double d) { return f >= d; }
	double df_du() { return fu; }
	double df_dv() { return fv; }
	double d2f_dudv() { return fuv; }
	void operator += (TwoJet x)
	{ f += x.f; fu += x.fu; fv += x.fv; fuv += x.fuv; }
	void operator += (double d)
	{ f += d; }
	void operator *= (TwoJet x)
	{
		fuv = f*x.fuv + fu*x.fv + fv*x.fu + fuv*x.f;
		fu = f*x.fu + fu*x.f;
		fv = f*x.fv + fv*x.f;
		f *= x.f;
	}
	void operator *=(double d)
	{ f *= d; fu *= d; fv *= d; fuv *= d; }
	void operator %=(double d)
	{ f = fmod(f, d); if(f < 0) f += d; }
	void operator ^=(double n)
	{
		if(f > 0)
			{
				double x0 = pow(f, n);
				double x1 = n * x0/f;
				double x2 = (n-1)*x1/f;
				fuv = x1*fuv + x2*fu*fv;
				fu = x1*fu;
				fv = x1*fv;
				f = x0;
			}
	}
	void Annihilate(int index)
	{
		if(index == 0)
			fu = 0;
		else if(index == 1)
			fv = 0;

		fuv = 0;
	}
	void TakeSin()
	{
		*this *= 2*M_PI;
		double s = sin(f), c = cos(f);
		f = s; fu = fu*c; fv = fv*c; fuv = c*fuv - s*fu*fv;
	}
	void TakeCos()
	{
		*this *= 2*M_PI;
		double s = cos(f), c = -sin(f);
		f = s; fu = fu*c; fv = fv*c; fuv = c*fuv - s*fu*fv;
	}

	friend TwoJet operator+(const TwoJet x, const TwoJet y);
	friend TwoJet operator*(const TwoJet x, const TwoJet y);
	friend TwoJet operator+(const TwoJet x, double d);
	friend TwoJet operator*(const TwoJet x, double d);
	friend TwoJet Sin(const TwoJet x);
	friend TwoJet Cos(const TwoJet x);
	friend TwoJet operator^(const TwoJet x, double n);
	friend TwoJet Annihilate(const TwoJet x, int index);
	friend TwoJet Interpolate(const TwoJet v1, const TwoJet v2, const TwoJet weight);
	friend class TwoJet D(const class ThreeJet x, int index);
	friend class ThreeJet;
};

// ----------------------------------------

TwoJet operator+(const TwoJet x, const TwoJet y)
{ return TwoJet(x.f+y.f, x.fu+y.fu, x.fv+y.fv, x.fuv + y.fuv); }

TwoJet operator*(const TwoJet x, const TwoJet y)
{
	return TwoJet(
		x.f*y.f,
		x.f*y.fu + x.fu*y.f,
		x.f*y.fv + x.fv*y.f,
		x.f*y.fuv + x.fu*y.fv + x.fv*y.fu + x.fuv*y.f
	);
}

TwoJet operator+(const TwoJet x, double d)
{ return TwoJet(x.f + d, x.fu, x.fv, x.fuv); }

TwoJet operator*(const TwoJet x, double d)
{ return TwoJet(d*x.f, d*x.fu, d*x.fv, d*x.fuv); }

TwoJet Sin(const TwoJet x)
{
	TwoJet t = x*(2*M_PI);
	double s = sin(t.f);
	double c = cos(t.f);
	return TwoJet(s, c*t.fu, c*t.fv, c*t.fuv - s*t.fu*t.fv);
}

TwoJet Cos(const TwoJet x)
{
	TwoJet t = x*(2*M_PI);
	double s = cos(t.f);
	double c = -sin(t.f);
	return TwoJet(s, c*t.fu, c*t.fv, c*t.fuv - s*t.fu*t.fv);
}

TwoJet operator^(const TwoJet x, double n)
{
	double x0 = pow(x.f, n);
	double x1 = (x.f == 0) ? 0 : n * x0/x.f;
	double x2 = (x.f == 0) ? 0 : (n-1)*x1/x.f;
	return TwoJet(x0, x1*x.fu, x1*x.fv, x1*x.fuv + x2*x.fu*x.fv);
}

TwoJet Annihilate(const TwoJet x, int index)
{ return TwoJet(x.f, index == 1 ? x.fu : 0, index == 0 ? x.fv : 0, 0); }

TwoJet Interpolate(const TwoJet v1, const TwoJet v2, const TwoJet weight)
{ return (v1) * ((weight) * (-1) + 1) + v2*weight; }

// ----------------------------------------

class ThreeJet
{
public:
	double f;

private:
	double fu, fv;
	double fuu, fuv, fvv;
	double fuuv, fuvv;

	ThreeJet(double d, double du, double dv, double duu, double duv, double dvv,
		double duuv, double duvv)
	{ f = d; fu = du; fv = dv; fuu = duu; fuv = duv; fvv = dvv;
		fuuv = duuv; fuvv = duvv; }

public:
	ThreeJet() {}
	ThreeJet(double d, double du, double dv)
	{ f = d; fu = du; fv = dv; fuu = fuv = fvv = fuuv = fuvv = 0;}
	operator TwoJet() { return TwoJet(f, fu, fv, fuv); }

	bool operator<(double d) { return f < d; }
	bool operator>(double d) { return f > d; }
	bool operator<=(double d) { return f <= d; }
	bool operator>=(double d) { return f >= d; }
	void operator%=(double d) { f = fmod(f, d); if(f < 0) f += d; }
	friend ThreeJet operator+(const ThreeJet x, const ThreeJet y);
	friend ThreeJet operator*(const ThreeJet x, const ThreeJet y);
	friend ThreeJet operator+(const ThreeJet x, double d);
	friend ThreeJet operator*(const ThreeJet x, double d);
	friend ThreeJet Sin(const ThreeJet x);
	friend ThreeJet Cos(const ThreeJet x);
	friend ThreeJet operator^(const ThreeJet x, double n);
	friend ThreeJet Annihilate(const ThreeJet x, int index);
	friend ThreeJet Interpolate(const ThreeJet v1, const ThreeJet v2, const ThreeJet weight);
	friend class TwoJet D(const class ThreeJet x, int index);
};

// ----------------------------------------

ThreeJet operator+(const ThreeJet x, const ThreeJet y)
{
	ThreeJet result;
	result.f = x.f + y.f;
	result.fu = x.fu + y.fu;
	result.fv = x.fv + y.fv;
	result.fuu = x.fuu + y.fuu;
	result.fuv = x.fuv + y.fuv;
	result.fvv = x.fvv + y.fvv;
	result.fuuv = x.fuuv + y.fuuv;
	result.fuvv = x.fuvv + y.fuvv;

	return result;
}

ThreeJet operator*(const ThreeJet x, const ThreeJet y)
{
	ThreeJet result;
	result.f = x.f*y.f;
	result.fu = x.f*y.fu + x.fu*y.f;
	result.fv = x.f*y.fv + x.fv*y.f;
	result.fuu = x.f*y.fuu + 2*x.fu*y.fu + x.fuu*y.f;
	result.fuv = x.f*y.fuv + x.fu*y.fv + x.fv*y.fu + x.fuv*y.f;
	result.fvv = x.f*y.fvv + 2*x.fv*y.fv + x.fvv*y.f;
	result.fuuv = x.f*y.fuuv + 2*x.fu*y.fuv + x.fv*y.fuu + 2*x.fuv*y.fu + x.fuu*y.fv + x.fuuv*y.f;
	result.fuvv = x.f*y.fuvv + 2*x.fv*y.fuv + x.fu*y.fvv + 2*x.fuv*y.fv + x.fvv*y.fu + x.fuvv*y.f;
	return result;
}

ThreeJet operator+(const ThreeJet x, double d)
{
	ThreeJet result;
	result = x;
	result.f += d;
	return result;
}

ThreeJet operator*(const ThreeJet x, double d)
{
	ThreeJet result;
	result.f = d*x.f;
	result.fu = d*x.fu;
	result.fv = d*x.fv;
	result.fuu = d*x.fuu;
	result.fuv = d*x.fuv;
	result.fvv = d*x.fvv;
	result.fuuv = d*x.fuuv;
	result.fuvv = d*x.fuvv;
	return result;
}

ThreeJet Sin(const ThreeJet x)
{
	ThreeJet result;
	ThreeJet t = x*(2*M_PI);
	double s = sin(t.f);
	double c = cos(t.f);
	result.f = s;
	result.fu = c*t.fu;
	result.fv = c*t.fv;
	result.fuu = c*t.fuu - s*t.fu*t.fu;
	result.fuv = c*t.fuv - s*t.fu*t.fv;
	result.fvv = c*t.fvv - s*t.fv*t.fv;
	result.fuuv = c*t.fuuv - s*(2*t.fu*t.fuv + t.fv*t.fuu) - c*t.fu*t.fu*t.fv;
	result.fuvv = c*t.fuvv - s*(2*t.fv*t.fuv + t.fu*t.fvv) - c*t.fu*t.fv*t.fv;
	return result;
}

ThreeJet Cos(const ThreeJet x)
{
	ThreeJet result;
	ThreeJet t = x*(2*M_PI);
	double s = cos(t.f);
	double c = -sin(t.f);
	result.f = s;
	result.fu = c*t.fu;
	result.fv = c*t.fv;
	result.fuu = c*t.fuu - s*t.fu*t.fu;
	result.fuv = c*t.fuv - s*t.fu*t.fv;
	result.fvv = c*t.fvv - s*t.fv*t.fv;
	result.fuuv = c*t.fuuv - s*(2*t.fu*t.fuv + t.fv*t.fuu) - c*t.fu*t.fu*t.fv;
	result.fuvv = c*t.fuvv - s*(2*t.fv*t.fuv + t.fu*t.fvv) - c*t.fu*t.fv*t.fv;
	return result;
}

ThreeJet operator^(const ThreeJet x, double n)
{
	double x0 = pow(x.f, n);
	double x1 = (x.f == 0) ? 0 : n * x0/x.f;
	double x2 = (x.f == 0) ? 0 : (n-1) * x1/x.f;
	double x3 = (x.f == 0) ? 0 : (n-2) * x2/x.f;
	ThreeJet result;
	result.f = x0;
	result.fu = x1*x.fu;
	result.fv = x1*x.fv;
	result.fuu = x1*x.fuu + x2*x.fu*x.fu;
	result.fuv = x1*x.fuv + x2*x.fu*x.fv;
	result.fvv = x1*x.fvv + x2*x.fv*x.fv;
	result.fuuv = x1*x.fuuv + x2*(2*x.fu*x.fuv + x.fv*x.fuu) + x3*x.fu*x.fu*x.fv;
	result.fuvv = x1*x.fuvv + x2*(2*x.fv*x.fuv + x.fu*x.fvv) + x3*x.fu*x.fv*x.fv;
	return result;
}

TwoJet D(const ThreeJet x, int index)
{
	TwoJet result;
	if(index == 0)
		{
			result.f = x.fu;
			result.fu = x.fuu;
			result.fv = x.fuv;
			result.fuv = x.fuuv;
		}
	else if(index == 1)
		{
			result.f = x.fv;
			result.fu = x.fuv;
			result.fv = x.fvv;
			result.fuv = x.fuvv;
		}
	else
		{
			result.f = result.fu = result.fv =
			result.fuv = 0;
		}

	return result;
}

ThreeJet Annihilate(const ThreeJet x, int index)
{
	ThreeJet result = ThreeJet(x.f,0,0);
	if(index == 0)
		{
			result.fv = x.fv;
			result.fvv = x.fvv;
		}
	else if(index == 1)
		{
			result.fu = x.fu;
			result.fuu = x.fuu;
		}

	return result;
}

ThreeJet Interpolate(const ThreeJet v1, const ThreeJet v2, const ThreeJet weight)
{ return (v1) * ((weight) * (-1) + 1) + v2*weight; }

// ----------------------------------------

struct TwoJetVec
{
	TwoJet x;
	TwoJet y;
	TwoJet z;
	TwoJetVec() {}
	TwoJetVec(TwoJet a, TwoJet b, TwoJet c) { x = a; y = b; z = c; }
};

TwoJetVec operator+(TwoJetVec v, TwoJetVec w);
TwoJetVec operator*(TwoJetVec v, TwoJet a);
TwoJetVec operator*(TwoJetVec v, double a);
TwoJetVec AnnihilateVec(TwoJetVec v, int index);
TwoJetVec Cross(TwoJetVec v, TwoJetVec w);
TwoJet Dot(TwoJetVec v, TwoJetVec w);
TwoJetVec Normalize(TwoJetVec v);
TwoJetVec RotateZ(TwoJetVec v, TwoJet angle);
TwoJetVec RotateY(TwoJetVec v, TwoJet angle);
TwoJetVec RotateX(TwoJetVec v, TwoJet angle);
TwoJetVec InterpolateVec(TwoJetVec v1, TwoJetVec v2, TwoJet weight);
TwoJet Length(TwoJetVec v);

// ----------------------------------------

TwoJetVec operator+(TwoJetVec v, TwoJetVec w)
{
	TwoJetVec result;
	result.x = v.x + w.x;
	result.y = v.y + w.y;
	result.z = v.z + w.z;
	return result;
}

TwoJetVec operator*(TwoJetVec v, TwoJet a)
{
	TwoJetVec result;
	result.x = v.x*a;
	result.y = v.y*a;
	result.z = v.z*a;
	return result;
}

TwoJetVec operator*(TwoJetVec v, double a)
{
	TwoJetVec result;
	result.x = v.x*a;
	result.y = v.y*a;
	result.z = v.z*a;
	return result;
}

TwoJetVec AnnihilateVec(TwoJetVec v, int index)
{
	TwoJetVec result;
	result.x = Annihilate(v.x, index);
	result.y = Annihilate(v.y, index);
	result.z = Annihilate(v.z, index);
	return result;
}

TwoJetVec Cross(TwoJetVec v, TwoJetVec w)
{
	TwoJetVec result;
	result.x = v.y*w.z + v.z*w.y*-1;
	result.y = v.z*w.x + v.x*w.z*-1;
	result.z = v.x*w.y + v.y*w.x*-1;
	return result;
}

TwoJet Dot(TwoJetVec v, TwoJetVec w)
{ return v.x*w.x + v.y*w.y + v.z*w.z; }

TwoJetVec Normalize(TwoJetVec v)
{
	TwoJet a;
	a = Dot(v,v);
	if(a > 0)
		a = a^-0.5;
	else
		a = TwoJet(0, 0, 0);

	return v*a;
}

TwoJetVec RotateZ(TwoJetVec v, TwoJet angle)
{
	TwoJetVec result;
	TwoJet s,c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x*c + v.y*s;
	result.y = v.x*s*-1 + v.y*c;
	result.z = v.z;
	return result;
}

TwoJetVec RotateY(TwoJetVec v, TwoJet angle)
{
	TwoJetVec result;
	TwoJet s, c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x*c + v.z*s*-1;
	result.y = v.y;
	result.z = v.x*s + v.z*c;
	return result;
}

TwoJetVec RotateX(TwoJetVec v, TwoJet angle)
{
	TwoJetVec result;
	TwoJet s,c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x;
	result.y = v.y*c + v.z*s;
	result.z = v.y*s*-1 + v.z*c;
	return result;
}

TwoJetVec InterpolateVec(TwoJetVec v1, TwoJetVec v2, TwoJet weight)
{
	return (v1) * (weight*-1 + 1) + v2*weight;
}

TwoJet Length(TwoJetVec v)
{
	return (TwoJet(v.x^2) + TwoJet(v.y^2)) ^ (.5);
}

// ----------------------------------------

struct ThreeJetVec
{
	ThreeJet x;
	ThreeJet y;
	ThreeJet z;
	operator TwoJetVec() { return TwoJetVec(x,y,z); }
};

ThreeJetVec operator+(ThreeJetVec v, ThreeJetVec w);
ThreeJetVec operator*(ThreeJetVec v, ThreeJet a);
ThreeJetVec operator*(ThreeJetVec v, double a);
ThreeJetVec AnnihilateVec(ThreeJetVec v, int index);
ThreeJetVec Cross(ThreeJetVec v, ThreeJetVec w);
ThreeJet Dot(ThreeJetVec v, ThreeJetVec w);
TwoJetVec D(ThreeJetVec x, int index);
ThreeJetVec Normalize(ThreeJetVec v);
ThreeJetVec RotateZ(ThreeJetVec v, ThreeJet angle);
ThreeJetVec RotateY(ThreeJetVec v, ThreeJet angle);
ThreeJetVec RotateX(ThreeJetVec v, ThreeJet angle);
ThreeJetVec InterpolateVec(ThreeJetVec v1, ThreeJetVec v2, ThreeJet weight);
ThreeJet Length(ThreeJetVec v);

// ----------------------------------------

ThreeJetVec operator+(ThreeJetVec v, ThreeJetVec w)
{
	ThreeJetVec result;
	result.x = v.x + w.x;
	result.y = v.y + w.y;
	result.z = v.z + w.z;
	return result;
}

ThreeJetVec operator*(ThreeJetVec v, ThreeJet a)
{
	ThreeJetVec result;
	result.x = v.x*a;
	result.y = v.y*a;
	result.z = v.z*a;
	return result;
}

ThreeJetVec operator*(ThreeJetVec v, double a)
{
	ThreeJetVec result;
	result.x = v.x*a;
	result.y = v.y*a;
	result.z = v.z*a;
	return result;
}

ThreeJetVec AnnihilateVec(ThreeJetVec v, int index)
{
	ThreeJetVec result;
	result.x = Annihilate(v.x, index);
	result.y = Annihilate(v.y, index);
	result.z = Annihilate(v.z, index);
	return result;
}

TwoJetVec D(ThreeJetVec x, int index)
{
	TwoJetVec result;
	result.x = D(x.x, index);
	result.y = D(x.y, index);
	result.z = D(x.z, index);
	return result;
}

ThreeJetVec Cross(ThreeJetVec v, ThreeJetVec w)
{
	ThreeJetVec result;
	result.x = v.y*w.z + v.z*w.y*-1;
	result.y = v.z*w.x + v.x*w.z*-1;
	result.z = v.x*w.y + v.y*w.x*-1;
	return result;
}

ThreeJet Dot(ThreeJetVec v, ThreeJetVec w)
{ return v.x*w.x + v.y*w.y + v.z*w.z; }

ThreeJetVec Normalize(ThreeJetVec v)
{
	ThreeJet a;
	a = Dot(v,v);
	if(a > 0)
		a = a^-0.5;
	else
		a = ThreeJet(0, 0, 0);

	return v*a;
}

ThreeJetVec RotateZ(ThreeJetVec v, ThreeJet angle)
{
	ThreeJetVec result;
	ThreeJet s,c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x*c + v.y*s;
	result.y = v.x*s*-1 + v.y*c;
	result.z = v.z;
	return result;
}

ThreeJetVec RotateY(ThreeJetVec v, ThreeJet angle)
{
	ThreeJetVec result;
	ThreeJet s, c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x*c + v.z*s*-1;
	result.y = v.y;
	result.z = v.x*s + v.z*c;
	return result;
}

ThreeJetVec RotateX(ThreeJetVec v, ThreeJet angle)
{
	ThreeJetVec result;
	ThreeJet s,c;
	s = Sin (angle);
	c = Cos (angle);
	result.x = v.x;
	result.y = v.y*c + v.z*s;
	result.z = v.y*s*-1 + v.z*c;
	return result;
}

ThreeJetVec InterpolateVec(ThreeJetVec v1, ThreeJetVec v2, ThreeJet weight)
{
	return (v1) * (weight*-1 + 1) + v2*weight;
}

ThreeJet Length(ThreeJetVec v)
{
	return (ThreeJet(v.x^2) + ThreeJet(v.y^2)) ^ (.5);
}

// ----------------------------------------

TwoJetVec FigureEight(TwoJetVec w, TwoJetVec h, TwoJetVec bend, TwoJet form, TwoJet v)
{
	TwoJet height;
	v %= 1;
	height = (Cos (v*2) + -1) * (-1);
	if(v > 0.25 && v < 0.75)
		height = height*-1 + 4;

	height = height*0.6;
	h = h + bend*(height*height*(1/64.0));

	return w*Sin (v*2) + (h) * (Interpolate((Cos (v) + -1) * (-2), height, form)) ;
}

TwoJetVec AddFigureEight(ThreeJetVec p, ThreeJet u, TwoJet v, ThreeJet form, ThreeJet scale, int numStrips)
{
	ThreeJet size = form * scale;
	form = form*2 + form*form*-1;
	TwoJetVec dv = AnnihilateVec(D(p, 1), 1);
	p = AnnihilateVec(p, 1);
	TwoJetVec du = Normalize(D(p, 0));
	TwoJetVec h = Normalize(Cross(du, dv))*TwoJet(size);
	TwoJetVec w = Normalize(Cross(h, du))*(TwoJet(size)*1.1);

	return RotateZ(
	TwoJetVec(p) +
	FigureEight(w, h, du*D(size, 0)*(D(u, 0)^(-1)), form, v),
	v*(1.0/numStrips)
	);
}

// ----------------------------------------

ThreeJetVec Arc(ThreeJet u, ThreeJet v, double xsize, double ysize, double zsize)
{
	ThreeJetVec result;
	u = u*0.25;
	result.x = Sin (u) * Sin (v) * xsize;
	result.y = Sin (u) * Cos (v) * ysize;
	result.z = Cos (u) * zsize;
	return result;
}

ThreeJetVec Straight(ThreeJet u, ThreeJet v, double xsize, double ysize, double zsize)
{
	ThreeJetVec result;
	u = u*0.25;

	result.x = Sin (v) * xsize;
	result.y = Cos (v) * ysize;
	result.z = Cos (u) * zsize;
	return result;
}

ThreeJet Param1(ThreeJet x)
{
	double offset = 0;
	x %= 4;
	if(x > 2)
	{
		x = x+(-2);
		offset = 2;
	}

	if(x <= 1)
		return x*2 + (x^2)*(-1) + offset;
	else
		return (x^2) + x*(-2) + (2 + offset);
}

ThreeJet Param2(ThreeJet x)
{
	double offset = 0;
	x %= 4;
	if(x > 2)
	{
		x = x+(-2);
		offset = 2;
	}

	if(x <= 1)
		return (x^2) + offset;
	else
		return (x^2)*(-1) + x*4 + (-2 + offset);
}

static inline ThreeJet TInterp(double x)
{ return ThreeJet(x,0,0); }

ThreeJet UInterp(ThreeJet x)
{
	x %= 2;
	if(x > 1)
		x = x*(-1) + 2;
	return (x^2)*3 + (x^3) * (-2);
}

const double FFPOW = 3;
ThreeJet FFInterp(ThreeJet x)
{
	x %= 2;
	if(x > 1)
		x = x*(-1) + 2;
	x = x*1.06 + -0.05;
	if(x < 0)
		return ThreeJet(0, 0, 0);
	else if(x > 1)
		return ThreeJet(0, 0, 0) + 1;
	else
		return (x ^ (FFPOW-1)) * (FFPOW) + (x^FFPOW) * (-FFPOW+1);
}

const double FSPOW = 3;
ThreeJet FSInterp(ThreeJet x)
{
	x %= 2;
	if(x > 1)
		x = x*(-1) + 2;

	return ((x ^ (FSPOW-1)) * (FSPOW) + (x^FSPOW) * (-FSPOW+1)) * (-0.2);
}

ThreeJetVec Stage0(ThreeJet u, ThreeJet v)
{ return Straight(u, v, 1, 1, 1); }

ThreeJetVec Stage1(ThreeJet u, ThreeJet v)
{ return Arc(u, v, 1, 1, 1); }

ThreeJetVec Stage2(ThreeJet u, ThreeJet v)
{
	return InterpolateVec(
		Arc(Param1(u), v, 0.9, 0.9, -1),
		Arc(Param2(u), v, 1, 1, 0.5),
		UInterp(u)
	);
}

ThreeJetVec Stage3(ThreeJet u, ThreeJet v)
{
	return InterpolateVec(
		Arc(Param1(u), v,-0.9,-0.9,-1),
		Arc(Param2(u), v,-1, 1,-0.5),
		UInterp(u)
	);
}

ThreeJetVec Stage4(ThreeJet u, ThreeJet v)
{ return Arc(u, v, -1,-1, -1); }

ThreeJetVec Scene01(ThreeJet u, ThreeJet v, double t)
{ return InterpolateVec(Stage0(u,v), Stage1(u,v), TInterp(t)); }

ThreeJetVec Scene12(ThreeJet u, ThreeJet v, double t)
{ return InterpolateVec(Stage1(u,v), Stage2(u,v), TInterp(t)); }

ThreeJetVec Scene23(ThreeJet u, ThreeJet v, double t)
{
	ThreeJet tmp = TInterp(t);
	t = tmp.f * 0.5;
	double tt = (u <= 1) ? t : -t;

	return InterpolateVec(
		RotateZ(Arc(Param1(u), v, 0.9, 0.9,-1), ThreeJet(tt,0,0)),
		RotateY(Arc(Param2(u), v, 1, 1, 0.5), ThreeJet(t,0,0)),
		UInterp(u)
	);
}

ThreeJetVec Scene34(ThreeJet u, ThreeJet v, double t)
{ return InterpolateVec(Stage3(u,v), Stage4(u,v), TInterp(t)); }

TwoJetVec BendIn(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	ThreeJet tmp = TInterp(t);
	t = tmp.f;

	return AddFigureEight(
		Scene01(u, ThreeJet(0, 0, 1), t),
		u, v, ThreeJet(0, 0, 0), FSInterp(u),
		numStrips
	);
}

TwoJetVec Corrugate(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	ThreeJet tmp = TInterp(t);
	t = tmp.f;

	return AddFigureEight(
		Stage1(u, ThreeJet(0, 0, 1)),
		u, v, FFInterp(u) * ThreeJet(t,0,0), FSInterp(u),
		numStrips
	);
}

TwoJetVec PushThrough(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	return AddFigureEight(
		Scene12(u,ThreeJet(0, 0, 1),t),
		u, v, FFInterp(u), FSInterp(u),
		numStrips
	);
}

TwoJetVec Twist(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	return AddFigureEight(
		Scene23(u,ThreeJet(0, 0, 1),t),
		u, v, FFInterp(u), FSInterp(u),
		numStrips
	);
}

TwoJetVec UnPush(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	return AddFigureEight(
		Scene34(u,ThreeJet(0, 0, 1),t),
		u, v, FFInterp(u), FSInterp(u),
		numStrips
	);
}

TwoJetVec UnCorrugate(ThreeJet u, ThreeJet v, double t, int numStrips)
{
	ThreeJet tmp;
	tmp = TInterp((t) * (-1) + 1);
	t = tmp.f;

	return AddFigureEight(
		Stage4(u,ThreeJet(0, 0, 1)),
		u, v, FFInterp(u) * ThreeJet(t,0,0), FSInterp(u),
		numStrips
	);
}

// ----------------------------------------

typedef std::map<k3d::vector3, k3d::point*> point_map_t;
point_map_t point_map;
k3d::imaterial* current_material;
double current_radius;

k3d::point* get_spline_point(k3d::mesh* mesh, TwoJetVec p, double ps, double pus, double pvs, double puvs)
{
	k3d::vector3 v;
	v[0] = (p.x.f)*ps + (p.x.df_du()*pus + p.x.df_dv()*pvs)/3.0 + p.x.d2f_dudv()*puvs/9.0;
	v[1] = (p.y.f)*ps + (p.y.df_du()*pus + p.y.df_dv()*pvs)/3.0 + p.y.d2f_dudv()*puvs/9.0;
	v[2] = (p.z.f)*ps + (p.z.df_du()*pus + p.z.df_dv()*pvs)/3.0 + p.z.d2f_dudv()*puvs/9.0;
	// Handedness change
	v[0] = -v[0];

	k3d::point* point;
	point_map_t::const_iterator mapped_point = point_map.find(v);
	if(mapped_point != point_map.end())
		{
			point = mapped_point->second;
		}
	else
		{
			point = new k3d::point(v * current_radius);
			point_map[v] = point;
			mesh->points.push_back(point);
		}

	return point;
}

void get_patch(k3d::mesh* mesh, TwoJetVec v00, TwoJetVec v01, TwoJetVec v10, TwoJetVec v11, double us, double vs, double s0, double s1, double t0, double t1)
{
	k3d::bicubic_patch* const patch = new k3d::bicubic_patch();

	patch->control_points[0] = get_spline_point(mesh, v00, 1, 0, 0, 0);
	patch->control_points[1] = get_spline_point(mesh, v00, 1, us, 0, 0);
	patch->control_points[2] = get_spline_point(mesh, v10, 1,-us, 0, 0);
	patch->control_points[3] = get_spline_point(mesh, v10, 1, 0, 0, 0);

	patch->control_points[4] = get_spline_point(mesh, v00, 1, 0, vs, 0);
	patch->control_points[5] = get_spline_point(mesh, v00, 1, us, vs, us*vs);
	patch->control_points[6] = get_spline_point(mesh, v10, 1,-us, vs,-us*vs);
	patch->control_points[7] = get_spline_point(mesh, v10, 1, 0, vs, 0);

	patch->control_points[8] = get_spline_point(mesh, v01, 1, 0,-vs, 0);
	patch->control_points[9] = get_spline_point(mesh, v01, 1, us,-vs,-us*vs);
	patch->control_points[10] = get_spline_point(mesh, v11, 1,-us,-vs, us*vs);
	patch->control_points[11] = get_spline_point(mesh, v11, 1, 0,-vs, 0);

	patch->control_points[12] = get_spline_point(mesh, v01, 1, 0, 0, 0);
	patch->control_points[13] = get_spline_point(mesh, v01, 1, us, 0, 0);
	patch->control_points[14] = get_spline_point(mesh, v11, 1,-us, 0, 0);
	patch->control_points[15] = get_spline_point(mesh, v11, 1, 0, 0, 0);

	// Texture coordinates?
	//s0,t0, s1,t0, s0,t1, s1,t1;

	mesh->bicubic_patches.push_back(patch);
	patch->material = current_material;
}

void get_polygon_point(TwoJetVec p, eversion_point& point)
{
	double x = p.x.f ;
	double y = p.y.f ;
	double z = p.z.f ;
	double nx = p.y.df_du() * p.z.df_dv() - p.z.df_du() * p.y.df_dv();
	double ny = p.z.df_du() * p.x.df_dv() - p.x.df_du() * p.z.df_dv();
	double nz = p.x.df_du() * p.y.df_dv() - p.y.df_du() * p.x.df_dv();

	double s = nx*nx + ny*ny + nz*nz;
	if(s > 0)
		s = sqrt(1/s);

	point.vertex[0] = x;
	point.vertex[1] = y;
	point.vertex[2] = z;
	point.normal[0] = -nx*s;
	point.normal[1] = -ny*s;
	point.normal[2] = -nz*s;
}

// ----------------------------------------

typedef TwoJetVec SurfaceTimeFunction(ThreeJet u, ThreeJet v, double t, int numStrips);

inline double sqr(double x)
{
	return x*x;
}
inline double calcSpeedV(TwoJetVec v)
{
	return sqrt(sqr(v.x.df_dv()) + sqr(v.y.df_dv()) + sqr(v.z.df_dv()));
}
inline double calcSpeedU(TwoJetVec v)
{
	return sqrt(sqr(v.x.df_du()) + sqr(v.y.df_du()) + sqr(v.z.df_du()));
}

void get_eversion(k3d::mesh* mesh, bool polygonal, SurfaceTimeFunction* func, double umin, double umax, int ucount, double vmin, double vmax, int vcount, double t, geometry_matrix_t& points, int numStrips)
{
	static TwoJetVec** values;

	if(ucount <= 0 || vcount <= 0)
		return;

	double delta_u = (umax-umin) / ucount;
	double delta_v = (vmax-vmin) / vcount;
	values = (TwoJetVec**)calloc(ucount+1, sizeof(TwoJetVec*));
	//double* speedv = (double*)calloc(ucount+1, sizeof(double));
	//double** speedu = (double**)calloc(ucount+1, sizeof(double*));
	for(int j = 0; j <= ucount; j++)
		{
			double u = umin + j*delta_u;
			values[j] = (TwoJetVec*)calloc(vcount+1, sizeof(TwoJetVec));
			//speedu[j] = (double*)calloc(vcount+1, sizeof(double));
			//speedv[j] = calcSpeedV((*func)(ThreeJet(u, 1, 0), ThreeJet(0, 0, 1), t, numStrips));
			double speed_v = calcSpeedV((*func)(ThreeJet(u, 1, 0), ThreeJet(0, 0, 1), t, numStrips));
			if(speed_v == 0)
			//if(speedv[j] == 0)
				{
					// Perturb a bit, hoping to avoid degeneracy
					u += (u < 1) ? 1e-9 : -1e-9;
					//speedv[j] = calcSpeedV((*func)(ThreeJet(u, 1, 0), ThreeJet(0, 0, 1), t, numStrips));
				}

			for(int k = 0; k <= vcount; k++)
				{
					double v = vmin + k*delta_v;
					values[j][k] = (*func)(ThreeJet(u, 1, 0), ThreeJet(v, 0, 1), t, numStrips);
					//speedu[j][k] = calcSpeedU(values[j][k]);
				}
		}

	if(!polygonal)
		{
			// Bezier patches
			for(int j = 0; j < ucount; ++j)
				//double u = umin + j*delta_u;
				for(int k = 0; k < vcount; ++k)
					//double v = vmin + k*delta_v;
					get_patch(mesh, values[j][k], values[j][k+1], values[j+1][k], values[j+1][k+1], delta_u, delta_v, umin + j*delta_u, umin + (j+1)*delta_u, vmin + k*delta_v, vmin + (k+1)*delta_v);
		}
	else
		{
			// Quadrilaterals
			for(int j = 0; j <= ucount; ++j)
				{
					geometry_vector_t point_vector;
					for(int k = 0; k <= vcount; ++k)
						{
							eversion_point new_point;
							get_polygon_point(values[j][k], new_point);
							point_vector.push_back(new_point);
						}

					points.push_back(point_vector);
				}
		}

	// Clean up
	for(int j = 0; j <= ucount; j++)
		{
			free(values[j]);
			//free(speedu[j]);
		}

	free(values);
	//free(speedu);
	//free(speedv);
}

// Generate an eversion strip
void generate_strip(
	// Mesh
	k3d::mesh* mesh,
	// Return Bezier patches (false) or polygons (true)
	bool polygonal,
	// An array of (1 + u_count) arrays of (1 + v_count) elements
	geometry_matrix_t& points,
	// Eversion time, must be between 0.0 and 1.0
	double time = 0,
	// If this is too small, there will be pinches in the eversion
	int numStrips = 8,
	// The u parameter corresponds to latitude, v to longitude:
	// - parallels [90 degrees north, 90 degrees south] are mapped to [0.0, 2.0]
	double u_min = 0.0,	// set to 0.0 to start at north pole
	int u_count = 12,	// recommended value: 12*(u_max-u_min)
	double u_max = 1.0,	// set to 1.0 to stop at equator, or 2.0 to stop at south pole
	// - meridians [0 degrees, 180 degrees west] are mapped to [0, numStrips]
	double v_min = 0.0,
	int v_count = 12,	// recommended value: 12*(v_max-v_min)
	double v_max = 1.0,

	// -1 means don't do bendtime at all
	double bendtime = -1.0,
	double corrStart = 0.00,	// start of corrugation
	double pushStart = 0.10,	// start of push (poles are pushed through each other)
	double twistStart = 0.23,	// start of twist (poles rotate in opposite directions)
	double unpushStart = 0.60,	// start of unpush (poles held fixed while corrugations pushed through center)
	double uncorrStart = 0.93)	// start of uncorrugation
{
	point_map.clear();
	if(bendtime >= 0.0)
		{
			get_eversion(mesh, polygonal, BendIn, u_min, u_max, u_count, v_min, v_max, v_count, bendtime, points, numStrips);
			return;
		}

	// time = (time - howfar) / chunk

	if(time >= uncorrStart)
		get_eversion(mesh, polygonal, UnCorrugate, u_min, u_max, u_count, v_min, v_max, v_count,
			(time - uncorrStart) / (1.0 - uncorrStart), points, numStrips);
	else if(time >= unpushStart)
		get_eversion(mesh, polygonal, UnPush, u_min, u_max, u_count, v_min, v_max, v_count,
			(time - unpushStart) / (uncorrStart - unpushStart), points, numStrips);
	else if(time >= twistStart)
		get_eversion(mesh, polygonal, Twist, u_min, u_max, u_count, v_min, v_max, v_count,
			(time - twistStart) / (unpushStart - twistStart), points, numStrips);
	else if(time >= pushStart)
		get_eversion(mesh, polygonal, PushThrough, u_min, u_max, u_count, v_min, v_max, v_count,
			(time - pushStart) / (twistStart - pushStart), points, numStrips);
	else if(time >= corrStart)
		get_eversion(mesh, polygonal, Corrugate, u_min, u_max, u_count, v_min, v_max, v_count,
			(time - corrStart) / (pushStart - corrStart), points, numStrips);
}

} // namespace detail

/////////////////////////////////////////////////////////////////////////////
// sphere_eversion_implementation

class sphere_eversion_implementation :
	public k3d::material_collection<k3d::mesh_source<k3d::persistent<k3d::object> > >
{
	typedef k3d::material_collection<k3d::mesh_source<k3d::persistent<k3d::object> > > base;

public:
	sphere_eversion_implementation(k3d::idocument& Document) :
		base(Document),
		m_document(Document),
		m_time(k3d::init_name("time") + k3d::init_description("Eversion Time [number]") + k3d::init_value(0.0) + k3d::init_constraint(k3d::constraint::minimum(0.0, k3d::constraint::maximum(1.0))) + k3d::init_precision(3) + k3d::init_step_increment(0.01) + k3d::init_units(typeid(k3d::measurement::scalar)) + k3d::init_document(Document)),
		m_radius(k3d::init_name("radius") + k3d::init_description("Radius [number]") + k3d::init_value(5.0) + k3d::init_precision(2) + k3d::init_step_increment(0.1) + k3d::init_units(typeid(k3d::measurement::scalar)) + k3d::init_document(Document)),
		m_latitudinal_patches(k3d::init_name("latitudinal_patches") + k3d::init_description("Latitudinal Patches Per Hemisphere [number]") + k3d::init_document(Document) + k3d::init_value(6) + k3d::init_precision(0) + k3d::init_constraint(k3d::constraint::minimum(2UL)) + k3d::init_step_increment(1) + k3d::init_units(typeid(k3d::measurement::scalar))),
		m_longitudinal_patches(k3d::init_name("longitudinal_patches") + k3d::init_description("Longitudinal Patches Per Strip [number]") + k3d::init_document(Document) + k3d::init_value(4) + k3d::init_precision(0) + k3d::init_constraint(k3d::constraint::minimum(2UL)) + k3d::init_step_increment(1) + k3d::init_units(typeid(k3d::measurement::scalar))),
		m_strips(k3d::init_name("strips") + k3d::init_description("Strip Number [number]") + k3d::init_document(Document) + k3d::init_value(8) + k3d::init_precision(0) + k3d::init_constraint(k3d::constraint::minimum(1UL)) + k3d::init_step_increment(1) + k3d::init_units(typeid(k3d::measurement::scalar))),
		m_display_strips(k3d::init_name("display_strips") + k3d::init_description("Display Strip Number [number]") + k3d::init_document(Document) + k3d::init_value(8) + k3d::init_precision(0) + k3d::init_constraint(k3d::constraint::minimum(1UL)) + k3d::init_step_increment(1) + k3d::init_units(typeid(k3d::measurement::scalar))),
		m_half_strips(k3d::init_name("half_strips") + k3d::init_description("Half Strips") + k3d::init_value(false) + k3d::init_document(Document)),
		m_lower_hemisphere(k3d::init_name("lower_hemisphere") + k3d::init_description("Show lower hemisphere") + k3d::init_value(true) + k3d::init_document(Document)),
		m_cylinder_test(k3d::init_name("cylinder_test") + k3d::init_description("Cylinder bending into sphere") + k3d::init_value(false) + k3d::init_document(Document)),
		m_polygonal(k3d::init_name("polygonal") + k3d::init_description("Draw polygons (true) or patches (false) [boolean]") + k3d::init_value(false) + k3d::init_document(Document))
	{
		enable_serialization(k3d::persistence::proxy(m_time));
		enable_serialization(k3d::persistence::proxy(m_radius));
		enable_serialization(k3d::persistence::proxy(m_latitudinal_patches));
		enable_serialization(k3d::persistence::proxy(m_longitudinal_patches));
		enable_serialization(k3d::persistence::proxy(m_strips));
		enable_serialization(k3d::persistence::proxy(m_display_strips));
		enable_serialization(k3d::persistence::proxy(m_half_strips));
		enable_serialization(k3d::persistence::proxy(m_lower_hemisphere));
		enable_serialization(k3d::persistence::proxy(m_cylinder_test));
		enable_serialization(k3d::persistence::proxy(m_polygonal));

		register_property(m_time);
		register_property(m_radius);
		register_property(m_latitudinal_patches);
		register_property(m_longitudinal_patches);
		register_property(m_strips);
		register_property(m_display_strips);
		register_property(m_half_strips);
		register_property(m_lower_hemisphere);
		register_property(m_cylinder_test);
		register_property(m_polygonal);

		m_material.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));

		m_time.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_radius.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_latitudinal_patches.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_longitudinal_patches.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_strips.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_display_strips.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_half_strips.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_lower_hemisphere.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_cylinder_test.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));
		m_polygonal.changed_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_reset_geometry));

		m_output_mesh.need_data_signal().connect(SigC::slot(*this, &sphere_eversion_implementation::on_create_geometry));
	}

	void on_reset_geometry()
	{
		m_output_mesh.reset();
	}

	k3d::mesh* on_create_geometry()
	{
		const double time = m_time.property_value();
		const double radius = m_radius.property_value();
		const unsigned long latitudinal = m_latitudinal_patches.property_value();
		const unsigned long longitudinal = m_longitudinal_patches.property_value();
		const unsigned long strip_number = m_strips.property_value();
		const unsigned long display_strips_number = m_display_strips.property_value();
		const bool half_strips = m_half_strips.property_value();
		const bool lower_hemisphere = m_lower_hemisphere.property_value();
		const bool cylinder_test = m_cylinder_test.property_value();
		const bool polygonal = m_polygonal.property_value();

		k3d::mesh* const strip_mesh = new k3d::mesh();
		detail::current_material = m_material.interface();
		detail::current_radius = radius;
		detail::geometry_matrix_t sphere_points;
		detail::generate_strip(strip_mesh, polygonal, sphere_points, time, strip_number,
			0.0,
			latitudinal,
			1.0,
			0.0,
			longitudinal,
			half_strips ? 0.5 : 1.0,
			cylinder_test ? time : -1.0
		);

		if(polygonal)
			{
				// Create polygonal strip mesh ...
				std::vector< std::vector<k3d::point*> > points;
				for(unsigned long j = 0; j <= latitudinal; ++j)
					{
						std::vector<k3d::point*> vpoints;
						for(unsigned long k = 0; k <= longitudinal; ++k)
							{
								double* p = sphere_points[j][k].vertex;
								double x = p[0] * radius;
								double y = p[1] * radius;
								double z = p[2] * radius;

								// Fix handedness
								k3d::point* point = new k3d::point(x, y, -z);
								return_val_if_fail(point, 0);

								// Set normal
								p = sphere_points[j][k].normal;
								x = p[0];
								y = p[1];
								z = p[2];
								point->vertex_data["N"] = k3d::ri::normal(k3d::vector3(x, y, -z));

								strip_mesh->points.push_back(point);
								vpoints.push_back(point);
							}

						points.push_back(vpoints);
					}

				strip_mesh->polyhedra.push_back(new k3d::polyhedron());
				k3d::polyhedron& polyhedron = *strip_mesh->polyhedra.back();
				polyhedron.material = detail::current_material;

				for(unsigned long j = 0; j < latitudinal; ++j)
					{
						k3d::point* p1 = 0;
						k3d::point* p2 = 0;
						k3d::point* p3 = 0;
						bool strip_toggle = false;

						for(unsigned long k = 0; k <= longitudinal; ++k)
						for(unsigned long n = 0; n < 2; n++)
						{
							p1 = p2;
							p2 = p3;
							p3 = points[j+n][k];
							if(!p1)
								continue;

							k3d::polyhedron::edges_t edges;
							edges.push_back(new k3d::split_edge(p1));
							edges.push_back(new k3d::split_edge(p2));
							edges.push_back(new k3d::split_edge(p3));
							if(strip_toggle)
								std::reverse(edges.begin(), edges.end());
							strip_toggle = !strip_toggle;

							k3d::loop_edges(edges.begin(), edges.end());

							polyhedron.edges.insert(polyhedron.edges.end(), edges.begin(), edges.end());

							k3d::face* new_face = new k3d::face(edges.front());
							return_val_if_fail(new_face, 0);
							polyhedron.faces.push_back(new_face);
						}
					}

				return_val_if_fail(is_valid(polyhedron), 0);
				//assert_warning(is_solid(polyhedron));
			}

		// Duplicate strip to strip number
		unsigned long original_points = strip_mesh->points.size();

		k3d::mesh* const output = new k3d::mesh();

		double strip_rotation = 2*k3d::pi() / static_cast<double>(strip_number);
		for(unsigned long n = 0; n < display_strips_number && n < strip_number; n++)
			{
				k3d::deep_copy(*strip_mesh, *output);

				k3d::matrix4 rotation(k3d::rotation3D(k3d::vector3(0, 0, static_cast<double>(n) * strip_rotation)));
				if(n > 0)
					for(k3d::mesh::points_t::iterator p = output->points.end() - original_points; p != output->points.end(); p++)
						{
							(*p)->position = rotation * (*p)->position;
							if(polygonal)
								(*p)->vertex_data["N"] = k3d::ri::normal(rotation * boost::any_cast<k3d::ri::normal>((*p)->vertex_data["N"]));
						}

				if(!lower_hemisphere)
					continue;

				if(!polygonal)
					{
						k3d::deep_copy(*strip_mesh, *output);
						k3d::matrix4 rotation2(k3d::rotation3D(k3d::vector3(k3d::pi(), 0, 0)));
						if(n > 0)
							rotation = rotation2 * rotation;
						else
							rotation = rotation2;

						for(k3d::mesh::points_t::iterator p = output->points.end() - original_points; p != output->points.end(); p++)
							(*p)->position = rotation * (*p)->position;
					}
				else
					{
						k3d::deep_copy(*strip_mesh, *output);
						for(k3d::mesh::points_t::iterator p = output->points.end() - original_points; p != output->points.end(); p++)
							{
								k3d::vector3 position = rotation * (*p)->position;
								k3d::ri::normal normal = k3d::ri::normal(rotation * boost::any_cast<k3d::ri::normal>((*p)->vertex_data["N"]));

								position.n[2] = -position.n[2];
								normal.n[2] = -normal.n[2];

								(*p)->position = position;
								(*p)->vertex_data["N"] = k3d::ri::normal(normal);
							}

						// Invert triangles
						k3d::polyhedron& polyhedron = *output->polyhedra.back();
						for(k3d::polyhedron::faces_t::iterator face = polyhedron.faces.begin(); face != polyhedron.faces.end(); face++)
							{
								k3d::split_edge* first = (*face)->first_edge;
								k3d::split_edge* second = first->face_clockwise;
								k3d::split_edge* third = second->face_clockwise;

								second->face_clockwise = first;
								third->face_clockwise = second;
								first->face_clockwise = third;
							}
					}

			}

		delete strip_mesh;

		return output;
	}

	k3d::iplugin_factory& factory()
	{
		return get_factory();
	}

	static k3d::iplugin_factory& get_factory()
	{
		static k3d::plugin_factory<k3d::document_plugin<sphere_eversion_implementation>, k3d::interface_list<k3d::imesh_source > > factory(
			k3d::uuid(0x2e885acd, 0x9e914078, 0xb4fcb4a7, 0x415231c5),
			"SphereEversion",
			"Generates sphere eversion at any stage",
			"Objects",
			k3d::iplugin_factory::STABLE);

		return factory;
	}

private:
	k3d::idocument& m_document;
	k3d_measurement_property(double, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_time;
	k3d_measurement_property(double, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_radius;
	k3d_measurement_property(unsigned long, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_latitudinal_patches;
	k3d_measurement_property(unsigned long, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_longitudinal_patches;
	k3d_measurement_property(unsigned long, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_strips;
	k3d_measurement_property(unsigned long, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::with_constraint) m_display_strips;
	k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_half_strips;
	k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_lower_hemisphere;
	k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_cylinder_test;
	k3d_data_property(bool, k3d::immutable_name, k3d::change_signal, k3d::with_undo, k3d::local_storage, k3d::no_constraint) m_polygonal;
};

/////////////////////////////////////////////////////////////////////////////
// sphere_eversion_factory

k3d::iplugin_factory& sphere_eversion_factory()
{
	return sphere_eversion_implementation::get_factory();
}

} // namespace libk3dmesh