// Copyright (c) 2000, 2001, 2002, 2003 by David Scherer and others.
// See the file license.txt for complete license terms.
// See the file authors.txt for a complete list of contributors.

#include <boost/python/detail/wrap_python.hpp>
#include <boost/crc.hpp>

#include "util/errors.hpp"
#include "util/gl_enable.hpp"

#include "python/slice.hpp"
#include "python/curve.hpp"

#include <stdexcept>
#include <cassert>
#include <sstream>
#include <iostream>

// Recall that the default constructor for object() is a reference to None.

namespace cvisual { namespace python {

curve::curve()
	: antialias( true), radius(0.0), sides(4)
{
	for (size_t i=0; i<sides; i++) {
		curve_sc[i]  = (float) std::cos(i * 2 * M_PI / sides);
		curve_sc[i+sides] = (float) std::sin(i * 2 * M_PI / sides);
	}

	// curve_slice is a list of indices for picking out the correct vertices from
	// a list of vertices representing one side of a thick-line curve. The lower
	// indices (0-255) are used for all but one of the sides. The upper indices
	// (256-511) are used for the final side.
	for (int i=0; i<128; i++) {
		curve_slice[i*2]       = i*sides;
		curve_slice[i*2+1]     = i*sides + 1;
		curve_slice[i*2 + 256] = i*sides + (sides - 1);
		curve_slice[i*2 + 257] = i*sides;
	}
}

void
curve::set_radius( const double& radius)
{
	this->radius = radius;
}

void
curve::set_antialias( bool aa)
{
	this->antialias = aa;
}

bool
curve::degenerate() const
{
	return count < 2;
}

bool
curve::monochrome(float* tcolor, size_t pcount)
{
	rgb first_color( tcolor[0], tcolor[1], tcolor[2]);
	size_t nn;

	for(nn=0; nn<pcount; nn++)  {
		if (tcolor[nn*3] != first_color.red)
			return false;
		if (tcolor[nn*3+1] != first_color.green)
			return false;
		if (tcolor[nn*3+2] != first_color.blue)
			return false;
	}

	return true;
}

namespace {
// Determines if two values differ by more than frac of either one.
bool
eq_frac( double lhs, double rhs, double frac)
{
	if (lhs == rhs)
		return true;
	double diff = fabs(lhs - rhs);
	lhs = fabs(lhs);
	rhs = fabs(rhs);
	return lhs*frac > diff && rhs*frac > diff;
}
} // !namespace (unnamed)

vector
curve::get_center() const
{
	// TODO: Optimize this by only recomputing the center when the checksum of
	// the entire object has changed.
	// TODO: Only add the "optimization" if the checksum is actually faster than
	// computing the average value every time...
	if (degenerate())
		return vector();
	vector ret;
	const double* pos_i = pos.data();
	const double* pos_end = pos.end();
	while (pos_i < pos_end) {
		ret.x += pos_i[0];
		ret.y += pos_i[1];
		ret.z += pos_i[2];
		pos_i += 3;
	}
	ret /= count;
	return ret;
}

void
curve::gl_pick_render( const view& scene)
{
	// Aack, I can't think of any obvious optimizations here.
	// But since Visual 3 didn't permit picking of curves, omit for now.
	// We can't afford it; serious impact on performance.
	//gl_render( scene);
}

void
curve::grow_extent( extent& world)
{
	if (degenerate())
		return;
	const double* pos_i = pos.data();
	const double* pos_end = pos.end();
	if (radius == 0.0)
		for ( ; pos_i < pos_end; pos_i += 3)
			world.add_point( vector(pos_i));
	else
		for ( ; pos_i < pos_end; pos_i += 3)
			world.add_sphere( vector(pos_i), radius);
	world.add_body();
}

bool
curve::adjust_colors( const view& scene, float* tcolor, size_t pcount)
{
	rgb rendered_color;
	bool mono = monochrome(tcolor, pcount);
	if (mono) {
		// We can get away without using a color array.
		rendered_color = rgb( tcolor[0], tcolor[1], tcolor[2]);
		if (scene.anaglyph) {
			if (scene.coloranaglyph)
				rendered_color.desaturate().gl_set(opacity);
			else
				rendered_color.grayscale().gl_set(opacity);
		}
		else
			rendered_color.gl_set(opacity);
	}
	else {
		glEnableClientState( GL_COLOR_ARRAY);
		if (scene.anaglyph) {
			// Must desaturate or grayscale the color.

			for (size_t i = 0; i < pcount; ++i) {
				rendered_color = rgb( tcolor[3*i], tcolor[3*i+1], tcolor[3*i+2]);
				if (scene.coloranaglyph)
					rendered_color = rendered_color.desaturate();
				else
					rendered_color = rendered_color.grayscale();
				tcolor[3*i] = rendered_color.red;
				tcolor[3*i+1] = rendered_color.green;
				tcolor[3*i+2] = rendered_color.blue;
			}
		}
	}
	return mono;
}

namespace {
template <typename T>
struct converter
{
	T data[3];
};
} // !namespace (anonymous)

void
curve::thickline( const view& scene, double* spos, float* tcolor, size_t pcount, double scaled_radius)
{
	float *cost = curve_sc;
	float *sint = cost + sides;

	vector lastA; // unit vector of previous segment

	if (pcount < 2) return;

	bool closed = vector(&spos[0]) == vector(&spos[(pcount-1)*3]);

	size_t vcount = pcount*2 - closed;  // The number of vertices along each edge of the curve
	std::vector<vector> projected( vcount*sides );
	std::vector<vector> normals( vcount*sides );
	std::vector<rgb> light( vcount*sides );

	// pos and color iterators
	const double* v_i = spos;
	const float* c_i = tcolor;
	size_t i = closed ? 0 : sides;
	bool mono = adjust_colors( scene, tcolor, pcount);

	for (size_t corner=0; corner < pcount; ++corner, v_i += 3, c_i += 3) {
		vector current( &v_i[0] );

		vector next, A, bisecting_plane_normal;
		double sectheta;
		if (corner != pcount-1) {
			next = vector( &v_i[3] ); // The next vector in spos
			A = (next - current).norm();
			if (!A) {
				if (corner == 0) {
					const double* tv_i = v_i;
					for (size_t tcorner=0; tcorner < pcount; ++tcorner, tv_i += 3) {
						A = (vector( &tv_i[3] ) - current).norm();
						if (!A) continue;
					}
					if (!A) { // all the points of this curve are at the same location; abort
						return;
					}
					lastA = A;
				} else {
					A = lastA;
				}
			}
			bisecting_plane_normal = (A + lastA).norm();
			if (!bisecting_plane_normal) {  //< Exactly 180 degree bend
				bisecting_plane_normal = vector(0,0,1).cross(A);
				if (!bisecting_plane_normal)
					bisecting_plane_normal = vector(0,1,0).cross(A);
			}
			sectheta = bisecting_plane_normal.dot( lastA );
			if (sectheta) sectheta = 1.0 / sectheta;
		}

		if (corner == 0) {
			vector y = vector(0,1,0);
			vector x = A.cross(y).norm();
			if (!x) {
				x = A.cross( vector(0, 0, 1)).norm();
			}
			y = x.cross(A).norm();

			if (!x || !y || x == y) {
				std::ostringstream msg;
				msg << "Degenerate curve case! please report the following "
					"information to visualpython-users@lists.sourceforge.net: ";
				msg << "current:" << current << " next:" << next
			 		<< " A:" << A << " x:" << x << " y:" << y
			 		<< std::endl;
			 	VPYTHON_WARNING( msg.str());
			}

			// scale radii
			x *= scaled_radius;
			y *= scaled_radius;

			for (size_t a=0; a < sides; a++) {
				vector rel = x*sint[a] + y*cost[a]; // first point is "up"

				normals[a+i] = rel.norm();
				projected[a+i] = current + rel;
				if (!mono) light[a+i] = rgb( c_i );

				if (!closed) {
					// Cap start of curve
					projected[a] = current;
					normals[a] = -A;
					if (!mono) light[a] = light[a+i];
				}
			}

			i += sides;
		} else {
			double Adot = A.dot(next - current);
			for (size_t a=0; a < sides; a++) {
				vector prev_start = projected[i+a-sides];
				vector rel = current - prev_start;
				double t = rel.dot(lastA);
				if (corner != pcount-1 && sectheta > 0.0) {
					double t1 = (rel.dot(bisecting_plane_normal)) * sectheta;
					t1 = std::max( t1, t - Adot );
					t = std::max( 0.0, std::min( t, t1 ) );
				}
				vector prev_end = prev_start + t*lastA;

				projected[i+a] = prev_end;
				normals[i+a] = normals[i+a-sides];
				if (!mono) light[i+a] = rgb( c_i );

				if (corner != pcount-1) {
					vector next_start = prev_end - 2*(prev_end-current).dot(bisecting_plane_normal)*bisecting_plane_normal;

					rel = next_start - current;

					projected[i+a+sides] = next_start;
					normals[i+a+sides] = (rel - A.dot(next_start-current)*A).norm();
					if (!mono) light[i+a+sides] = light[i+a];
				} else if (!closed) {
					// Cap end of curve
					for (size_t a=0; a < sides; a++) {
						projected[i+a+sides] = current;
						normals[i+a+sides] = lastA;
						if (!mono) light[i+a+sides] = light[a+i];
					}
				}
			}
			i += 2*sides;
		}
		lastA = A;
	}

	if (closed) {
		// Connect the end of the curve to the start... can be ugly because the basis has gotten
		//   twisted around!
		size_t i = (vcount - 1)*sides;
		for(size_t a=0; a<sides; a++) {
			projected[i+a] = projected[a];
			normals[i+a] = normals[a];
			if (!mono) light[i+a] = light[a];
		}
	}

	// Thick lines are often used to represent smooth curves, so we want
	// to smooth the normals at the joints.  But that can make a sharp corner
	// do odd things, so we smoothly disable the smoothing when the joint angle
	// is too big.  This is somewhat arbitrary but seems to work well.
	size_t prev_i = closed ? (vcount-1)*sides : 0;
	for( i = closed ? 0 : sides; i < vcount*sides; i += 2*sides ) {
		for(size_t a=0; a<sides; a++) {
			vector& n1 = normals[i+a];
			vector& n2 = normals[prev_i+a];
			double smooth_amount = (n1.dot(n2) - .65) * 4.0;
			smooth_amount = std::min(1.0, std::max(0.0, smooth_amount));
			if (smooth_amount) {
				vector n_smooth = (n1+n2).norm() * smooth_amount;
				n1 = n1 * (1-smooth_amount) + n_smooth;
				n2 = n2 * (1-smooth_amount) + n_smooth;
			}
		}
		prev_i = i + sides;
	}

	gl_enable_client vertex_arrays( GL_VERTEX_ARRAY);
	gl_enable_client normal_arrays( GL_NORMAL_ARRAY);
	if (!mono) {
		glEnableClientState( GL_COLOR_ARRAY);
	}

	int *ind = curve_slice;
	for (size_t a=0; a < sides; a++) {
		size_t ai = a;
		if (a == sides-1) {
			ind += 256; // upper portion of curve_slice indices, for the last side
			ai = 0;
		}

		// List all the vertices for the ai-th side of the thick line:
		for (size_t i = 0; i < vcount; i += 127u) {
			glVertexPointer(3, GL_DOUBLE, sizeof( vector), &projected[i*sides + ai].x);
			if (!mono)
				glColorPointer(3, GL_FLOAT, sizeof( rgb), &light[(i*sides + ai)].red );
			glNormalPointer( GL_DOUBLE, sizeof(vector), &normals[i*sides + ai].x);
			if (vcount-i < 128)
				glDrawElements(GL_TRIANGLE_STRIP, 2*(vcount-i), GL_UNSIGNED_INT, ind);
			else
				glDrawElements(GL_TRIANGLE_STRIP, 256u, GL_UNSIGNED_INT, ind);
		}
	}
	if (!mono)
		glDisableClientState( GL_COLOR_ARRAY);
}

void
curve::gl_render( const view& scene)
{
	if (degenerate())
		return;
	const size_t true_size = count;
	// Set up the leading and trailing points for the joins.  See
	// glePolyCylinder() for details.  The intent is to create joins that are
	// perpendicular to the path at the last segment.  When the path appears
	// to be closed, it should be rendered that way on-screen.
	// The maximum number of points to display.
	const int LINE_LENGTH = 1000;
	// Data storage for the position and color data (plus room for 3 extra points)
	double spos[3*(LINE_LENGTH+3)];
	float tcolor[3*(LINE_LENGTH+3)]; // opacity not yet implemented for curves
	float fstep = (float)(count-1)/(float)(LINE_LENGTH-1);
	if (fstep < 1.0F) fstep = 1.0F;
	size_t iptr=0, iptr3, pcount=0;

	const double* p_i = pos.data();
	const double* c_i = color.data();

	// Choose which points to display
	for (float fptr=0.0; iptr < count && pcount < LINE_LENGTH; fptr += fstep, iptr = (int) (fptr+.5), ++pcount) {
		iptr3 = 3*iptr;
		spos[3*pcount] = p_i[iptr3];
		spos[3*pcount+1] = p_i[iptr3+1];
		spos[3*pcount+2] = p_i[iptr3+2];
		tcolor[3*pcount] = c_i[iptr3];
		tcolor[3*pcount+1] = c_i[iptr3+1];
		tcolor[3*pcount+2] = c_i[iptr3+2];
	}

	// Do scaling if necessary
	double scaled_radius = radius;
	if (scene.gcf != 1.0 || (scene.gcfvec[0] != scene.gcfvec[1])) {
		scaled_radius = radius*scene.gcfvec[0];
		for (size_t i = 0; i < pcount; ++i) {
			spos[3*i] *= scene.gcfvec[0];
			spos[3*i+1] *= scene.gcfvec[1];
			spos[3*i+2] *= scene.gcfvec[2];
		}
	}

	clear_gl_error();

	if (radius == 0.0) {
		glEnableClientState( GL_VERTEX_ARRAY);
		glDisable( GL_LIGHTING);
		// Assume monochrome.
		if (antialias) {
			glEnable( GL_LINE_SMOOTH);
		}

		glVertexPointer( 3, GL_DOUBLE, 0, spos);
		bool mono = adjust_colors( scene, tcolor, pcount);
		if (!mono) glColorPointer( 3, GL_FLOAT, 0, tcolor);
		glDrawArrays( GL_LINE_STRIP, 0, pcount);
		glDisableClientState( GL_VERTEX_ARRAY);
		glDisableClientState( GL_COLOR_ARRAY);
		glEnable( GL_LIGHTING);
		if (antialias) {
			glDisable( GL_LINE_SMOOTH);
		}
	}
	else {
		thickline( scene, spos, tcolor, pcount, scaled_radius);
	}

	check_gl_error();
}

void
curve::outer_render( const view& v ) {
	if (radius)
		arrayprim::outer_render(v);
	else
		gl_render(v);  //< no materials
}

void
curve::get_material_matrix( const view& v, tmatrix& out ) {
	if (degenerate()) return;

	// TODO: note this code is identical to faces::get_material_matrix, except for considering radius

	// TODO: Add some caching for extent with grow_extent etc
	vector min_extent, max_extent;
	const double* pos_i = pos.data();
	const double* pos_end = pos.end();
	min_extent = max_extent = vector( pos_i ); pos_i += 3;
	while (pos_i < pos_end)
		for(int j=0; j<3; j++) {
			if (*pos_i < min_extent[j]) min_extent[j] = *pos_i;
			else if (*pos_i > max_extent[j]) max_extent[j] = *pos_i;
			pos_i++;
		}

	min_extent -= vector(radius,radius,radius);
	max_extent += vector(radius,radius,radius);

	out.translate( vector(.5,.5,.5) );
	out.scale( vector(1,1,1) * (.999 / (v.gcf * std::max(max_extent.x-min_extent.x, std::max(max_extent.y-min_extent.y, max_extent.z-min_extent.z)))) );
	out.translate( -.5 * v.gcf * (min_extent + max_extent) );
}

} } // !namespace cvisual::python