Feat: Add connection angle calculation

CityDoctorParent/CityDoctorModel/src/main/java/de/hft/stuttgart/citydoctor2/math/Triangle3d.java

@@ -20,6 +20,8 @@ package de.hft.stuttgart.citydoctor2.math;
 
 import java.io.Serial;
 import java.io.Serializable;
+import java.util.ArrayList;
+import java.util.List;
 
 import de.hft.stuttgart.citydoctor2.tesselation.TesselatedPolygon;
 
@@ -208,6 +210,32 @@ public class Triangle3d implements Serializable {
 		return intersection1 != null;
 	}
 
+	/**
+	 * Computes the points where this triangle's edges cross the other triangle and vice versa.
+	 * For two generically intersecting (non-coplanar) triangles this yields the two endpoints
+	 * of their intersection segment.
+	 *
+	 * @param other the other triangle
+	 * @param eps   epsilon used for the underlying edge/triangle intersection test
+	 * @return the distinct intersection points found, usually 0 (no intersection) or 2
+	 */
+	public List<Vector3d> intersectionPoints(Triangle3d other, double eps) {
+		List<Vector3d> points = new ArrayList<>();
+		addDistinct(points, new Segment3d(p1, p2).intersection(other, eps));
+		addDistinct(points, new Segment3d(p1, p3).intersection(other, eps));
+		addDistinct(points, new Segment3d(p2, p3).intersection(other, eps));
+		addDistinct(points, new Segment3d(other.p1, other.p2).intersection(this, eps));
+		addDistinct(points, new Segment3d(other.p1, other.p3).intersection(this, eps));
+		addDistinct(points, new Segment3d(other.p2, other.p3).intersection(this, eps));
+		return points;
+	}
+
+	private static void addDistinct(List<Vector3d> points, Vector3d point) {
+		if (point != null && !points.contains(point)) {
+			points.add(point);
+		}
+	}
+
 	public boolean containsPoint(Vector3d point) {
 		Vector3d a = p1;
 		Vector3d b = p2;

CityDoctorParent/CityDoctorValidation/src/main/java/de/hft/stuttgart/citydoctor2/checks/topology/MisalignedUtilityConnectionCheck.java

@@ -30,6 +30,8 @@ import de.hft.stuttgart.citydoctor2.datastructure.Polygon;
 import de.hft.stuttgart.citydoctor2.datastructure.Vertex;
 import de.hft.stuttgart.citydoctor2.math.OrthogonalRegressionPlane;
 import de.hft.stuttgart.citydoctor2.math.Plane;
+import de.hft.stuttgart.citydoctor2.math.Ray;
+import de.hft.stuttgart.citydoctor2.math.Triangle3d;
 import de.hft.stuttgart.citydoctor2.math.Vector3d;
 import de.hft.stuttgart.citydoctor2.parser.ParserConfiguration;
 import de.hft.stuttgart.citydoctor2.utils.Pair;
@@ -165,6 +167,11 @@ public class MisalignedUtilityConnectionCheck extends Check {
                             .toList();
                     List<Pair<Vertex, Vertex>> pairs = getVertexPairs(sGeo, cGeo);
                     ConnectionAlignmentType alignment = analyseVertexPairConnectionAlignment(sGeo, cGeo, pairs, intersections);
+                    if (alignment == ConnectionAlignmentType.UNKNOWN) {
+                        alignment = !intersections.isEmpty()
+                                ? analyseIntersectionConnectionAlignment(intersections)
+                                : analyseContainmentAlignment(sGeo, cGeo);
+                    }
                     switch(alignment) {
                         case ALIGNED:
                             connected = true;
@@ -183,46 +190,19 @@ public class MisalignedUtilityConnectionCheck extends Check {
                             connected = true;
                             misalignment = true;
                             break;
+                        case DISCONNECTED:
+                            // sGeo does not connect to this candidate's geometry at all; handled
+                            // globally by the freeFloating check once every candidate was tried
+                            break;
                         case UNKNOWN:
+                            subject.addCheckResult(new CheckResult(this,
+                                    ResultStatus.WARNING, new MisalignedConnectionError(sGeo, subject, candidate.getGmlId(),
+                                    MisalignmentType.UNKNOWN, ErrorType.WARNING)));
+                            misalignment = true;
                             break;
                         default:
                             throw new IllegalStateException("Vertex pairs connection analysis returned illegal alignment type");
                     }
-                    if(alignment == ConnectionAlignmentType.UNKNOWN && !intersections.isEmpty()){
-
-                        List<Pair<Double, Double>> dihedralAngles = new ArrayList<>();
-                        double intersectionLengthMax = 0d;
-                        for(PolygonIntersection intersection : intersections){
-
-                           Polygon polyA = intersection.getP1();
-                           Polygon polyB = intersection.getP2();
-                           Vector3d normalA = polyA.calculateNormal();
-                           Vector3d normalB = polyB.calculateNormal();
-                           // TODO: Möller-Trumbore intersection of polyA and polyB
-                           double intersectionLength = 1d;
-                           if (intersectionLengthMax < intersectionLength) {
-                               intersectionLengthMax = intersectionLength;
-                           }
-                           double thetaN = Math.atan2(normalA.cross(normalB).getLength(), normalA.dot(normalB));
-                           dihedralAngles.add(new Pair<>((180d-thetaN),intersectionLength));
-                        }
-                        double finalIntersectionLengthMax = intersectionLengthMax;
-                        double weightedSum = dihedralAngles.stream()
-                                .mapToDouble(
-                                        pair->
-                                                (pair.getValue0() * pair.getValue1()) / finalIntersectionLengthMax)
-                                .sum();
-                        double mean = weightedSum / dihedralAngles.size();
-                        double std = 0.0;
-                        for (Pair<Double, Double> pair : dihedralAngles) {
-                            std += Math.pow(pair.getValue0() - mean, 2);
-                        }
-                    } else if (alignment == ConnectionAlignmentType.UNKNOWN) {
-                        //TODO: Check if sGeo is contained within the interior volumetric space of  cGeo
-
-                    }
-
-
                 }
                 if (misalignment) {
                     misalignedConnections = true;
@@ -307,11 +287,68 @@ public class MisalignedUtilityConnectionCheck extends Check {
         return centroidsAligned ? ConnectionAlignmentType.ALIGNED : ConnectionAlignmentType.SHEARED;
     }
 
-    private ConnectionAlignmentType analyseIntersectionConnectionAlignment(Geometry subject, Geometry candidate,
-                                                                           List<PolygonIntersection> intersections) {
+    private ConnectionAlignmentType analyseIntersectionConnectionAlignment(List<PolygonIntersection> intersections) {
+        List<Pair<Double, Double>> dihedralAngles = new ArrayList<>();
+        double totalLength = 0d;
+        for (PolygonIntersection intersection : intersections) {
+            List<Vector3d> points = intersection.getT1().intersectionPoints(intersection.getT2(), epsilon);
+            double intersectionLength = maxPairwiseDistance(points);
+            if (intersectionLength <= 0d) {
+                continue;
+            }
+            Vector3d normalA = intersection.getP1().calculateNormal();
+            Vector3d normalB = intersection.getP2().calculateNormal();
+            double thetaN = Math.toDegrees(Math.atan2(normalA.cross(normalB).getLength(), normalA.dot(normalB)));
+            dihedralAngles.add(new Pair<>(180d - thetaN, intersectionLength));
+            totalLength += intersectionLength;
+        }
+        if (dihedralAngles.isEmpty()) {
+            return ConnectionAlignmentType.UNKNOWN;
+        }
+
+        double mean = dihedralAngles.stream()
+                .mapToDouble(pair -> pair.getValue0() * pair.getValue1())
+                .sum() / totalLength;
+
+        if (mean <= angleTolerance) {
+            // utilities are essentially in a line; only the fine tolerance applies
+            return ConnectionAlignmentType.ALIGNED;
+        }
+        if (mean < shallownessThreshold) {
+            return ConnectionAlignmentType.BENT;
+        }
+        // bend is too severe to confidently call this a misaligned-but-real connection
         return ConnectionAlignmentType.UNKNOWN;
     }
 
+    private double maxPairwiseDistance(List<Vector3d> points) {
+        double max = 0d;
+        for (int i = 0; i < points.size(); i++) {
+            for (int j = i + 1; j < points.size(); j++) {
+                max = Math.max(max, points.get(i).getDistance(points.get(j)));
+            }
+        }
+        return max;
+    }
+
+    private ConnectionAlignmentType analyseContainmentAlignment(Geometry subject, Geometry candidate) {
+        Vertex point = UtilityNetworkUtils.getCentroid(subject.getVertices());
+        return isPointInsideSolid(point, candidate) ? ConnectionAlignmentType.ALIGNED : ConnectionAlignmentType.DISCONNECTED;
+    }
+
+    private boolean isPointInsideSolid(Vertex point, Geometry solid) {
+        Ray ray = new Ray(point, new Vector3d(1, 1, 1));
+        int crossings = 0;
+        for (Polygon p : solid.getPolygons()) {
+            for (Triangle3d t : p.tesselate().getTriangles()) {
+                if (ray.intersectTriangle(t) != null) {
+                    crossings++;
+                }
+            }
+        }
+        return crossings % 2 == 1;
+    }
+
     private Plane deriveConnectionPlane(List<Vertex> sPairVerts, List<PolygonIntersection> intersections) {
         // Calculate intersection plane
         if (!intersections.isEmpty()) {

CityDoctorParent/CityDoctorValidation/src/main/java/de/hft/stuttgart/citydoctor2/checks/util/UtilityNetworkUtils.java

@@ -4,6 +4,7 @@ import de.hft.stuttgart.citydoctor2.datastructure.Edge;
 import de.hft.stuttgart.citydoctor2.datastructure.Geometry;
 import de.hft.stuttgart.citydoctor2.datastructure.Vertex;
 import de.hft.stuttgart.citydoctor2.math.Plane;
+import de.hft.stuttgart.citydoctor2.math.Vector3d;
 
 import java.util.ArrayList;
 import java.util.Collection;
@@ -58,9 +59,10 @@ public class UtilityNetworkUtils {
     }
 
     public static Vertex getCentroid(Collection<Vertex> vertices) {
-        Vertex centroid = new Vertex(0,0,0);
-        vertices.forEach(centroid::plus);
-        centroid.mult(((double) 1 /vertices.size()));
-        return centroid;
+        Vector3d sum = new Vector3d(0, 0, 0);
+        for (Vertex v : vertices) {
+            sum = sum.plus(v);
+        }
+        return new Vertex(sum.mult(1d / vertices.size()));
     }
 }