/* Copyright (C) 1999-2004 by Peter Eastman

   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. */

package artofillusion;

import artofillusion.math.*;
import artofillusion.object.*;
import artofillusion.object.TriangleMesh.*;
import artofillusion.ui.*;
import buoy.event.*;
import buoy.widget.*;

/** TriMeshSimplifier defines methods for simplifying a TriangleMesh.  The algorithm
   is based on
   <p>
   R. Ronfard and J. Rossignac, "Full-range Approximation of Triangulated Polyhedra."
   Computer Graphics Forum, vol. 15, no. 3, 1996.
   <p>
   I have made the following modifications to their algorithm:
   <p>
   1. When calculating the cost for an edge collapse, they arbitrarily choose which direction
   to collapse it in.  This results in meshes not being simplified as far as they could be.
   Instead, I always test both collapse directions to see which one has the lower cost.
   <p>
   2. After collapsing vertex v1 to v2, they state that it is only necessary to recalculate
   costs for edges which have been modified, that is, those edges in the star of v1 or v2.
   This is incorrect.  It is also necessary to recalculate the cost for any edge which shares
   a vertex with an edge in v1's star.
   <p>
   3. They define the local tesselation error for a face as being proportional to its angle
   of rotation.  I replace the angle by 1-cos(angle).  This has the same effect, and saves
   the need to evalate an acos.
   <p>
   4. Rather than always simplifying the entire mesh, I allow the user to select only a
   subset of the mesh to be simplified.
   <p>
   5. Their algorithm only works for closed meshes (meshes without a boundary).  I have
   modified it to also work for open meshes.
   <p>
   6. I omit the relaxation step, since this would produce bad results for smoothed
   meshes. */

import artofillusion.ui.*;
public class TriMeshSimplifier extends Thread
{
  class Constraint
  {
    double a, b, c, d;
  }

  class Zone
  {
    int constraints;
    Constraint con[];
    Zone next;
  }

  class MeshEdge
  {
    int v1, v2;
    MeshFace f1, f2;
    float smoothness;
    double cost;
    boolean selected;
  }

  class MeshFace
  {
    int v1, v2, v3, index;
    MeshEdge e1, e2, e3;
    Vec3 normal;
  }

  class VertexInfo
  {
    Vertex theVert;
    Vec3 pos;
    int faces, edges;
    float smoothness;
    MeshEdge star[];
    MeshFace crown[];
    Zone zone;
  }

  private TriangleMesh mesh;
  private VertexInfo vertex[];
  private MeshEdge edge[];
  private MeshFace face[];
  private Vec3 temp1, temp2, temp3;
  private int faces;
  private BDialog dial;
  private BButton cancelButton;
  private BLabel status, numLabel;
  private boolean cancel;
  private double tol;

  public TriMeshSimplifier(TriangleMesh theMesh, boolean selection[], double tolerance, BFrame fr)
  {
    mesh = theMesh;
    tol = tolerance*tolerance;
    try
    {
      buildDataStructures(selection);
    }
    catch(ArrayIndexOutOfBoundsException e)
    {
      if (fr != null)
      {
        new MessageDialog(fr, new String [] {
          "Sorry, the mesh simplifier could not handle this mesh."}, new String [] {Translate.text("button.ok")});
      }
      return;
    }

    if (fr == null)
      start();
    else
    {
      createDialog(fr);
      start();
      Thread update = new Thread() {
        public void run()
        {
          try
          {
            while (!cancel)
            {
              sleep(200);
              status.setText(String.valueOf(faces));
            }
          }
          catch (InterruptedException ex)
          {
          }
        }
      };
      update.start();
      dial.setVisible(true);
    }
  }

  public void run()
  {
    doSimplification();
    if (!cancel)
      buildMesh();
    cancel = true;
    if (numLabel != null)
    {
      numLabel.setText(Translate.text("Final")+":");
      status.setText(String.valueOf(faces));
      cancelButton.setText(Translate.text("Done"));
    }
  }

  /* Put up a status dialog. */

  void createDialog(BFrame fr)
  {
    dial = new BDialog(fr, true);
    BorderContainer content = new BorderContainer();
    dial.setContent(BOutline.createEmptyBorder(content, ModellingApp.standardDialogInsets));
    content.add(Translate.label("numTrianglesInMesh"), BorderContainer.NORTH);
    GridContainer grid = new GridContainer(2, 2);
    grid.add(new BLabel(Translate.text("Original")+":"), 0, 0);
    grid.add(numLabel = new BLabel(Translate.text("Current")+":"), 0, 1);
    grid.add(new BLabel(String.valueOf(faces)), 1, 0);
    grid.add(status = new BLabel(String.valueOf(faces)), 1, 1);
    content.add(grid, BorderContainer.CENTER);
    content.add(cancelButton = Translate.button("cancel", new Object() {
      void processEvent()
      {
        cancel = true;
        dial.dispose();
      }
    }, "processEvent"), BorderContainer.SOUTH, new LayoutInfo());
    dial.pack();
    dial.setResizable(false);
    UIUtilities.centerWindow(dial);
  }

  /* Build the new TriangleMesh. */

  void buildMesh()
  {
    int i, j, k, index[] = new int [vertex.length], f[][] = new int [faces][];
    Vertex vert[];
    Edge ed[];

    // Find the indices of all unique vertices.

    for (i = 0, k = 0; i < vertex.length; i++)
      {
	index[i] = k++;
	for (j = 0; j < i; j++)
	  if (vertex[i] == vertex[j])
	    {
	      index[i] = index[j];
	      k--;
	      break;
	    }
      }

    // Update the mesh shape.

    vert = new Vertex [k];
    for (i = 0; i < vertex.length; i++)
      {
	vert[index[i]] = new TriangleMesh.Vertex(vertex[i].theVert);
        vert[index[i]].r = new Vec3(vertex[i].pos);
	vert[index[i]].smoothness = vertex[i].smoothness;
	vert[index[i]].edges = 0;
	vert[index[i]].firstEdge = -1;
      }
    for (i = 0; i < faces; i++)
      f[i] = new int [] {index[face[i].v1], index[face[i].v2], index[face[i].v3]};
    mesh.setShape(vert, f);

    // Set the smoothness values for edges.

    ed = mesh.getEdges();
    for (i = 0; i < ed.length; i++)
      for (j = 0; j < edge.length; j++)
	if ((index[edge[j].v1] == ed[i].v1 && index[edge[j].v2] == ed[i].v2) ||
		(index[edge[j].v2] == ed[i].v1 && index[edge[j].v1] == ed[i].v2))
	  {
	    ed[i].smoothness = edge[j].smoothness;
	    break;
	  }
  }

  /* Build the data structures necessary for simplifying the mesh. */

  void buildDataStructures(boolean selection[])
  {
    TriangleMesh.Vertex v[] = (Vertex []) mesh.getVertices();
    TriangleMesh.Edge e[] = mesh.getEdges();
    TriangleMesh.Face f[] = mesh.getFaces();
    int i, j, edgeCount[], faceCount[];
    MeshEdge tempEdge;
    VertexInfo v1, v2;
    Constraint con;
    double cost;

    vertex = new VertexInfo [v.length];
    edge = new MeshEdge [e.length];
    face = new MeshFace [f.length];
    edgeCount = new int [v.length];
    faceCount = new int [v.length];
    faces = f.length;

    // Create temporary vectors for later use.

    temp1 = new Vec3();
    temp2 = new Vec3();
    temp3 = new Vec3();

    // Count the edges and faces touching each vertex, and build the list of vertices.

    for (i = 0; i < e.length; i++)
      {
	edgeCount[e[i].v1]++;
	edgeCount[e[i].v2]++;
      }
    for (i = 0; i < f.length; i++)
      {
	faceCount[f[i].v1]++;
	faceCount[f[i].v2]++;
	faceCount[f[i].v3]++;
      }
    for (i = 0; i < v.length; i++)
      {
	vertex[i] = new VertexInfo();
	vertex[i].theVert = v[i];
	vertex[i].pos = v[i].r;
	vertex[i].smoothness = v[i].smoothness;
	vertex[i].star = new MeshEdge [edgeCount[i]];
	vertex[i].crown = new MeshFace [faceCount[i]];
	vertex[i].zone = new Zone();
	if (faceCount[i] == edgeCount[i])
	  vertex[i].zone.con = new Constraint [faceCount[i]];
	else
	  vertex[i].zone.con = new Constraint [faceCount[i]+2];
      }

    // Build the list of edges.

    for (i = 0; i < e.length; i++)
      {
	edge[i] = new MeshEdge();
	edge[i].v1 = e[i].v1;
	edge[i].v2 = e[i].v2;
	edge[i].smoothness = e[i].smoothness;
	v1 = vertex[e[i].v1];
	v2 = vertex[e[i].v2];
        v1.star[v1.edges++] = edge[i];
	v2.star[v2.edges++] = edge[i];
        edge[i].selected = selection[i];
      }

    // Build the list of faces, and record the constraint associated with each one.

    for (i = 0; i < f.length; i++)
      {
	face[i] = new MeshFace();
	face[i].v1 = f[i].v1;
	face[i].v2 = f[i].v2;
	face[i].v3 = f[i].v3;
	face[i].e1 = edge[f[i].e1];
	face[i].e2 = edge[f[i].e2];
	face[i].e3 = edge[f[i].e3];
	face[i].index = i;
	face[i].normal = new Vec3();
	findNormal(vertex[f[i].v1].pos, vertex[f[i].v2].pos, vertex[f[i].v3].pos, face[i].normal);
	vertex[f[i].v1].crown[vertex[f[i].v1].faces++] = face[i];
	vertex[f[i].v2].crown[vertex[f[i].v2].faces++] = face[i];
	vertex[f[i].v3].crown[vertex[f[i].v3].faces++] = face[i];
	con = new Constraint();
	con.a = face[i].normal.x;
	con.b = face[i].normal.y;
	con.c = face[i].normal.z;
	con.d = -(con.a*v[f[i].v1].r.x + con.b*v[f[i].v1].r.y + con.c*v[f[i].v1].r.z);
	vertex[f[i].v1].zone.con[vertex[f[i].v1].zone.constraints++] = con;
	vertex[f[i].v2].zone.con[vertex[f[i].v2].zone.constraints++] = con;
	vertex[f[i].v3].zone.con[vertex[f[i].v3].zone.constraints++] = con;
	if (edge[f[i].e1].f1 == null)
	  edge[f[i].e1].f1 = face[i];
	else
	  edge[f[i].e1].f2 = face[i];
	if (edge[f[i].e2].f1 == null)
	  edge[f[i].e2].f1 = face[i];
	else
	  edge[f[i].e2].f2 = face[i];
	if (edge[f[i].e3].f1 == null)
	  edge[f[i].e3].f1 = face[i];
	else
	  edge[f[i].e3].f2 = face[i];
      }

    // Record additional constraints for boundary edges.

    for (i = 0; i < e.length; i++)
      if (edge[i].f2 == null)
	{
	  temp1.set(v[edge[i].v1].r);
	  temp1.subtract(v[edge[i].v2].r);
	  temp2.set(edge[i].f1.normal);
	  con = new Constraint();
	  con.a = temp1.y*temp2.z - temp1.z*temp2.y;
	  con.b = temp1.z*temp2.x - temp1.x*temp2.z;
	  con.c = temp1.x*temp2.y - temp1.y*temp2.x;
	  con.d = -(con.a*v[edge[i].v1].r.x + con.b*v[edge[i].v1].r.y + con.c*v[edge[i].v1].r.z);
	  vertex[edge[i].v1].zone.con[vertex[edge[i].v1].zone.constraints++] = con;
	  vertex[edge[i].v2].zone.con[vertex[edge[i].v2].zone.constraints++] = con;
	}

    // Record the initial cost associated with each edge, and bring the lowest cost edge
    // to the front of the list.

    for (i = 0; i < e.length; i++)
      updateCost(edge[i]);
    j = 0;
    for (i = 0; i < e.length; i++)
      if (edge[i].cost < edge[j].cost)
	j = i;
    tempEdge = edge[j];
    edge[j] = edge[0];
    edge[0] = tempEdge;
  }

  /* Simplify the mesh! */

  void doSimplification()
  {
    int i, j, k, skip;
    MeshEdge rem1, rem2, rep1, rep2;
    MeshEdge e, star[], tempEdge;
    MeshFace tempCrown[], face1, face2;
    VertexInfo v1, v2, v3, v4, tempVert;
    Zone zone1, zone2;
    boolean boundary;

    for (i = 0; i < edge.length; i += skip)
      {
	e = edge[i];
	if (e.cost >= tol || cancel)
	  return;
	v1 = vertex[e.v1];
	v2 = vertex[e.v2];
	boundary = (e.f2 == null);

	// Move v1 to v2.

        v2.pos = v1.pos;

	// Flag all the edges whose costs must be updated.  This includes all edges in v1's
	// star or v2's star, and all edges which share a vertex with an edge in v1's star.

	star = v1.star;
	for (j = 0; j < v1.edges; j++)
	  {
	    if (vertex[star[j].v1] == v1)
	      tempVert = vertex[star[j].v2];
	    else
	      tempVert = vertex[star[j].v1];
	    for (k = 0; k < tempVert.edges; k++)
	      tempVert.star[k].cost = -1.0;
	  }
	star = v2.star;
	for (j = 0; j < v2.edges; j++)
	  star[j].cost = -1.0;

	// Find the two extra edges to remove, and the ones to replace them with.

	if (e.f1.e1 == e)
	  {
	    if (vertex[e.f1.e2.v1] == v2 || vertex[e.f1.e2.v2] == v2)
	      {
		rem1 = e.f1.e2;
		rep1 = e.f1.e3;
	      }
	    else
	      {
		rem1 = e.f1.e3;
		rep1 = e.f1.e2;
	      }
	  }
	else if (e.f1.e2 == e)
	  {
	    if (vertex[e.f1.e3.v1] == v2 || vertex[e.f1.e3.v2] == v2)
	      {
		rem1 = e.f1.e3;
		rep1 = e.f1.e1;
	      }
	    else
	      {
		rem1 = e.f1.e1;
		rep1 = e.f1.e3;
	      }
	  }
	else
	  {
	    if (vertex[e.f1.e1.v1] == v2 || vertex[e.f1.e1.v2] == v2)
	      {
		rem1 = e.f1.e1;
		rep1 = e.f1.e2;
	      }
	    else
	      {
		rem1 = e.f1.e2;
		rep1 = e.f1.e1;
	      }
	  }
	if (boundary)
	  rem2 = rep2 = null;
	else
	  {
	    if (e.f2.e1 == e)
	      {
		if (vertex[e.f2.e2.v1] == v2 || vertex[e.f2.e2.v2] == v2)
		  {
		    rem2 = e.f2.e2;
		    rep2 = e.f2.e3;
		  }
		else
		  {
		    rem2 = e.f2.e3;
		    rep2 = e.f2.e2;
		  }
	      }
	    else if (e.f2.e2 == e)
	      {
		if (vertex[e.f2.e3.v1] == v2 || vertex[e.f2.e3.v2] == v2)
		  {
		    rem2 = e.f2.e3;
		    rep2 = e.f2.e1;
		  }
		else
		  {
		    rem2 = e.f2.e1;
		    rep2 = e.f2.e3;
		  }
	      }
	    else
	      {
		if (vertex[e.f2.e1.v1] == v2 || vertex[e.f2.e1.v2] == v2)
		  {
		    rem2 = e.f2.e1;
		    rep2 = e.f2.e2;
		  }
		else
		  {
		    rem2 = e.f2.e2;
		    rep2 = e.f2.e1;
		  }
	      }
	  }

	// Identify the other two points whose stars and crowns must be modified.

	if (vertex[rep1.v1] == v1)
	  v3 = vertex[rep1.v2];
	else
	  v3 = vertex[rep1.v1];
	if (boundary)
	  v4 = null;
	else
	  {
	    if (vertex[rep2.v1] == v1)
	      v4 = vertex[rep2.v2];
	    else
	      v4 = vertex[rep2.v1];
	  }

	// Update v1's crown.

	if (boundary && v1.faces+v2.faces-2 > v1.crown.length)
	  tempCrown = new MeshFace [v1.faces+v2.faces-2];
	else if (v1.faces+v2.faces-4 > v1.crown.length)
	  tempCrown = new MeshFace [v1.faces+v2.faces-4];
	else
	  tempCrown = v1.crown;
	for (j = 0, k = 0; k < v1.faces; k++)
	  if (v1.crown[k] != e.f1 && v1.crown[k] != e.f2)
	    tempCrown[j++] = v1.crown[k];
	for (k = 0; k < v2.faces; k++)
	  if (v2.crown[k] != e.f1 && v2.crown[k] != e.f2)
	    tempCrown[j++] = v2.crown[k];
	v1.crown = tempCrown;
	if (boundary)
	  v1.faces = v1.faces+v2.faces-2;
	else
	  v1.faces = v1.faces+v2.faces-4;

	// Update v1's star.

	if (v1.edges+v2.edges-2 > v1.star.length)
	  star = new MeshEdge [v1.edges+v2.edges-2];
	else
	  star = v1.star;
	for (j = 0, k = 0; k < v1.edges; k++)
	  if (v1.star[k] != e)
	    star[j++] = v1.star[k];
	for (k = 0; k < v2.edges; k++)
	  if (v2.star[k] != e)
	    star[j++] = v2.star[k];
	v1.star = star;
	v1.edges = v1.edges+v2.edges-2;

	// Update v3 and v4's crowns.

	for (j = 0, k = 0; k < v3.faces; k++)
	  if (v3.crown[k] != e.f1)
	    v3.crown[j++] = v3.crown[k];
	v3.faces--;
	if (!boundary)
	  {
	    for (j = 0, k = 0; k < v4.faces; k++)
	      if (v4.crown[k] != e.f2)
		v4.crown[j++] = v4.crown[k];
	    v4.faces--;
	  }

	// Update v3 and v4's stars.

	for (j = 0, k = 0; k < v3.edges; k++)
	  if (v3.star[k] != rem1)
	    v3.star[j++] = v3.star[k];
	v3.edges--;
	if (!boundary)
	  {
	    for (j = 0, k = 0; k < v4.edges; k++)
	      if (v4.star[k] != rem2)
		v4.star[j++] = v4.star[k];
	    v4.edges--;
	  }

	// Eliminate the faces which share this edge.

	if (boundary)
	  {
	    j = e.f1.index;
	    face[j] = face[faces-1];
	    face[j].index = j;
	    faces--;
	  }
	else
	  {
	    j = e.f1.index;
	    k = e.f2.index;
	    if (j < faces-2 && k < faces-2)
	      {
		face[j] = face[faces-1];
		face[k] = face[faces-2];
	      }
	    else
	      {
		if (j != faces-1 && k != faces-1)
		  face[Math.min(j,k)] = face[faces-1];
		else
		  face[Math.min(j,k)] = face[faces-2];
	      }
	    face[j].index = j;
	    face[k].index = k;
	    faces -= 2;
	  }

	// Add the constraints from v2 to v1.  We do this by removing any duplicate
	// constraints, then appending v2's zone to v1's zone.

	zone2 = v2.zone;
	while (zone2 != null)
	  {
	    zone1 = v1.zone;
	    while (zone1 != null)
	      {
		for (j = 0; j < zone2.constraints; j++)
		  for (k = 0; k < zone1.constraints; k++)
		    if (zone2.con[j] == zone1.con[k])
		      zone2.con[j] = zone2.con[--zone2.constraints];
		zone1 = zone1.next;
	      }
	    zone2 = zone2.next;
	  }

	// Eliminate any zones which are completely free of constraints.

	zone2 = v2.zone;
	while (zone2.next != null)
	  {
	    if (zone2.next.constraints == 0)
	      zone2.next = zone2.next.next;
	    else
	      zone2 = zone2.next;
	  }

	// Add v2's zone to v1's zone.

	zone2 = v2.zone;
	zone1 = v1.zone;
	while (zone1.next != null)
	  zone1 = zone1.next;
	if (zone2.constraints == 0)
	  zone1.next = zone2.next;
	else
	  zone1.next = zone2;

	// Go through the faces in v2's crown, and replace references to rem1 and rem2
	// by rep1 and rep2.

	for (j = 0; j < v2.faces; j++)
	  {
	    if (v2.crown[j].e1 == rem1)
	      v2.crown[j].e1 = rep1;
	    else if (v2.crown[j].e2 == rem1)
	      v2.crown[j].e2 = rep1;
	    else if (v2.crown[j].e3 == rem1)
	      v2.crown[j].e3 = rep1;
	  }
	if (!boundary)
	  for (j = 0; j < v2.faces; j++)
	    {
	      if (v2.crown[j].e1 == rem2)
		v2.crown[j].e1 = rep2;
	      else if (v2.crown[j].e2 == rem2)
	 	v2.crown[j].e2 = rep2;
	      else if (v2.crown[j].e3 == rem2)
		v2.crown[j].e3 = rep2;
	    }

	// Similarly, go through the edges adjacent to the faces being removed, and
	// replace references.

	if (rem1.f1 == e.f1)
	  {
	    if (rep1.f1 == e.f1)
	      rep1.f1 = rem1.f2;
	    else
	      rep1.f2 = rem1.f2;
	  }
	else
	  {
	    if (rep1.f1 == e.f1)
	      rep1.f1 = rem1.f1;
	    else
	      rep1.f2 = rem1.f1;
	  }
	if (rep1.f1 == null)
	  {
	    rep1.f1 = rep1.f2;
	    rep1.f2 = null;
	  }
	if (!boundary)
	  {
	    if (rem2.f1 == e.f2)
	      {
		if (rep2.f1 == e.f2)
		  rep2.f1 = rem2.f2;
		else
		  rep2.f2 = rem2.f2;
	      }
	    else
	      {
		if (rep2.f1 == e.f2)
		  rep2.f1 = rem2.f1;
		else
		  rep2.f2 = rem2.f1;
	      }
	    if (rep2.f1 == null)
	      {
		rep2.f1 = rep2.f2;
		rep2.f2 = null;
	      }
	  }

	// Replace all references to v2 by v1.

	for (j = 0; j < vertex.length; j++)
	  if (vertex[j] == v2)
	    vertex[j] = v1;

	// Update smoothness values for vertices and edges.

	v1.smoothness = Math.min(v1.smoothness, v2.smoothness);
	rep1.smoothness = Math.min(rep1.smoothness, rem1.smoothness);
	if (!boundary)
	  rep2.smoothness = Math.min(rep2.smoothness, rem2.smoothness);

	// Remove rem1 and rem2 from the list, then resort it.

	for (j = edge.length-1, k = edge.length-1; j > i; j--)
	  if (edge[j] != rem1 && edge[j] != rem2)
	    edge[k--] = edge[j];

	// Update the cost for edges which need it.

	if (boundary)
	  skip = 2;
	else
	  skip = 3;
	k = i+skip;
	for (j = i+skip; j < edge.length; j++)
	  {
	    if (edge[j].cost == -1.0)
	      updateCost(edge[j]);
	    if (edge[j].cost < edge[k].cost)
	      k = j;
	  }
	tempEdge = edge[i+skip];
	edge[i+skip] = edge[k];
	edge[k] = tempEdge;
      }
  }

  /* Update the cost for an edge.  Try both directions to see which gives a lower cost. */

  void updateCost(MeshEdge ed)
  {
    double cost = findCost(ed);
    int i = ed.v1;

    ed.v1 = ed.v2;
    ed.v2 = i;
    ed.cost = findCost(ed);
    if (ed.cost > cost)
      {
  	ed.v2 = ed.v1;
	ed.v1 = i;
  	ed.cost = cost;
      }
  }

  /* Find the cost associated with contracting a given edge.  If the cost is determined
     to be greater than tol, the method simply returns tol.  This allows findCost() to
     return as soon as it determines the edge will never be contracted, and also speeds
     up the sorting of edges (since many edges will have identical costs). */

  double findCost(MeshEdge ed)
  {
    VertexInfo v1 = vertex[ed.v1], v2 = vertex[ed.v2];
    MeshFace f;
    Zone zone = v2.zone;
    double cost, max = 0.0;
    Vec3 p;
    int i;

    if (!ed.selected)
      return tol;

    // First calculate the local tesselation error.

    for (i = 0; i < v2.faces; i++)
      if (v2.crown[i] != ed.f1 && v2.crown[i] != ed.f2)
	{
	  f = v2.crown[i];
	  if (vertex[f.v1] == v2)
	    findNormal(v1.pos, vertex[f.v2].pos, vertex[f.v3].pos, temp3);
	  else if (vertex[f.v2] == v2)
	    findNormal(vertex[f.v1].pos, v1.pos, vertex[f.v3].pos, temp3);
	  else
	    findNormal(vertex[f.v1].pos, vertex[f.v2].pos, v1.pos, temp3);
	  cost = tol*(1.0-temp3.dot(f.normal));
	  if (cost >= tol || Double.isNaN(cost))
	    return tol;
	  if (cost > max)
	    max = cost;
	}

    // Now calculate the local geometric error.

    while (zone != null)
      {
	for (i = 0; i < zone.constraints; i++)
	  {
	    cost = zone.con[i].a*v1.pos.x + zone.con[i].b*v1.pos.y + zone.con[i].c*v1.pos.z + zone.con[i].d;
	    cost *= cost;
	    if (cost >= tol)
	      return tol;
	    if (cost > max)
	      max = cost;
	  }
	zone = zone.next;
      }

    // If this edge is adjacent to a face whose other edges are both boundaries, then
    // collapsing this edge would leave a dangling edge.  This should not be allowed.

    f = ed.f1;
    if ((f.e1 == ed || f.e1.f2 == null) && (f.e2 == ed || f.e2.f2 == null) && (f.e3 == ed || f.e3.f2 == null))
      return tol;
    f = ed.f2;
    if (f != null)
      if ((f.e1 == ed || f.e1.f2 == null) && (f.e2 == ed || f.e2.f2 == null) && (f.e3 == ed || f.e3.f2 == null))
	return tol;
    return max;
  }

  /* Given three points which define a face, find the normal vector. */

  void findNormal(Vec3 v1, Vec3 v2, Vec3 v3, Vec3 normal)
  {
    temp1.set(v2.x, v2.y, v2.z);
    temp1.subtract(v1);
    temp2.set(v3.x, v3.y, v3.z);
    temp2.subtract(v1);
    normal.set(temp1.y*temp2.z-temp1.z*temp2.y, temp1.z*temp2.x-temp1.x*temp2.z, temp1.x*temp2.y-temp1.y*temp2.x);
    normal.normalize();
  }
}