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

some refectoring before merging the master in.

parent 31b50e65
......@@ -23,16 +23,13 @@ import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.shapes.IPoint;
import org.vadere.util.geometry.shapes.VPoint;
import org.vadere.util.geometry.shapes.VShape;
import org.vadere.util.geometry.shapes.VTriangle;
import org.vadere.util.math.IDistanceFunction;
import org.vadere.util.math.InterpolationUtil;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Optional;
import java.util.function.Function;
import java.util.stream.Collectors;
......@@ -40,14 +37,17 @@ import java.util.stream.Collectors;
/**
* @author Benedikt Zoennchen
*
* Abstract class that solves the eikonal equation on a Mesh i.e. 2-D triangular mesh.
* It is the basis for the implementation of the FMM, FIM, FSM and IFIM on triangular meshes.
*
* @param <V> the type of the vertices of the triangulation
* @param <E> the type of the half-edges of the triangulation
* @param <F> the type of the faces of the triangulation
*/
public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge, F extends IFace> implements MeshEikonalSolver<V, E, F>, ITriEventListener<V, E, F> {
@Nullable private ITimeCostFunction timeCostFunction;
@Nullable private ITimeCostFunctionMesh<V> meshTimeCostFunction;
private ITimeCostFunctionMesh<V> meshTimeCostFunction;
@Nullable IDistanceFunction distanceFunction;
private final IIncrementalTriangulation<V, E, F> triangulation;
private Collection<V> initialVertices;
......@@ -79,6 +79,14 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
private MeshEikonalSolver.LocalSover localSover = MeshEikonalSolver.LocalSover.SETHIAN;
/**
* Default constructor.
*
* @param identifier the identifier is used for all container names. A container is basically a bijective function, that maps
* a data point (Object or primitive data type) to a vertex, halfe-edge or face.
* @param triangulation the 2-D triangular mesh i.e. the discretization the solution is based on.
* @param timeCostFunction the timeCostFunction of the eikonal equation, i.e. 1/F(x).
*/
public AMeshEikonalSolver(@NotNull final String identifier,
@NotNull final IIncrementalTriangulation<V, E, F> triangulation,
@NotNull final ITimeCostFunction timeCostFunction) {
......@@ -87,7 +95,7 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
// unsave cast
this.meshTimeCostFunction = (ITimeCostFunctionMesh<V>)timeCostFunction;
} else {
this.timeCostFunction = timeCostFunction;
this.meshTimeCostFunction = ITimeCostFunctionMesh.convert(timeCostFunction);
}
this.triangulation = triangulation;
......@@ -124,6 +132,13 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
}
}
/**
* Marks and returns the initial vertices, i.e. vertices at which the wavefront propagation starts.
*
* @param targetShapes the agent destination/target of the eikonal equation
*
* @return
*/
protected List<V> findInitialVertices(final Collection<VShape> targetShapes) {
//TODO a more clever init!
List<V> initialVertices = new ArrayList<>();
......@@ -142,6 +157,10 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
return initialVertices;
}
/**
* Resets all computed travel times but keeps the computed values which are mesh dependent
* (and therefore, do not change if the mesh does not change), for example, the angles inside a triangle.
*/
protected void unsolve() {
triangulation.getMesh().streamVerticesParallel().filter(v -> !isInitialVertex(v)).forEach(v -> {
setUndefined(v);
......@@ -152,6 +171,9 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
//solved = false;
}
/**
* Computes all mesh dependent measures such as angles and virtual supports.
*/
protected void prepareMesh() {
triangulation.getMesh().streamEdgesParallel().forEach(e -> setVirtualSupport(e, Collections.EMPTY_LIST));
triangulation.getMesh().streamEdgesParallel().forEach(e -> setAccuteEdge(e));
......@@ -173,7 +195,7 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
@Override
public ITimeCostFunction getTimeCostFunction() {
return timeCostFunction != null ? timeCostFunction : meshTimeCostFunction;
return meshTimeCostFunction;
}
@Override
......@@ -306,7 +328,8 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
V v2 = getMesh().getVertex(prev);
double potential = Double.MAX_VALUE;
//TODO: the code below should be used to enable virtual simplexes, however this might currently fail because their is no neighboring connectifity
// for those simplices defined by the mesh and the FIM might stall because it tests the wrong connectivities not part of the DAG!
/*if(isNonAcute(edge)) {
V v = getMesh().getVertex(edge);
List<Pair<V, V>> list = getVirtualSupport(edge);
......@@ -441,6 +464,17 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
return getTimeCostFunction().needsUpdate();
}
/**
* Returns barycentric interpolated value at (x,y) based on the {@link IVertexContainerDouble} containerDouble.
*
* @param triangulation the triangular mesh (this is most of the time {@link this#triangulation})
* @param containerDouble some vertex container
* @param x x-coordinate of the request point
* @param y y-coordinate of the request point
* @param caller the caller object to optimize the triangle walk
*
* @return the barycentric interpolated value at (x,y)
*/
protected double getInterpolatedPotential(
@NotNull final IIncrementalTriangulation<V, E, F> triangulation,
@NotNull final IVertexContainerDouble<V, E, F> containerDouble,
......@@ -613,12 +647,7 @@ public abstract class AMeshEikonalSolver<V extends IVertex, E extends IHalfEdge,
}
protected void setTimeCost(@NotNull final V v) {
if(meshTimeCostFunction != null) {
setTimeCost(v, meshTimeCostFunction.costAt(v));
} else {
setTimeCost(v, timeCostFunction.costAt(getMesh().toPoint(v), v));
}
setTimeCost(v, meshTimeCostFunction.costAt(v));
}
protected void setTimeCost(@NotNull final V v, final double value) {
......
......@@ -22,10 +22,9 @@ public abstract class AMeshEikonalSolverFMM<V extends IVertex, E extends IHalfEd
Comparator<V> pointComparator = (v1, v2) -> {
double alpha = 0.3;
double beta = 1;
double key1 = alpha * getPotential(v1)/* + (1-alpha)*beta*distToDest.apply(getMesh().toPoint(v1))*/;
double key2 = alpha * getPotential(v2)/* + (1-alpha)*beta*distToDest.apply(getMesh().toPoint(v2))*/;
double key1 = alpha * getPotential(v1);
double key2 = alpha * getPotential(v2);
if (key1 < key2) {
return -1;
......
......@@ -18,6 +18,18 @@ import java.util.Optional;
import java.util.Set;
import java.util.function.Predicate;
/**
* This callback class updates the so called agent/pedestrian density every time an agent moves,
* i.e. whenever this listener is called. Therefore, the listener has to be added in the first place.
* {@link DensityUpdater#nameAgentDensity} is the name of the container of the mesh where the agent
* density (a {@link Double}) value is saved.
*
* @param <V> the type of the vertices of the triangulation
* @param <E> the type of the half-edges of the triangulation
* @param <F> the type of the faces of the triangulation
*
* @author Benedikt Zoennchen
*/
public class DensityUpdater<V extends IVertex, E extends IHalfEdge, F extends IFace> implements IMoveDynamicElementListener {
public static final String nameAgentDensity = "agent_density";
......
......@@ -12,6 +12,9 @@ import org.vadere.util.geometry.shapes.VPoint;
public interface MeshEikonalSolver<V extends IVertex, E extends IHalfEdge, F extends IFace> extends EikonalSolver {
/**
* Types of gradient approximations.
*/
enum LocalSover {
SETHIAN, MATRIX;
}
......
......@@ -83,12 +83,12 @@ public class MeshEikonalSolverFIM<V extends IVertex, E extends IHalfEdge, F exte
this.identifier = identifier;
this.activeList = new LinkedList<>();
File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
/*File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
try {
bufferedWriter = IOUtils.getWriter("updates_fim.csv", dir);
} catch (IOException e) {
e.printStackTrace();
}
}*/
//TODO a more clever init!
List<V> initialVertices = new ArrayList<>();
......@@ -113,7 +113,7 @@ public class MeshEikonalSolverFIM<V extends IVertex, E extends IHalfEdge, F exte
double ms = System.currentTimeMillis();
getTriangulation().enableCache();
nUpdates = 0;
narrowBandSizes.add(new ArrayList<>());
//narrowBandSizes.add(new ArrayList<>());
if(!solved || needsUpdate()) {
if(!solved) {
......@@ -129,17 +129,17 @@ public class MeshEikonalSolverFIM<V extends IVertex, E extends IHalfEdge, F exte
}
solved = true;
updates.add(nUpdates);
//updates.add(nUpdates);
double runTime = (System.currentTimeMillis() - ms);
logger.debug("fim run time = " + runTime);
logger.debug("#nUpdates = " + nUpdates);
logger.debug("#nVertices = " + (getMesh().getNumberOfVertices() - (int)getMesh().streamVertices().filter(v -> isInitialVertex(v)).count()));
if(iteration % 100 == 0) {
/*if(iteration % 100 == 0) {
writeNarrowBandSize();
}
if(iteration == 3354) {
writeUpdates();
}
}*/
iteration++;
//logger.debug(getMesh().toPythonTriangulation(v -> getPotential(v)));
}
......@@ -163,9 +163,9 @@ public class MeshEikonalSolverFIM<V extends IVertex, E extends IHalfEdge, F exte
private void march() {
ArrayList<Integer> narrowBandSize=null;
if(iteration % 100 == 0) {
/*if(iteration % 100 == 0) {
narrowBandSize = narrowBandSizes.get(narrowBandSizes.size()-1);
}
}*/
while(!activeList.isEmpty()) {
//logger.debug("#activeList = " + activeList.size());
......
......@@ -25,7 +25,7 @@ import java.util.stream.IntStream;
/**
* This class computes the traveling time T using the fast iterative method for arbitrary triangulated meshes.
* This class computes the traveling time T using a lock free variant of the fast iterative method for arbitrary triangulated meshes.
* The quality of the result depends on the quality of the triangulation. For a high accuracy the triangulation
* should not contain too many non-acute triangles.
*
......
......@@ -57,12 +57,12 @@ public class MeshEikonalSolverFMM<V extends IVertex, E extends IHalfEdge, F exte
super(solver.identifier, triangulation, solver.getTimeCostFunction());
this.identifier = solver.identifier;
setInitialVertices(initialVertices, p -> 0.0);
File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
/*File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
try {
bufferedWriter = IOUtils.getWriter("floorfields.csv", dir);
} catch (IOException e) {
e.printStackTrace();
}
}*/
}
/**
......@@ -80,12 +80,12 @@ public class MeshEikonalSolverFMM<V extends IVertex, E extends IHalfEdge, F exte
super(identifier, triangulation, timeCostFunction);
this.identifier = identifier;
File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
/*File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
try {
bufferedWriter = IOUtils.getWriter("floorfields.csv", dir);
} catch (IOException e) {
e.printStackTrace();
}
}*/
HashSet<V> targetVertices = new HashSet<>();
IDistanceFunction distFunc = p -> IDistanceFunction.createToTargetPoints(targetPoints).apply(p);
......
......@@ -29,6 +29,15 @@ import java.util.function.Predicate;
// TODO this is experimental code!
/**
* This code solves the eikonal equation by iteratively refine the mesh as described in the PhD thesis of B. Zoennchen.
* The implementation is experimental and not jet ready for general usage, i.e. the code is still not very well structured
* and rather hard to read.
*
* @param <V> the type of the vertices of the triangulation
* @param <E> the type of the half-edges of the triangulation
* @param <F> the type of the faces of the triangulation
*/
public class MeshEikonalSolverFMMIterative<V extends IVertex, E extends IHalfEdge, F extends IFace> implements EikonalSolver, ITriEventListener<V, E, F> {
private static Logger logger = Logger.getLogger(MeshEikonalSolverFMMIterative.class);
......
......@@ -33,7 +33,7 @@ import java.util.stream.Collectors;
/**
* This class computes the traveling time T using the fast iterative method for arbitrary triangulated meshes.
* This class computes the traveling time T using the informed fast iterative method (IFIM) for arbitrary triangulated meshes.
* The quality of the result depends on the quality of the triangulation. For a high accuracy the triangulation
* should not contain too many non-acute triangles.
*
......@@ -115,12 +115,12 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
this.speedChange = getMesh().getBooleanVertexContainer(identifier + "_" + nameSpeedChanged);
setInitialVertices(findInitialVertices(targetShapes), IDistanceFunction.createToTargets(targetShapes));
File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
/*File dir = new File("/Users/bzoennchen/Development/workspaces/hmRepo/PersZoennchen/PhD/trash/generated/floorFieldPlot/");
try {
bufferedWriter = IOUtils.getWriter("updates_ifim.csv", dir);
} catch (IOException e) {
e.printStackTrace();
}
}*/
}
......@@ -129,7 +129,7 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
double ms = System.currentTimeMillis();
getTriangulation().enableCache();
nUpdates = 0;
narrowBandSizes.add(new ArrayList<>());
//narrowBandSizes.add(new ArrayList<>());
if(!solved || needsUpdate()) {
if(!solved) {
......@@ -145,17 +145,17 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
}
solved = true;
updates.add(nUpdates);
//updates.add(nUpdates);
double runTime = (System.currentTimeMillis() - ms);
logger.debug("fim run time = " + runTime);
logger.debug("#nUpdates = " + nUpdates);
logger.debug("#nVertices = " + (getMesh().getNumberOfVertices() - (int)getMesh().streamVertices().filter(v -> isInitialVertex(v)).count()));
if(iteration % 100 == 0) {
/*if(iteration % 100 == 0) {
writeNarrowBandSize();
}
if(iteration == 3354) {
writeUpdates();
}
}*/
iteration++;
//logger.debug("#nVertices = " + getMesh().getNumberOfVertices());
//logger.debug(getMesh().toPythonTriangulation(v -> getPotential(v)));
......@@ -227,10 +227,10 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
}
private void march() {
ArrayList<Integer> narrowBandSize=null;
/*ArrayList<Integer> narrowBandSize=null;
if(iteration % 100 == 0) {
narrowBandSize = narrowBandSizes.get(narrowBandSizes.size()-1);
}
}*/
while(!activeList.isEmpty()) {
ListIterator<V> listIterator = activeList.listIterator();
//logger.debug("#activeList = " + activeList.size());
......@@ -271,9 +271,9 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
this.definingSimplex.setValue(xn, Pair.of(triple2.getMiddle(), triple2.getRight()));
setPotential(xn, qq);
newActiveList.add(xn);
if(iteration % 100 == 0) {
/*if(iteration % 100 == 0) {
narrowBandSize.add(newActiveList.size()+activeList.size());
}
}*/
setBurning(xn);
setUnburned(xn);
......@@ -281,9 +281,9 @@ public class MeshEikonalSolverIFIM<V extends IVertex, E extends IHalfEdge, F ext
}
}
listIterator.remove();
if(iteration % 100 == 0) {
/*if(iteration % 100 == 0) {
narrowBandSize.add(newActiveList.size()+activeList.size());
}
}*/
if(updated) {
nUpdates++;
}
......
......@@ -30,7 +30,7 @@ import java.util.stream.IntStream;
/**
* This class computes the traveling time T using the fast iterative method for arbitrary triangulated meshes.
* This class computes the traveling time T using the lock free variant of the fast iterative method for arbitrary triangulated meshes.
* The quality of the result depends on the quality of the triangulation. For a high accuracy the triangulation
* should not contain too many non-acute triangles.
*
......
package org.vadere.simulator.models.potential.solver.calculators.mesh;
package org.vadere.simulator.models.potential.solver.calculators.mesh.examples;
import org.vadere.meshing.examples.MeshExamples;
import org.vadere.meshing.mesh.gen.PFace;
......@@ -12,6 +12,8 @@ import org.vadere.meshing.mesh.triangulation.improver.eikmesh.impl.PEikMesh;
import org.vadere.meshing.utils.io.poly.PSLGGenerator;
import org.vadere.meshing.utils.io.tex.TexGraphGenerator;
import org.vadere.simulator.models.potential.solver.calculators.EikonalSolver;
import org.vadere.simulator.models.potential.solver.calculators.mesh.MeshEikonalSolverFMM;
import org.vadere.simulator.models.potential.solver.calculators.mesh.PotentialPoint;
import org.vadere.simulator.models.potential.solver.timecost.UnitTimeCostFunction;
import org.vadere.util.geometry.GeometryUtils;
import org.vadere.util.geometry.shapes.VPoint;
......@@ -27,22 +29,38 @@ import java.util.Arrays;
import java.util.Collection;
import java.util.function.Function;
/**
* Just an example of how to use the FMM on a triangular mesh using EikMesh.
*
* The following tasks are executed:
* (1) load the planar straight-line graph from a file (the definition of the spatial domain)
* (2) generate a mesh using EikMesh (if {@link FMMTriangulationExamples#visualize} is true this is displayed)
* (3) define a target area (polygon)
* (4) apply the FMM on a triangular mesh
*/
public class FMMTriangulationExamples {
private static Logger log = Logger.getLogger(FMMTriangulationExamples.class);
private static boolean visualize = true;
private static boolean systemprint = true;
public static void main(String... args) throws IOException, InterruptedException {
// (1) + (2)
var triangulation = bridge();
VPolygon targetShape = GeometryUtils.toPolygon(new VPoint(33.90000000002328, 0.20000000018626451),
new VPoint(29.900000000023283, 0.5),
new VPoint(32.300000000046566, 6.0),
new VPoint(36.40000000002328, 4.900000000372529));
EikonalSolver solver = new MeshEikonalSolverFMM(Arrays.asList(targetShape), new UnitTimeCostFunction(), triangulation);
MeshEikonalSolverFMM solver = new MeshEikonalSolverFMM(Arrays.asList(targetShape), new UnitTimeCostFunction(), triangulation);
log.info("start FFM");
solver.solve();
log.info("FFM finished");
System.out.println(triangulation.getMesh().toPythonTriangulation(p -> triangulation.getMesh().getDoubleData(p, MeshEikonalSolverFMM.namePotential)));
if(systemprint) {
System.out.println(triangulation.getMesh().toPythonTriangulation(vertex -> solver.getPotential(vertex)));
}
}
public static IIncrementalTriangulation<
......@@ -69,6 +87,8 @@ public class FMMTriangulationExamples {
);
var mesh = meshImprover.getMesh();
Color green = new Color(85, 168, 104);
Color red = new Color(196,78,82);
Color blue = new Color(76,114,202);
......@@ -84,18 +104,26 @@ public class FMMTriangulationExamples {
}
};
var meshPanel = new PMeshPanel(meshImprover.getMesh(), 1000, 800, colorFunction);
meshPanel.display("Combined distance functions " + h0);
var meshPanel = visualize ? new PMeshPanel(meshImprover.getMesh(), 1000, 800, colorFunction) : null;
if(visualize) {
meshPanel.display("Combined distance functions " + h0);
}
while (!meshImprover.isFinished()) {
meshImprover.improve();
Thread.sleep(20);
meshPanel.repaint();
if(visualize) {
meshPanel.repaint();
}
}
//meshImprover.generate();
// display the mesh
meshPanel.display("Combined distance functions " + h0);
System.out.println(TexGraphGenerator.toTikz(mesh, colorFunction, null, 1.0f, true));
if(systemprint) {
System.out.println(TexGraphGenerator.toTikz(mesh, colorFunction, null, 1.0f, true));
}
return meshImprover.getTriangulation();
}
}
package org.vadere.simulator.models.potential.solver.timecost;
import org.jetbrains.annotations.NotNull;
import org.vadere.meshing.mesh.inter.IVertex;
import org.vadere.util.geometry.shapes.IPoint;
/**
* @author Benedikt Zoennchen
......@@ -13,4 +15,28 @@ public interface ITimeCostFunctionMesh<V extends IVertex> extends ITimeCostFunct
default double costAt(V p, Object caller) {
return costAt(p);
}
static <V extends IVertex> ITimeCostFunctionMesh convert(@NotNull final ITimeCostFunction timeCostFunction) {
return new ITimeCostFunctionMesh<V>() {
@Override
public double costAt(V p) {
return timeCostFunction.costAt(p);
}
@Override
public double costAt(V p, Object obj) {
return timeCostFunction.costAt(p, obj);
}
@Override
public double costAt(IPoint p) {
return timeCostFunction.costAt(p);
}
public double costAt(IPoint p, Object obj) {
return timeCostFunction.costAt(p, obj);
}
};
}
}
......@@ -27,7 +27,12 @@ import java.util.Set;
import java.util.function.Predicate;
/**
* A TimeCostFunction which reduces the travelling speed decreases with the obstacle density.
* A TimeCostFunction which reduces the travelling speed decreases with the pedestrian density.
* At the moment the code to dynamically refine the mesh is disabled
* (compare the unused method {@link TimeCostGaussianPedestrianDensityMesh#refineMesh}).
*
* Enabling the dynamic refinement is experimental code! It is not clear how often one should refine
* and how large the refinement cutoff should be used.
*
* @param <V>
* @param <E>
......@@ -46,6 +51,8 @@ public class TimeCostGaussianPedestrianDensityMesh<V extends IVertex, E extends
private GenRegularRefinement<V, E, F> refiner;
private final double R = 0.7;
// the radius in meter for which a pedestrian can contribute to the density, i.e. the cutoff radius.
private final int influenceRadius = 5;
private final double a;
private final double Sp;
......@@ -115,7 +122,8 @@ public class TimeCostGaussianPedestrianDensityMesh<V extends IVertex, E extends
VTriangle triangle = triangulation.getMesh().toTriangle(triangulation.getMesh().getFace(e));
if(/*!refiner.isGreen(e) || */triangulation.getMesh().toLine(e).length() > h_max) {
for(Pedestrian pedestrian : topography.getPedestrianDynamicElements().getElements()) {
if(pedestrian.getPosition().distanceSq(triangle.midPoint()) < influenceRadius * influenceRadius) {
// influenceRadius is cutoff radius for the mesh refinement is 4 times the density c