/*
 * Copyright 2006, NVIDIA Corporation Ignacio Castano <icastano@nvidia.com>
 *
 * Modifications copyright (c) 2011, Blender Foundation. All rights reserved.
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include "subd_build.h"
#include "subd_edge.h"
#include "subd_face.h"
#include "subd_ring.h"
#include "subd_mesh.h"
#include "subd_patch.h"
#include "subd_stencil.h"
#include "subd_vert.h"

#include "util_algorithm.h"
#include "util_debug.h"
#include "util_math.h"
#include "util_string.h"

CCL_NAMESPACE_BEGIN

/* Subd Builder */

SubdBuilder *SubdBuilder::create(bool linear)
{
	if(linear)
		return new SubdLinearBuilder();
	else
		return new SubdAccBuilder();
}

/* Gregory ACC Stencil */

class GregoryAccStencil {
public:
	SubdFaceRing *ring;
	StencilMask stencil[20];

	GregoryAccStencil(SubdFaceRing *ring_)
	{
		ring = ring_;

		for(int i = 0; i < 20; i++)
			stencil[i].resize(ring->num_verts());
	}

	StencilMask& get(int i)
	{
		assert(i < 20);
		return stencil[i];
	}

	float& get(int i, SubdVert *vert)
	{
		assert(i < 20);
		return stencil[i][ring->vert_index(vert)];
	}
};

static float pseudoValence(SubdVert *vert)
{
	float valence = (float)vert->valence();

	if(vert->is_boundary()) {
		/* we treat boundary verts as being half a closed mesh. corners are
		 * special case. n = 4 for corners and n = 2*(n-1) for boundaries. */
		if(valence == 2) return 4;
		return (valence - 1)*2;
	}

	return valence;
}

/* Subd ACC Builder */

SubdAccBuilder::SubdAccBuilder()
{
}

SubdAccBuilder::~SubdAccBuilder()
{
}

Patch *SubdAccBuilder::run(SubdFace *face)
{
	SubdFaceRing ring(face, face->edge);
	GregoryAccStencil stencil(&ring);
	float3 position[20];

	computeCornerStencil(&ring, &stencil);
	computeEdgeStencil(&ring, &stencil);
	computeInteriorStencil(&ring, &stencil);

	ring.evaluate_stencils(position, stencil.stencil, 20);

	if(face->num_edges() == 3) {
		GregoryTrianglePatch *patch = new GregoryTrianglePatch();
		memcpy(patch->hull, position, sizeof(float3)*20);
		return patch;
	}
	else if(face->num_edges() == 4) {
		GregoryQuadPatch *patch = new GregoryQuadPatch();
		memcpy(patch->hull, position, sizeof(float3)*20);
		return patch;
	}

	assert(0); /* n-gons should have been split already */
	return NULL;
}

/* Gregory Patch */

void SubdAccBuilder::computeCornerStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
{
	const int cornerIndices[7] = {8, 11, 19, 16,   6, 9, 12};
	int primitiveOffset = ring->is_quad()? 0: 4;

	SubdEdge *firstEdge = ring->firstEdge();

	/* compute corner control points */
	int v = 0;

	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
		SubdVert *vert = it.current()->from();
		int valence = vert->valence();
		int cid = cornerIndices[primitiveOffset+v];

		if(vert->is_boundary()) {
			/* compute vertex limit position */
			SubdEdge *edge0 = vert->edge;
			SubdEdge *edge1 = vert->edge->prev;

			assert(edge0->face == NULL);
			assert(edge0->to() != vert);
			assert(edge1->face == NULL);
			assert(edge1->from() != vert);

			stencil->get(cid, vert) = 2.0f/3.0f;
			stencil->get(cid, edge0->to()) = 1.0f/6.0f;
			stencil->get(cid, edge1->from()) = 1.0f/6.0f;

			assert(stencil->get(cid).is_normalized());
		}
		else {
			stencil->get(cid, vert) = 3.0f*valence*valence;

			for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance()) {
				SubdEdge *edge = eit.current();
				assert(vert->co == edge->from()->co);

				stencil->get(cid, edge->to()) = 12.0f;

				if(SubdFaceRing::is_triangle(edge->face)) {
					/* distribute weight to all verts */
					stencil->get(cid, vert) += 1.0f;
					stencil->get(cid, edge->to()) += 1.0f;
					stencil->get(cid, edge->next->to()) += 1.0f;
				}
				else
					stencil->get(cid, edge->next->to()) = 3.0f;
			}

			/* normalize stencil. */
			stencil->get(cid).normalize();
		}
	}
}

void SubdAccBuilder::computeEdgeStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
{
	const int cornerIndices[7] = {8, 11, 19, 16,    6,  9, 12};
	const int edge1Indices[7] = {9, 13, 18, 14,    7, 10, 13};
	const int edge2Indices[7] = {12, 10, 15, 17,    14,  8, 11};
	int primitiveOffset = ring->is_quad()? 0: 4;

	float tangentScales[14] = {
		0.0f, 0.0f, 0.0f, 0.667791f, 1.0f,
		1.11268f, 1.1284f, 1.10289f, 1.06062f,
		1.01262f, 0.963949f, 0.916926f, 0.872541f, 0.831134f
	};

	SubdEdge *firstEdge = ring->firstEdge();

	/* compute corner / edge control points */
	int v = 0;

	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
		SubdVert *vert = it.current()->from();
		int valence = vert->valence();
		int cid = cornerIndices[primitiveOffset+v];

		int i1 = 0, i2 = 0, j = 0;

		for(SubdVert::EdgeIterator eit(vert->edge); !eit.isDone(); eit.advance(), j++) {
			SubdEdge *edge = eit.current();

			/* find index of "our" edge for edge control points */
			if(edge == it.current())
				i1 = j;
			if(edge == it.current()->prev->pair)
				i2 = j;
		}

		if(vert->is_boundary()) {
			int num_verts = ring->num_verts();
			StencilMask	r0(num_verts);
			StencilMask	r1(num_verts);

			computeBoundaryTangentStencils(ring, vert, r0, r1);

			int k = valence - 1;
			float omega = M_PI_F / k;

			int eid1 = edge1Indices[primitiveOffset + v];
			int eid2 = edge2Indices[primitiveOffset + v];

			if(it.current()->is_boundary()) {
				assert(it.current()->from() == vert);

				stencil->get(eid1, vert) = 2.0f / 3.0f;
				stencil->get(eid1, it.current()->to()) = 1.0f / 3.0f;

				assert(stencil->get(eid1).is_normalized());

				if(valence == 2) {
					for(int i = 0; i < num_verts; i++)
						stencil->get(eid1)[i] += r0[i] * 0.0001f;
				}
			}
			else {
				stencil->get(eid1) = stencil->get(cid);

				/* compute index of it.current() around vert */
				int idx = 0;

				for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
					if(eit.current() == it.current())
						break;

				assert(idx != valence);

				float c = cosf(idx * omega);
				float s = sinf(idx * omega);

				for(int i = 0; i < num_verts; i++)
					stencil->get(eid1)[i] += (r0[i] * s + r1[i] * c) / 3.0f;
			}

			if(it.current()->prev->is_boundary()) {
				assert(it.current()->prev->pair->from() == vert);

				stencil->get(eid2, vert) = 2.0f / 3.0f;
				stencil->get(eid2, it.current()->prev->pair->to()) = 1.0f / 3.0f;

				assert(stencil->get(eid2).is_normalized());

				if(valence == 2) {
					for(int i = 0; i < num_verts; i++)
						stencil->get(eid2)[i] += r0[i] * 0.0001f;
				}
			}
			else {
				stencil->get(eid2) = stencil->get(cid);

				/* compute index of it.current() around vert */
				int idx = 0;

				for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), idx++)
					if(eit.current() == it.current()->prev->pair)
						break;

				assert(idx != valence);

				float c = cosf(idx * omega);
				float s = sinf(idx * omega);

				for(int i = 0; i < num_verts; i++)
					stencil->get(eid2)[i] += (r0[i] * s + r1[i] * c) / 3;
			}
		}
		else {
			float costerm = cosf(M_PI_F / valence);
			float sqrtterm = sqrtf(4.0f + costerm*costerm);

			/* float tangentScale = 1.0f; */
			float tangentScale = tangentScales[min(valence, 13U)];

			float alpha = (1.0f +  costerm / sqrtterm) / (3.0f * valence) * tangentScale;
			float beta  = 1.0f / (3.0f * valence * sqrtterm) * tangentScale;


			int eid1 = edge1Indices[primitiveOffset + v];
			int eid2 = edge2Indices[primitiveOffset + v];

			stencil->get(eid1) = stencil->get(cid);
			stencil->get(eid2) = stencil->get(cid);

			j = 0;
			for(SubdVert::EdgeIterator eit(vert->edges()); !eit.isDone(); eit.advance(), j++) {
				SubdEdge *edge = eit.current();
				assert(vert->co == edge->from()->co);

				float costerm1_a = cosf(M_PI_F * 2 * (j-i1) / valence);
				float costerm1_b = cosf(M_PI_F * (2 * (j-i1)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */

				float costerm2_a = cosf(M_PI_F * 2 * (j-i2) / valence);
				float costerm2_b = cosf(M_PI_F * (2 * (j-i2)-1) / valence); /* -1 instead of +1 b/c of edge->next->to() */


				stencil->get(eid1, edge->to()) += alpha * costerm1_a;
				stencil->get(eid2, edge->to()) += alpha * costerm2_a;

				if(SubdFaceRing::is_triangle(edge->face)) {
					/* @@ this probably does not provide watertight results!! (1/3 + 1/3 + 1/3 != 1) */

					/* distribute weight to all verts */
					stencil->get(eid1, vert) += beta * costerm1_b / 3.0f;
					stencil->get(eid1, edge->to()) += beta * costerm1_b / 3.0f;
					stencil->get(eid1, edge->next->to()) += beta * costerm1_b / 3.0f;

					stencil->get(eid2, vert) += beta * costerm2_b / 3.0f;
					stencil->get(eid2, edge->to()) += beta * costerm2_b / 3.0f;
					stencil->get(eid2, edge->next->to()) += beta * costerm2_b / 3.0f;
				}
				else {
					stencil->get(eid1, edge->next->to()) += beta * costerm1_b;
					stencil->get(eid2, edge->next->to()) += beta * costerm2_b;
				}
			}
		}
	}
}

void SubdAccBuilder::computeInteriorStencil(SubdFaceRing *ring, GregoryAccStencil *stencil)
{
	static int corner1Indices[7] = {8, 11, 19, 16,    6, 9, 12};
	static int corner2Indices[7] = {11, 19, 16, 8,    9, 12, 6};
	static int edge1Indices[7] = {9, 13, 18, 14,    7, 10, 13};
	static int edge2Indices[7] = {10, 15, 17, 12,    8, 11, 14};
	static int interior1Indices[7] = {1, 3, 6, 4,    1, 3, 5};
	static int interior2Indices[7] = {2, 7, 5, 0,    2, 4, 0};

	int primitiveOffset = ring->is_quad()? 0: 4;

	SubdFace * face = ring->face();
	SubdEdge *firstEdge = ring->firstEdge();

	/* interior control points */
	int v = 0;
	for(SubdFace::EdgeIterator it(firstEdge); !it.isDone(); it.advance(), v++) {
		SubdEdge *edge = it.current();

		if(edge->is_boundary()) {
			float valence1 = pseudoValence(edge->from());
			float valence2 = pseudoValence(edge->to());

			float weights1[4];
			float weights2[4];

			if(ring->is_quad()) {
				weights1[0] = 3 * valence1;
				weights1[1] = 6;
				weights1[2] = 3;
				weights1[3] = 6;

				weights2[0] = 6;
				weights2[1] = 3 * valence2;
				weights2[2] = 6;
				weights2[3] = 3;
			}
			else {
				assert(ring->is_triangle());
				weights1[0] = 3 * valence1 + 1;
				weights1[1] = 7;
				weights1[2] = 7;

				weights2[0] = 7;
				weights2[1] = 3 * valence2 + 1;
				weights2[2] = 7;
			}

			int idx1 = interior1Indices[primitiveOffset+v];
			int idx2 = interior2Indices[primitiveOffset+v];

			int i = 0;
			for(SubdFace::EdgeIterator it_sub(face->edges(edge)); !it_sub.isDone(); it_sub.advance(), i++) {
				SubdVert *vert = it_sub.current()->from();
				stencil->get(idx1, vert) += weights1[i];
				stencil->get(idx2, vert) += weights2[i];
			}

			stencil->get(idx1).normalize();
			stencil->get(idx2).normalize();
		}
		else {
			SubdVert *e0 = edge->from();
			float costerm0 = cosf(2.0f * M_PI_F / pseudoValence(e0));

			SubdVert *f0 = edge->to();
			float costerm1 = cosf(2.0f * M_PI_F / pseudoValence(f0));

			/*  p0 +------+ q0
			 *	 |	  |
			 *  f0 +======+ e0 <=== current edge
			 *	 |	  |
			 *  p1 +------+ q1
			 */

			SubdVert *q0 = edge->next->to();
			SubdVert *p0 = edge->prev->from();

			SubdVert *p1 = edge->pair->next->to();
			SubdVert *q1 = edge->pair->prev->from();


			StencilMask	x(ring->num_verts());
			StencilMask	y(ring->num_verts());

			for(int i = 0; i < ring->num_verts(); i++) {
				x[i] =
					(costerm1 * stencil->get(corner1Indices[primitiveOffset+v])[i] -
					(2*costerm0 + costerm1) * stencil->get(edge1Indices[primitiveOffset+v])[i] +
					2*costerm0 * stencil->get(edge2Indices[primitiveOffset+v])[i]) / 3.0f;
			}

			/* y = (2*( midedgeA1 - midedgeB1) + 4*(centroidA - centroidB))/18.0f; */
			y[ring->vert_index(p0)] = 1;
			y[ring->vert_index(p1)] = -1;

			/* add centroidA */
			if(ring->is_triangle()) {
				y[ring->vert_index(p0)] += 4.0f / 3.0f;
				y[ring->vert_index(e0)] += 4.0f / 3.0f;
				y[ring->vert_index(f0)] += 4.0f / 3.0f;
			}
			else {
				y[ring->vert_index(p0)] += 1;
				y[ring->vert_index(q0)] += 1;
				y[ring->vert_index(e0)] += 1;
				y[ring->vert_index(f0)] += 1;
			}

			/* sub centroidB */
			if(SubdFaceRing::is_triangle(edge->pair->face)) {
				y[ring->vert_index(p1)] -= 4.0f / 3.0f;
				y[ring->vert_index(e0)] -= 4.0f / 3.0f;
				y[ring->vert_index(f0)] -= 4.0f / 3.0f;

			}
			else {
				y[ring->vert_index(p1)] -= 1;
				y[ring->vert_index(q1)] -= 1;
				y[ring->vert_index(e0)] -= 1;
				y[ring->vert_index(f0)] -= 1;
			}

			y /= 18.0f;

			if(ring->is_triangle()) {
				x *= 3.0f / 4.0f;
				y *= 3.0f / 4.0f;
			}

			/* this change makes the triangle boundaries smoother, but distorts the quads next to them */
#if 0
			if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face)) {
				y *= 4.0f / 3.0f;
			}
#endif

			stencil->get(interior1Indices[primitiveOffset+v]) = stencil->get(edge1Indices[primitiveOffset+v]);
			stencil->get(interior1Indices[primitiveOffset+v]) += x;
			stencil->get(interior1Indices[primitiveOffset+v]) += y;

			for(int i = 0; i < ring->num_verts(); i++) {
				x[i] =
					(costerm0 * stencil->get(corner2Indices[primitiveOffset+v])[i] -
					(2*costerm1 + costerm0) * stencil->get(edge2Indices[primitiveOffset+v])[i] + 
					2*costerm1 * stencil->get(edge1Indices[primitiveOffset+v])[i]) / 3.0f;
			}

			/* y = (2*( midedgeA2 - midedgeB2) + 4*(centroidA - centroidB))/18.0f; */
			y = 0.0f;

			/* (2*( midedgeA2 - midedgeB2) */
			y[ring->vert_index(q0)] = 1;
			y[ring->vert_index(q1)] = -1;

			/* add centroidA */
			if(ring->is_triangle()) {
				y[ring->vert_index(p0)] += 4.0f / 3.0f;
				y[ring->vert_index(e0)] += 4.0f / 3.0f;
				y[ring->vert_index(f0)] += 4.0f / 3.0f;
			}
			else {
				y[ring->vert_index(p0)] += 1;
				y[ring->vert_index(q0)] += 1;
				y[ring->vert_index(e0)] += 1;
				y[ring->vert_index(f0)] += 1;
			}

			/* sub centroidB */
			if(SubdFaceRing::is_triangle(edge->pair->face)) {
				y[ring->vert_index(p1)] -= 4.0f / 3.0f;
				y[ring->vert_index(e0)] -= 4.0f / 3.0f;
				y[ring->vert_index(f0)] -= 4.0f / 3.0f;

			}
			else {
				y[ring->vert_index(p1)] -= 1;
				y[ring->vert_index(q1)] -= 1;
				y[ring->vert_index(e0)] -= 1;
				y[ring->vert_index(f0)] -= 1;
			}

			y /= 18.0f;

			if(ring->is_triangle()) {
				x *= 3.0f / 4.0f;
				y *= 3.0f / 4.0f;
			}

			/* this change makes the triangle boundaries smoother, but distorts the quads next to them. */
#if 0
			if(ring->is_triangle() || SubdFaceRing::is_triangle(edge->pair->face))
				y *= 4.0f / 3.0f;
#endif

			stencil->get(interior2Indices[primitiveOffset+v]) = stencil->get(edge2Indices[primitiveOffset+v]);
			stencil->get(interior2Indices[primitiveOffset+v]) += x;
			stencil->get(interior2Indices[primitiveOffset+v]) += y;
		}
	}
}

void SubdAccBuilder::computeBoundaryTangentStencils(SubdFaceRing *ring, SubdVert *vert, StencilMask & r0, StencilMask & r1)
{
	assert(vert->is_boundary());
	assert(r0.size() == ring->num_verts());
	assert(r1.size() == ring->num_verts());

	SubdEdge *edge0 = vert->edge;
	assert(edge0->face == NULL);
	assert(edge0->to() != vert);

	SubdEdge *edgek = vert->edge->prev;
	assert(edgek->face == NULL);
	assert(edgek->from() != vert);

	int valence = vert->valence();

	int k = valence - 1;
	float omega = M_PI_F / k;
	float s = sinf(omega);
	float c = cosf(omega);

	float factor = 1.0f / (3 * k + c);

	float gamma = -4 * s * factor;
	r0[ring->vert_index(vert)] = gamma;
	/* r1[ring->vert_index(vert)] = 0; */

	float salpha0 = -((1 + 2 * c) * sqrtf(1 + c)) * factor / sqrtf(1 - c);
	float calpha0 = 1.0f / 2.0f;

	r0[ring->vert_index(edge0->to())] = salpha0;
	r1[ring->vert_index(edge0->to())] = calpha0;

	float salphak = salpha0;
	float calphak = -1.0f / 2.0f;

	r0[ring->vert_index(edgek->from())] = salphak;
	r1[ring->vert_index(edgek->from())] = calphak;

	int j = 0;
	for(SubdVert::EdgeIterator it(vert->edges()); !it.isDone(); it.advance(), j++) {
		SubdEdge *edge = it.current();

		if(j == k) break;

		SubdVert *p = edge->to();
		SubdVert *q = edge->pair->prev->from();

		float alphaj = 4 * sinf(j * omega) * factor;
		float betaj = (sinf(j * omega) + sinf((j + 1) * omega)) * factor;

		if(j != 0)
			r0[ring->vert_index(p)] += alphaj;

		if(edge->pair->prev->prev->prev == edge->pair) {
			r0[ring->vert_index(vert)] += betaj / 3.0f;
			r0[ring->vert_index(edge->pair->from())] += betaj / 3.0f;
			r0[ring->vert_index(q)] += betaj / 3.0f;
		}
		else
			r0[ring->vert_index(q)] += betaj;
	}

	if(valence == 2) {
		/* r0 perpendicular to r1 */
		r0[ring->vert_index(vert)] = -4.0f / 3.0f;
		r0[ring->vert_index(edgek->from())] = 1.0f / 2.0f;
		r0[ring->vert_index(edge0->to())] = 1.0f / 2.0f;
		r0[ring->vert_index(edge0->next->to())] = 1.0f / 3.0f;
	}
}

/* Subd Linear Builder */

SubdLinearBuilder::SubdLinearBuilder()
{
}

SubdLinearBuilder::~SubdLinearBuilder()
{
}

Patch *SubdLinearBuilder::run(SubdFace *face)
{
	Patch *patch;
	float3 *hull;

	if(face->num_edges() == 3) {
		LinearTrianglePatch *lpatch = new LinearTrianglePatch();
		hull = lpatch->hull;
		patch = lpatch;
	}
	else if(face->num_edges() == 4) {
		LinearQuadPatch *lpatch = new LinearQuadPatch();
		hull = lpatch->hull;
		patch = lpatch;
	}
	else {
		assert(0); /* n-gons should have been split already */
		return NULL;
	}

	int i = 0;

	for(SubdFace::EdgeIterator it(face->edge); !it.isDone(); it.advance())
		hull[i++] = it.current()->from()->co;

	if(face->num_edges() == 4)
		swap(hull[2], hull[3]);

	return patch;
}

CCL_NAMESPACE_END