Commit a138df4b authored by Matthias Betz's avatar Matthias Betz
Browse files

changing math api to be more robust

parent f74c2f52
Pipeline #4354 passed with stage
in 2 minutes and 21 seconds
......@@ -22,6 +22,9 @@ import java.util.ArrayList;
import java.util.List;
import de.hft.stuttgart.citydoctor2.check.Check;
import de.hft.stuttgart.citydoctor2.math.ProjectionAxis;
import de.hft.stuttgart.citydoctor2.math.Ring2d;
import de.hft.stuttgart.citydoctor2.math.UnitVector3d;
import de.hft.stuttgart.citydoctor2.math.Vector2d;
import de.hft.stuttgart.citydoctor2.math.Vector3d;
......@@ -57,37 +60,13 @@ public class LinearRing extends GmlElement {
*/
public boolean isPointInside(Vector3d v) {
// project to 2d ring
List<Vector2d> projectedRing = new ArrayList<>();
Vector2d point;
Vector3d normal = calculateNormalNormalized();
double x = Math.abs(normal.getX());
double y = Math.abs(normal.getY());
double z = Math.abs(normal.getZ());
if (x > y && x > z) {
for (Vertex vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getY(), vert.getZ());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getY(), v.getZ());
} else if (y > x && y > z) {
for (Vertex vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getX(), vert.getZ());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getX(), v.getZ());
} else {
for (Vertex vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getX(), vert.getY());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getX(), v.getY());
}
ProjectionAxis axis = ProjectionAxis.of(this);
Vector2d point = axis.project(v);
Ring2d ring = Ring2d.of(this, axis);
int t = -1;
for (int i = 0; i < projectedRing.size() - 1; i++) {
t = t * crossProdTest(point, projectedRing.get(i), projectedRing.get(i + 1));
for (int i = 0; i < ring.getVertices().size() - 1; i++) {
t = t * crossProdTest(point, ring.getVertices().get(i), ring.getVertices().get(i + 1));
if (t == 0) {
return true;
}
......@@ -135,31 +114,8 @@ public class LinearRing extends GmlElement {
*
* @return the normal as a normalized vector
*/
public Vector3d calculateNormalNormalized() {
double[] coords = new double[3];
for (int i = 0; i < vertices.size() - 1; i++) {
Vertex current = vertices.get(i + 0);
Vertex next = vertices.get(i + 1);
coords[0] += (current.getZ() + next.getZ()) * (current.getY() - next.getY());
coords[1] += (current.getX() + next.getX()) * (current.getZ() - next.getZ());
coords[2] += (current.getY() + next.getY()) * (current.getX() - next.getX());
}
if (coords[0] == 0 && coords[1] == 0 && coords[2] == 0) {
// no valid normal vector found
if (vertices.size() < 3) {
// no three points, return x-axis
return new Vector3d(1, 0, 0);
}
Vertex v1 = vertices.get(0);
Vertex v2 = vertices.get(1);
Vertex v3 = vertices.get(2);
return calculateNormalWithCross(v1, v2, v3);
}
Vector3d v = new Vector3d(coords);
v.normalize();
return v;
public UnitVector3d calculateNormalNormalized() {
return calculateNormal().normalize();
}
/**
......@@ -183,7 +139,7 @@ public class LinearRing extends GmlElement {
// no valid normal vector found
if (vertices.size() < 3) {
// no three points, return x-axis
return new Vector3d(1, 0, 0);
return UnitVector3d.X_AXIS;
}
Vertex v1 = vertices.get(0);
......@@ -194,7 +150,7 @@ public class LinearRing extends GmlElement {
return new Vector3d(coords);
}
private Vector3d calculateNormalWithCross(Vertex v1, Vertex v2, Vertex v3) {
private UnitVector3d calculateNormalWithCross(Vertex v1, Vertex v2, Vertex v3) {
Vector3d dir1 = v2.minus(v1);
Vector3d dir2 = v3.minus(v1);
Vector3d cross = dir1.cross(dir2);
......
......@@ -33,22 +33,18 @@ public class CovarianceMatrix {
}
/**
* Calculates the covariance matrix of the given points, with the given expected
* values.
* see {@link CovarianceMatrix#calculateCovarianceMatrix(List, double[])}
*
* @param vertices the vertices for which the matrix is calculated
* @param expected the expected values
* @param expected the expected values as a vector
* @return the covariance 3 x 3 matrix
*/
public static double[][] calculateCovarianceMatrix(List<? extends Vector3d> vertices, double[] expected) {
if (expected.length != 3) {
throw new IllegalArgumentException("for 3D points 3 expected values have to be provided");
}
public static double[][] calculateCovarianceMatrix(List<? extends Vector3d> vertices, Vector3d expected) {
double[][] covValues = new double[3][3];
for (Vector3d v : vertices) {
double xdiff = v.getX() - expected[0];
double ydiff = v.getY() - expected[1];
double zdiff = v.getZ() - expected[2];
double xdiff = v.getX() - expected.getX();
double ydiff = v.getY() - expected.getY();
double zdiff = v.getZ() - expected.getZ();
covValues[0][0] += xdiff * xdiff;
covValues[0][1] += xdiff * ydiff;
covValues[0][2] += xdiff * zdiff;
......@@ -69,17 +65,6 @@ public class CovarianceMatrix {
return covValues;
}
/**
* see {@link CovarianceMatrix#calculateCovarianceMatrix(List, double[])}
*
* @param vertices the vertices for which the matrix is calculated
* @param expected the expected values as a vector
* @return the covariance 3 x 3 matrix
*/
public static double[][] calculateCovarianceMatrix(List<? extends Vector3d> vertices, Vector3d expected) {
return calculateCovarianceMatrix(vertices, expected.getCoordinates());
}
public static Vector3d getCentroid(List<? extends Vector3d> vertices) {
if (vertices.isEmpty()) {
throw new IllegalArgumentException("Vertices may not be empty to calculate centroid");
......
......@@ -45,16 +45,16 @@ public class MovedRing {
}
return indRing;
}
private MovedRing() {
vertices = new ArrayList<>();
}
public LinearRing getOriginal() {
return original;
}
public MovedRing() {
vertices = new ArrayList<>();
}
public void addVertex(Vector3d v) {
private void addVertex(Vector3d v) {
vertices.add(v);
}
......@@ -71,75 +71,10 @@ public class MovedRing {
*/
public boolean isPointInside(Vector3d v) {
// project to 2d ring
List<Vector2d> projectedRing = new ArrayList<>();
Vector2d point;
Vector3d normal = calculateNormal();
double x = Math.abs(normal.getX());
double y = Math.abs(normal.getY());
double z = Math.abs(normal.getZ());
if (x >= y && x >= z) {
for (Vector3d vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getY(), vert.getZ());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getY(), v.getZ());
} else if (y >= x && y >= z) {
for (Vector3d vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getX(), vert.getZ());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getX(), v.getZ());
} else {
for (Vector3d vert : vertices) {
Vector2d projCoords = new Vector2d(vert.getX(), vert.getY());
projectedRing.add(projCoords);
}
point = new Vector2d(v.getX(), v.getY());
}
int t = -1;
for (int i = 0; i < projectedRing.size() - 1; i++) {
t = t * crossProdTest(point, projectedRing.get(i), projectedRing.get(i + 1));
if (t == 0) {
return true;
}
}
return t >= 0;
}
private int crossProdTest(Vector2d a, Vector2d b, Vector2d c) {
if (a.getY() == b.getY() && a.getY() == c.getY()) {
if ((b.getX() <= a.getX() && a.getX() <= c.getX()) || (c.getX() <= a.getX() && a.getX() <= b.getX())) {
return 0;
} else {
return 1;
}
}
if (a.getY() == b.getY() && a.getX() == b.getX()) {
return 0;
}
if (b.getY() > c.getY()) {
Vector2d temp = b;
b = c;
c = temp;
}
if (a.getY() <= b.getY() || a.getY() > c.getY()) {
return 1;
}
return calculateDelta(a, b, c);
}
private int calculateDelta(Vector2d a, Vector2d b, Vector2d c) {
double delta = (b.getX() - a.getX()) * (c.getY() - a.getY()) - (b.getY() - a.getY()) * (c.getX() - a.getX());
if (delta > 0) {
return -1;
} else if (delta < 0) {
return 1;
} else {
return 0;
}
ProjectionAxis axis = ProjectionAxis.of(this);
Ring2d projectedRing = Ring2d.of(this, axis);
Vector2d point = axis.project(v);
return projectedRing.containsPoint(point);
}
/**
......@@ -149,37 +84,7 @@ public class MovedRing {
*
* @return the normal as a normalized vector
*/
public Vector3d calculateNormal() {
double[] coords = new double[3];
for (int i = 0; i < vertices.size() - 1; i++) {
Vector3d current = vertices.get(i + 0);
Vector3d next = vertices.get(i + 1);
coords[0] += (current.getZ() + next.getZ()) * (current.getY() - next.getY());
coords[1] += (current.getX() + next.getX()) * (current.getZ() - next.getZ());
coords[2] += (current.getY() + next.getY()) * (current.getX() - next.getX());
}
if (coords[0] == 0 && coords[1] == 0 && coords[2] == 0) {
// no valid normal vector found
if (vertices.size() < 3) {
// no three points, return x-axis
return new Vector3d(1, 0, 0);
}
Vector3d v1 = vertices.get(0);
Vector3d v2 = vertices.get(1);
Vector3d v3 = vertices.get(2);
return calculateNormalWithCross(v1, v2, v3);
}
Vector3d v = new Vector3d(coords);
v.normalize();
return v;
}
private Vector3d calculateNormalWithCross(Vector3d v1, Vector3d v2, Vector3d v3) {
Vector3d dir1 = v2.minus(v1);
Vector3d dir2 = v3.minus(v1);
Vector3d cross = dir1.cross(dir2);
return cross.normalize();
public UnitVector3d calculateNormalNormalized() {
return original.calculateNormalNormalized();
}
}
......@@ -39,7 +39,7 @@ public class Plane implements Serializable {
private static final double EPSILON = 0.0001;
private Vector3d normal;
private UnitVector3d normal;
private Vector3d point;
private double d;
......@@ -51,11 +51,10 @@ public class Plane implements Serializable {
*/
public Plane(Vector3d normal, Vector3d p) {
point = p;
this.normal = normal.copy();
this.normal.normalize();
this.normal = normal.normalize();
d = this.normal.dot(point);
if (d < 0) {
this.normal = this.normal.mult(-1);
this.normal = this.normal.invert();
d = -d;
}
}
......@@ -138,7 +137,7 @@ public class Plane implements Serializable {
*
* @return the normal vector of the plane
*/
public Vector3d getNormal() {
public UnitVector3d getNormal() {
return normal;
}
......
......@@ -24,7 +24,6 @@ import java.util.List;
import de.hft.stuttgart.citydoctor2.datastructure.LinearRing;
import de.hft.stuttgart.citydoctor2.datastructure.Polygon;
import de.hft.stuttgart.citydoctor2.datastructure.Vertex;
/**
* A two dimensional polygon. Has an 2d exterior ring and interior rings
......@@ -38,81 +37,36 @@ public class Polygon2d {
private List<Ring2d> innerRings;
public static Polygon2d withProjection(Polygon poly) {
Vector3d normal = poly.calculateNormalNormalized();
int[] axis = getProjectionAxis(normal);
ProjectionAxis axis = ProjectionAxis.of(poly);
return projectTo2D(poly, axis);
}
public static Polygon2d withProjection(MovedPolygon poly, int[] projectionAxis) {
public static Polygon2d withProjection(MovedPolygon poly, ProjectionAxis projectionAxis) {
return projectTo2D(poly, projectionAxis);
}
public static Polygon2d withProjection(Polygon poly, int[] projectionAxis) {
public static Polygon2d withProjection(Polygon poly, ProjectionAxis projectionAxis) {
return projectTo2D(poly, projectionAxis);
}
private static Polygon2d projectTo2D(Polygon p, int[] axis) {
private static Polygon2d projectTo2D(Polygon p, ProjectionAxis axis) {
List<Ring2d> interior = new ArrayList<>();
Ring2d exterior = projectRing(p.getExteriorRing(), axis);
Ring2d exterior = Ring2d.of(p.getExteriorRing(), axis);
for (LinearRing innerRing : p.getInnerRings()) {
interior.add(projectRing(innerRing, axis));
interior.add(Ring2d.of(innerRing, axis));
}
return new Polygon2d(exterior, interior);
}
private static Polygon2d projectTo2D(MovedPolygon p, int[] axis) {
private static Polygon2d projectTo2D(MovedPolygon p, ProjectionAxis axis) {
List<Ring2d> interior = new ArrayList<>();
Ring2d exterior = projectRing(p.getExteriorRing(), axis);
Ring2d exterior = Ring2d.of(p.getExteriorRing(), axis);
for (MovedRing innerRing : p.getInnerRings()) {
interior.add(projectRing(innerRing, axis));
interior.add(Ring2d.of(innerRing, axis));
}
return new Polygon2d(exterior, interior);
}
public static int[] getProjectionAxis(Vector3d normal) {
double nx = Math.abs(normal.getX());
double ny = Math.abs(normal.getY());
double nz = Math.abs(normal.getZ());
int[] axis;
if (nx >= ny && nx >= nz) {
axis = new int[] { 1, 2 };
} else if (ny >= nx && ny >= nz) {
axis = new int[] { 0, 2 };
} else {
axis = new int[] { 0, 1 };
}
return axis;
}
private static Ring2d projectRing(LinearRing r, int[] axis) {
List<Vector2d> projectedVertices = new ArrayList<>();
for (Vertex v : r.getVertices()) {
projectedVertices.add(new Vector2d(v.getCoordinate(axis[0]), v.getCoordinate(axis[1])));
}
return new Ring2d(projectedVertices, r);
}
private static Ring2d projectRing(MovedRing r, int[] axis) {
List<Vector2d> projectedVertices = new ArrayList<>();
for (Vector3d v : r.getVertices()) {
projectedVertices.add(new Vector2d(v.getCoordinate(axis[0]), v.getCoordinate(axis[1])));
}
return new Ring2d(projectedVertices, r.getOriginal());
}
public static Polygon2d withRotationMatrix(Polygon poly) {
Matrix3x3d rotMatrix = PolygonUtils.calculateRotationMatrix(poly);
LinearRing lr = poly.getExteriorRing();
Ring2d exterior = createRing2d(rotMatrix, lr);
List<Ring2d> innerRings = new ArrayList<>();
for (LinearRing innerRing : poly.getInnerRings()) {
innerRings.add(createRing2d(rotMatrix, innerRing));
}
return new Polygon2d(exterior, innerRings);
}
/**
* Converts a Polygon to a 2d polygon. This will change the area of the polygon.
*
......@@ -126,16 +80,6 @@ public class Polygon2d {
this.innerRings = interior;
}
private static Ring2d createRing2d(Matrix3x3d rotMatrix, LinearRing lr) {
List<Vector2d> ringVertices = new ArrayList<>();
for (Vertex v : lr.getVertices()) {
Vector3d rotated = rotMatrix.mult(v);
Vector2d v2d = new Vector2d(rotated.getX(), rotated.getY());
ringVertices.add(v2d);
}
return new Ring2d(ringVertices, lr);
}
public Ring2d getExterior() {
return exterior;
}
......
/*-
* Copyright 2020 Beuth Hochschule für Technik Berlin, Hochschule für Technik Stuttgart
*
* This file is part of CityDoctor2.
*
* CityDoctor2 is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* CityDoctor2 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with CityDoctor2. If not, see <https://www.gnu.org/licenses/>.
*/
package de.hft.stuttgart.citydoctor2.math;
import java.util.Arrays;
import de.hft.stuttgart.citydoctor2.datastructure.LinearRing;
import de.hft.stuttgart.citydoctor2.datastructure.Polygon;
import de.hft.stuttgart.citydoctor2.datastructure.Vertex;
/**
* This class is needed for projecting rings or polygons to a 2D axis plane. It
* uses the longest component of the normal vector to determine which plane is
* the best plane to project the ring or polygon on while conserving the most
* area of the original polygon.
*
* @author Matthias Betz
*
*/
public class ProjectionAxis {
private static final String DIVISOR_IS_0 = "Divisor is 0";
private int[] axis;
public static ProjectionAxis of(Polygon p) {
return getProjectionAxis(p.calculateNormal());
}
public static ProjectionAxis of(MovedRing movedRing) {
// the normal should be the same even if the ring has moved
return of(movedRing.getOriginal());
}
public static ProjectionAxis of(LinearRing lr) {
return getProjectionAxis(lr.calculateNormal());
}
public static ProjectionAxis of(Plane plane) {
return getProjectionAxis(plane.getNormal());
}
/**
* Calculates the indices of the two lowest coordinate values of this vector.
* This is used to determine on which axis plane a ring or polygon should be
* projected.
*
* @param normal the normal vector, does not need to be normalized
* @return an array with the indices of the coordinate array which should be
* used in a 2d environment.
*/
private static ProjectionAxis getProjectionAxis(Vector3d normal) {
double nx = Math.abs(normal.getX());
double ny = Math.abs(normal.getY());
double nz = Math.abs(normal.getZ());
double biggestValue = nx;
int[] axis = new int[] { 1, 2 };
if (ny > biggestValue) {
biggestValue = ny;
axis[0] = 0;
axis[1] = 2;
}
if (nz > biggestValue) {
axis[0] = 0;
axis[1] = 1;
}
return new ProjectionAxis(axis);
}
private ProjectionAxis(int[] axis) {
this.axis = axis;
}
public Vector2d project(Vector3d v) {
return new Vector2d(v.getCoordinate(axis[0]), v.getCoordinate(axis[1]));
}
/**
* calculates the missing coordinate for 3d vector from the plane and this axis.
* @return the projected 3d point.
*/
public Vector3d projectToPlane(Plane plane, Vector2d v) {
return projectToPlane(plane, v.getX(), v.getY());
}
/**
* calculates the missing coordinate for 3d vector from the plane and this axis.
* @return the projected 3d point.
*/
public Vector3d projectToPlane(Plane plane, double vectorX, double vectorY) {
UnitVector3d normal = plane.getNormal();
double a = normal.getX();
double b = normal.getY();
double c = normal.getZ();
double d = plane.getD();
if (axis[0] == 0 && axis[1] == 1) {
double x = vectorX;
double y = vectorY;
if (c == 0) {
throw new IllegalStateException(DIVISOR_IS_0);
}
double z = (d - a * x - b * y) / c;
return new Vertex(x, y, z);
} else if (axis[0] == 0 && axis[1] == 2) {
double x = vectorX;
double z = vectorY;
if (b == 0) {
throw new IllegalStateException(DIVISOR_IS_0);
}
double y = (d - a * x - c * z) / b;
return new Vertex(x, y, z);
} else if (axis[0] == 1 && axis[1] == 2) {
double y = vectorX;
double z = vectorY;
if (a == 0) {
throw new IllegalStateException(DIVISOR_IS_0);
}
double x = (d - b * y - c * z) / a;
return new Vertex(x, y, z);
} else {
throw new IllegalStateException("Unknown axis: " + axis);
}
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
result = prime * result + Arrays.hashCode(axis);
return result;
}
@Override
public boolean equals(Object obj) {
if (this == obj)
return true;
if (obj == null)
return false;
if (getClass() != obj.getClass())
return false;
ProjectionAxis other = (ProjectionAxis) obj;
return Arrays.equals(axis, other.axis);
}
@Override
public String toString() {
return "ProjectionAxis [axis=" + Arrays.toString(axis) + "]";
}
}
......@@ -34,7 +34,7 @@ public class Ray {
public static final double EPSILON = 0.0001;
private Vector3d origin;
private Vector3d direction;
private UnitVector3d direction;
public Ray(Vector3d origin, Vector3d direction) {
this.origin = origin;
......
......@@ -18,9 +18,11 @@
*/
package de.hft.stuttgart.citydoctor2.math;
import java.util.ArrayList;
import java.util.List;
import de.hft.stuttgart.citydoctor2.datastructure.LinearRing;
import de.hft.stuttgart.citydoctor2.datastructure.Vertex;
/**
* A two dimensional ring containing the 2d vertices
......@@ -32,6 +34,27 @@ public class Ring2d {
private List<Vector2d> ringVertices;
private LinearRing original;
public static Ring2d of(MovedRing movedRing) {
return of(movedRing, ProjectionAxis.of(movedRing));
}
public static Ring2d of(LinearRing ring, ProjectionAxis axis) {
List<Vector2d> projectedVertices = new ArrayList<>();
for (Vertex v : ring.getVertices()) {
projectedVertices.add(axis.project(v));
}
return new Ring2d(projectedVertices, ring);
}
public static Ring2d of(MovedRing ring, ProjectionAxis axis) {
List<Vector2d> projectedVertices = new ArrayList<>();
for (Vector3d v : ring.getVertices()) {
projectedVertices.add(axis.project(v));
}
return new Ring2d(projectedVertices, ring.getOriginal());
}
public Ring2d(List<Vector2d> ringVertices, LinearRing original) {
this.ringVertices = ringVertices;
......
......@@ -18,25 +18,83 @@
*/
package de.hft.stuttgart.citydoctor2.math;
import de.hft.stuttgart.citydoctor2.datastructure.Polygon;
import java.util.Arrays;
/**
* Utility class for calculation concerning polygons
*
* A unit vector3d which always has the length 1.
* @author Matthias Betz
*
*/
public class PolygonUtils {
public class UnitVector3d extends Vector3d {
private static final String UNIT_VECTOR_IS_IMMUTABLE = "Unit vector is immutable";
private PolygonUtils() {
// only static utility
private static final long serialVersionUID = -374685263673211587L;
public static final UnitVector3d X_AXIS = new UnitVector3d(1, 0, 0);
public static final UnitVector3d Y_AXIS = new UnitVector3d(0, 1, 0);
public static final UnitVector3d Z_AXIS = new UnitVector3d(0, 0, 1);
public static UnitVector3d of(Vector3d v) {
double length = v.getLength();
double x = v.getX() / length;
double y = v.getY() / length;
double z = v.getZ() / length;
return new UnitVector3d(x, y, z);
}
private UnitVector3d(double x, double y, double z) {
super(x, y, z);
}
@Override
public double getLength() {
return 1;
}
@Override
public double getSquaredLength() {
return 1;
}
@Override
public UnitVector3d normalize() {
return this;
}
public static Matrix3x3d calculateRotationMatrix(Polygon p) {
Vector3d normal = p.calculateNormalNormalized();
Vector3d zAxis = new Vector3d(0d, 0d, 1d);
Quaternion rot = Quaternion.fromToRotation(normal, zAxis);
return rot.toRotationMatrix();
/**
* creates a new unit vector in the other direction, still with length 1
* @return the normal vector in the opposite direction
*/
public UnitVector3d invert() {
return new UnitVector3d(-getX(), -getY(), -getZ());
}
@Override
public UnitVector3d copy() {
return new UnitVector3d(getX(), getY(), getZ());
}
@Override
public double[] getCoordinates() {
return Arrays.copyOf(super.getCoordinates(), 3);
}
@Override
public void setX(double x) {
throw new UnsupportedOperationException(UNIT_VECTOR_IS_IMMUTABLE);
}
@Override
public void setY(double y) {
throw new UnsupportedOperationException(UNIT_VECTOR_IS_IMMUTABLE);
}
@Override
public void setZ(double z) {
throw new UnsupportedOperationException(UNIT_VECTOR_IS_IMMUTABLE);