Commit cafd2cf8 authored by BZoennchen's avatar BZoennchen

save computation in fmm for triangulations.

parent ecb24671
......@@ -78,7 +78,7 @@
"finishTime" : 400.0,
"simTimeStepLength" : 0.4,
"realTimeSimTimeRatio" : 0.0,
"writeSimulationData" : true,
"writeSimulationData" : false,
"visualizationEnabled" : true,
"printFPS" : false,
"digitsPerCoordinate" : 2,
......
package org.vadere.simulator.models.potential.solver.calculators.mesh;
import org.apache.commons.lang3.tuple.Triple;
import org.apache.log4j.Level;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.vadere.meshing.utils.debug.DebugGui;
import org.vadere.meshing.utils.debug.SimpleTriCanvas;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.meshing.mesh.inter.IFace;
import org.vadere.meshing.mesh.inter.IHalfEdge;
......@@ -13,7 +12,6 @@ import org.vadere.meshing.mesh.inter.IMesh;
import org.vadere.meshing.mesh.inter.IIncrementalTriangulation;
import org.vadere.meshing.mesh.inter.IVertex;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VShape;
import org.vadere.util.math.InterpolationUtil;
......@@ -24,13 +22,16 @@ import org.vadere.util.data.cellgrid.IPotentialPoint;
import org.vadere.simulator.models.potential.solver.timecost.ITimeCostFunction;
import org.vadere.util.math.IDistanceFunction;
import java.awt.*;
import java.util.Collection;
import java.util.Comparator;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.Optional;
import java.util.PriorityQueue;
import java.util.Set;
import java.util.function.Function;
import java.util.function.Predicate;
......@@ -48,6 +49,10 @@ import java.util.function.Predicate;
public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends IVertex<P>, E extends IHalfEdge<P>, F extends IFace<P>> implements EikonalSolver {
private static Logger logger = LogManager.getLogger(EikonalSolverFMMTriangulation.class);
private Set<F> nonAccuteTris = new HashSet<>();
private Map<Triple<P, P, P>, Double> angles = new HashMap();
private Map<Triple<P, P, P>, Double> sinPhis = new HashMap();
private Map<Triple<P, P, P>, Double> cosPhis = new HashMap();
static {
logger.setLevel(Level.INFO);
......@@ -58,6 +63,10 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
*/
private ITimeCostFunction timeCostFunction;
private IDistanceFunction distFunc;
private Collection<V> targetVertices;
/**
* The triangulation the solver uses.
*/
......@@ -154,30 +163,12 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
this.calculationFinished = false;
this.timeCostFunction = timeCostFunction;
this.narrowBand = new PriorityQueue<>(pointComparator);
for(V vertex : targetVertices) {
P potentialPoint = getMesh().getPoint(vertex);
double distance = -distFunc.apply(potentialPoint);
if(potentialPoint.getPathFindingTag() != PathFindingTag.Undefined) {
narrowBand.remove(vertex);
}
potentialPoint.setPotential(Math.min(potentialPoint.getPotential(), distance / timeCostFunction.costAt(potentialPoint)));
potentialPoint.setPathFindingTag(PathFindingTag.Reached);
narrowBand.add(vertex);
for(V v : triangulation.getMesh().getAdjacentVertexIt(vertex)) {
P potentialP = getMesh().getPoint(v);
if(potentialP.getPathFindingTag() == PathFindingTag.Undefined) {
double dist = Math.max(-distFunc.apply(potentialP), 0);
logger.debug("T at " + potentialP + " = " + dist);
potentialP.setPotential(Math.min(potentialP.getPotential(), dist / timeCostFunction.costAt(potentialP)));
potentialP.setPathFindingTag(PathFindingTag.Reachable);
narrowBand.add(v);
}
}
this.distFunc = distFunc;
this.targetVertices = targetVertices;
for(F face : triangulation.getMesh().getFaces()) {
if(isNonAcute(face)) {
nonAccuteTris.add(face);
}
}
}
......@@ -218,6 +209,37 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
}
}
public void reset() {
triangulation.getMesh().streamPoints().forEach(p -> p.setPathFindingTag(PathFindingTag.Undefined));
triangulation.getMesh().streamPoints().forEach(p -> p.setPotential(Double.MAX_VALUE));
calculationFinished = false;
for(V vertex : targetVertices) {
P potentialPoint = getMesh().getPoint(vertex);
double distance = -distFunc.apply(potentialPoint);
if(potentialPoint.getPathFindingTag() != PathFindingTag.Undefined) {
narrowBand.remove(vertex);
}
potentialPoint.setPotential(Math.min(potentialPoint.getPotential(), distance / timeCostFunction.costAt(potentialPoint)));
potentialPoint.setPathFindingTag(PathFindingTag.Reached);
narrowBand.add(vertex);
for(V v : triangulation.getMesh().getAdjacentVertexIt(vertex)) {
P potentialP = getMesh().getPoint(v);
if(potentialP.getPathFindingTag() == PathFindingTag.Undefined) {
double dist = Math.max(-distFunc.apply(potentialP), 0);
logger.debug("T at " + potentialP + " = " + dist);
potentialP.setPotential(Math.min(potentialP.getPotential(), dist / timeCostFunction.costAt(potentialP)));
potentialP.setPathFindingTag(PathFindingTag.Reachable);
narrowBand.add(v);
}
}
}
}
// unknownPenalty is ignored.
@Override
public double getPotential(@NotNull final IPoint pos, final double unknownPenalty, final double weight) {
......@@ -358,6 +380,10 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
return angle1 > rightAngle + GeometryUtils.DOUBLE_EPS;
}
private boolean isNonAcute(@NotNull final F face) {
return isNonAcute(triangulation.getMesh().getEdge(face));
}
/**
* Updates a point given a triangle. The point can only be updated if the
* triangle triangleContains it and the other two points are in the frozen band.
......@@ -378,10 +404,10 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
if(isFeasibleForComputation(p1) && isFeasibleForComputation(p2)) {
if(!isNonAcute(halfEdge)) {
if(!nonAccuteTris.contains(face)) {
double potential = computePotential(point, p1, p2);
//logger.info("compute potential " + potential);
return computePotential(point, p1, p2);
return computePotential(point, p1, p2);
} // we only try to find a virtual vertex if both points are already frozen
else {
logger.debug("special case for non-acute triangle");
......@@ -520,7 +546,6 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
* @return the traveling time T at <tt>point</tt> by using the triangle (point, point1, point2) for the computation
*/
private double computePotential(final P point, final P point1, final P point2) {
// see: Sethian, Level Set Methods and Fast Marching Methods, page 124.
P p1; // A
P p2; // B
......@@ -543,15 +568,17 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
double b = p1.distance(point);
double c = p1.distance(p2);
double phi = GeometryUtils.angle(p1, point, p2);
double cosphi = Math.cos(phi);
//double phi = angle(p1, point, p2);
double cosphi = cosPhi(p1, point, p2);
double sinPhi = sinPhi(p1, point, p2);
double F = 1.0 / timeCostFunction.costAt(point);
// solve x2 t^2 + x1 t + x0 == 0
double x2 = a * a + b * b - 2 * a * b * cosphi;
double x1 = 2 * b * u * (a * cosphi - b);
double x0 = b * b * (u * u - F * F * a * a * Math.sin(phi) * Math.sin(phi));
double x0 = b * b * (u * u - F * F * a * a * sinPhi * sinPhi);
double t = solveQuadratic(x2, x1, x0);
double inTriangle = (b * (t - u) / t);
......@@ -562,6 +589,42 @@ public class EikonalSolverFMMTriangulation<P extends IPotentialPoint, V extends
}
}
private double sinPhi(P p1, P p, P p2) {
Double sinPhi = sinPhis.get(Triple.of(p, p1, p2));
if(sinPhi == null) {
sinPhi = Math.sin(angle(p1, p, p2));
sinPhis.put(Triple.of(p, p1, p2), sinPhi);
return sinPhi;
}
else {
return sinPhi;
}
}
private double cosPhi(P p1, P p, P p2) {
Double cosPhi = cosPhis.get(Triple.of(p, p1, p2));
if(cosPhi == null) {
cosPhi = Math.cos(angle(p1, p, p2));
cosPhis.put(Triple.of(p, p1, p2), cosPhi);
return cosPhi;
}
else {
return cosPhi;
}
}
private double angle(P p1, P p, P p2) {
Double angle = angles.get(Triple.of(p, p1, p2));
if(angle == null) {
angle = GeometryUtils.angle(p1, p, p2);
angles.put(Triple.of(p, p1, p2), angle);
return angle;
}
else {
return angle;
}
}
/**
* Solves the quadratic equation given by (a*x^2+b*x+c=0).
*
......
......@@ -49,7 +49,7 @@ public class PerformanceTriangleFMM {
private static Logger log = LogManager.getLogger(TestFFMNonUniformTriangulation.class);
private static final VRectangle bbox = new VRectangle(-12, -12, 24, 24);
private static final IDistanceFunction distanceFunc = p -> Math.abs(6 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 4;
private static final double initialEdgeLen = 0.1;
private static final double initialEdgeLen = 0.6;
static {
//LogManager.shutdown();
......@@ -61,31 +61,13 @@ public class PerformanceTriangleFMM {
return meshGenerator.generate();
}
private static void solve(IIncrementalTriangulation<PotentialPoint, PVertex<PotentialPoint>, PHalfEdge<PotentialPoint>, PFace<PotentialPoint>> triangulation) {
triangulation.getMesh().streamPoints().forEach(p -> p.setPathFindingTag(PathFindingTag.Undefined));
/**
* (2) define target points
*/
List<PVertex<PotentialPoint>> targetVertices = triangulation.getMesh().getBoundaryVertices().stream().collect(Collectors.toList());
/**
* (3) solve the eikonal equation on the mesh
*/
EikonalSolver solver = new EikonalSolverFMMTriangulation(
new UnitTimeCostFunction(),
triangulation,
targetVertices,
distanceFunc);
private static void solve(EikonalSolverFMMTriangulation<PotentialPoint, PVertex<PotentialPoint>, PHalfEdge<PotentialPoint>, PFace<PotentialPoint>> solver) {
solver.reset();
long ms = System.currentTimeMillis();
System.out.println("start FFM");
solver.initialize();
System.out.println("FFM finished");
System.out.println("time: " + (System.currentTimeMillis() - ms));
triangulation.getMesh().garbageCollection();
System.out.println("nPoints: " + (triangulation.getMesh().getNumberOfVertices()));
}
......@@ -95,11 +77,23 @@ public class PerformanceTriangleFMM {
*/
IIncrementalTriangulation<PotentialPoint, PVertex<PotentialPoint>, PHalfEdge<PotentialPoint>, PFace<PotentialPoint>> triangulation = createTriangulation();
/**
* (2) define target points
*/
/**
* (2) define target points
*/
List<PVertex<PotentialPoint>> targetVertices = triangulation.getMesh().getBoundaryVertices().stream().collect(Collectors.toList());
/**
* (3) solve the eikonal equation on the mesh
*/
EikonalSolverFMMTriangulation<PotentialPoint, PVertex<PotentialPoint>, PHalfEdge<PotentialPoint>, PFace<PotentialPoint>> solver = new EikonalSolverFMMTriangulation(
new UnitTimeCostFunction(),
triangulation,
targetVertices,
distanceFunc);
/**
* (3) solve the eikonal equation on the mesh
*/
......@@ -110,7 +104,8 @@ public class PerformanceTriangleFMM {
} catch (InterruptedException e) {
e.printStackTrace();
}*/
solve(triangulation);
solve(solver);
System.out.println("nPoints: " + (triangulation.getMesh().getNumberOfVertices()));
}
}
}
......@@ -4,6 +4,7 @@ import jdk.nashorn.internal.runtime.regexp.joni.exception.InternalException;
public class MPoint implements IPoint, Cloneable{
private VPoint point;
private int hashCode = -1;
public MPoint(final double x, final double y){
this.point = new VPoint(x, y);
......@@ -30,36 +31,42 @@ public class MPoint implements IPoint, Cloneable{
@Override
public MPoint add(final IPoint point) {
this.point = this.point.add(point);
hashCode = -1;
return this;
}
@Override
public MPoint addPrecise(final IPoint point) {
this.point = this.point.addPrecise(point);
hashCode = -1;
return this;
}
@Override
public MPoint subtract(IPoint point) {
this.point = this.point.subtract(point);
hashCode = -1;
return this;
}
@Override
public MPoint multiply(IPoint point) {
this.point = this.point.multiply(point);
hashCode = -1;
return this;
}
@Override
public MPoint scalarMultiply(double factor) {
this.point = this.point.scalarMultiply(factor);
hashCode = -1;
return this;
}
@Override
public MPoint rotate(double radAngle) {
this.point = this.point.rotate(radAngle);
hashCode = -1;
return this;
}
......@@ -71,6 +78,7 @@ public class MPoint implements IPoint, Cloneable{
@Override
public MPoint norm() {
this.point = this.point.norm();
hashCode = -1;
return this;
}
......@@ -82,15 +90,19 @@ public class MPoint implements IPoint, Cloneable{
@Override
public MPoint normZeroSafe() {
this.point = this.point.normZeroSafe();
hashCode = -1;
return this;
}
@Override
public int hashCode() {
// hashCode of java.awt.geom.Point2D
long bits = java.lang.Double.doubleToLongBits(getX());
bits ^= java.lang.Double.doubleToLongBits(getY()) * 31;
return (((int) bits) ^ ((int) (bits >> 32)));
if(hashCode == -1) {
// hashCode of java.awt.geom.Point2D
long bits = java.lang.Double.doubleToLongBits(getX());
bits ^= java.lang.Double.doubleToLongBits(getY()) * 31;
hashCode = (((int) bits) ^ ((int) (bits >> 32)));
}
return hashCode;
}
@Override
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment