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

before introducing IVertex for the mesh

parent f6178c6a
......@@ -91,6 +91,26 @@ public class GeometryUtils {
return orderedList;
}
/**
* Computes area (it maybe a negative area) of the parallelogram defined by p, q, r.
* Note: + This area is zero if p, q, r lie on one line.
* + The area is < 0 if the points p, q, r are aligned clockwise (order matters).
* This is equivalent to: r lies on the right side of the line defined by p, q.
* + The area is > 0 if the points p, q, r are aligned counter-clockwise (order matters).
* This is equivalent to: r lies on the left side of the line defined by p, q.
*
* @param pX x-coordinate of p
* @param pY y-coordinate of p
* @param qX x-coordinate of q
* @param qY y-coordinate of q
* @param rX x-coordinate of r
* @param rY y-coordinate of r
* @return
*/
public static double ccw(final double qX, final double qY, final double pX, final double pY, final double rX, final double rY) {
return (qX - pX) * (rY - pY) - (rX - pX) * (qY - pY);
}
/**
* Calculate the counter clockwise result for the three given points.<br>
* ccw(p1,p2,p3) < 0 if p3 is left of Line(p1,p2)<br>
......@@ -105,8 +125,93 @@ public class GeometryUtils {
* third point
* @return ccw(p1 p2 p3)
*/
public static double ccw(VPoint p1, VPoint p2, VPoint p3) {
return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x);
public static double ccw(final IPoint p1, final IPoint p2, final IPoint p3) {
return (p2.getX() - p1.getX()) * (p3.getY() - p1.getY()) - (p2.getY() - p1.getY()) * (p3.getX() - p1.getX());
}
public static boolean isCCW(final double qX, final double qY, final double pX, final double pY, final double rX, final double rY) {
return ccw(qX, qY, pX, pY, rX, rY) > 0;
}
public static boolean isCCW(final IPoint p1, final IPoint p2, final IPoint p3) {
return ccw(p1, p2, p3) > 0;
}
public static boolean isCW(final double qX, final double qY, final double pX, final double pY, final double rX, final double rY) {
return ccw(qX, qY, pX, pY, rX, rY) < 0;
}
public static boolean isCW(final IPoint p1, final IPoint p2, final IPoint p3) {
return ccw(p1, p2, p3) < 0;
}
/**
* Tests if the line-segment (p1, p2) intersects the line defined by p, q.
* @param p point defining the line
* @param q point defining the line
* @param p1 point defining the line-segment
* @param p2 point defining the line-segment
* @return true if the line-segment intersects the line defined, otherwise false.
*/
public static boolean intersectLine(final IPoint p, final IPoint q, final IPoint p1, final IPoint p2) {
double ccw1 = ccw(p, q, p1);
double ccw2 = ccw(p, q, p2);
return (ccw1 < 0 && ccw2 > 0) || (ccw1 > 0 && ccw2 < 0);
}
/**
* Tests if the half-line-segment starting at p in the direction (q-p) intersects the line-segment (p1,p2).
* @param p the starting point of the half-line-segment
* @param q the point defining the direction (q-p) of the half-line-segment
* @param p1 point defining the line-segment
* @param p2 point defining the line-segment
* @return true if the line-segment intersects the half-line-segment defined, otherwise false.
*/
public static boolean intersectHalfLineSegment(final IPoint p, final IPoint q, final IPoint p1, final IPoint p2) {
double ccw1 = ccw(p, q, p1);
double ccw2 = ccw(p, q, p2);
if((ccw1 < 0 && ccw2 > 0)) {
return isCCW(p, p2, p1);
}
else if((ccw1 > 0 && ccw2 < 0)) {
return isCCW(p, p1, p2);
}
else {
return false;
}
}
/**
* Tests if the first line-segment (p,q) intersects the second line-segment (p1,p2).
* @param p point defining the first line-segment
* @param q point defining the first line-segment
* @param p1 point defining the second line-segment
* @param p2 point defining the second line-segment
* @return true if the first line-segment intersects the second line-segment, otherwise false.
*/
public static boolean intersectLineSegment(final IPoint p, final IPoint q, final IPoint p1, final IPoint p2) {
return intersectLine(p, q, p1, p2) && intersectLine(p1, p2, p, q);
}
/**
* Tests if the triangle (p1,p2,p3) contains the point r.
*
* @param p1 point of the triangle
* @param p2 point of the triangle
* @param p3 point of the triangle
* @param r point which the triangle might contain.
* @return if the triangle (p1,p2,p3) contains the point r, otherwise false.
*/
public static boolean contains(final IPoint p1, final IPoint p2, final IPoint p3, final IPoint r) {
boolean b1, b2, b3;
double d1 = GeometryUtils.ccw(r, p1, p2);
double d2 = GeometryUtils.ccw(r, p2, p3);
double d3 = GeometryUtils.ccw(r, p3, p1);
b1 = d1 < 0.0;
b2 = d2 < 0.0;
b3 = d3 < 0.0;
return ((b1 == b2) && (b2 == b3));
}
/**
......@@ -262,13 +367,6 @@ public class GeometryUtils {
return (angle1-angle2) < 0 ? (angle1-angle2) + 2*Math.PI :(angle1-angle2);
}
public static double sign(final IPoint p1, final IPoint p2, final IPoint p3) {
return (p1.getX() - p3.getX()) * (p2.getY() - p3.getY()) - (p2.getX() - p3.getX()) * (p1.getY() - p3.getY());
}
public static double sign(final double x1, final double y1, final double x2, final double y2, final double x3, final double y3) {
return (x1 - x3) * (y2 - y3) - (x2 -x3) * (y1 - y3);
}
public static <P extends IPoint> VRectangle bound(final Collection<P> points) {
......
package org.vadere.util.geometry.mesh.inter;
import org.vadere.util.geometry.shapes.IPoint;
/**
......
......@@ -3,6 +3,7 @@ package org.vadere.util.geometry.mesh.inter;
import org.apache.commons.collections.IteratorUtils;
import org.apache.commons.lang3.tuple.Triple;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.mesh.impl.PFace;
import org.vadere.util.geometry.mesh.impl.PHalfEdge;
import org.vadere.util.geometry.mesh.impl.PMesh;
......@@ -368,6 +369,19 @@ public interface IMesh<P extends IPoint, E extends IHalfEdge<P>, F extends IFace
return streamEdges(face).filter(e -> getVertex(e).distance(x, y) <= epsilon).findAny();
}
// TODO: rename?
default Optional<E> getEdgeCloseToVertex(F face, double x, double y, double epsilon) {
for(E halfEdge : getEdgeIt(face)) {
P p1 = getVertex(halfEdge);
P p2 = getVertex(getPrev(halfEdge));
if(Math.abs(GeometryUtils.ccw(p1.getX(), p1.getY(), p2.getX(), p2.getY(), x, y)) < epsilon) {
return Optional.of(halfEdge);
}
}
return Optional.empty();
}
Collection<P> getVertices();
int getNumberOfVertices();
......
......@@ -241,7 +241,7 @@ public interface ITriConnectivity<P extends IPoint, E extends IHalfEdge<P>, F ex
* @param p the point which splits the triangle
* @param face the triangle face we split
*
* returns a half-edge which has p as its end vertex
* returns a list of all newly created face.
*/
default List<F> splitTriangle(@NotNull F face, P p, boolean legalize) {
assert isTriangle(face);
......@@ -552,7 +552,7 @@ public interface ITriConnectivity<P extends IPoint, E extends IHalfEdge<P>, F ex
P v2 = getMesh().getVertex(getMesh().getPrev(halfEdge));
// TODO: think about the epsilon, absolute value seems to be a really bad idea!
if(!getMesh().isBoundary(getMesh().getTwinFace(halfEdge)) && Math.abs(GeometryUtils.sign(x, y, v1.getX(), v1.getY(), v2.getX(), v2.getY())) == 0.0) {
if(!getMesh().isBoundary(getMesh().getTwinFace(halfEdge)) && Math.abs(GeometryUtils.ccw(x, y, v1.getX(), v1.getY(), v2.getX(), v2.getY())) == 0.0) {
faces.add(getMesh().getTwinFace(halfEdge));
break;
}
......
package org.vadere.util.geometry.mesh.inter;
import org.vadere.util.geometry.shapes.IPoint;
/**
* @author Benedikt Zoennchen
*/
public interface IVertex<P extends IPoint> extends IPoint {
P getPoint();
@Override
default double getX() {
return getPoint().getX();
}
@Override
default double getY() {
return getPoint().getY();
}
@Override
default IPoint add(final IPoint point) {
return getPoint().add(point);
}
@Override
default IPoint addPrecise(IPoint point) {
return getPoint().addPrecise(point);
}
@Override
default IPoint subtract(IPoint point) {
return getPoint().subtract(point);
}
@Override
default IPoint multiply(IPoint point) {
return getPoint().multiply(point);
}
@Override
default IPoint scalarMultiply(double factor) {
return getPoint().scalarMultiply(factor);
}
@Override
default IPoint rotate(double radAngle) {
return getPoint().rotate(radAngle);
}
@Override
default double scalarProduct(IPoint point) {
return getPoint().scalarProduct(point);
}
@Override
default IPoint norm() {
return getPoint().norm();
}
@Override
default IPoint normZeroSafe() {
return getPoint().normZeroSafe();
}
@Override
default double distance(final IPoint other) {
return getPoint().distance(other);
}
@Override
default double distance(double x, double y) {
return getPoint().distance(x, y);
}
@Override
default double distanceToOrigin() {
return getPoint().distanceToOrigin();
}
}
package org.vadere.util.geometry.mesh.triangulations;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.mesh.inter.IFace;
import org.vadere.util.geometry.mesh.inter.IHalfEdge;
import org.vadere.util.geometry.mesh.inter.IMesh;
......@@ -9,7 +10,9 @@ import org.vadere.util.geometry.mesh.inter.ITriangulation;
import org.vadere.util.geometry.shapes.IPoint;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
......@@ -30,11 +33,14 @@ public class DelaunayHierarchy<P extends IPoint, E extends IHalfEdge<P>, F exten
private List<Map<E, E>> hierarchyConnector;
private ITriConnectivity<P, E, F> base;
private ITriangulation<P, E, F> base;
private Supplier<ITriangulation<P, E, F>> triangulationSupplier;
private double alpha;
// see delaunay-hierarchy paper!
private double alpha = 40;
private double epsilon = 0.00001;
private Random random;
......@@ -43,6 +49,8 @@ public class DelaunayHierarchy<P extends IPoint, E extends IHalfEdge<P>, F exten
this.hierarchyConnector = new ArrayList<>();
this.random = new Random();
this.triangulationSupplier = triangulationSupplier;
this.base = base;
this.base.init();
hierarchySets.add(base);
hierarchyConnector.add(new HashMap<>());
......@@ -55,15 +63,18 @@ public class DelaunayHierarchy<P extends IPoint, E extends IHalfEdge<P>, F exten
public void flipEdgeEvent(F f1, F f2) {}
@Override
public void insertEvent(final E halfEdge) {
public void insertEvent(@NotNull final E halfEdge) {
P vertex = base.getMesh().getVertex(halfEdge);
E lastEdge = halfEdge;
for(int i = 1; i < hierarchySets.size(); ++i) {
int limit = hierarchySets.size();
for(int i = 1; i <= limit; ++i) {
if(random.nextDouble() < alpha) {
if(hierarchySets.size() <= i) {
hierarchySets.add(triangulationSupplier.get());
ITriangulation<P, E, F> triangulation = triangulationSupplier.get();
triangulation.init();
hierarchySets.add(triangulation);
}
E edge = hierarchySets.get(i).insert(vertex);
......@@ -72,7 +83,7 @@ public class DelaunayHierarchy<P extends IPoint, E extends IHalfEdge<P>, F exten
hierarchyConnector.add(new HashMap<>());
}
hierarchyConnector.get(i-1).put(lastEdge, edge);
hierarchyConnector.get(i-1).put(edge, lastEdge);
lastEdge = edge;
}
......@@ -89,24 +100,56 @@ public class DelaunayHierarchy<P extends IPoint, E extends IHalfEdge<P>, F exten
@Override
public Collection<F> locatePoint(P point, boolean insertion) {
return null;
Optional<F> optFace = locate(point);
if(optFace.isPresent()) {
F face = optFace.get();
if(!insertion) {
return Collections.singleton(face);
}
else {
Optional<E> optEdge = base.getMesh().getMemberEdge(face, point.getX(), point.getY(), epsilon);
// ignore point
if(optEdge.isPresent()) {
return Collections.emptyList();
}
else {
optEdge = base.getMesh().getEdgeCloseToVertex(face, point.getX(), point.getY(), epsilon);
if(optEdge.isPresent()) {
return Arrays.asList(base.getMesh().getFace(optEdge.get()), base.getMesh().getTwinFace(optEdge.get()));
}
else {
return Collections.singleton(face);
}
}
}
}
else {
return Collections.emptyList();
}
}
@Override
public Optional<F> locate(final P point) {
Optional<F> optStartFace = Optional.empty();
for(int i = hierarchySets.size()-1; i >= 0; --i) {
if(!optStartFace.isPresent()) {
if(i == hierarchySets.size()-1) {
optStartFace = hierarchySets.get(i).locate(point.getX(), point.getY());
}
else {
E edge = getNearestPoint(hierarchySets.get(i-1), optStartFace.get(), point);
E newEdge = hierarchyConnector.get(i-1).get(edge);
optStartFace = hierarchySets.get(i).locate(point.getX(), point.getY(), hierarchySets.get(i).getMesh().getFace(newEdge));
if(!optStartFace.isPresent()) {
return Optional.empty();
if(i > 0) {
E edge = getNearestPoint(hierarchySets.get(i+1), optStartFace.get(), point);
E newEdge = hierarchyConnector.get(i).get(edge);
optStartFace = hierarchySets.get(i).locate(point.getX(), point.getY(), hierarchySets.get(i).getMesh().getFace(newEdge));
}
else {
E edge = getNearestPoint(hierarchySets.get(i+1), optStartFace.get(), point);
E newEdge = hierarchyConnector.get(i).get(edge);
return base.locate(point.getX(), point.getY(), base.getMesh().getFace(newEdge));
}
}
}
......
......@@ -44,8 +44,7 @@ public class DelaunayTree<P extends IPoint, E extends IHalfEdge<P>, F extends IF
@Override
public Collection<F> locatePoint(final P point, final boolean insertion) {
checkRoot();
Set<DAG<DAGElement<P, F>>> leafs = new HashSet<>();
LinkedList<DAG<DAGElement<P, F>>> nodesToVisit = new LinkedList<>();
nodesToVisit.add(dag);
......
......@@ -3,6 +3,7 @@ package org.vadere.util.geometry.mesh.triangulations;
import org.apache.commons.collections.IteratorUtils;
import org.apache.log4j.LogManager;
import org.apache.log4j.Logger;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.mesh.impl.PFace;
import org.vadere.util.geometry.mesh.iterators.FaceIterator;
......@@ -71,6 +72,7 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
this.points = points;
this.illegalPredicate = illegalPredicate;
this.bound = GeometryUtils.bound(points);
this.finalized = false;
}
public IncrementalTriangulation(final Set<P> points) {
......@@ -83,6 +85,7 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
this.points = new HashSet<>();
this.illegalPredicate = illegalPredicate;
this.bound = bound;
this.finalized = false;
}
public IncrementalTriangulation(final VRectangle bound) {
......@@ -97,10 +100,6 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
this.mesh = mesh;
}
public F getSuperTriangle() {
return superTriangle;
}
@Override
public void init() {
double gap = 1.0;
......@@ -117,7 +116,6 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
he1 = borderEdges.get(1);
he2 = borderEdges.get(2);
this.finalized = false;
this.initialized = true;
}
......@@ -138,6 +136,10 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
@Override
public E insert(P point) {
if(!initialized) {
init();
}
Collection<F> faces = this.pointLocator.locatePoint(point, true);
int numberOfFaces = faces.size();
......@@ -151,8 +153,8 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
else if(faces.size() == 1) {
log.info("splitTriangle:" + point);
F face = faces.iterator().next();
splitTriangle(face, point, true);
insertedEdge = mesh.getEdge(point);
List<F> createdFaces = splitTriangle(face, point, true);
insertedEdge = mesh.getMemberEdge(createdFaces.get(0), point.getX(), point.getY()).get();
insertEvent(insertedEdge);
} // point lies on an edge of 2 triangles
......@@ -389,7 +391,7 @@ public class IncrementalTriangulation<P extends IPoint, E extends IHalfEdge<P>,
}
@Override
public void insertEvent(final E halfEdge) {
public void insertEvent(@NotNull final E halfEdge) {
pointLocator.insertEvent(halfEdge);
}
......
......@@ -79,6 +79,12 @@ public class UniformRefinementTriangulation<P extends IPoint> {
}
}
removeTrianglesInsideObstacles();
logger.info("end computation");
}
private void removeTrianglesInsideObstacles() {
for(VShape shape : boundary) {
// 1. find a triangle inside the boundary
VPoint centroid = shape.getCentroid();
......@@ -102,7 +108,6 @@ public class UniformRefinementTriangulation<P extends IPoint> {
logger.warn("no face found");
}
}
logger.info("end computation");
}
public Set<VLine> getEdges() {
......
......@@ -58,22 +58,14 @@ public class VTriangle extends VPolygon {
@Override
public boolean contains(final IPoint point) {
boolean b1, b2, b3;
double d1 = GeometryUtils.sign(point, p1, p2);
double d2 = GeometryUtils.sign(point, p2, p3);
double d3 = GeometryUtils.sign(point, p3, p1);
b1 = d1 < 0.0;
b2 = d2 < 0.0;
b3 = d3 < 0.0;
return ((b1 == b2) && (b2 == b3));
return GeometryUtils.contains(p1, p2, p3, point);
}
// TODO: find better name
public boolean isPartOf(final IPoint point, final double eps) {
boolean b1, b2, b3;
double d1 = GeometryUtils.sign(point, p1, p2);
double d2 = GeometryUtils.sign(point, p2, p3);
double d3 = GeometryUtils.sign(point, p3, p1);
double d1 = GeometryUtils.ccw(point, p1, p2);
double d2 = GeometryUtils.ccw(point, p2, p3);
double d3 = GeometryUtils.ccw(point, p3, p1);
return (d1 <= eps && d2 <= eps && d3 <= eps) || (d1 >= -eps && d2 >= -eps && d3 >= -eps);
}
......
......@@ -137,7 +137,7 @@ public class TestBoyerWatson {
long ms = System.currentTimeMillis();
ITriangulation<VPoint, PHalfEdge<VPoint>, PFace<VPoint>> delaunayTriangulation = ITriangulation.createPTriangulation(IPointLocator.Type.DELAUNAY_TREE, points, (x, y) -> new VPoint(x, y));
delaunayTriangulation.finalize();
log.info("runtime of the BowyerWatson for " + numberOfPoints + " vertices =" + (System.currentTimeMillis() - ms));
log.info("runtime of the BowyerWatson for " + numberOfPoints + " vertices =" + (System.currentTimeMillis() - ms) + " using the delaunay-tree");
log.info("start checking the delaunay property, this can take some time");
Collection<VTriangle> triangles = delaunayTriangulation.streamTriangles().collect(Collectors.toList());
......@@ -151,6 +151,24 @@ public class TestBoyerWatson {
}
}
log.info("end checking the delaunay property");
ms = System.currentTimeMillis();
delaunayTriangulation = ITriangulation.createPTriangulation(IPointLocator.Type.DELAUNAY_HIERARCHY, points, (x, y) -> new VPoint(x, y));
delaunayTriangulation.finalize();
log.info("runtime of the BowyerWatson for " + numberOfPoints + " vertices =" + (System.currentTimeMillis() - ms) + " using the delaunay-hierarchy");
log.info("start checking the delaunay property, this can take some time");
triangles = delaunayTriangulation.streamTriangles().collect(Collectors.toList());
for(VTriangle triangle : triangles) {
List<VPoint> trianglePoints = triangle.getPoints();
for(VTriangle t : triangles) {
assertTrue(t.getPoints().stream().noneMatch(p -> !trianglePoints.contains(p) && triangle.isInCircumscribedCycle(p)));
}
}
log.info("end checking the delaunay property");
}