/*! gts_vertex_gaussian_curvature:
 *  @v: a #WVertex.
 *  @s: a #GtsSurface.
 *  @Kg: the Discrete Gaussian Curvature approximation at @v.
 *
 *  Computes the Discrete Gaussian Curvature approximation at @v.
 *
 *  This approximation is from the paper:
 *      Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
 *      Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr
 *      VisMath '02, Berlin (Germany)
 *      http://www-grail.usc.edu/pubs.html
 *
 *  Returns: %true if the operator could be evaluated, %false if the evaluation failed for some reason (@v is
 * boundary or is the endpoint of a non-manifold edge.)
 */
bool gts_vertex_gaussian_curvature(WVertex *v, real *Kg)
{
	real area = 0.0;
	real angle_sum = 0.0;

	if (!v)
		return false;
	if (!Kg)
		return false;

	/* this operator is not defined for boundary edges */
	if (v->isBoundary()) {
		*Kg = 0.0;
		return false;
	}

	WVertex::incoming_edge_iterator itE;
	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
		area += (*itE)->GetaFace()->getArea();
	}

	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
		WOEdge *e = (*itE)->getPrevOnFace();
		WVertex *v1 = e->GetaVertex();
		WVertex *v2 = e->GetbVertex();
		angle_sum += angle_from_cotan(v, v1, v2);
	}

	*Kg = (2.0 * M_PI - angle_sum) / area;

	return true;
}

bool AdvanceToEdge(WFace * face, Vec3r currentPoint, int edgeIndex, Vec3r & nextPoint)
{
  Vec3r viewpoint = SilhouetteGeomEngine::GetViewpoint();
  
  WOEdge * edge = face->GetOEdge(edgeIndex);

  Vec3r A = edge->GetaVertex()->GetVertex();
  Vec3r B = edge->GetbVertex()->GetVertex();

  Vec3r edgedir = B-A;
  Vec3r planeNormal = (viewpoint - currentPoint) ^ face->GetNormal();
  real t;
  GeomUtils::intersection_test result = GeomUtils::intersectLinePlanePN(A, edgedir, planeNormal, currentPoint, t);

  //  printf("\te = %d, result = %d, want: %d, t = %f\n", edgeIndex, result, GeomUtils::DO_INTERSECT, t);
  if (result != GeomUtils::DO_INTERSECT || t < 0 || t > 1)
    return false;

  nextPoint = A + t*edgedir;

  return true;
  /*
  real dp = (nextPoint - currentPoint) * (viewpoint - currentPoint) ;

  printf("\tdp = %f\n", dp);

  return dp > 0;
  */
}

static bool angle_obtuse(WVertex *v, WFace *f)
{
	WOEdge *e;
	f->getOppositeEdge(v, e);

	Vec3r vec1(e->GetaVertex()->GetVertex() - v->GetVertex());
	Vec3r vec2(e->GetbVertex()->GetVertex() - v->GetVertex());
	return ((vec1 * vec2) < 0);
}

WVertex::incoming_edge_iterator WVertex::incoming_edges_end()
{
	WOEdge *begin;
	WOEdge *aOEdge = _EdgeList.front()->GetaOEdge();
	if (aOEdge->GetbVertex() == this)
		begin = aOEdge;
	else
		begin = _EdgeList.front()->GetbOEdge();
	return incoming_edge_iterator(this, begin, 0);
}

void WVertex::incoming_edge_iterator::increment()
{
	WOEdge *twin = _current->twin();
	if (!twin) {
		// we reached a hole
		_current = 0;
		return;
	}
	WOEdge *next = twin->getPrevOnFace();
	if (next == _begin) {
		next = NULL;
	}
	_current = next;
}

/*! gts_vertex_mean_curvature_normal:
 *  @v: a #WVertex.
 *  @s: a #GtsSurface.
 *  @Kh: the Mean Curvature Normal at @v.
 *
 *  Computes the Discrete Mean Curvature Normal approximation at @v.
 *  The mean curvature at @v is half the magnitude of the vector @Kh.
 *
 *  Note: the normal computed is not unit length, and may point either into or out of the surface, depending on
 *  the curvature at @v. It is the responsibility of the caller of the function to use the mean curvature normal
 *  appropriately.
 *
 *  This approximation is from the paper:
 *      Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
 *      Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr
 *      VisMath '02, Berlin (Germany)
 *      http://www-grail.usc.edu/pubs.html
 *
 *  Returns: %true if the operator could be evaluated, %false if the evaluation failed for some reason (@v is
 *  boundary or is the endpoint of a non-manifold edge.)
 */
bool gts_vertex_mean_curvature_normal(WVertex *v, Vec3r &Kh)
{
	real area = 0.0;

	if (!v)
		return false;

	/* this operator is not defined for boundary edges */
	if (v->isBoundary())
		return false;

	WVertex::incoming_edge_iterator itE;

	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
		area += (*itE)->GetaFace()->getArea();

	Kh = Vec3r(0.0, 0.0, 0.0);

	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
		WOEdge *e = (*itE)->getPrevOnFace();
#if 0
		if ((e->GetaVertex() == v) || (e->GetbVertex() == v))
			cerr<< "BUG ";
#endif
		WVertex *v1 = e->GetaVertex();
		WVertex *v2 = e->GetbVertex();
		real temp;

		temp = cotan(v1, v, v2);
		Kh = Vec3r(Kh + temp * (v2->GetVertex() - v->GetVertex()));

		temp = cotan(v2, v, v1);
		Kh = Vec3r(Kh + temp * (v1->GetVertex() - v->GetVertex()));
	}
	if (area > 0.0) {
		Kh[0] /= 2 * area;
		Kh[1] /= 2 * area;
		Kh[2] /= 2 * area;
	}
	else {
		return false;
	}

	return true;
}

void compute_curvature_tensor_one_ring(WVertex *start, NormalCycle& nc)
{
	// in case we have a non-manifold vertex, skip it...
	if (start->isBoundary())
		return;

	WVertex::incoming_edge_iterator woeit = start->incoming_edges_begin();
	WVertex::incoming_edge_iterator woeitend = start->incoming_edges_end();
	for (; woeit != woeitend; ++woeit) {
		WOEdge *h = (*woeit)->twin();
		nc.accumulate_dihedral_angle(h->GetVec(), h->GetAngle());
		WOEdge *hprev = h->getPrevOnFace();
		nc.accumulate_dihedral_angle(hprev->GetVec(), hprev->GetAngle());
	}
}

// TODO: check optimizations:
// use marking ? (measure *timings* ...)
void compute_curvature_tensor(WVertex *start, real radius, NormalCycle& nc)
{
	// in case we have a non-manifold vertex, skip it...
	if (start->isBoundary())
		return;

	std::set<WVertex*> vertices;
	const Vec3r& O = start->GetVertex();
	std::stack<WVertex*> S;
	S.push(start);
	vertices.insert(start);
	while (!S.empty()) {
		WVertex *v = S.top();
		S.pop();
		if (v->isBoundary())
			continue;
		const Vec3r& P = v->GetVertex();
		WVertex::incoming_edge_iterator woeit = v->incoming_edges_begin();
		WVertex::incoming_edge_iterator woeitend = v->incoming_edges_end();
		for (; woeit != woeitend; ++woeit) {
			WOEdge *h = *woeit;
			if ((v == start) || h->GetVec() * (O - P) > 0.0) {
				Vec3r V(-1 * h->GetVec());
				bool isect = sphere_clip_vector(O, radius, P, V);
				assert (h->GetOwner()->GetNumberOfOEdges() == 2); // Because otherwise v->isBoundary() would be true
				nc.accumulate_dihedral_angle(V, h->GetAngle());

				if (!isect) {
					WVertex *w = h->GetaVertex();
					if (vertices.find(w) == vertices.end()) {
						vertices.insert(w);
						S.push(w);
					}
				}
			}
		}
	}
}

WFace *WShape::MakeFace(vector<WVertex *>& iVertexList, vector<bool>& iFaceEdgeMarksList, unsigned iMaterial,
                        WFace *face)
{
	int id = _FaceList.size();

	face->setFrsMaterialIndex(iMaterial);

	// Check whether we have a degenerated face:

	// LET'S HACK IT FOR THE TRIANGLE CASE:

	if (3 == iVertexList.size()) {
		if ((iVertexList[0] == iVertexList[1]) ||
		    (iVertexList[0] == iVertexList[2]) ||
		    (iVertexList[2] == iVertexList[1]))
		{
			cerr << "Warning: degenerated triangle detected, correcting" << endl;
			return NULL;
		}
	}

	vector<WVertex *>::iterator it;

	// compute the face normal (v1v2 ^ v1v3)
	WVertex *v1, *v2, *v3;
	it = iVertexList.begin();
	v1 = *it;
	it++;
	v2 = *it;
	it++;
	v3 = *it;

	Vec3r vector1(v2->GetVertex() - v1->GetVertex());
	Vec3r vector2(v3->GetVertex() - v1->GetVertex());

	Vec3r normal(vector1 ^ vector2);
	normal.normalize();
	face->setNormal(normal);

	vector<bool>::iterator mit = iFaceEdgeMarksList.begin();
	face->setMark(*mit);
	mit++;

	// vertex pointers used to build each edge
	vector<WVertex *>::iterator va, vb;

	va = iVertexList.begin();
	vb = va;
	for (; va != iVertexList.end(); va = vb) {
		++vb;
		// Adds va to the vertex list:
		//face->AddVertex(*va);

		WOEdge *oedge;
		if (*va == iVertexList.back())
			oedge = face->MakeEdge(*va, iVertexList.front()); //for the last (closing) edge
		else
			oedge = face->MakeEdge(*va, *vb);

		if (!oedge)
			return NULL;

		WEdge *edge = oedge->GetOwner();
		if (1 == edge->GetNumberOfOEdges()) {
			// means that we just created a new edge and that we must add it to the shape's edges list
			edge->setId(_EdgeList.size());
			AddEdge(edge);
			// compute the mean edge value:
			_meanEdgeSize += edge->GetaOEdge()->GetVec().norm();
		}

		edge->setMark(*mit);
		++mit;
	}

	// Add the face to the shape's faces list:
	face->setId(id);
	AddFace(face);

	return face;
}

WShape::WShape(WShape& iBrother)
{
	_Id = iBrother.GetId();
	_Name = iBrother._Name;
	_FrsMaterials = iBrother._FrsMaterials;
	_meanEdgeSize = iBrother._meanEdgeSize;
	iBrother.bbox(_min, _max);
	vector<WVertex *>& vertexList = iBrother.getVertexList();
	vector<WVertex *>::iterator v = vertexList.begin(), vend = vertexList.end();
	for (; v != vend; ++v) {
		//WVertex *newVertex = new WVertex(*(*v));
		WVertex *newVertex = (*v)->duplicate();

		newVertex->setShape(this);
		AddVertex(newVertex);
	}

	vector<WEdge *>& edgeList = iBrother.getEdgeList();
	vector<WEdge *>::iterator e = edgeList.begin(), eend = edgeList.end();
	for (; e != eend; ++e) {
		//WEdge *newEdge = new WEdge(*(*e));
		WEdge *newEdge = (*e)->duplicate();
		AddEdge(newEdge);
	}

	vector<WFace *>& faceList = iBrother.GetFaceList();
	vector<WFace *>::iterator f = faceList.begin(), fend = faceList.end();
	for (; f != fend; ++f) {
		//WFace *newFace = new WFace(*(*f));
		WFace *newFace = (*f)->duplicate();
		AddFace(newFace);
	}

	// update all pointed addresses thanks to the newly created objects:
	vend = _VertexList.end();
	for (v = _VertexList.begin(); v != vend; ++v) {
		const vector<WEdge *>& vedgeList = (*v)->GetEdges();
		vector<WEdge *> newvedgelist;
		unsigned int i;
		for (i = 0; i < vedgeList.size(); i++) {
			WEdge *current = vedgeList[i];
			edgedata *currentvedata = (edgedata *)current->userdata;
			newvedgelist.push_back(currentvedata->_copy);
		}
		(*v)->setEdges(newvedgelist);
	}

	eend = _EdgeList.end();
	for (e = _EdgeList.begin(); e != eend; ++e) {
		// update aOedge:
		WOEdge *aoEdge = (*e)->GetaOEdge();
		aoEdge->setaVertex(((vertexdata *)(aoEdge->GetaVertex()->userdata))->_copy);
		aoEdge->setbVertex(((vertexdata *)(aoEdge->GetbVertex()->userdata))->_copy);
		if (aoEdge->GetaFace())
			aoEdge->setaFace(((facedata *)(aoEdge->GetaFace()->userdata))->_copy);
		aoEdge->setbFace(((facedata *)(aoEdge->GetbFace()->userdata))->_copy);
		aoEdge->setOwner(((edgedata *)(aoEdge->GetOwner()->userdata))->_copy);

		// update bOedge:
		WOEdge *boEdge = (*e)->GetbOEdge();
		if (boEdge) {
			boEdge->setaVertex(((vertexdata *)(boEdge->GetaVertex()->userdata))->_copy);
			boEdge->setbVertex(((vertexdata *)(boEdge->GetbVertex()->userdata))->_copy);
			if (boEdge->GetaFace())
				boEdge->setaFace(((facedata *)(boEdge->GetaFace()->userdata))->_copy);
			boEdge->setbFace(((facedata *)(boEdge->GetbFace()->userdata))->_copy);
			boEdge->setOwner(((edgedata *)(boEdge->GetOwner()->userdata))->_copy);
		}
	}

	fend = _FaceList.end();
	for (f = _FaceList.begin(); f != fend; ++f) {
		unsigned int i;
		const vector<WOEdge *>& oedgeList = (*f)->getEdgeList();
		vector<WOEdge *> newoedgelist;

		unsigned int n = oedgeList.size();
		for (i = 0; i < n; i++) {
			WOEdge *current = oedgeList[i];
			oedgedata *currentoedata = (oedgedata *)current->userdata;
			newoedgelist.push_back(currentoedata->_copy);
			//oedgeList[i] = currentoedata->_copy;
			//oedgeList[i] = ((oedgedata *)(oedgeList[i]->userdata))->_copy;
		}
		(*f)->setEdgeList(newoedgelist);
	}

	// Free all memory (arghh!)
	// Vertex
	vend = iBrother.getVertexList().end();
	for (v = iBrother.getVertexList().begin(); v != vend; ++v) {
		delete (vertexdata *)((*v)->userdata);
		(*v)->userdata = NULL;
	}

	// Edges and OEdges:
	eend = iBrother.getEdgeList().end();
	for (e = iBrother.getEdgeList().begin(); e != eend; ++e) {
		delete (edgedata *)((*e)->userdata);
		(*e)->userdata = NULL;
//.........这里部分代码省略.........

WOEdge *WFace::MakeEdge(WVertex *v1, WVertex *v2)
{
	// First check whether the same oriented edge already exists or not:
	vector<WEdge *>& v1Edges = v1->GetEdges();
	for (vector<WEdge*>::iterator it1 = v1Edges.begin(), end = v1Edges.end(); it1 != end; it1++) {
		WEdge *we = (*it1);
		WOEdge *woea = we->GetaOEdge();

		//if ((*it1)->GetbVertex() == v2) {
		if ((woea->GetaVertex() == v1) && (woea->GetbVertex() == v2)) {
			// The oriented edge already exists
			cerr << "Warning: edge " << v1->GetId() << " - " << v2->GetId() << " appears twice, correcting" << endl;
			// Adds the edge to the face
			//AddEdge((*it1)->GetaOEdge());
			AddEdge(woea);
			(*it1)->setNumberOfOEdges((*it1)->GetNumberOfOEdges() + 1);
			//sets these vertices as border:
			v1->setBorder(true);
			v2->setBorder(true);
			//return (*it1)->GetaOEdge();
			return woea;
		}

		WOEdge *woeb = we->GetbOEdge();
		//if ((*it1)->GetbVertex() == v2)
		if (woeb && (woeb->GetaVertex() == v1) && (woeb->GetbVertex() == v2)) {
			// The oriented edge already exists
			cerr << "Warning: edge " << v1->GetId() << " - " << v2->GetId() << " appears twice, correcting" << endl;
			// Adds the edge to the face
			//AddEdge((*it1)->GetaOEdge());
			AddEdge(woeb);
			(*it1)->setNumberOfOEdges((*it1)->GetNumberOfOEdges() + 1);
			//sets these vertices as border:
			v1->setBorder(true);
			v2->setBorder(true);
			//return (*it1)->GetaOEdge();
			return woeb;
		}
	}

	// the oriented edge we're about to build
	WOEdge *pOEdge = new WOEdge;
	// The edge containing the oriented edge.
	WEdge *edge;

	// checks whether this edge already exists or not
	// If it exists, it points outward v2
	bool exist = false;
	WOEdge *pInvertEdge = NULL; // The inverted edge if it exists
	vector<WEdge *>& v2Edges = v2->GetEdges();
	vector<WEdge *>::iterator it;
	for (it = v2Edges.begin(); it != v2Edges.end(); it++) {
		if ((*it)->GetbVertex() == v1) {
			// The invert edge already exists
			exist = true;
			pInvertEdge = (*it)->GetaOEdge();
			break;
		}
	}

	//DEBUG:
	if (true == exist) { // The invert edge already exists
		// Retrieves the corresponding edge
		edge = pInvertEdge->GetOwner();

		// Sets the a Face (retrieved from pInvertEdge
		pOEdge->setaFace(pInvertEdge->GetbFace());

		// Updates the invert edge:
		pInvertEdge->setaFace(this);
	}
	else { // The invert edge does not exist yet
		// we must create a new edge
		//edge = new WEdge;
		edge = instanciateEdge();

		// updates the a,b vertex edges list:
		v1->AddEdge(edge);
		v2->AddEdge(edge);
	}

	pOEdge->setOwner(edge);
	// Add the vertices:
	pOEdge->setaVertex(v1);
	pOEdge->setbVertex(v2);

	// Debug:
	if (v1->GetId() == v2->GetId())
		cerr << "Warning: edge " << this << " null with vertex " << v1->GetId() << endl;

	edge->AddOEdge(pOEdge);
	//edge->setNumberOfOEdges(edge->GetNumberOfOEdges() + 1);

	// Add this face (the b face)
	pOEdge->setbFace(this);

	// Adds the edge to the face
	AddEdge(pOEdge);

	return pOEdge;
//.........这里部分代码省略.........

void FEdgeXDetector::ProcessRidgeFace(WXFace *iFace)
{
	// RIDGE LAYER
	// Compute the RidgeFunction, that is the derivative of the ppal curvature along e1 at each vertex of the face
	WVertex *v;
	Vec3r v1v2;
	real t;
	vector<WXFaceLayer*> SmoothLayers;
	WXFaceLayer *faceLayer;
	Face_Curvature_Info *layer_info;
	real K1_a(0), K1_b(0);
	Vec3r Inter_a, Inter_b;

	// find the ridge layer of the face
	iFace->retrieveSmoothLayers(Nature::RIDGE, SmoothLayers);
	if ( SmoothLayers.size()!=1 )
		return;
	faceLayer = SmoothLayers[0];
	// retrieve the curvature info of this layer
	layer_info = (Face_Curvature_Info *)faceLayer->userdata;

	int numVertices = iFace->numberOfVertices();
	for (int i = 0; i < numVertices; i++) {
		v = iFace->GetVertex(i);
		// vec_curvature_info[i] contains the curvature info of this vertex
		Vec3r e2 = layer_info->vec_curvature_info[i]->K2*layer_info->vec_curvature_info[i]->e2;
		Vec3r e1 = layer_info->vec_curvature_info[i]->K1*layer_info->vec_curvature_info[i]->e1;
		e2.normalize();

		WVertex::face_iterator fit = v->faces_begin();
		WVertex::face_iterator fitend = v->faces_end();
		for (; fit != fitend; ++fit) {
			WXFace *wxf = dynamic_cast<WXFace*>(*fit);
			WOEdge *oppositeEdge;
			if (!(wxf->getOppositeEdge(v, oppositeEdge)))
				continue;
			v1v2 = oppositeEdge->GetbVertex()->GetVertex() - oppositeEdge->GetaVertex()->GetVertex();
			GeomUtils::intersection_test res;
			res = GeomUtils::intersectRayPlane(oppositeEdge->GetaVertex()->GetVertex(), v1v2, e2, -(v->GetVertex()*e2),
			                                   t, 1.0e-06);
			if ((res == GeomUtils::DO_INTERSECT) && (t >= 0.0) && (t <= 1.0)) {
				vector<WXFaceLayer*> second_ridge_layer;
				wxf->retrieveSmoothLayers(Nature::RIDGE, second_ridge_layer);
				if (second_ridge_layer.size() != 1)
					continue;
				Face_Curvature_Info *second_layer_info = (Face_Curvature_Info*)second_ridge_layer[0]->userdata;

				unsigned index1 = wxf->GetIndex(oppositeEdge->GetaVertex());
				unsigned index2 = wxf->GetIndex(oppositeEdge->GetbVertex());
				real K1_1 = second_layer_info->vec_curvature_info[index1]->K1;
				real K1_2 = second_layer_info->vec_curvature_info[index2]->K1;
				real K1 = (1.0 - t) * K1_1 + t * K1_2;
				Vec3r inter((1.0 - t) * oppositeEdge->GetaVertex()->GetVertex() +
				            t * oppositeEdge->GetbVertex()->GetVertex());
				Vec3r vtmp(inter - v->GetVertex());
				// is it K1_a or K1_b ?
				if (vtmp * e1 > 0) {
					K1_b = K1;
					Inter_b = inter;
				}
				else {
					K1_a = K1;
					Inter_a = inter;
				}
			}
		}
		// Once we have K1 along the ppal direction compute the derivative : K1b - K1a put it in DotP
		//real d = fabs(K1_b) - fabs(K1_a);
		real d = 0;
		real threshold = _meanK1 + (_maxK1 - _meanK1) / 7.0;
		//real threshold = _meanK1;
		//if ((fabs(K1_b) > threshold) || ((fabs(K1_a) > threshold)))
		d = (K1_b) - (K1_a) / (Inter_b - Inter_a).norm();
		faceLayer->PushDotP(d);
		//faceLayer->PushDotP(layer_info->vec_curvature_info[i]->K1);
	}

	// Make the values relevant by checking whether all principal directions have the "same" direction:
	Vec3r e0((layer_info->vec_curvature_info[0]->K1 * layer_info->vec_curvature_info[0]->e1));
	e0.normalize();
	Vec3r e1((layer_info->vec_curvature_info[1]->K1 * layer_info->vec_curvature_info[1]->e1));
	e1.normalize();
	Vec3r e2((layer_info->vec_curvature_info[2]->K1 * layer_info->vec_curvature_info[2]->e1));
	e2.normalize();
	if (e0 * e1 < 0)
		// invert dotP[1]
		faceLayer->ReplaceDotP(1, -faceLayer->dotP(1));
	if (e0 * e2 < 0)
		// invert dotP[2]
		faceLayer->ReplaceDotP(2, -faceLayer->dotP(2));

#if 0 // remove the weakest values;
	real minDiff = (_maxK1 - _minK1) / 10.0;
	real minDiff = _meanK1;
	if ((faceLayer->dotP(0) < minDiff) && (faceLayer->dotP(1) < minDiff) && (faceLayer->dotP(2) < minDiff)) {
		faceLayer->ReplaceDotP(0, 0);
		faceLayer->ReplaceDotP(1, 0);
		faceLayer->ReplaceDotP(2, 0);
	}
#endif
//.........这里部分代码省略.........

/*! gts_vertex_principal_directions:
 *  @v: a #WVertex.
 *  @s: a #GtsSurface.
 *  @Kh: mean curvature normal (a #Vec3r).
 *  @Kg: Gaussian curvature (a real).
 *  @e1: first principal curvature direction (direction of largest curvature).
 *  @e2: second principal curvature direction.
 *
 *  Computes the principal curvature directions at a point given @Kh and @Kg, the mean curvature normal and
 *  Gaussian curvatures at that point, computed with gts_vertex_mean_curvature_normal() and
 *  gts_vertex_gaussian_curvature(), respectively.
 *
 *  Note that this computation is very approximate and tends to be unstable. Smoothing of the surface or the principal
 *  directions may be necessary to achieve reasonable results.
 */
void gts_vertex_principal_directions(WVertex *v, Vec3r Kh, real Kg, Vec3r &e1, Vec3r &e2)
{
	Vec3r N;
	real normKh;

	Vec3r basis1, basis2, d, eig;
	real ve2, vdotN;
	real aterm_da, bterm_da, cterm_da, const_da;
	real aterm_db, bterm_db, cterm_db, const_db;
	real a, b, c;
	real K1, K2;
	real *weights, *kappas, *d1s, *d2s;
	int edge_count;
	real err_e1, err_e2;
	int e;
	WVertex::incoming_edge_iterator itE;

	/* compute unit normal */
	normKh = Kh.norm();

	if (normKh > 0.0) {
		Kh.normalize();
	}
	else {
		/* This vertex is a point of zero mean curvature (flat or saddle point). Compute a normal by averaging
		 * the adjacent triangles
		 */
		N[0] = N[1] = N[2] = 0.0;

		for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++)
			N = Vec3r(N + (*itE)->GetaFace()->GetNormal());
		real normN = N.norm();
		if (normN <= 0.0)
			return;
		N.normalize();
	}

	/* construct a basis from N: */
	/* set basis1 to any component not the largest of N */
	basis1[0] =  basis1[1] =  basis1[2] = 0.0;
	if (fabs (N[0]) > fabs (N[1]))
		basis1[1] = 1.0;
	else
		basis1[0] = 1.0;

	/* make basis2 orthogonal to N */
	basis2 = (N ^ basis1);
	basis2.normalize();

	/* make basis1 orthogonal to N and basis2 */
	basis1 = (N ^ basis2);
	basis1.normalize();

	aterm_da = bterm_da = cterm_da = const_da = 0.0;
	aterm_db = bterm_db = cterm_db = const_db = 0.0;
	int nb_edges = v->GetEdges().size();

	weights = (real *)malloc(sizeof(real) * nb_edges);
	kappas = (real *)malloc(sizeof(real) * nb_edges);
	d1s = (real *)malloc(sizeof(real) * nb_edges);
	d2s = (real *)malloc(sizeof(real) * nb_edges);
	edge_count = 0;

	for (itE = v->incoming_edges_begin(); itE != v->incoming_edges_end(); itE++) {
		WOEdge *e;
		WFace *f1, *f2;
		real weight, kappa, d1, d2;
		Vec3r vec_edge;
		if (!*itE)
			continue;
		e = *itE;

		/* since this vertex passed the tests in gts_vertex_mean_curvature_normal(), this should be true. */
		//g_assert(gts_edge_face_number (e, s) == 2);

		/* identify the two triangles bordering e in s */
		f1 = e->GetaFace();
		f2 = e->GetbFace();

		/* We are solving for the values of the curvature tensor
		 *     B = [ a b ; b c ].
		 *  The computations here are from section 5 of [Meyer et al 2002].
		 *
		 *  The first step is to calculate the linear equations governing the values of (a,b,c). These can be computed
		 *  by setting the derivatives of the error E to zero (section 5.3).
//.........这里部分代码省略.........

OWXFaceLayer ViewEdgeXBuilder::FindPreviousFaceLayer(const OWXFaceLayer& iFaceLayer)
{
	WXFace *previousFace = NULL;
	WOEdge *woebegin;
	real tend;
	if (iFaceLayer.order) {
		woebegin = iFaceLayer.fl->getSmoothEdge()->woea();
		tend = iFaceLayer.fl->getSmoothEdge()->ta();
	}
	else {
		woebegin = iFaceLayer.fl->getSmoothEdge()->woeb();
		tend = iFaceLayer.fl->getSmoothEdge()->tb();
	}

	// special case of EDGE_VERTEX config:
	if ((tend == 0.0) || (tend == 1.0)) {
		WVertex *previousVertex;
		if (tend == 0.0)
			previousVertex = woebegin->GetaVertex();
		else
			previousVertex = woebegin->GetbVertex();
		if (previousVertex->isBoundary()) // if it's a non-manifold vertex -> ignore
			return OWXFaceLayer(NULL, true);
		bool found = false;
		WVertex::face_iterator f = previousVertex->faces_begin();
		WVertex::face_iterator fend = previousVertex->faces_end();
		for (; (!found) && (f != fend); ++f) {
			previousFace = dynamic_cast<WXFace*>(*f);
			if ((0 != previousFace) && (previousFace != iFaceLayer.fl->getFace())) {
				vector<WXFaceLayer*> sameNatureLayers;
				previousFace->retrieveSmoothEdgesLayers(iFaceLayer.fl->nature(), sameNatureLayers);
				// don't know... Maybe should test whether this face has also a vertex_edge configuration
				if (sameNatureLayers.size() == 1) {
					WXFaceLayer *winner = sameNatureLayers[0];
					// check face mark continuity
					if (winner->getFace()->GetMark() != iFaceLayer.fl->getFace()->GetMark())
						return OWXFaceLayer(NULL, true);
					if (woebegin == winner->getSmoothEdge()->woeb()->twin())
						return OWXFaceLayer(winner, true);
					else
						return OWXFaceLayer(winner, false);
				}
			}
		}
	}
	else {
		previousFace = dynamic_cast<WXFace*>(iFaceLayer.fl->getFace()->GetBordingFace(woebegin));
		if (0 == previousFace)
			return OWXFaceLayer(NULL, true);
		// if the next face layer has either no smooth edge or no smooth edge of same nature, no next face
		if (!previousFace->hasSmoothEdges())
			return OWXFaceLayer(NULL, true);
		vector<WXFaceLayer*> sameNatureLayers;
		previousFace->retrieveSmoothEdgesLayers(iFaceLayer.fl->nature(), sameNatureLayers);
		// don't know how to deal with several edges of same nature on a single face
		if ((sameNatureLayers.empty()) || (sameNatureLayers.size() != 1)) {
			return OWXFaceLayer(NULL, true);
		}
		else {
			WXFaceLayer *winner = sameNatureLayers[0];
			// check face mark continuity
			if (winner->getFace()->GetMark() != iFaceLayer.fl->getFace()->GetMark())
				return OWXFaceLayer(NULL, true);
			if (woebegin == winner->getSmoothEdge()->woeb()->twin())
				return OWXFaceLayer(winner, true);
			else
				return OWXFaceLayer(winner, false);
		}
	}
	return OWXFaceLayer(NULL, true);
}