// 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 Bart Janssens <bart.janssens@lid.kviv.be>
*/

#include "subdiv_algorithms.h"

#include <k3dsdk/result.h>
#include <k3dsdk/irenderman.h>

#include <cmath>
#include <map>
#include <set>

#define DEBUG_MESSAGE(Message)
//#define DEBUG_MESSAGE(Message) Message

namespace subdiv
{

// Helper class to associate a vertex to an edge.
// Used to i.e. associate face vertices to edges
class vertex_edge_container
{
	typedef std::map<k3d::split_edge*, k3d::point*> vertex_map;
	typedef vertex_map::value_type vertex_map_t;
	typedef vertex_map::iterator vertex_map_it;

public:
	// Add a new edge-vertex pair, if the edge is not registered already
	void add(k3d::split_edge* edge, k3d::point* vertex)
	{
		if(!((m_vertices.insert(vertex_map_t(edge, vertex))).second))
			DEBUG_MESSAGE(std::cerr << error << "Edge already has a vertex associated" << std::endl;);
	}

	// Return the vertex associated with edge, if there is one.
	k3d::point* get(k3d::split_edge* edge)
	{
		vertex_map_it it = m_vertices.find(edge);
		if(it != m_vertices.end())
			return it->second;

		return 0;
	}

private:
	vertex_map m_vertices;
};

// Helper class to associate 2 edges
class edge_container
{
	typedef std::map<k3d::split_edge*, k3d::split_edge*> edge_map;
	typedef edge_map::value_type edge_map_t;
	typedef edge_map::iterator edge_map_it;

public:
	// Add a new pair
	void add(k3d::split_edge* key_edge, k3d::split_edge* edge)
	{
		if(!((m_edges.insert(edge_map_t(key_edge, edge))).second))
			DEBUG_MESSAGE(std::cerr << error << "Edge already has an edge associated" << std::endl;);
	}

	// Return the edge associated to <edge>
	k3d::split_edge* get(k3d::split_edge* edge)
	{
		edge_map_it it = m_edges.find(edge);
		if(it != m_edges.end())
			return it->second;

		return 0;
	}

private:
	edge_map m_edges;
};

// Extra data stored by the catmull-clark algorithm
class catmull_clark_data
{
	typedef std::set<k3d::point*> point_set;
	typedef std::set<k3d::point*>::value_type point_set_t;
	typedef std::set<k3d::point*>::iterator point_set_it;

public:
	catmull_clark_data(k3d::polyhedron* polyhedron, k3d::mesh::points_t* output_points, bool ignore_selection) : m_ignore_selection(ignore_selection)
	{
		m_polyhedron = polyhedron;
		m_original_faces = m_polyhedron->faces;
		m_original_edges = m_polyhedron->edges;
		m_next_face = m_original_faces.begin();
		m_output_points = output_points;
	}

	// Return the companion of an edge
	k3d::split_edge* companion(k3d::split_edge* edge)
	{
		// return the old companion, if there is one
		k3d::split_edge* old_companion = m_saved_companions.get(edge);
		if(old_companion)
			return old_companion;

		// otherwise return the real companion
		return edge->companion;
	}

	// Return the face_clockwise of an edge
	k3d::split_edge* clockwise(k3d::split_edge* edge)
	{
		// return the old clockwise, if there is one
		k3d::split_edge* old_clockwise = m_saved_clockwise.get(edge);
		if(old_clockwise)
			return old_clockwise;

		// otherwise return the real face_clockwise
		return edge->face_clockwise;
	}

	// Helper function to return the face vertex to the right of an edge
	k3d::point* face_vertex(k3d::split_edge* edge)
	{
		// If this face vertex was calculated earlier, return the pointer to it
		return_val_if_fail(edge, 0);
		k3d::point* existing_face_vertex = m_face_vertices.get(edge);
		if(existing_face_vertex)
			return existing_face_vertex;

		// Else calculate the new face vertex and add it to the list of calculated vertices.
		k3d::point* vf = new k3d::point(edge->vertex->position);
		// Number of vertices in the face
		int n = 1;
		k3d::split_edge* ej = edge->face_clockwise;
		while(ej != edge)
			{
				// Associate the new face vertex with all surrounding edges
				m_face_vertices.add(ej, vf);
				vf->position += ej->vertex->position;
				ej = ej->face_clockwise;
				n++;
			}

		m_face_vertices.add(edge, vf);
		vf->position /= n;
		m_output_points->push_back(vf);
		vf->selected = !m_ignore_selection;

		return vf;
	}

	// Return the edge vertex associated with a given edge and face vertex
	k3d::point* edge_vertex(k3d::split_edge* edge, k3d::point* vf)
	{
		// If this edge vertex exists, return a pointer to it
		k3d::point* existing_edge_vertex = m_edge_vertices.get(edge);
		if(existing_edge_vertex)
			return existing_edge_vertex;

		// Else calculate a new one
		k3d::split_edge* edge_companion = companion(edge);
		if(!(edge_companion) || (!m_ignore_selection && edge_companion->selected == false))
		{
			// if edge is a boundary edge or at the edge of the selection
			k3d::point* ev = new k3d::point((edge->vertex->position
					+ edge->face_clockwise->vertex->position)/2);
			m_output_points->push_back(ev);
			if(edge_companion)
				m_edge_vertices.add(edge_companion, ev);

			ev->selected = !m_ignore_selection;

			return ev;
		}

		// For a normal edge:
		k3d::point* ev = new k3d::point((edge->vertex->position
			+ edge->face_clockwise->vertex->position
			+ vf->position
			+ face_vertex(edge_companion)->position)/4);

		// The companion edge has the same edge vertex, so store it:
		m_edge_vertices.add(edge_companion, ev);
		m_output_points->push_back(ev);
		ev->selected = !m_ignore_selection;

		return ev;
	}

	// Return the next edge to be treated, or 0 if there is none
	k3d::split_edge* next_edge()
	{
		if(m_next_face != m_original_faces.end())
		{
			while(!m_ignore_selection)
			{
				if((*m_next_face)->selected)
					return (*(m_next_face++))->first_edge;
				else
					++m_next_face;

				if(m_next_face == m_original_faces.end())
					return 0;
			}

			return (*(m_next_face++))->first_edge;
		}

		return 0;
	}

	// Join 2 edges, and store the old companion of the first edge
	void save_join(k3d::split_edge* first, k3d::split_edge* second)
	{
		if(first == 0 || second == 0 || first->companion == 0)
			return;

		m_saved_companions.add(first, first->companion);
		k3d::join_edges(*first, *second);
	}

	// Set the clockwise on an edge, saving the old one
	void save_clockwise(k3d::split_edge* edge, k3d::split_edge* clockwise)
	{
		if(edge == 0 || clockwise == 0 || edge->face_clockwise == 0)
			return;

		m_saved_clockwise.add(edge, edge->face_clockwise);
		edge->face_clockwise = clockwise;
	}

	// Reposition a point
	void reposition_point(k3d::split_edge* edge)
	{
		// Valence of the vertex
		int n = 1;

		if(!(edge->companion))
		{
			// Boundary edge
			k3d::split_edge* ecl = edge->face_clockwise->face_clockwise->face_clockwise;
			if(!m_ignore_selection && ecl->selected == false)
				return;
			if(!m_ignore_selection && edge->face_clockwise->selected == false)
				return;

			k3d::vector3 v = 0.25*edge->face_clockwise->vertex->position;
			if(!(ecl->companion))
				v += 0.25*ecl->vertex->position;
			else
				v += 0.25*ecl->companion->face_clockwise->face_clockwise->face_clockwise->vertex->position;

			edge->vertex->position = v + 0.5*edge->vertex->position;

			return;
		}

		k3d::split_edge* ecl1 = edge->face_clockwise;
		if(!m_ignore_selection && ecl1->selected == false)
			return;
		k3d::split_edge* ecl2 = ecl1->face_clockwise;
		if(!m_ignore_selection && ecl2->selected == false)
			return;
		k3d::split_edge* ecl3 = ecl2->face_clockwise;
		if(!m_ignore_selection && ecl3->selected == false)
			return;

		k3d::vector3 v = ecl1->vertex->position; // first edge vertex
		v += ecl2->vertex->position; // first face vertex
		k3d::split_edge* ei = ecl3->companion;
		if(!ei)
			{
				// we started on the wrong edge around a valence 3 vertex
				reposition_point(edge->companion->face_clockwise);
				return;
			}

		while(ei != edge)
		{
			n++;
			ecl1 = ei->face_clockwise;
			if(!m_ignore_selection && ecl1->selected == false)
				return;
			ecl2 = ecl1->face_clockwise;
			if(!m_ignore_selection && ecl2->selected == false)
				return;
			ecl3 = ecl2->face_clockwise;
			if(!m_ignore_selection && ecl3->selected == false)
				return;

			v += ecl1->vertex->position;
			v += ecl2->vertex->position;
			ei = ecl3->companion;
			if(!ei)
				{
					// we started on the wrong edge around a valence 3 vertex
					reposition_point(edge->companion->face_clockwise);
					return;
				}
		}

		edge->vertex->position = (n-2.0)/n * edge->vertex->position + 1.0/(n*n) * v;
	}

	// Reposition all the original points
	void reposition_points()
	{
		point_set visited_points;
		k3d::polyhedron::edges_t::iterator it = m_original_edges.begin();
		while(it != m_original_edges.end())
		{
			if(*it && (m_ignore_selection || (!m_ignore_selection && (*it)->selected == true)) && visited_points.find((*it)->vertex) == visited_points.end())
				{
					reposition_point(*it);
					visited_points.insert(point_set_t((*it)->vertex));
				}

			++it;
		}
	}

private:
	vertex_edge_container m_face_vertices;
	vertex_edge_container m_edge_vertices;
	edge_container m_saved_companions;
	edge_container m_saved_clockwise;
	k3d::polyhedron* m_polyhedron;
	k3d::polyhedron::faces_t m_original_faces;
	k3d::polyhedron::faces_t::iterator m_next_face;
	k3d::mesh::points_t* m_output_points;
	k3d::polyhedron::edges_t m_original_edges;
	const bool m_ignore_selection;
};

void catmull_clark(const unsigned long Levels, const k3d::mesh& Input, k3d::mesh& Output, const bool ignore_selection)
{
	// clear the output
	Output.polyhedra.clear();
	k3d::deep_copy(Input, Output);
	for(k3d::mesh::polyhedra_t::iterator it = Output.polyhedra.begin(); it != Output.polyhedra.end(); ++it) {
		if((*it)->type == k3d::polyhedron::CATMULL_CLARK_SUBDIVISION_MESH || !ignore_selection) {
			k3d::polyhedron* output_polyhedron = *(it);
			for(unsigned long i = 1; i <= Levels; i++) {
				assert(output_polyhedron);
				return_if_fail(k3d::is_valid(*output_polyhedron));
				// extra containers and functions we need:
				catmull_clark_data data(output_polyhedron, &(Output.points), ignore_selection);
				DEBUG_MESSAGE(std::cerr << debug << "Going to SDS level " << i << std::endl;);

				// clear the faces of the output:
				if(ignore_selection) {
					output_polyhedron->faces.clear();
				} else {
					k3d::polyhedron::faces_t original_faces = output_polyhedron->faces;
					output_polyhedron->faces.clear();
					for(k3d::polyhedron::faces_t::iterator f = original_faces.begin(); f != original_faces.end(); ++f) {
						if(!((*f)->selected)) {
							output_polyhedron->faces.push_back(*f);
						}
					}
				}
				k3d::split_edge* e_start = 0;
				k3d::split_edge* ei = 0;
				k3d::split_edge* ei_old_clockwise = 0;
				k3d::split_edge* ev_edge = 0;
				k3d::split_edge* ev_center = 0;
				k3d::split_edge* ev_edge_last = 0;
				k3d::split_edge* vf_center = 0;
				k3d::split_edge* vf_center_last = 0;
				k3d::split_edge* vf_center_first = 0;
				k3d::point* vf = 0;
				k3d::point* ev_last = 0;
				k3d::point* ev = 0;
				k3d::point* ev_first = 0;
				// loop through all original faces:
				e_start = data.next_edge();
				while(e_start) {
					// face loop initialisation:
					ei = data.clockwise(e_start);
					vf = data.face_vertex(e_start);
					ev_last = data.edge_vertex(e_start, vf);
					// fix the neighbour face, to allow local subdivision
					k3d::split_edge* e_start_companion = data.companion(e_start);
					if(e_start_companion && (ignore_selection || (!ignore_selection && e_start_companion->selected == true))) {
						data.face_vertex(e_start_companion);
					}
					if(!(e_start->face_clockwise->vertex == ev_last)) {
						ev_edge_last = new k3d::split_edge(ev_last, ei, 0);
						ev_edge_last->selected = !ignore_selection;
						output_polyhedron->edges.push_back(ev_edge_last);
						data.save_join(e_start_companion, ev_edge_last);
						if(e_start_companion) {
							k3d::split_edge* ev_comp_edge = new k3d::split_edge(ev_last, (e_start_companion)->face_clockwise);
							ev_comp_edge->selected = (!ignore_selection && e_start_companion->selected);
							output_polyhedron->edges.push_back(ev_comp_edge);
							data.save_clockwise(e_start_companion,ev_comp_edge);
							data.save_join(e_start, ev_comp_edge);
						}
					} else {
						ev_edge_last = e_start->face_clockwise;
					}
					vf_center = new k3d::split_edge(vf, ev_edge_last, 0);
					vf_center->selected = !ignore_selection;
					output_polyhedron->edges.push_back(vf_center);
					vf_center_first = vf_center;
					vf_center_last = vf_center;
					ev_first = ev_last;
					// the loop for one face:
					while(ei != e_start) {
						ei_old_clockwise = data.clockwise(ei);
						ev = data.edge_vertex(ei, vf);
						// fix the neighbour face, to allow local subdivision
						k3d::split_edge* ei_companion = data.companion(ei);
						if(ei_companion && (ignore_selection || (!ignore_selection && ei_companion->selected == true))) {
							data.face_vertex(ei_companion);
						}
						if(!(ei->face_clockwise->vertex == ev)) {
							ev_edge = new k3d::split_edge(ev, ei->face_clockwise, 0);
							ev_edge->selected = !ignore_selection;
							output_polyhedron->edges.push_back(ev_edge);
							data.save_join(ei_companion, ev_edge);
							if(ei_companion) {
								k3d::split_edge* ev_comp_edge = new k3d::split_edge(ev, ei_companion->face_clockwise);
								ev_comp_edge->selected = (!ignore_selection && ei_companion->selected);
								output_polyhedron->edges.push_back(ev_comp_edge);
								data.save_clockwise(ei_companion, ev_comp_edge);
								data.save_join(ei, ev_comp_edge);
							}
						} else {
							ev_edge = ei->face_clockwise;
						}
						vf_center = new k3d::split_edge(vf, ev_edge, 0);
						vf_center->selected = !ignore_selection;
						output_polyhedron->edges.push_back(vf_center);
						ev_center = new k3d::split_edge(ev, vf_center_last, 0);
						ev_center->selected = !ignore_selection;
						output_polyhedron->edges.push_back(ev_center);
						k3d::join_edges(*vf_center, *ev_center);
						ei->face_clockwise = ev_center;
						vf_center_last = vf_center;
						output_polyhedron->faces.push_back(new k3d::face(ev_center));
						output_polyhedron->faces.back()->selected = !ignore_selection;
						ei = ei_old_clockwise;
					}
					// finish up the loose edges
					ev_center = new k3d::split_edge(ev_first, vf_center_last, 0);
					ev_center->selected = !ignore_selection;
					output_polyhedron->edges.push_back(ev_center);
					e_start->face_clockwise = ev_center;
					k3d::join_edges(*vf_center_first, *ev_center);
					output_polyhedron->faces.push_back(new k3d::face(ev_center));
					output_polyhedron->faces.back()->selected = !ignore_selection;
					e_start = data.next_edge();
				}
				data.reposition_points();
			}
			return_if_fail(k3d::is_valid(*output_polyhedron));
		}
	}
}
/// class splitter
k3d::split_edge* splitter::split_edge(k3d::split_edge& Edge, const double factor)
{
	double sharpness = get_sharpness(Edge);
	double comp_sharpness = Edge.companion ? get_sharpness(*(Edge.companion)) : 0.0;

	// first point of the edge:
	k3d::point* this_point = Edge.vertex;
	// second point of the edge:
	k3d::point* other_point = Edge.face_clockwise->vertex;
	std::cerr << debug << "Starting split_edge with [" << this_point->position[0] << ", " << this_point->position[1] << ", " << this_point->position[2] << "] and [" << other_point->position[0] << ", " << other_point->position[1] << ", " << other_point->position[2] << "]" << std::endl;
	double f = factor;
	// new point position:
	k3d::vector3 new_pos;
	k3d::split_edge* ret_edge = 0;
	point_map::iterator it = m_other_points.find(this_point);
	if(it != m_other_points.end()) {
		this_point = (other_point == it->second.first) ? it->second.second : it->second.first;
		ret_edge = &Edge;
	} else {
		it = m_other_points.find(other_point);
		if(it != m_other_points.end()) {
			other_point = (this_point == it->second.first) ? it->second.second : it->second.first;
			ret_edge = Edge.face_clockwise;
		}
	}
	std::cerr << debug << "Finishing split_edge with [" << this_point->position[0] << ", " << this_point->position[1] << ", " << this_point->position[2] << "] and [" << other_point->position[0] << ", " << other_point->position[1] << ", " << other_point->position[2] << "]" << std::endl;
	new_pos = k3d::mix<k3d::vector3>(this_point->position, other_point->position, f);
	if(Edge.companion && m_split_edges.find(Edge.companion) != m_split_edges.end()) {
		std::cerr << debug << "Mixing old and new split-point." << std::endl;
		it->first->position = k3d::mix<k3d::vector3>(new_pos, it->first->position, 0.5);
		return ret_edge;
	}
	std::cerr << debug << "Split edge [" << Edge.vertex->position[0] << ", " << Edge.vertex->position[1] << ", " << Edge.vertex->position[2] << "]->[" << other_point->position[0] << ", " << other_point->position[1] << ", " << other_point->position[2] << "] at [" << new_pos[0] << ", " << new_pos[1] << ", " << new_pos[2] << "]" << std::endl;
	// new point:
	std::cerr << debug << "Adding new point [" << new_pos[0] << ", " << new_pos[1] << ", " << new_pos[2] << "]" << std::endl;
	k3d::point* new_point = new k3d::point(new_pos);
	m_Points.push_back(new_point);
	// new edge:
	k3d::split_edge* new_edge = new k3d::split_edge(new_point, Edge.face_clockwise);
	Edge.face_clockwise = new_edge;
	new_edge->tags["crease"] = sharpness;
	// companion edge
	if(Edge.companion) {
		k3d::split_edge* new_comp_edge = new k3d::split_edge(new_point, Edge.companion->face_clockwise);
		Edge.companion->face_clockwise = new_comp_edge;
		k3d::join_edges(*Edge.companion, *new_edge);
		k3d::join_edges(Edge, *new_comp_edge);
		m_Polyhedron.edges.push_back(new_comp_edge);
		new_comp_edge->tags["crease"] = comp_sharpness;
	}
	m_Polyhedron.edges.push_back(new_edge);
	m_other_points[new_edge->vertex] = std::pair<k3d::point*, k3d::point*>(this_point, other_point);
	//m_split_edges.insert(Edge.companion);
	return new_edge;
}

k3d::split_edge* splitter::detach_edge_novertex(k3d::split_edge& Edge, const k3d::vector3& centroid)
{
	std::cerr << debug << "detach_edge_novertex" << std::endl;
	k3d::split_edge* clockwise = Edge.face_clockwise;
	//double f = adapt_factor(*clockwise, m_f);
	double f = 2*sqrt(m_f * m_parallel_edges[clockwise]);
	f = f > 1.0 ? 1.0 : f;
	std::cerr << debug << "Moving [" << clockwise->vertex->position[0] << ", " << clockwise->vertex->position[1] << ", " << clockwise->vertex->position[2] << "] to [";
	//clockwise->vertex->position = k3d::mix(clockwise->vertex->position, clockwise->face_clockwise->vertex->position, f);
	clockwise->vertex->position = k3d::mix(Edge.vertex->position, centroid, f);
	m_other_points.erase(clockwise->vertex);
	std::cerr << debug << clockwise->vertex->position[0] << ", " << clockwise->vertex->position[1] << ", " << clockwise->vertex->position[2] << "]" << std::endl;
	m_last_factor[clockwise] = f;
	clockwise->companion->face_clockwise->vertex = Edge.vertex;
	m_last_factor.erase(clockwise->companion->face_clockwise);
	m_split_edges.erase(clockwise->companion->face_clockwise);
	if(clockwise->companion->face_clockwise->companion) {
		clockwise->companion->face_clockwise->companion->face_clockwise = clockwise->companion->face_clockwise->companion->face_clockwise->face_clockwise;
		m_last_factor.erase(clockwise->companion->face_clockwise->companion);
		m_split_edges.erase(clockwise->companion->face_clockwise->companion);
	}
	k3d::split_edge* companion = Edge.companion;
	if(!companion) {
		companion = new k3d::split_edge(clockwise->vertex, clockwise->companion->face_clockwise);
		k3d::join_edges(*companion, Edge);
		m_Polyhedron.edges.push_back(companion);
	}
	companion->face_clockwise = clockwise->companion->face_clockwise;
	clockwise->companion->face_clockwise = companion;
	edge_face_map::iterator it = m_face_starts.find(Edge.companion);
	if(it != m_face_starts.end()) {
		k3d::face* removed_face = it->second;
		removed_face->first_edge = companion->face_clockwise->companion;
		m_face_starts.erase(it);
		m_face_starts[companion->face_clockwise->companion] = removed_face;
	}

	return companion;
}

// Detach the starting point of a parallel split face, where the detached vertex belongs to Edge.
k3d::split_edge* splitter::detach_edge_vertex(k3d::split_edge& Edge, const k3d::vector3& centroid)
{
	std::cerr << debug << "detach_edge_vertex" << std::endl;
	k3d::split_edge* c_clockwise = counter_clockwise(Edge);
	k3d::split_edge* other_edge = counter_clockwise(*(c_clockwise->companion));
	//double f = adapt_factor(*c_clockwise, m_f);
	double f = 2*sqrt(m_f * m_parallel_edges[c_clockwise]);
	f = f > 1.0 ? 1.0 : f;
	std::cerr << debug << "Moving [" << Edge.vertex->position[0] << ", " << Edge.vertex->position[1] << ", " << Edge.vertex->position[2] << "] to [";
	//Edge.vertex->position = k3d::mix(Edge.vertex->position, c_clockwise->vertex->position, f);
	Edge.vertex->position = k3d::mix(Edge.face_clockwise->vertex->position, centroid, f);
	m_other_points.erase(Edge.vertex);
	std::cerr << debug << Edge.vertex->position[0] << ", " << Edge.vertex->position[1] << ", " << Edge.vertex->position[2] << "]" << std::endl;
	m_last_factor[c_clockwise] = f;
	k3d::split_edge* companion = Edge.companion;
	if(companion) {
		k3d::split_edge* companion_ccl = counter_clockwise(*companion);
		companion->face_clockwise->vertex = Edge.face_clockwise->vertex;
		m_last_factor.erase(companion->face_clockwise);
		m_split_edges.erase(companion->face_clockwise);
		companion_ccl->face_clockwise = companion->face_clockwise;
	} else {
		companion = new k3d::split_edge(Edge.face_clockwise->vertex, c_clockwise->companion);
		k3d::join_edges(Edge, *companion);
		m_Polyhedron.edges.push_back(companion);
	}
	companion->face_clockwise = c_clockwise->companion;
	other_edge->face_clockwise = companion;
	m_last_factor.erase(other_edge);
	m_split_edges.erase(other_edge);
	edge_face_map::iterator it = m_face_starts.find(Edge.companion);
	if(it != m_face_starts.end()) {
		k3d::face* removed_face = it->second;
		removed_face->first_edge = other_edge->companion;
		m_face_starts.erase(it);
		m_face_starts[other_edge->companion] = removed_face;
	}
	return companion->face_clockwise;
}

k3d::split_edge* splitter::split_face_along_edge(k3d::split_edge& Edge, const double factor, const k3d::vector3& centroid)
{
	m_f = factor;
	// get the two neigbour edges
	k3d::split_edge* clockwise = Edge.face_clockwise;
	k3d::split_edge* c_clockwise = counter_clockwise(Edge);
	// Check if the clockwise edge is new
	bool is_new = false;
	bool is_new_companion = false;
	bool is_new_cc = false;
	k3d::split_edge* other_end = 0;
	if(m_parallel_edges.find(clockwise) != m_parallel_edges.end())
		is_new = true;
	if(m_parallel_edges.find(clockwise->companion) != m_parallel_edges.end())
		is_new_companion = true;
	if(Edge.face_clockwise->companion && Edge.face_clockwise->companion->face_clockwise->face_clockwise->companion &&  (m_parallel_edges.find(Edge.face_clockwise->companion->face_clockwise->face_clockwise->companion) != m_parallel_edges.end())) {
		is_new_cc = true;
		other_end = Edge.face_clockwise->companion->face_clockwise->face_clockwise->companion->face_clockwise;
	}
	k3d::split_edge* new_ccl_edge = 0;
	k3d::split_edge* new_cl_edge = 0;
	if(is_new || is_new_companion || is_new_cc) {
		if(is_new) {
			std::cerr << debug << "is_new" << std::endl;
			new_ccl_edge = detach_edge_novertex(Edge, centroid);
			if(!is_new_cc) {
				new_cl_edge = this->split_edge(*(Edge.companion->face_clockwise->face_clockwise), factor);
				m_last_factor[new_cl_edge] = factor;
				m_split_edges.insert(Edge.companion->face_clockwise->face_clockwise);
			}
		}
		if(is_new_companion) {
			std::cerr << debug << "is_new_companion" << std::endl;
			new_cl_edge = detach_edge_vertex(*(Edge.face_clockwise->companion->face_clockwise), centroid);
			if(!is_new_cc) {
				new_ccl_edge = this->split_edge(*c_clockwise, 1-factor);
				m_last_factor[c_clockwise] = factor;
				m_split_edges.insert(new_ccl_edge);
				if(c_clockwise == new_ccl_edge)
					c_clockwise = counter_clockwise(*c_clockwise);
			}
		}
		if(is_new_cc) {
			std::cerr << debug << "is_new_cc" << std::endl;
			new_cl_edge = detach_edge_vertex(*other_end, centroid);
		}
		clockwise = counter_clockwise(*new_cl_edge);
		c_clockwise = counter_clockwise(*new_ccl_edge);
	} else {
		// split neighbour edges at the desired length
		new_ccl_edge = this->split_edge(*c_clockwise, 1-factor);
		m_last_factor[c_clockwise] = factor;
		m_split_edges.insert(new_ccl_edge);
		if(c_clockwise == new_ccl_edge)
			c_clockwise = counter_clockwise(*c_clockwise);
		new_cl_edge = this->split_edge(*clockwise, factor);
		m_last_factor[new_cl_edge] = factor;
		m_split_edges.insert(clockwise);
	}
	// create new edges
	k3d::split_edge* new_edge = new k3d::split_edge(new_ccl_edge->vertex, new_cl_edge);
	k3d::split_edge* new_edge_companion = new k3d::split_edge(new_cl_edge->vertex, new_ccl_edge);
	// add them to the list
	m_Polyhedron.edges.push_back(new_edge);
	m_Polyhedron.edges.push_back(new_edge_companion);
	// join them
	k3d::join_edges(*new_edge, *new_edge_companion);
	// fix the other face_clockwise pointers
	clockwise->face_clockwise = new_edge_companion;
	c_clockwise->face_clockwise = new_edge;
	// add new faces
	if(!is_face(*new_edge)) {
		add_face(*new_edge);
	} else {
		add_face(*new_edge_companion);
	}
	m_parallel_edges[new_edge_companion] = factor;
	if(is_new)
		m_last_factor[new_edge_companion] = m_parallel_edges[(counter_clockwise(*new_edge))->companion];
	if(is_new_companion)
		m_last_factor[new_edge_companion] = m_parallel_edges[new_edge->face_clockwise->companion];
	std::cerr << debug << "Created parallel edge [" << new_edge_companion->vertex->position[0] << ", " << new_edge_companion->vertex->position[1] << ", " << new_edge_companion->vertex->position[2] << "]->[" << new_edge->vertex->position[0] << ", " << new_edge->vertex->position[1] << ", " << new_edge->vertex->position[2] << "]" << std::endl;
	return new_edge_companion;
}

k3d::split_edge* splitter::counter_clockwise(k3d::split_edge& Edge)
{
	k3d::split_edge* last_edge = &Edge;
	k3d::split_edge* this_edge = Edge.face_clockwise;

	while(this_edge != &Edge)
		{
			last_edge = this_edge;
			this_edge = this_edge->face_clockwise;
		}

	return last_edge;
}

bool splitter::is_face(k3d::split_edge& Edge)
{
	k3d::split_edge* last_edge = &Edge;
	k3d::split_edge* this_edge = Edge.face_clockwise;

	while(this_edge != &Edge)
		{
			if(m_face_starts.find(last_edge) != m_face_starts.end())
				return true;

			last_edge = this_edge;
			this_edge = this_edge->face_clockwise;
		}

	return (m_face_starts.find(last_edge) != m_face_starts.end());
}

void splitter::add_face(k3d::split_edge& Edge)
{
	k3d::face* new_face = new k3d::face(&Edge);
	m_Polyhedron.faces.push_back(new_face);
	m_face_starts[&Edge] = new_face;
}

double splitter::adapt_factor(k3d::split_edge& Edge, const double factor)
{
	double f = factor;
	double f2 = -1.0;
	edge_double_map::iterator it = m_last_factor.find(&Edge);

	if(it != m_last_factor.end())
		f2 = it->second;
	else
		{
			it = m_last_factor.find(Edge.companion);
			if(it != m_last_factor.end())
				f2 = (it->second);
		}

	if(f2 != -1.0)
		{
			if(f < 0.5)
				f /= 1-f2;
			else
				f = 1 - f2/f;
		}

	if(f > 1.0)
		f = 1.0;
	std::cerr << debug << "Factor " << f << " instead of " << factor << "." << std::endl;

	return f;
}

double get_sharpness(const k3d::split_edge& Edge)
{
	for(k3d::parameters_t::const_iterator tag = Edge.tags.begin(); tag != Edge.tags.end(); ++tag)
	{
		if(tag->first == "crease" && (tag->second.type() == typeid(k3d::ri::real)))
		{
			std::cerr << debug << "Found sharpness " << boost::any_cast<double>(tag->second) << " on [" << Edge.vertex->position[0] << ", " << Edge.vertex->position[1] << ", " << Edge.vertex->position[2] << "]->[" << Edge.face_clockwise->vertex->position[0] << ", " << Edge.face_clockwise->vertex->position[1] << ", " << Edge.face_clockwise->vertex->position[2] << "]" << std::endl;
			return boost::any_cast<double>(tag->second);
		}
	}

	std::cerr << debug << "Found no sharpness on [" << Edge.vertex->position[0] << ", " << Edge.vertex->position[1] << ", " << Edge.vertex->position[2] << "]->[" << Edge.face_clockwise->vertex->position[0] << ", " << Edge.face_clockwise->vertex->position[1] << ", " << Edge.face_clockwise->vertex->position[2] << "]" << std::endl;
	return 0.0;
}

k3d::vector3 get_centroid(const k3d::split_edge& Edge)
{
	k3d::vector3 centroid = Edge.vertex->position;
	unsigned long edges = 1;
	const k3d::split_edge* last_edge = &Edge;
	const k3d::split_edge* this_edge = Edge.face_clockwise;
	while(this_edge != &Edge)
		{
			++edges;
			centroid += this_edge->vertex->position;
			last_edge = this_edge;
			this_edge = this_edge->face_clockwise;
		}

	return centroid / edges;
}

void crease(const k3d::mesh& Input, k3d::mesh& Output)
{
	k3d::deep_copy(Input, Output);
	for(k3d::mesh::polyhedra_t::iterator it = Output.polyhedra.begin(); it != Output.polyhedra.end(); ++it)
	{
		//if((*it)->type == k3d::polyhedron::CATMULL_CLARK_SUBDIVISION_MESH)
		k3d::polyhedron& Polyhedron = **it;
		return_if_fail(k3d::is_valid(Polyhedron));
		splitter Splitter(Polyhedron, Output.points);
		unsigned long len = Polyhedron.faces.size();
		for(unsigned long i = 0; i < len; ++i)
		{
			std::vector<k3d::split_edge*> edges;
			k3d::split_edge* first = Polyhedron.faces[i]->first_edge;
			edges.push_back(first);
			k3d::split_edge* next = first->face_clockwise;
			while(next != first)
				{
					edges.push_back(next);
					next = next->face_clockwise;
				}

			std::cerr << debug << "Face has " << edges.size() << " edges." << std::endl;
			unsigned long elen = edges.size();
			k3d::vector3 centroid(0,0,0);

			for(unsigned long j = 0; j < elen; ++j)
				centroid += edges[j]->vertex->position;

			centroid /= elen;

			for(unsigned long j = 0; j < elen; ++j)
			{
				k3d::split_edge& edge = *(edges[j]);
				double s = get_sharpness(edge);
				if(s > 0)
				{
					Splitter.split_face_along_edge(edge,splitter::sharpness_to_factor(s), centroid);
					edge.tags["crease"] = 0.0;
				}
			}
		}

		for(unsigned long j = 0; j < Polyhedron.edges.size(); ++j)
			Polyhedron.edges[j]->tags.erase("crease");

		return_if_fail(k3d::is_valid(Polyhedron));
	}
}

} // namespace subdiv