Commit 09abb91f authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

much faster implementation of the BowyerWatson algorithm and a splitting...

much faster implementation of the BowyerWatson algorithm and a splitting algorithm for obscuse triangles
parent 210a40b1
......@@ -64,7 +64,7 @@ public abstract class DefaultModel<T extends DefaultConfig> extends Observable i
public T config;
public List<VTriangle> triangulation;
public Collection<? extends VTriangle> triangulation;
public DefaultModel(final T config) {
this.config = config;
......@@ -504,7 +504,7 @@ public abstract class DefaultModel<T extends DefaultConfig> extends Observable i
/*
* returns the adaptive triangulation (see persson-2004 'A Simple Mesh Generator in MATLAB.')
*/
public Collection<VTriangle> getTriangulation() {
public Collection<? extends VTriangle> getTriangulation() {
if(triangulation.isEmpty()) {
PerssonStrangDistmesh psd = new PerssonStrangDistmesh(
new VRectangle(getTopographyBound()),
......@@ -513,7 +513,7 @@ public abstract class DefaultModel<T extends DefaultConfig> extends Observable i
false,
l -> 0.0,
"Distmesh");
triangulation = psd.getTriangulation().getVTriangles();
triangulation = psd.getTriangles();
}
return triangulation;
}
......
package org.vadere.util.data;
import org.jetbrains.annotations.NotNull;
import java.util.ArrayList;
import java.util.Collection;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.function.Predicate;
public class DAG<E> {
private final E root;
private final List<DAG<E>> children;
public DAG(@NotNull final E element) {
root = element;
children = new ArrayList<>();
}
public List<DAG<E>> getChildren() {
return children;
}
public void addChild(E child) {
this.children.add(new DAG<E>(child));
}
public void addChild(DAG<E> child) {
this.children.add(child);
}
public E getRoot() {
return root;
}
public Collection<E> collectLeafs() {
Collection<E> leafs = new ArrayList<E>();
LinkedList<DAG<E>> nodesToVisit = new LinkedList<>();
nodesToVisit.add(this);
while (!nodesToVisit.isEmpty()) {
DAG<E> currentNode = nodesToVisit.removeFirst();
nodesToVisit.addAll(currentNode.children);
if(currentNode.isLeaf())
leafs.add(currentNode.getRoot());
}
return leafs;
}
public boolean isLeaf(){
return children.isEmpty();
}
/**
* Finds the first DAG-node element in a dept first fashion.
* @param test the predicate the element of the DAG-node has to fullfill.
* @return
*/
public Optional<E> findFirstElement(final Predicate<E> test){
Optional<DAG<E>> optDag = findFirst(test);
if(optDag.isPresent()) {
return Optional.of(optDag.get().getRoot());
}
else {
return Optional.empty();
}
}
/**
* Finds the first DAG-node in a dept first fashion.
* @param test the predicate the element of the DAG-node has to fullfill.
* @return
*/
public Optional<DAG<E>> findFirst(final Predicate<E> test){
if(test.test(root)) {
return Optional.of(this);
}
else {
return children.stream().map(child -> child.findFirst(test)).filter(opt -> opt.isPresent()).map(opt -> opt.get()).findFirst();
}
}
public Optional<DAG<E>> matchAll(final Predicate<E> test) {
if(test.test(root)) {
if(isLeaf()) {
return Optional.of(this);
}
else {
return children.stream().map(child -> child.matchAll(test)).filter(opt -> opt.isPresent()).map(opt -> opt.get()).findFirst();
}
}
return Optional.empty();
}
}
package org.vadere.util.delaunay;
import org.apache.commons.lang3.tuple.ImmutableTriple;
import org.apache.commons.lang3.tuple.Triple;
import org.vadere.util.geometry.LinkedCellsGrid;
import org.jetbrains.annotations.NotNull;
import org.vadere.util.data.DAG;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.shapes.VCircle;
import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VRectangle;
import org.vadere.util.geometry.shapes.VTriangle;
import javax.swing.*;
import java.awt.*;
import java.util.*;
import java.util.Arrays;
import java.util.Collection;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.function.BiFunction;
import java.util.function.Predicate;
import java.util.Optional;
import java.util.Random;
import java.util.Set;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
import java.util.stream.Stream;
/**
* @author Benedikt Zoennchen
*
* This class is for computing the DelaunayTriangulation using the BowyerWatson-Algorithm. In average the algorithm should perfom in O(n log(n)) but
* in degenerated cases its runtime can be in O(n^2) where n is the number of points.
*/
public class BowyerWatson<P extends VPoint> {
private List<Triple<P, P, P>> triangles;
private Collection<P> points;
private List<P> initPoints;
private final BiFunction<Double, Double, P> pointConstructor;
public BowyerWatson(final Collection<P> points, final BiFunction<Double, Double, P> pointConstructor) {
this.points = points;
this.pointConstructor = pointConstructor;
}
public void execute() {
P max = points.parallelStream().reduce(pointConstructor.apply(Double.MIN_VALUE, Double.MIN_VALUE), (a, b) -> pointConstructor.apply(Math.max(a.getX(), b.getX()), Math.max(a.getY(), b.getY())));
P min = points.parallelStream().reduce(pointConstructor.apply(Double.MIN_VALUE, Double.MIN_VALUE), (a, b) -> pointConstructor.apply(Math.min(a.getX(), b.getX()), Math.min(a.getY(), b.getY())));
VRectangle bound = new VRectangle(min.getX(), min.getY(), max.getX()-min.getX(), max.getY()- min.getY());
init(bound);
points.stream().forEach(point -> handle(point));
cleanUp();
}
public List<Triple<P, P, P>> getTriangles() {
return triangles;
}
public List<VTriangle> getVTriangles() {
return triangles.stream().map(this::pointsToTriangle).collect(Collectors.toList());
}
public Set<VLine> getEdges() {
return triangles.parallelStream().flatMap(triangle -> Stream.of(new VTriangle(triangle.getLeft(), triangle.getMiddle(), triangle.getRight()).getLines())).collect(Collectors.toSet());
}
private void init(final VRectangle bound) {
triangles = new ArrayList<>();
initPoints = new ArrayList<>();
Triple<P, P, P> superTriangle = getSuperTriangle(bound);
triangles.add(superTriangle);
initPoints.add(superTriangle.getLeft());
initPoints.add(superTriangle.getMiddle());
initPoints.add(superTriangle.getRight());
}
private Triple<P, P, P> getSuperTriangle(final VRectangle bound) {
double gap = 1.0;
double max = Math.max(bound.getWidth(), bound.getHeight());
P p1 = pointConstructor.apply(bound.getX() - max - gap, bound.getY() - gap);
P p2 = pointConstructor.apply(bound.getX() + 2 * max + gap, bound.getY() - gap);
P p3 = pointConstructor.apply(bound.getX() + (max+2*gap)/2, bound.getY() + 2 * max+ gap);
return ImmutableTriple.of(p1, p2, p3);
}
private void handle(final P point) {
HashSet<Line> edges = new HashSet<>();
Map<Boolean, List<Triple<P, P, P>>> partition = triangles.parallelStream().collect(Collectors.partitioningBy(t -> pointsToTriangle(t).isInCircumscribedCycle(point)));
List<Triple<P, P, P>> badTriangles = partition.get(true);
triangles = partition.get(false);
IntStream s;
HashSet<Line> toRemove = new HashSet<>();
// duplicated edges
badTriangles.stream().flatMap(t -> getEdges(t).stream()).forEach(line -> {
if(!edges.add(line)) {
toRemove.add(line);
}
});
toRemove.stream().forEach(removeEdge -> edges.remove(removeEdge));
// identifier ?
edges.stream().forEach(edge -> triangles.add(Triple.of(edge.p1, edge.p2, point)));
}
private List<Line> getEdges(Triple<P, P, P> triangle) {
List<Line> list = new ArrayList<>();
list.add(new Line(triangle.getLeft(), triangle.getMiddle()));
list.add(new Line(triangle.getMiddle(), triangle.getRight()));
list.add(new Line(triangle.getRight(), triangle.getLeft()));
return list;
}
private void cleanUp() {
triangles = triangles.stream().filter(triangle -> !isTriangleConnectedToInitialPoints(triangle)).collect(Collectors.toList());
}
private boolean isTriangleConnectedToInitialPoints(final Triple<P, P, P> trianglePoints) {
return Stream.of(pointsToTriangle(trianglePoints).getLines()).anyMatch(edge -> {
VPoint p1 = new VPoint(edge.getP1().getX(), edge.getP1().getY());
VPoint p2 = new VPoint(edge.getP2().getX(), edge.getP2().getY());
return initPoints.stream().anyMatch(initPoint -> p1.equals(initPoint) || p2.equals(initPoint));
});
}
private VTriangle pointsToTriangle(Triple<P, P, P> points) {
return new VTriangle(points.getLeft(), points.getMiddle(), points.getRight());
}
private class Line {
final P p1;
final P p2;
private Line(P p1, P p2) {
this.p1 = p1;
this.p2 = p2;
}
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
Line line = (Line) o;
return (p1.equals(line.p1) && p2.equals(line.p2)) || (p2.equals(line.p1) && p1.equals(line.p2));
}
@Override
public int hashCode() {
return p1.hashCode() * p2.hashCode();
}
}
import javax.swing.*;
public class BowyerWatson<P extends VPoint, T extends VTriangle> {
private final Collection<P> points;
private final PointConstructor<P> pointConstructor;
private final TriangleConstructor<P, T> triangleConstructor;
private P p0;
private P p1;
private P p2;
private DAG<DAGElement<P, T>> dag;
private final HashMap<Face<P>, DAG<DAGElement<P, T>>> map;
public BowyerWatson(
final Collection<P> points,
final PointConstructor<P> pointConstructor,
final TriangleConstructor<P, T> triangleConstructor) {
this.points = points;
this.pointConstructor = pointConstructor;
this.triangleConstructor = triangleConstructor;
this.map = new HashMap<>();
}
public void init() {
P p_max = points.parallelStream().reduce(pointConstructor.create(Double.MIN_VALUE, Double.MIN_VALUE), (a, b) -> pointConstructor.create(Math.max(a.getX(), b.getX()), Math.max(a.getY(), b.getY())));
P p_min = points.parallelStream().reduce(pointConstructor.create(Double.MIN_VALUE, Double.MIN_VALUE), (a, b) -> pointConstructor.create(Math.min(a.getX(), b.getX()), Math.min(a.getY(), b.getY())));
VRectangle bound = new VRectangle(p_min.getX(), p_min.getY(), p_max.getX()-p_min.getX(), p_max.getY()- p_min.getY());
double gap = 1.0;
double max = Math.max(bound.getWidth(), bound.getHeight());
p0 = pointConstructor.create(bound.getX() - max - gap, bound.getY() - gap);
p1 = pointConstructor.create(bound.getX() + 2 * max + gap, bound.getY() - gap);
p2 = pointConstructor.create(bound.getX() + (max+2*gap)/2, bound.getY() + 2 * max+ gap);
Face superTriangle = Face.of(p0, p1, p2);
this.dag = new DAG<>(new DAGElement<>(superTriangle, Triple.of(p0, p1, p2), triangleConstructor));
this.map.put(superTriangle, dag);
}
public void execude() {
// 1. compute a random permutation of the point set/collection.
for(P p : points) {
// find triangle containing p using the DAG.
Optional<DAG<DAGElement<P, T>>> opt = dag.matchAll(dagElement -> dagElement.getTriangle().contains(p));
assert opt.isPresent();
DAG<DAGElement<P, T>> trianglePair = opt.get();
T vTriangle = trianglePair.getRoot().getTriangle();
Face<P> fTriangle = trianglePair.getRoot().getFace();
// p lies in the interior of the triangle
if(vTriangle.contains(p)){
// split the triangle into 3 new triangles.
split(p, trianglePair);
}
else {
throw new IllegalStateException("this should never happen!");
}
}
}
public Collection<T> getTriangles() {
return map.values().stream().map(dagElement -> dagElement.getRoot().getTriangle()).collect(Collectors.toList());
}
public Set<VLine> getEdges() {
return map.values().stream().flatMap(dagElement -> dagElement.getRoot().getTriangle().getLineStream()).collect(Collectors.toSet());
}
/**
* Splits the triangle xyz into three new triangles xyp, yzp and zxp.
*
* @param p
* @param xyzDag
*/
public DAG<DAGElement<P, T>> split(@NotNull P p, @NotNull DAG<DAGElement<P, T>> xyzDag) {
VTriangle xyzTriangle = xyzDag.getRoot().getTriangle();
Face<P> xyzFace = xyzDag.getRoot().getFace();
List<HalfEdge<P>> edges = xyzFace.getEdges();
Face xyp = new Face();
Face yzp = new Face();
Face zxp = new Face();
HalfEdge<P> zx = edges.get(0);
HalfEdge<P> xy = edges.get(1);
HalfEdge<P> yz = edges.get(2);
P x = zx.getEnd();
P y = xy.getEnd();
P z = yz.getEnd();
HalfEdge<P> yp = new HalfEdge<P>(p, xyp);
HalfEdge<P> py = new HalfEdge<P>(y, yzp);
yp.setTwin(py);
HalfEdge<P> xp = new HalfEdge<P>(p, zxp);
HalfEdge<P> px = new HalfEdge<P>(x, xyp);
xp.setTwin(px);
HalfEdge<P> zp = new HalfEdge<P>(p, yzp);
HalfEdge<P> pz = new HalfEdge<P>(z, zxp);
zp.setTwin(pz);
zx.setNext(xp);
xp.setNext(pz);
pz.setNext(zx);
xy.setNext(yp);
yp.setNext(px);
px.setNext(xy);
yz.setNext(zp);
zp.setNext(py);
py.setNext(yz);
xyp.setEdge(yp);
yzp.setEdge(py);
zxp.setEdge(xp);
xy.setFace(xyp);
zx.setFace(zxp);
yz.setFace(yzp);
DAG<DAGElement<P, T>> xypDag = new DAG<>(new DAGElement(xyp, Triple.of(x, y , p), triangleConstructor));
DAG<DAGElement<P, T>> yzpDag = new DAG<>(new DAGElement(yzp, Triple.of(y, z, p), triangleConstructor));
DAG<DAGElement<P, T>> zxpDag = new DAG<>(new DAGElement(zxp, Triple.of(z, x, p), triangleConstructor));
xyzDag.addChild(xypDag);
xyzDag.addChild(yzpDag);
xyzDag.addChild(zxpDag);
map.remove(xyzDag.getRoot().getFace());
map.put(xyp, xypDag);
map.put(yzp, yzpDag);
map.put(zxp, zxpDag);
legalize(p, zx);
legalize(p, xy);
legalize(p, yz);
return xyzDag;
}
public DAG<DAGElement<P, T>> getDag() {
return dag;
}
/**
* Checks if the edge xy of the triangle xyz is illegal, which is the case if:
* There is a point p and a triangle yxp and p is in the circumscribed cycle of xyz.
*
* @param p
* @param edge
* @return
*/
private boolean isIllegalEdge(P p, HalfEdge<P> edge){
if(edge.hasTwin()) {
P x = edge.getTwin().get().getEnd();
P y = edge.getTwin().get().getNext().getEnd();
P z = edge.getTwin().get().getNext().getNext().getEnd();
VTriangle triangle = new VTriangle(x,y,z);
return triangle.isInCircumscribedCycle(p);
}
return false;
}
/**
* Legalizes an edge xy of a triangle xyz where xy is the twin of yx of the triangle
* yxp and p is the point on the other side of xy. So legalizing means to flip xy with
* zp.
*
* @param p p is the point on the other side of xy
* @param edge an edge xy of a triangle xyz
*/
private void legalize(final P p, final HalfEdge<P> edge) {
if(isIllegalEdge(p, edge)) {
// read all informations
HalfEdge<P> zx = edge;
HalfEdge<P> xp = zx.getNext();
HalfEdge<P> pz = xp.getNext();
Face<P> f2 = zx.getFace();
HalfEdge<P> xz = zx.getTwin().get();
HalfEdge<P> zy = xz.getNext();
HalfEdge<P> yx = zy.getNext();
P x = zx.getEnd();
P y = zy.getEnd();
P z = pz.getEnd();
Face<P> f1 = xz.getFace();
DAG<DAGElement<P, T>> xpzDag = map.get(f2);
DAG<DAGElement<P, T>> xzyDag = map.get(f1);
// update data structure
HalfEdge<P> py = new HalfEdge<P>(y, f2);
HalfEdge<P> yp = new HalfEdge<P>(p, f1);
py.setTwin(yp);
xp.setNext(py);
py.setNext(yx);
yx.setNext(xp);
xp.setFace(f2);
py.setFace(f2);
yx.setFace(f2);
f2.setEdge(py);
pz.setNext(zy);
zy.setNext(yp);
yp.setNext(pz);
f1.setEdge(yp);
pz.setFace(f1);
zy.setFace(f1);
yp.setFace(f1);
// update DAG
DAG<DAGElement<P, T>> yzpDag = new DAG<>(new DAGElement(f1, Triple.of(y, z, p), triangleConstructor));
DAG<DAGElement<P, T>> xypDag = new DAG<>(new DAGElement(f2, Triple.of(x, y, p), triangleConstructor));
xpzDag.addChild(yzpDag);
xpzDag.addChild(xypDag);
xzyDag.addChild(yzpDag);
xzyDag.addChild(xypDag);
map.put(f1, yzpDag);
map.put(f2, xypDag);
legalize(p, py.getNext());
legalize(p, yp.getPrevious());
}
}
/*public class DAGElement {
private Face<P> face;
private Triple<P, P, P> vertices;
private T triangle;
public DAGElement(final Face<P> face, final Triple<P, P, P> vertices) {
this.face = face;
this.vertices = vertices;
this.triangle = triangleConstructor.create(vertices.getLeft(), vertices.getMiddle(), vertices.getRight());
}
public Face<P> getFace() {
return face;
}
public T getTriangle() {
return triangle;
}
public Triple<P, P, P> getVertices() {
return vertices;
}
}*/
// TODO: the following code can be deleted, this is only for visual checks
public static void main(String[] args) {
// TODO Auto-generated method stub
int height = 1000;
int width = 1000;
int max = Math.max(height, width);
// TODO: the following code can be deleted, this is only for visual checks
public static void main(String[] args) {
// TODO Auto-generated method stub
int height = 700;
int width = 700;
int max = Math.max(height, width);
Set<VPoint> points = new HashSet<>();
Set<VPoint> points = new HashSet<>();
/*points.add(new VPoint(20,20));
points.add(new VPoint(20,40));
points.add(new VPoint(75,53));
points.add(new VPoint(80,70));*/
Random r = new Random();
for(int i=0; i<100; i++) {
VPoint point = new VPoint(width*r.nextDouble(), height*r.nextDouble());