Commit 2db821e0 authored by Benedikt Zoennchen's avatar Benedikt Zoennchen
Browse files

add adpative triangulation implementation

parent 43e77c48
package org.vadere.util.delaunay;
import org.apache.commons.lang3.tuple.ImmutablePair;
import org.apache.commons.lang3.tuple.ImmutableTriple;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.commons.lang3.tuple.Triple;
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.List;
import java.util.function.BiFunction;
import java.util.function.Predicate;
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();
}
}
// 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);
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<10000; i++) {
VPoint point = new VPoint(width*r.nextDouble(), height*r.nextDouble());
points.add(point);
}
BowyerWatson<VPoint> bw = new BowyerWatson<VPoint>(points, (x, y) -> new VPoint(x, y));
bw.execute();
Set<VLine> edges = bw.getEdges();
JFrame window = new JFrame();
window.setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE);
window.setBounds(0, 0, max, max);
window.getContentPane().add(new Lines(edges, points, max));
window.setVisible(true);
}
private static class Lines extends JComponent{
private Set<VLine> edges;
private Set<VPoint> points;
private final int max;
public Lines(final Set<VLine> edges, final Set<VPoint> points, final int max){
this.edges = edges;
this.points = points;
this.max = max;
}
public void paint(Graphics g) {
Graphics2D g2 = (Graphics2D) g;
g2.setBackground(Color.white);
g2.setStroke(new BasicStroke(1.0f));
g2.setColor(Color.gray);
edges.stream().forEach(edge -> {
Shape k = new VLine(edge.getP1().getX(), edge.getP1().getY(), edge.getP2().getX(), edge.getP2().getY());
g2.draw(k);
});
points.stream().forEach(point -> {
VCircle k = new VCircle(point.getX(), point.getY(), 1.0);
g2.draw(k);
});
}
}
public void removeTriangleIf(final Predicate<Triple<P, P, P>> predicate) {
triangles.removeIf(predicate);
}
}
......@@ -171,6 +171,10 @@ public class VTriangle extends VPolygon {
return getCircumcenter().distance(point) < getCircumscribedRadius();
}
public boolean isInCircumscribedCycle(final VPoint point) {
return getCircumcenter().distance(point) < getCircumscribedRadius();
}
/**
* Computes the inward facing normal vector for the given points of the
* triangle.
......@@ -196,10 +200,6 @@ public class VTriangle extends VPolygon {
return this.p1.equals(point) || this.p2.equals(point) || this.p3.equals(point);
}
public VLine[] getLines() {
return lines;
}
public Stream<VLine> getLineStream() {
return Arrays.stream(getLines());
}
......@@ -215,4 +215,8 @@ public class VTriangle extends VPolygon {
public String toString() {
return p1 + "-" + p2 + "-" + p3;
}
public VLine[] getLines() {
return lines;
}
}
package org.vadere.util.triangulation.adaptive;
import org.apache.commons.lang3.tuple.Pair;
import org.vadere.util.geometry.shapes.VLine;
public class IndexedVLine extends VLine {
private final IndexedVPoint p1;
private final IndexedVPoint p2;
public IndexedVLine(IndexedVPoint p1, IndexedVPoint p2) {
super(p1, p2);
this.p1 = p1;
this.p2 = p2;
}
@Override
public boolean equals(Object obj) {
if(obj == null) return false;
if(obj == this) return true;
if(obj.getClass() != getClass()) return false;
IndexedVLine line = (IndexedVLine) obj;
return line.getP1().equals(getP1()) && line.getP2().equals(getP2()) || line.getP2().equals(getP1()) && line.getP1().equals(getP2());
}
@Override
public int hashCode() {
return p1.hashCode() * p2.hashCode();
}
public Pair<Integer, Integer> getId() {
return Pair.of(p1.getId(), p2.getId());
}
}
package org.vadere.util.triangulation.adaptive;
import org.vadere.util.geometry.shapes.VPoint;
/**
* Created by bzoennchen on 10.11.16.
*/
public class IndexedVPoint extends VPoint {
private int id;
public IndexedVPoint(final VPoint point, final int id){
this(point.getX(), point.getY(), id);
}
public IndexedVPoint(final double x, final double y, final int id){
super(x, y);
this.id = id;
}
public IndexedVPoint subtract(VPoint point) {
return new IndexedVPoint(super.subtract(point), id);
}
public int getId() {
return id;
}
public void setId(int id) {
this.id = id;
}
}
......@@ -96,7 +96,7 @@ public class PSDistmesh {
}
public boolean hasMaximalSteps() {
return steps >= Parameters.MAX_NUMBER_OF_STEPS;
return steps >= 1000;
}
/**
......
......@@ -5,17 +5,16 @@ package org.vadere.util.triangulation.adaptive;
*/
public class Parameters {
final static double TOL = .1;
final static double FSCALE = 1.0;
final static double FSCALE = 1.2;
final static double DELTAT = 0.2;
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.95;
public final static String method = "Density"; // "Distmesh" or "Density"
final static double qualityMeasurement = 0.875;
final static double MINIMUM = 0.25;
final static double DENSITYWEIGHT = 2;
final static int NPOINTS = 100000;
final static int SAMPLENUMBER = 10;
final static int SAMPLEDIVISION = 10;
static final int SEGMENTDIVISION = 0;
final static int MAX_NUMBER_OF_STEPS = 20;
}
package org.vadere.util.triangulation.adaptive;
import org.vadere.util.triangulation.BowyerWatson3;
import org.apache.commons.lang3.tuple.Pair;
import org.vadere.util.delaunay.BowyerWatson;
import org.vadere.util.geometry.shapes.VLine;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VRectangle;
......@@ -18,9 +18,9 @@ import java.util.stream.IntStream;
import static java.lang.Math.pow;
public class PerssonStrangDistmesh {
private List<VPoint> points = new ArrayList<>();
private List<VPoint> oldPoints = new ArrayList<>();
private BowyerWatson3 triangulation;
private List<IndexedVPoint> points = new ArrayList<>();
private List<IndexedVPoint> oldPoints = new ArrayList<>();
private BowyerWatson<IndexedVPoint> triangulation;
private Function<VPoint, Double> fd;
private Function<VPoint, Double> fh;
......@@ -38,6 +38,7 @@ public class PerssonStrangDistmesh {
boolean uniform,
Function<VPoint, Double> density,
String method) {
long now = System.currentTimeMillis();
this.h0 = h0;
this.geps = .001*h0;
this.deps = 1.4901e-8*h0;
......@@ -87,10 +88,6 @@ public class PerssonStrangDistmesh {
addAll(obstacles);
add(box);
}});
}
public void execude() {
long now = System.currentTimeMillis();
setOldPointsToInf();
work();
Date date = new Date(System.currentTimeMillis() - now);
......@@ -128,36 +125,22 @@ public class PerssonStrangDistmesh {
{
while(true)
{
//double maxMove = largeMovement();
if(largeMovement())
{
copyPoints();
long now = System.currentTimeMillis();
triangulation = new BowyerWatson3(points);
triangulation = new BowyerWatson<>(points, (x, y) -> new IndexedVPoint(x, y, -1));
triangulation.execute();
//System.out.println("different edges:" + triangulation.getTriangles().stream().flatMap(t -> t.getLineStream()).collect(Collectors.toSet()).size());
Set<VLine> setLine = new HashSet<>(triangulation.getTriangles().stream().flatMap(t -> t.getLineStream()).collect(Collectors.toList()));
System.out.println("different edges:" + setLine.size());
//triangulation.setTriangles(triangulation.getTriangles().parallelStream().filter(p ->
// fd.apply(p.midPoint()) < -geps).collect(Collectors.toList()));
int size = triangulation.getTriangles().size();
/*triangulation.removeTriangleIf(triple -> fd.apply(
new VTriangle(triple.getLeft(), triple.getMiddle(), triple.getRight()).midPoint()) >= -geps);*/
triangulation.setTriangles(triangulation.getTriangles().parallelStream().filter(p ->
fd.apply(p.midPoint()) < -geps).collect(Collectors.toList()));
System.out.println("deleted:" + (size - triangulation.getTriangles().size()));
triangulation.removeTriangleIf(t -> fd.apply(new VTriangle(t.getLeft(), t.getMiddle(), t.getRight()).midPoint()) < -geps);
Date date = new Date(System.currentTimeMillis() - now);
System.out.println(new SimpleDateFormat("mm:ss:SSS").format(date) + " Bowyer-Watson-Algo");
}
long now = System.currentTimeMillis();
HashMap<VLine, Integer> reversedLines = createBars();
HashMap<IndexedVLine, Integer> reversedLines = createBars();
Date date = new Date(System.currentTimeMillis() - now);
System.out.println(new SimpleDateFormat("mm:ss:SSS").format(date) + " CreateBars");
now = System.currentTimeMillis();
List<VLine> lines = new ArrayList<>(reversedLines.keySet());
List<IndexedVLine> lines = new ArrayList<>(reversedLines.keySet());
date = new Date(System.currentTimeMillis() - now);
System.out.println(new SimpleDateFormat("mm:ss:SSS").format(date) + " CompleteLineList");
now = System.currentTimeMillis();
......@@ -196,12 +179,10 @@ public class PerssonStrangDistmesh {
private void doEuler(double[][] array)
{
points = IntStream.range(0,points.size()).parallel().mapToObj(i -> {
if(points.get(i).identifier != -1) {
return new VPoint(array[i][0] * Parameters.DELTAT + points.get(i).getX(), Parameters.DELTAT * array[i][1] + points.get(i).getY(), i);
} else {
if(points.get(i).getId() != -1)
return new IndexedVPoint(array[i][0] * Parameters.DELTAT + points.get(i).getX(), Parameters.DELTAT * array[i][1] + points.get(i).getY(), i);
else
return points.get(i);
}
}).collect(Collectors.toList());
}
......@@ -221,7 +202,7 @@ public class PerssonStrangDistmesh {
path.next();
path.currentSegment(coordinates);
VPoint[] points = divLine(coordinates[0], coordinates[1], tempCoords[0], tempCoords[1], Parameters.SEGMENTDIVISION);
IndexedVPoint[] points = divLine(coordinates[0], coordinates[1], tempCoords[0], tempCoords[1], Parameters.SEGMENTDIVISION);
if (coordinates[0] == tempCoords[0] && coordinates[1] == tempCoords[1]) {
break;
}
......@@ -234,15 +215,15 @@ public class PerssonStrangDistmesh {
/*
Unterteilt eine Linie eines Objekts in Fixpunkte
*/
private VPoint[] divLine(double x1, double y1, double x2, double y2, int segments)
private IndexedVPoint[] divLine(double x1, double y1, double x2, double y2, int segments)
{
VPoint[] points = new VPoint[segments+1];
IndexedVPoint[] points = new IndexedVPoint[segments+1];
double dX = (x1-x2)/segments;
double dY = (y1-y2)/segments;
for(int i = 1; i < points.length; i++)
{
points[i] = new VPoint(x2 + i*dX, y2 + i*dY, -1);
points[i] = new IndexedVPoint(x2 + i*dX, y2 + i*dY, -1);
}
return points;
}
......@@ -254,7 +235,7 @@ public class PerssonStrangDistmesh {
{
double ave = 0;
int i = 0;
for(VTriangle v : triangulation.getTriangles()) {
for(VTriangle v : triangulation.getVTriangles()) {
VLine[] line = v.getLines();
double a = Math.sqrt(Math.pow(line[0].getX1() - line[0].getX2(), 2) + Math.pow(line[0].getY1() - line[0].getY2(), 2));
double b = Math.sqrt(Math.pow(line[1].getX1() - line[1].getX2(), 2) + Math.pow(line[1].getY1() - line[1].getY2(), 2));
......@@ -299,11 +280,11 @@ public class PerssonStrangDistmesh {
int j = 0;
for(int i : positions)
{
int save = points.get(i).identifier;
int save = points.get(i).getId();
if(save != -1)
{
points.set(i, points.get(i).subtract(new VPoint(fd.apply(points.get(i)) * dGradPX[j], fd.apply(points.get(i)) * dGradPY[j])));
points.get(i).identifier = save;
points.get(i).setId(save);
}
j++;
}
......@@ -312,60 +293,60 @@ public class PerssonStrangDistmesh {
/*
Erzeugt alle Kräfte der Kanten (Teilkräfte von Ftot)
*/
private List<VPoint> generateFVec(List<VLine> lines)
private List<VPoint> generateFVec(List<IndexedVLine> lines)
{
List<Double> L = lines.parallelStream().map(l -> l.getP1().distance(l.getP2())).collect(Collectors.toList());
List<Double> L0 = calculateDesiredLengths(lines, L);
List<Double> F = IntStream.range(0, L0.size()).parallel().mapToDouble(i -> Math.max(L0.get(i) - L.get(i),0)).boxed().collect(Collectors.toList());
return IntStream.range(0, L.size()).parallel().mapToObj(i ->
new VPoint((F.get(i) / L.get(i)) * (lines.get(i).getX1()-lines.get(i).getX2()),
new IndexedVPoint((F.get(i) / L.get(i)) * (lines.get(i).getX1()-lines.get(i).getX2()),
(F.get(i) / L.get(i)) * (lines.get(i).getY1() - lines.get(i).getY2()), i)).collect(Collectors.toList());
}
/*
Berechnet die gewünschten Kantenlängen
*/
private List<Double> calculateDesiredLengths(List<VLine> allLines, List<Double> L)
private List<Double> calculateDesiredLengths(List<IndexedVLine> allLines, List<Double> L)
{
List<Double> hbars = allLines.parallelStream().map(p -> fh.apply(new VPoint((p.x1+p.x2)/2, (p.y1+p.y2)/2))).collect(Collectors.toList());
double sumOfL = L.parallelStream().mapToDouble(s -> Math.pow(s, 2)).sum();
double sumOfhbars = hbars.parallelStream().mapToDouble(s -> Math.pow(s, 2)).sum();
System.out.println("scale factor" + Math.sqrt(sumOfL/sumOfhbars));
return hbars.parallelStream().mapToDouble(h -> h*Parameters.FSCALE*Math.sqrt(sumOfL/sumOfhbars)).boxed().collect(Collectors.toList());
}
/*
Erzeugt die insgesamt wirkenden Kräfte im Szenario
*/
private HashMap<String, Double> generateFtot(HashMap<VLine, Integer> filteredLines, List<VLine> lines, List<VPoint> Fvec)
private HashMap<String, Double> generateFtot(HashMap<IndexedVLine, Integer> filteredLines, List<IndexedVLine> lines, List<VPoint> Fvec)
{
HashMap<String, Double> fTot = new HashMap<>();
for(VLine bar : lines)
{
String[] id = bar.getIdentifier().split(":");
for(int i = 0; i < 2; i++)
for(IndexedVLine bar : lines)
{
Pair<Integer, Integer> ids = bar.getId();
if(!(fTot.containsKey(ids.getLeft()+":0") && fTot.containsKey(ids.getLeft()+":1"))) {
VPoint f = Fvec.get(filteredLines.get(bar));
if(!(fTot.containsKey(id[i]+":0") && fTot.containsKey(id[i]+":1")))
{
if(i==0) {
fTot.put(id[i]+":0", f.getX());
fTot.put(id[i]+":1", f.getY());
} else {
fTot.put(id[i]+":0", -f.getX());
fTot.put(id[i]+":1", -f.getY());
fTot.put(ids.getLeft()+":0", f.getX());
fTot.put