Commit f88ed04d authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

triangulation updates

parent 7e38204f
......@@ -110,7 +110,8 @@
"treadCount" : 15,
"upwardDirection" : {
"x" : 1.0,
"y" : 0.0
"y" : 0.0,
"identifier" : -1
}
} ],
"targets" : [ {
......
......@@ -146,7 +146,8 @@
"treadCount" : 13,
"upwardDirection" : {
"x" : 0.0,
"y" : 1.0
"y" : 1.0,
"identifier" : -1
}
}, {
"shape" : {
......@@ -160,7 +161,8 @@
"treadCount" : 13,
"upwardDirection" : {
"x" : 0.0,
"y" : -1.0
"y" : -1.0,
"identifier" : -1
}
} ],
"targets" : [ {
......
......@@ -288,6 +288,13 @@ public interface IPolyConnectivity<P extends IPoint, V extends IVertex<P>, E ext
return GeometryUtils.isRightOf(p1, p2, x1, y1);
}
/**
* Tests if the line-segment edge intersects the line defined by p1 and p2.
* @param p1
* @param p2
* @param edge
* @return
*/
default boolean intersects(final IPoint p1, final IPoint p2, E edge) {
VPoint q1 = getMesh().toPoint(getMesh().getVertex(getMesh().getPrev(edge)));
VPoint q2 = getMesh().toPoint(getMesh().getVertex(edge));
......
......@@ -751,30 +751,17 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
* @param stopCondition
* @return
*/
default F straightWalk2D(final E startVertex, final VPoint direction, final Predicate<E> stopCondition) {
default F straightWalk2D(final E startVertex, final F face, final VPoint direction, final Predicate<E> stopCondition) {
E startEdge = getMesh().getPrev(startVertex);
E edge = getMesh().getPrev(startVertex);
VPoint p1 = getMesh().toPoint(startVertex);
VPoint p2 = p1.add(direction);
VPoint p3 = getMesh().toPoint(startEdge);
VPoint p4 = getMesh().toPoint(getMesh().getPrev(startEdge));
VPoint q = p1;
VPoint p = GeometryUtils.intersectionPoint(p1.getX(), p1.getY(), p2.getX(), p2.getY(), p3.getX(), p3.getY(), p4.getX(), p4.getY());
F face;
E edge;
if(isRightOf(p.getX(), p.getY(), startEdge)) {
edge = startEdge;
face = getMesh().getFace(edge);
}
else {
edge = getMesh().getTwin(startEdge);
face = getMesh().getFace(edge);
}
return straightWalk2D(q, p, face, edge, stopCondition);
assert intersects(p1, p1.add(direction), edge);
assert getMesh().getEdges(face).contains(startVertex);
// TODO: quick solution!
VPoint q = p1.add(direction.scalarMultiply(Double.MAX_VALUE * 0.5));
return straightWalk2D(p1, q, face, edge, e -> (stopCondition.test(e) || !isRightOf(q.x, q.y, e)));
}
default F straightWalk2D(final double x1, final double y1, final F startFace, final Predicate<E> stopCondition) {
......@@ -789,6 +776,19 @@ public interface ITriConnectivity<P extends IPoint, V extends IVertex<P>, E exte
return straightWalk2D(q, p, face, edge, stopCondition);
}
/**
* Walks along the line defined by q and p. The walking direction should be controlled by the stopCondition e.g.
* e -> !isRightOf(x1, y1, e) stops the walk if (x1, y1) is on the left side of each edge which is the case if the
* point is inside the reached face. The walk starts at the startFace and continues in the direction of line defined
* by q and p using the any edge which does not fulfill the stopCondition.
*
* @param q
* @param p
* @param startFace
* @param startEdge
* @param stopCondition
* @return
*/
default F straightWalk2D(final VPoint q, final VPoint p, final F startFace, final E startEdge, final Predicate<E> stopCondition) {
Optional<E> optEdge;
F face = startFace;
......
......@@ -33,6 +33,10 @@ public class CellGrid {
/** Number of points along y axis. */
protected final int numPointsY;
protected final double xMin;
protected final double yMin;
protected CellState[][] values;
/**
......@@ -40,10 +44,12 @@ public class CellGrid {
* point values are initialized with 'value'.
*/
public CellGrid(double width, double height, double resolution,
CellState value) {
CellState value, double xMin, double yMin) {
this.width = width;
this.height = height;
this.resolution = resolution;
this.xMin = xMin;
this.yMin = yMin;
/* 0.001 avoids that numPointsX/Y are too small due to numerical errors. */
numPointsX = (int) Math.floor(width / resolution + 0.001) + 1;
......@@ -54,6 +60,14 @@ public class CellGrid {
reset(value);
}
/**
* Creates an grid with the given width, height and resolution. All grid
* point values are initialized with 'value'.
*/
public CellGrid(double width, double height, double resolution, CellState value) {
this(width, height, resolution, value, 0, 0);
}
/**
* Creates a deep copy of the given grid.
*/
......@@ -64,6 +78,8 @@ public class CellGrid {
numPointsX = grid.numPointsX;
numPointsY = grid.numPointsY;
values = new CellState[numPointsX][numPointsY];
xMin = grid.xMin;
yMin = grid.yMin;
for (int row = 0; row < numPointsY; row++) {
for (int col = 0; col < numPointsX; col++) {
......@@ -138,7 +154,7 @@ public class CellGrid {
* Converts the matrix indices to coordinates.
*/
public VPoint pointToCoord(int pointX, int pointY) {
return new VPoint(pointX * resolution, pointY * resolution);
return new VPoint(xMin + pointX * resolution, yMin + pointY * resolution);
}
/**
......@@ -183,20 +199,20 @@ public class CellGrid {
* Returns the closest grid point (matrix index) to the given coordinates.
*/
public Point getNearestPoint(double x, double y) {
if (x < 0) {
x = 0;
if (x < xMin) {
x = xMin;
}
if (y < 0) {
y = 0;
if (y < yMin) {
y = yMin;
}
if (y > getHeight()) {
y = getHeight();
if (y > getHeight() + yMin) {
y = yMin+getHeight();
}
if (x > getWidth()) {
x = getWidth();
if (x > getWidth() + xMin) {
x = getWidth() + xMin;
}
return new Point((int) (x / resolution + 0.5),
(int) (y / resolution + 0.5));
return new Point((int) ((x - xMin) / resolution + 0.5),
(int) ((y - yMin) / resolution + 0.5));
}
/**
......@@ -204,19 +220,19 @@ public class CellGrid {
* towards origin.
*/
public Point getNearestPointTowardsOrigin(double x, double y) {
if (x < 0) {
x = 0;
if (x < xMin) {
x = xMin;
}
if (y < 0) {
y = 0;
if (y < yMin) {
y = yMin;
}
if (y > getHeight()) {
y = getHeight();
if (y > getHeight() + yMin) {
y = yMin+getHeight();
}
if (x > getWidth()) {
x = getWidth();
if (x > getWidth() + xMin) {
x = getWidth() + xMin;
}
return new Point((int) (x / resolution), (int) (y / resolution));
return new Point((int) ((x - xMin) / resolution), (int) ((y - yMin) / resolution));
}
/**
......@@ -288,6 +304,10 @@ public class CellGrid {
.flatMap(stream -> stream);
}
public double getMinX() { return xMin; }
public double getMinY() { return yMin; }
public boolean isValidPoint(Point point) {
Point p = (point);
......
......@@ -27,11 +27,11 @@ public abstract class AbstractGridEikonalSolver implements EikonalSolver {
int incX = 1, incY = 1;
double gridPotentials[];
if (x >= potentialField.getWidth()) {
if (x >= potentialField.getWidth() + potentialField.getMinY()) {
incX = 0;
}
if (y >= potentialField.getHeight()) {
if (y >= potentialField.getHeight() + potentialField.getMinY()) {
incY = 0;
}
......
......@@ -14,6 +14,7 @@ import org.vadere.util.potential.CellGrid;
import org.vadere.util.potential.CellState;
import org.vadere.util.potential.PathFindingTag;
import org.vadere.util.potential.timecost.ITimeCostFunction;
import org.vadere.util.triangulation.adaptive.IDistanceFunction;
/**
* EikonalSolverFMM initializes a potential field on basis
......@@ -29,6 +30,8 @@ public class EikonalSolverFMM implements EikonalSolver {
protected CellGrid cellGrid;
protected List<Point> targetPoints;
protected List<VShape> targetShapes;
protected IDistanceFunction distFunc;
boolean isHighAccuracy = false;
/** only for logging */
......@@ -68,13 +71,41 @@ public class EikonalSolverFMM implements EikonalSolver {
}
}
/**
* Initializes the FM potential calculator with a time cost function F > 0.
*/
public EikonalSolverFMM(
final CellGrid potentialField,
final IDistanceFunction distFunc,
final boolean isHighAccuracy,
final ITimeCostFunction timeCostFunction,
final double unknownPenalty) {
super(potentialField, unknownPenalty);
this.cellGrid = potentialField;
this.targetPoints = cellGrid.pointStream().filter(p -> cellGrid.getValue(p).tag == PathFindingTag.Target)
.collect(Collectors.toList());
this.distFunc = distFunc;
this.isHighAccuracy = isHighAccuracy;
ComparatorPotentialFieldValue comparator = new ComparatorPotentialFieldValue(
potentialField);
this.narrowBand = new PriorityQueue<>(50, comparator);
this.timeCostFunction = timeCostFunction;
if (targetPoints.size() == 0) {
logger.error("PotentialFieldInitializerFastMarching::Run(): "
+ "Warning, no target points given. Target missing or grid resolution too low.");
return;
}
}
@Override
public void initialize() {
for (Point point : targetPoints) {
if(!targetShapes.isEmpty()) {
//if(!targetShapes.isEmpty()) {
setTargetNeighborsDistances(point);
}
//}
narrowBand.add(point);
}
......@@ -195,7 +226,7 @@ public class EikonalSolverFMM implements EikonalSolver {
cellGrid.setValue(
neighbor,
new CellState(minDistanceToTarget(cellGrid.pointToCoord(neighbor)),
new CellState(minDistanceToTarget(cellGrid.pointToCoord(neighbor)) / timeCostFunction.costAt(cellGrid.pointToCoord(neighbor)),
PathFindingTag.Reachable));
narrowBand.add(neighbor);
}
......@@ -209,6 +240,8 @@ public class EikonalSolverFMM implements EikonalSolver {
// computed starting there.
VPoint dp = new VPoint(cellGrid.getWidth() / (cellGrid.getNumPointsX() - 1) / 2.0,
cellGrid.getHeight() / (cellGrid.getNumPointsY() - 1) / 2.0);
if(targetShapes != null && !targetShapes.isEmpty()) {
for (VShape targetShape : targetShapes) {
// negative distances are possible when point is inside the target
tmp = Math.max(0, targetShape.distance(point.add(dp)));
......@@ -216,6 +249,14 @@ public class EikonalSolverFMM implements EikonalSolver {
minDistance = tmp;
}
}
}
else if(distFunc != null) {
tmp = Math.max(0, -distFunc.apply(point));
if (tmp < minDistance) {
minDistance = tmp;
}
}
return minDistance;
}
}
......@@ -129,7 +129,10 @@ public class PSMeshing implements IPSMeshing {
quality = getQuality();
log.info("quality = " + quality);
}
finalize();
}
public void finalize() {
computeScalingFactor();
computeForces();
computeDelta();
......@@ -176,10 +179,10 @@ public class PSMeshing implements IPSMeshing {
//log.info(scalingFactor);
//long ms = System.currentTimeMillis();
computeForces();
//ms = System.currentTimeMillis() - ms;
//log.info("ms: " + ms);
computeScalingFactor();
computeForces();
//computeDelta();
updateVertices();
......@@ -189,9 +192,10 @@ public class PSMeshing implements IPSMeshing {
numberOfIllegalMovementTests++;
}
if(illegalMovement) {
//retriangulate();
while (flipEdges());
if(illegalMovement) {
retriangulate();
//while (flipEdges());
numberOfRetriangulations++;
}
......@@ -199,15 +203,15 @@ public class PSMeshing implements IPSMeshing {
flipEdges();
}
if(minDeltaTravelDistance < 0) {
if(minDeltaTravelDistance <= 0) {
computeMaxLegalMovements();
}
numberOfIterations++;
/*log.info("#illegalMovementTests: " + numberOfIllegalMovementTests);
log.info("#illegalMovementTests: " + numberOfIllegalMovementTests);
log.info("#retriangulations: " + numberOfRetriangulations);
log.info("#steps: " + numberOfIterations);
log.info("#points: " + getMesh().getVertices().size());*/
log.info("#points: " + getMesh().getVertices().size());
}
public boolean isMovementIllegal() {
......
......@@ -11,7 +11,7 @@ public class Parameters {
public final static double h0 = 0.15;
public final static boolean uniform = false;
public final static String method = "Distmesh"; // "Distmesh" or "Density"
final static double qualityMeasurement = 0.928;
final static double qualityMeasurement = 0.95;
final static double MINIMUM = 0.25;
final static double DENSITYWEIGHT = 2;
final static int NPOINTS = 100000;
......
......@@ -2,6 +2,7 @@ package org.vadere.util.triangulation.adaptive;
import org.apache.log4j.Logger;
import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.tex.TexGraphGenerator;
import java.util.ArrayList;
......@@ -21,7 +22,10 @@ public class TestEnhancedVersion3 extends JFrame {
IDistanceFunction distanceFunc = p -> Math.abs(6 - Math.sqrt(p.getX() * p.getX() + p.getY() * p.getY())) - 4;
//IDistanceFunction distanceFunc4 = p -> Math.max(Math.abs(p.getY()) - 4, Math.abs(p.getX()) - 25);
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + p.distanceToOrigin();
IEdgeLengthFunction edgeLengthFunc = p -> 0.5;
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0;
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.min(Math.abs(distanceFunc.apply(p) + 4), Math.abs(distanceFunc.apply(p)));
IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.abs(distanceFunc.apply(p)*0.5);
//IDistanceFunction distanceFunc = p -> Math.max(Math.max(Math.max(distanceFunc1.apply(p), distanceFunc2.apply(p)), distanceFunc3.apply(p)), distanceFunc4.apply(p));
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.abs(distanceFunc.apply(p))/2;
......@@ -29,7 +33,7 @@ public class TestEnhancedVersion3 extends JFrame {
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.min(Math.abs(distanceFunc.apply(p) + 4), Math.abs(distanceFunc.apply(p)));
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0;
VRectangle bbox = new VRectangle(-11, -11, 22, 22);
PSMeshing meshGenerator = new PSMeshing(distanceFunc, edgeLengthFunc, 6.0, bbox, new ArrayList<>());
PSMeshing meshGenerator = new PSMeshing(distanceFunc, edgeLengthFunc, 0.6, bbox, new ArrayList<>());
meshGenerator.initialize();
PSMeshingPanel distmeshPanel = new PSMeshingPanel(meshGenerator, 1000, 800);
......@@ -44,7 +48,7 @@ public class TestEnhancedVersion3 extends JFrame {
int counter = 0;
long time = 0;
while (counter <= 560) {
while (counter <= 1000) {
//obscuteTriangles = meshGenerator.getTriangles().stream().filter(tri -> tri.isNonAcute()).count();
//PriorityQueue<PFace<MeshPoint>> priorityQueue = meshGenerator.getQuailties();
//avgQuality = priorityQueue.stream().reduce(0.0, (aDouble, meshPointPFace) -> aDouble + meshGenerator.faceToQuality(meshPointPFace), (d1, d2) -> d1 + d2) / priorityQueue.size();
......@@ -64,7 +68,9 @@ public class TestEnhancedVersion3 extends JFrame {
System.out.println("Quality: " + meshGenerator.getQuality());
System.out.println("Step-Time: " + ms);
}
meshGenerator.finalize();
System.out.print("overall time: " + time);
System.out.print(TexGraphGenerator.meshToGraph(meshGenerator.getMesh()));
//System.out.print("finished:" + meshGenerator.getMesh().getVertices().stream().filter(v -> !meshGenerator.getMesh().isDestroyed(v)).count());
//System.out.print("finished:" + avgQuality);
......
......@@ -2,12 +2,14 @@ package org.vadere.util.potential;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.junit.Before;
import org.junit.Test;
import org.vadere.util.geometry.mesh.gen.PFace;
import org.vadere.util.geometry.mesh.gen.PHalfEdge;
import org.vadere.util.geometry.mesh.gen.PVertex;
import org.vadere.util.geometry.mesh.inter.ITriangulation;
import org.vadere.util.geometry.mesh.inter.IVertex;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VCircle;
import org.vadere.util.geometry.shapes.VPoint;
......@@ -25,12 +27,15 @@ import org.vadere.util.triangulation.adaptive.PSDistmesh;
import org.vadere.util.triangulation.adaptive.PSDistmeshPanel;
import org.vadere.util.triangulation.adaptive.PSMeshing;
import java.awt.*;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Comparator;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.Stream;
import javax.swing.*;
......@@ -63,8 +68,9 @@ public class TestFFMNonUniformTriangulation {
@Test
public void testTriangulationFMM() {
IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.min(Math.abs(distanceFunc.apply(p) + 4), Math.abs(distanceFunc.apply(p)));
IEdgeLengthFunction unifromEdgeLengthFunc = p -> 1.0;
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + 0.5*Math.min(Math.abs(distanceFunc.apply(p) + 4), Math.abs(distanceFunc.apply(p)));
//IEdgeLengthFunction unifromEdgeLengthFunc = p -> 1.0;
IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.abs(distanceFunc.apply(p)) * 0.5;
List<VRectangle> targetAreas = new ArrayList<>();
List<IPoint> targetPoints = new ArrayList<>();
PSMeshing meshGenerator = new PSMeshing(distanceFunc, edgeLengthFunc, 0.6, bbox, new ArrayList<>());
......@@ -77,16 +83,22 @@ public class TestFFMNonUniformTriangulation {
VRectangle rect = new VRectangle(width / 2, height / 2, 100, 100);
targetAreas.add(rect);
EikonalSolver solver = new EikonalSolverFMMTriangulation(new UnitTimeCostFunction(), triangulation,
triangulation.getMesh().getBoundaryVertices().stream().collect(Collectors.toList()));
EikonalSolver solver = new EikonalSolverFMMTriangulation(
new UnitTimeCostFunction(),
triangulation,
triangulation.getMesh().getBoundaryVertices().stream().collect(Collectors.toList()),
distanceFunc);
log.info("start FFM");
solver.initialize();
log.info("FFM finished");
double maxError = 0;
double sum = 0;
int counter = 0;
try {
//System.out.println(getClass().getClassLoader().getResource("./potentialField.csv").getFile());
FileWriter writer = new FileWriter("./potentialField.csv");
FileWriter writer = new FileWriter("./potentialField_adapt_0_7.csv");
for(double y = bbox.getMinY()+2; y <= bbox.getMaxY()-2; y += 0.1) {
for(double x = bbox.getMinX()+2; x < bbox.getMaxX()-2; x += 0.1) {
double val = solver.getValue(x ,y);
......@@ -94,6 +106,87 @@ public class TestFFMNonUniformTriangulation {
double side = Math.min((new VPoint(x, y).distanceToOrigin()-2.0), (10 - new VPoint(x, y).distanceToOrigin()));
side = Math.max(side, 0.0);
maxError = Math.max(maxError, Math.abs(val - side));
sum += Math.abs(val - side) * Math.abs(val - side);
counter++;
}
writer.write(""+solver.getValue(x ,y) + " ");
}
writer.write("\n");
}
writer.flush();
} catch (FileNotFoundException e) {
e.printStackTrace();
} catch (IOException e) {
e.printStackTrace();
}
log.info(triangulation.getMesh().getVertices().size());
log.info("max edge length: " + triangulation.getMesh().streamEdges().map(e -> triangulation.getMesh().toLine(e).length()).max(Comparator.comparingDouble(d -> d)));
log.info("min edge length: " +triangulation.getMesh().streamEdges().map(e -> triangulation.getMesh().toLine(e).length()).min(Comparator.comparingDouble(d -> d)));
log.info("max distance to boundary: " + triangulation.getMesh().getBoundaryVertices().stream().map(p -> Math.abs(distanceFunc.apply(p))).max(Comparator.comparingDouble(d -> d)));
//log.info("L2-Error: " + computeL2Error(triangulation, distanceFunc));
log.info("max error: " + maxError);
log.info("max error-2: " + triangulation.getMesh().getVertices().stream().map(p -> Math.abs(Math.abs(p.getPoint().getPotential() + distanceFunc.apply(p)))).max(Comparator.comparingDouble(d -> d)));
log.info("L2-error: " + Math.sqrt(sum / counter));
log.info("L2-error-2: " + Math.sqrt(triangulation.getMesh().getVertices().stream()
.map(p -> Math.abs(Math.abs(p.getPoint().getPotential() + distanceFunc.apply(p))))
.map(val -> val * val)
.reduce(0.0, (d1, d2) -> d1 + d2) / triangulation.getMesh().getNumberOfVertices()));
//assertTrue(0.0 == solver.getValue(5, 5));
//assertTrue(0.0 < solver.getValue(1, 7));
}
@Test
public void testTriangulationFMMCase2() {
//IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.min(Math.abs(distanceFunc.apply(p) + 4) * 2, Math.abs(distanceFunc.apply(p)) * 2);
//IEdgeLengthFunction unifromEdgeLengthFunc = p -> 1.0;
IEdgeLengthFunction edgeLengthFunc = p -> 1.0 + Math.abs(distanceFunc.apply(p)*0.5);
List<VRectangle> targetAreas = new ArrayList<>();
PSMeshing meshGenerator = new PSMeshing(distanceFunc, edgeLengthFunc, 0.6, bbox, new ArrayList<>());
meshGenerator.execute();
triangulation = meshGenerator.getTriangulation();
//targetPoints.add(new MeshPoint(0, 0, false));
VRectangle rect = new VRectangle(width / 2, height / 2, 100, 100);
targetAreas.add(rect);
List<PVertex<MeshPoint>> targetPoints = triangulation.getMesh().getVertices().stream()
.filter(v -> triangulation.getMesh().isAtBoundary(v))
.filter(p-> (Math.abs(new VPoint(p.getX(), p.getY()).distanceToOrigin()-2.0)) < 2).collect(Collectors.toList());
log.info(targetPoints);
EikonalSolver solver = new EikonalSolverFMMTriangulation(
new