/* jCAE stand for Java Computer Aided Engineering. Features are : Small CAD
modeler, Finite element mesher, Plugin architecture.
Copyright (C) 2003,2005, by EADS CRC
Copyright (C) 2007,2008,2010, by EADS France
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
package org.jcae.mesh.amibe.algos3d;
import org.jcae.mesh.amibe.ds.Mesh;
import org.jcae.mesh.amibe.ds.Triangle;
import org.jcae.mesh.amibe.ds.AbstractHalfEdge;
import org.jcae.mesh.amibe.ds.Vertex;
import org.jcae.mesh.amibe.projection.QuadricProjection;
import org.jcae.mesh.amibe.projection.LocalSurfaceProjection;
import org.jcae.mesh.amibe.util.QSortedTree;
import org.jcae.mesh.amibe.util.PAVLSortedTree;
import org.jcae.mesh.xmldata.MeshReader;
import org.jcae.mesh.xmldata.MeshWriter;
import java.util.Map;
import java.util.HashMap;
import java.util.LinkedHashSet;
import java.util.Collection;
import java.util.Iterator;
import java.io.IOException;
import gnu.trove.map.hash.TObjectDoubleHashMap;
import java.util.logging.Level;
import java.util.logging.Logger;
/**
* Node smoothing. Triangle quality is computed for all triangles,
* and vertex quality is the lowest value of its incident triangles.
* Vertices are sorted according to their quality, and processed
* iteratively by beginning with worst vertex. A modified Laplacian
* smoothing is performed, as briefly explained in
* <a href="http://www.ann.jussieu.fr/~frey/publications/ijnme4198.pdf">Adaptive Triangular-Quadrilateral Mesh Generation</a>, by Houman Borouchaky and
* Pascal J. Frey.
* If final position improves vertex quality, point is moved.
*/
public class SmoothNodes3D
{
private static final Logger LOGGER=Logger.getLogger(SmoothNodes3D.class.getName());
private final Mesh mesh;
private double sizeTarget = -1.0;
private int nloop = 10;
private double tolerance = Double.MAX_VALUE / 2.0;
private boolean preserveBoundaries = false;
private boolean checkQuality = true;
private int progressBarStatus = 10000;
private static final double scaleFactor = 12.0 * Math.sqrt(3.0);
private double relaxation = 0.6;
private final Vertex c;
private final QSortedTree<Vertex> tree = new PAVLSortedTree<Vertex>();
private boolean refresh = false;
private int processed = 0;
private int notProcessed = 0;
private TObjectDoubleHashMap<Triangle> qualityMap;
private Map<Vertex, LocalSurfaceProjection> nodeProjection;
private Collection<Vertex> nodeset;
/**
* Creates a <code>SmoothNodes3D</code> instance.
*
* @param m the <code>Mesh</code> instance to refine.
*/
public SmoothNodes3D(Mesh m)
{
this(m, new HashMap<String, String>());
}
/**
* Creates a <code>SmoothNodes3D</code> instance.
*
* @param m the <code>Mesh</code> instance to refine.
* @param options map containing key-value pairs to modify algorithm
* behaviour. Valid keys are <code>size</code>,
* <code>iterations</code>, <code>boundaries</code>,
* <code>tolerance</code>, <code>refresh</code> and
* <code>relaxation</code>.
*/
public SmoothNodes3D(final Mesh m, final Map<String, String> options)
{
mesh = m;
c = mesh.createVertex(0.0, 0.0, 0.0);
for (final Map.Entry<String, String> opt: options.entrySet())
{
final String key = opt.getKey();
final String val = opt.getValue();
if (key.equals("size"))
sizeTarget = Double.valueOf(val).doubleValue();
else if (key.equals("iterations"))
nloop = Integer.valueOf(val).intValue();
else if (key.equals("boundaries"))
preserveBoundaries = Boolean.valueOf(val).booleanValue();
else if (key.equals("tolerance"))
tolerance = Double.valueOf(val).doubleValue();
else if (key.equals("refresh"))
refresh = true;
else if (key.equals("check"))
checkQuality = Boolean.valueOf(val).booleanValue();
else if (key.equals("relaxation"))
relaxation = Double.valueOf(val).doubleValue();
else
throw new RuntimeException("Unknown option: "+key);
}
if (LOGGER.isLoggable(Level.FINE))
{
if (sizeTarget > 0.0)
LOGGER.fine("Size: "+sizeTarget);
LOGGER.fine("Iterations: "+nloop);
LOGGER.fine("Refresh: "+refresh);
LOGGER.fine("Relaxation: "+relaxation);
LOGGER.fine("Tolerance: "+tolerance);
LOGGER.fine("Preserve boundaries: "+preserveBoundaries);
}
}
public void setProgressBarStatus(int n)
{
progressBarStatus = n;
}
/**
* Moves all nodes until all iterations are done.
*/
private void computeTriangleQuality()
{
AbstractHalfEdge ot = null;
for (Triangle f: mesh.getTriangles())
{
if (f.hasAttributes(AbstractHalfEdge.OUTER))
continue;
ot = f.getAbstractHalfEdge(ot);
double val = triangleQuality(ot);
qualityMap.put(f, val);
}
}
public final void compute()
{
LOGGER.info("Run "+getClass().getName());
if (nloop > 0)
{
// First compute triangle quality
qualityMap = new TObjectDoubleHashMap<Triangle>(mesh.getTriangles().size());
computeTriangleQuality();
nodeset = mesh.getNodes();
nodeProjection = new HashMap<Vertex, LocalSurfaceProjection>(mesh.getTriangles().size() / 2);
if (nodeset == null)
{
nodeset = new LinkedHashSet<Vertex>(mesh.getTriangles().size() / 2);
for (Triangle f: mesh.getTriangles())
{
if (f.hasAttributes(AbstractHalfEdge.OUTER))
continue;
f.addVertexTo(nodeset);
}
}
for (Vertex v: nodeset)
{
if (!v.isManifold() || !v.isMutable() || v.getRef() > 0)
continue;
LocalSurfaceProjection qP = new QuadricProjection(v, true);
if (!qP.canProject())
continue;
nodeProjection.put(v, qP);
}
for (int i = 0; i < nloop; i++)
{
processAllNodes();
postProcessIteration(mesh, i);
}
}
LOGGER.info("Number of moved points: "+processed);
LOGGER.info("Total number of points not moved during processing: "+notProcessed);
}
/*
* Moves all nodes using a modified Laplacian smoothing.
*/
private void processAllNodes()
{
AbstractHalfEdge ot = null;
// Compute vertex quality
tree.clear();
for (Vertex v: nodeset)
{
if (!v.isManifold() || !v.isMutable() || v.getRef() > 0)
continue;
Triangle f = (Triangle) v.getLink();
ot = f.getAbstractHalfEdge(ot);
if (ot.destination() == v)
ot = ot.next();
else if (ot.apex() == v)
ot = ot.prev();
assert ot.origin() == v;
double qv = vertexQuality(ot);
if (qv <= tolerance)
tree.insert(v, qv);
}
// Now smooth nodes iteratively
while (!tree.isEmpty())
{
Iterator<QSortedTree.Node<Vertex>> itt = tree.iterator();
QSortedTree.Node<Vertex> q = itt.next();
if (q.getValue() > tolerance)
break;
Vertex v = q.getData();
tree.remove(v);
if (v.getRef() != 0 && preserveBoundaries)
{
notProcessed++;
continue;
}
if (smoothNode(v, ot, q.getValue()))
{
processed++;
if (processed > 0 && (processed % progressBarStatus) == 0)
LOGGER.info("Vertices processed: "+processed);
if (!refresh)
continue;
assert ot != null;
// Update triangle quality
Vertex d = ot.destination();
do
{
ot = ot.nextOriginLoop();
if (ot.hasAttributes(AbstractHalfEdge.OUTER))
continue;
double qt = triangleQuality(ot);
qualityMap.put(ot.getTri(), qt);
}
while (ot.destination() != d);
// Update neighbor vertex quality
do
{
ot = ot.nextOriginLoop();
Vertex n = ot.destination();
if (n == mesh.outerVertex || !tree.contains(n))
continue;
if (ot.hasAttributes(AbstractHalfEdge.OUTER))
continue;
ot = ot.next();
double qv = vertexQuality(ot);
ot = ot.prev();
if (qv <= tolerance)
tree.update(n, qv);
else
{
tree.remove(n);
notProcessed++;
}
}
while (ot.destination() != d);
}
else
notProcessed++;
}
}
final void postProcessIteration(Mesh mesh, int i)
{
// Can be overridden
}
private boolean smoothNode(Vertex n, AbstractHalfEdge ot, double quality)
{
Triangle f = (Triangle) n.getLink();
ot = f.getAbstractHalfEdge(ot);
if (ot.destination() == n)
ot = ot.next();
else if (ot.apex() == n)
ot = ot.prev();
assert ot.origin() == n;
// Compute 3D coordinates centroid
int nn = 0;
c.moveTo(0, 0, 0);
assert n.isManifold();
Vertex d = ot.destination();
do
{
ot = ot.nextOriginLoop();
Vertex v = ot.destination();
if (v != mesh.outerVertex)
{
nn++;
double l = n.distance3D(v);
if (sizeTarget > 0.0)
{
if (l > 1.0)
{
// Find the point on this edge which has the
// desired length
l = sizeTarget / l;
c.moveTo(
c.getX() + v.getX() + l * (n.getX() - v.getX()),
c.getY() + v.getY() + l * (n.getY() - v.getY()),
c.getZ() + v.getZ() + l * (n.getZ() - v.getZ()));
}
else
c.add(v);
}
else
c.add(v);
}
}
while (ot.destination() != d);
assert (nn > 0);
c.scale(1.0/nn);
c.moveTo(
n.getX() + relaxation * (c.getX() - n.getX()),
n.getY() + relaxation * (c.getY() - n.getY()),
n.getZ() + relaxation * (c.getZ() - n.getZ()));
if (!mesh.checkNewRingNormals(ot, c))
{
LOGGER.finer("Point not moved, some triangles would become inverted");
return false;
}
LocalSurfaceProjection tr = nodeProjection.get(n);
if (tr == null || !tr.canProject())
{
LOGGER.finer("Point cannot be projected into surface");
return false;
}
tr.project(c);
double saveX = c.getX();
double saveY = c.getY();
double saveZ = c.getZ();
n.moveTo(c);
if (checkQuality)
{
// Check that quality has not been degraded
if (vertexQuality(ot) < quality)
{
n.moveTo(saveX, saveY, saveZ);
LOGGER.finer("Point not moved, quality decreases");
return false;
}
}
return true;
}
private double triangleQuality(AbstractHalfEdge edge)
{
Triangle f = edge.getTri();
assert f.getV0() != mesh.outerVertex && f.getV1() != mesh.outerVertex && f.getV2() != mesh.outerVertex : f;
double p = f.getV0().distance3D(f.getV1()) + f.getV1().distance3D(f.getV2()) + f.getV2().distance3D(f.getV0());
double area = edge.area(mesh);
double ret = scaleFactor * area / p / p;
assert ret >= 0.0 && ret <= 1.01;
return ret;
}
private double vertexQuality(AbstractHalfEdge edge)
{
Vertex d = edge.destination();
double ret = Double.MAX_VALUE;
do
{
edge = edge.nextOriginLoop();
if (edge.hasAttributes(AbstractHalfEdge.OUTER))
continue;
double qt = triangleQuality(edge);
if (qt < ret)
ret = qt;
}
while (edge.destination() != d);
return ret;
}
private static void usage(int rc)
{
System.out.println("Usage: SmoothNodes3D [options] xmlDir outDir");
System.out.println("Options:");
System.out.println(" -h, --help Display this message and exit");
System.out.println(" --iterations <n> Iterate <n> times over all nodes");
System.out.println(" --size <s> Set target size");
System.out.println(" --tolerance <t> Consider only nodes with quality lower than <t>");
System.out.println(" --relaxation <r> Set relaxation factor");
System.out.println(" --refresh Update vertex quality before each iteration");
System.exit(rc);
}
/**
*
* @param args [options] xmlDir outDir
*/
public static void main(String[] args)
{
org.jcae.mesh.amibe.traits.MeshTraitsBuilder mtb = org.jcae.mesh.amibe.traits.MeshTraitsBuilder.getDefault3D();
mtb.addNodeList();
Mesh mesh = new Mesh(mtb);
Map<String, String> opts = new HashMap<String, String>();
int argc = 0;
for (String arg: args)
if (arg.equals("--help") || arg.equals("-h"))
usage(0);
while (argc < args.length-1)
{
if (args[argc].length() < 2 || args[argc].charAt(0) != '-' || args[argc].charAt(1) != '-')
break;
if (args[argc].equals("--refresh") || args[argc].equals("--boundaries"))
{
opts.put(args[argc].substring(2), "true");
argc++;
}
else
{
opts.put(args[argc].substring(2), args[argc+1]);
argc += 2;
}
}
if (argc + 2 != args.length)
usage(1);
try
{
MeshReader.readObject3D(mesh, args[argc]);
}
catch (IOException ex)
{
ex.printStackTrace();
throw new RuntimeException(ex);
}
new SmoothNodes3D(mesh, opts).compute();
try
{
MeshWriter.writeObject3D(mesh, args[argc+1], null);
}
catch (IOException ex)
{
ex.printStackTrace();
throw new RuntimeException(ex);
}
}
}