Rounding errors and side effects with straight skeleton algorithm

[micycle1]

A user writes...

I have hit an issue that has to do with the skeleton algorithm. The algorithm is being used with path2D (Double) as input and output. In between I use several layers of transformations including your port of the excellent clipper2 library.

I get also various error messages when trying to use triangulations but lets focus on the skeleton algorithm that gives in my particular case the very best, aesthetically, results.

I think I am being hit by rounding errors on double to float conversions because of PShape/PVector design.
This is apparent on almost all cases but it is more pronounced when using holes. There is a partial workarounds by using layering of clipper2 offsets with very small (close to zero)offsets, so usually the calculation allows some refinement of the rounding errors, or at least space can be filled with different "skeleton faces". I have tried some other ways with no results.

I use single and Multiple (totally internal) holes and holes that can touch or intersect main polygon, I create holes by using Area class operations.

PGS 1.3.0 version skeleton faces does not support multi-holes. It gives many artifacts on multi-hole paths and it breaks XOR operations too. Mixed results.
PGS 1.4.0 version is much better, XOR operations work, INTERSECTION IS ALSO OK on multi hole paths (as I use it).

Now about the SUBTRACTION, it generally seems to work BUT it often produces some randomly exposed artifacts, and I am not sure if it is the float/double issue or a bug/inconsistency in skeleton or jts conversion code. Those artifacts are not deterministic. But when I use the clipper offset transfor layer they become deterministic (but still unpredictably random). By deterministic I mean when recalculating the result is the same.

I have also noted that hole polygons that touch the base polygon usually (but not always) get the skeleton right. Only totally inner polygons can create intense artifacts, so a guess is that rounding errors get more pronounced, exaggerated on internal only holes.
Those artifacts are being revealed when iterating thousands of polygons, you might not have the chance to hit this bugs otherwise as it can be 1-10 per 1000 polygons.

I also often see errors like:

java.lang.RuntimeException: Matrix is singular.
at Jama.LUDecomposition.solve(LUDecomposition.java:282)
at Jama.Matrix.solve(Matrix.java:815)
at org.twak.utils.Jama.solve(Jama.java:34)
at org.twak.utils.geom.LinearForm3D.collide(LinearForm3D.java:195)
at org.twak.camp.CoSitedCollision.validateChains(CoSitedCollision.java:180)
at org.twak.camp.HeightCollision.process(HeightCollision.java:105)
at org.twak.camp.Skeleton.skeleton(Skeleton.java:211)
at micycle.pgs.PGS_Contour.straightSkeletonTwak(PGS_Contour.java:301)
at micycle.pgs.PGS_Contour.straightSkeleton(PGS_Contour.java:272)
at micycle.pgs.PGS_Contour.straightSkeleton(PGS_Contour.java:253)

This error is not appearing on every artifact, but it could be related, in fact it mostly appears on very high resolution (very small, near zero length) structures that possibly need better accuracy to be processed and rendered, and can be a side effect of the rounding errors.

I cannot at the moment create a reproduction code, and I am not sure if I can reproduce this outside my workflow.
I am not asking for a solution or course, I m just reporting my experience.
I have also created a feature request in processing4 GitHub but I do not think it will never be touched.
processing/processing4#893

I have though, created a very quick n brutal (only geometry) version of processing-core for doubles, fairly quickly (by replacing of all "float" to "double" words in the codebase, n then trying to fix some errors in various places in order for the codebase be able to compile), but then pgs uses float internally so it breaks with errors. That was just late night code overdose joke, but it could potentially give a valuable result if pgs codebase was also modified accordingly. I didn't try it, I cannot go further, I am not that good.

Assuming that this is a rounding error issue only (not 100% sure on all levels, but it seems to be!), what do you think are the options?
Maybe try to re-create the skeleton algorithm for Path2D i/o?
Or try to avoid the usage of PShape class completely and rely only on jts classes code in a new method in PGS codebase? With the appropriate to/fromJava2D converters. would that be possible?

If you have any ideas or any other/additional experimental skeleton code, please let me know...

[micycle1]

@istinnstudio

Some thoughts...

Or try to avoid the usage of PShape class completely and rely only on jts classes code in a new method in PGS codebase? With the appropriate to/fromJava2D converters. would that be possible?

It's fairly possible and I'd try that first to check whether it is a problem with precision. For instance, the internal straightSkeletonTwak() method takes in JTS Polygon types, whose coordinates are defined with doubles. You could do an initial conversion of Java2D --> JTS, then work entirely with JTS Geometries. Many of PGS methods operate on JTS Geometries under the hood so you could lift the methods from source and change the inputs from PShapes.

An alternative is the skeleton algorithm itself cannot handle the "very small, near zero length" geometries you are attempting to process, even at float-like precision (assuming they are valid geometries). Another alternative is that after subtraction, some polygon's widths get collapsed to zero, so a skeleton is impossible (this could make sense given offsetting seems to help a bit).

Also, you shouldn't need to use clipper2 alongside PGS/JTS. Instead of using clipper2 to offset, you can use geom.buffer() (this is what PGS_Morphology.buffer() uses under the hood).

And if you can create a minimal reproducible example, that would be handy -- I could poke around with things too.

[istinnstudio]

Thanks for the suggestions and the informative reply. It helps a lot.

[istinnstudio]

I have created a rough version of processing4-core that supports double coordinates and an equivalent version of pgs, The issue remains so it looks like is not a float to double rounding error at this level (if my double version is ok). It could rely somewhere in the conversions or skeleton code itself, I cannot go that deep. So for now I am not able to do more or create a reproducible example, I will keep this version and edit artifacts in post-processing. At least the result is more than acceptable!

to sum up, I am posting the 2 most common console errors I get in very high resolution vector artworks, those are not directly related 1-1 messages with the issue, but could be close enough to what it is going on:

using clipper2 offset (on Path2D level before and after skeleton) helps a lot to eliminate them, and redirect the amount and positioning of the faces artifacts

(this one is on the double coords version)
java.lang.IllegalArgumentException: Points of LinearRing do not form a closed linestring
	at org.locationtech.jts.geom.LinearRing.validateConstruction(LinearRing.java:91)
	at org.locationtech.jts.geom.LinearRing.<init>(LinearRing.java:86)
	at org.locationtech.jts.geom.GeometryFactory.createLinearRing(GeometryFactory.java:380)
	at org.locationtech.jts.geom.GeometryFactory.createLinearRing(GeometryFactory.java:367)
	at org.locationtech.jts.awt.ShapeReader.read(ShapeReader.java:98)
	at org.locationtech.jts.awt.ShapeReader.read(ShapeReader.java:59)
	at micycle.pgs.PGS_Conversion.fromJava2D(PGS_Conversion.java:1277)

java.lang.RuntimeException: Matrix is singular.
	at Jama.LUDecomposition.solve(LUDecomposition.java:282)
	at Jama.Matrix.solve(Matrix.java:815)
	at org.twak.utils.Jama.solve(Jama.java:34)
	at org.twak.utils.geom.LinearForm3D.collide(LinearForm3D.java:195)
	at org.twak.camp.CoSitedCollision.validateChains(CoSitedCollision.java:180)
	at org.twak.camp.HeightCollision.process(HeightCollision.java:105)
	at org.twak.camp.Skeleton.skeleton(Skeleton.java:211)
	at micycle.pgs.PGS_Contour.straightSkeletonTwak(PGS_Contour.java:299)
	at micycle.pgs.PGS_Contour.straightSkeleton(PGS_Contour.java:268)
	at micycle.pgs.PGS_Contour.straightSkeleton(PGS_Contour.java:248)

[micycle1]

The first error means your Java2D shapes don't form polygons (they are an unclosed string of vertices).

[istinnstudio]

Yes apparently it does and it is about the skeleton faces and not the source polygon. The only pattern I could recognize is that more complex polygons are more prone to artifacts. Polygons with holes are such a category. In my geometries there are always exactly same polygons that skeletonize right. That's why I initially thought that this is a rounding error.  For my experience there are 2 kind "triggers" of artifacts. One appears on skeletonized polygons without holes and usually can be fixed with small offset values. The other one is most apparent on multi-hole source polygons (so complex that can scramble the skeletonized faces). This one might has to do with the complexity of the source polygon and is not easily fixed, as a result and moves around on offset width change, and one other same initial polygon inherits the artifact!  For what is worth I have created a github repository with the experimental double precision versions of processing4 and pgs. I cannot say if is useful for you:   GitHub - istinnstudio/processing4-double: This is a rough modification of processing4 core in order to support double precision coordinates

[istinnstudio]

As far as I can see and test, both kendzi and twak straight skeleton algorithms have issues on my polygons. I do not know exactly what causes those issues, probably some extreme precision calculations fail to deliver exact results. In fact using kendzi and if it fails then switch to twak actually give better (but still with artifacts) results than using them standalone. As an alternative, outside of PGS scope, I am thinking of trying to access CGAL straight skeleton implementation using the Java 21 Foreign Function and Memory API. Java swig bindings for CGAL does not support straight skeleton at the moment.

[istinnstudio]

I have used this class: https://github.com/next-mad-hatter/straight-skeleton-study/blob/master/triton/src/main/java/util/CoordinatesScaler.java
Inside PGS code, transforming coordinates to those arithmetic formats increases the calculation precision in Twak algorithm. Then reverting back to double coordinates for rendering. This one eliminates artifacts to a controllable state, especially on "very small" polygons. It is far from perfect but it is an improvement, artifacts now are controllable. I could not do this for Kendzi, though, but now Twak is at least usable in my specific case. So it looks like precision is a factor but I also think that both algorithms suffer from internal degenerated results too.

[micycle1]

That's good to know. What kind of scaling factor are you using?

…

________________________________ From: istinnstudio ***@***.***> Sent: 28 March 2025 20:36 To: micycle1/PGS ***@***.***> Cc: Michael Carleton ***@***.***>; Author ***@***.***> Subject: Re: [micycle1/PGS] Rounding errors and side effects with straight skeleton algorithm (Discussion #133) I have used this class: https://github.com/next-mad-hatter/straight-skeleton-study/blob/master/triton/src/main/java/util/CoordinatesScaler.java Inside PGS code, transforming coordinates to those arithmetic formats increases the calculation precision in Twak algorithm. Then reverting back to double coordinates for rendering. This one eliminates artifacts to a minimum state, especially on "very small" polygons. It is not perfect but it is a significant improvement, artifacts now are rare and controllable. I could not do this for Kendzi, though, but now Twak is at least usable in my specific case. — Reply to this email directly, view it on GitHub<#133 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACG7RKRLMU423XEOCH7NZ332WWQEZAVCNFSM6AAAAABUL5MC72VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTENRVHAYTAMQ>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[istinnstudio]

Big one, like a million give ok results, even Integer.MAX_VALUE can be used, it is a trial n error situation.

[istinnstudio]

Those algorithms behave unexpectedly if pushed hard. Before this "mod" there were frequent missing faces, now those are rare, but I get overlaid faces instead, most times in a symmetrical form, it is operation sensitive though (e.g. rotation, scaling). This is the code, it is AI in most parts, I paste it for reference:

  // --- Main Method with Scaling Integration ---
	private static PShape straightSkeletonTwakBig(Polygon polygon) {
		if (polygon == null || polygon.isEmpty()) {
            System.err.println("Input polygon is null or empty.");
            return new PShape(PConstants.GROUP);
        }
            
    polygon =   cleanPolygonWithBufferZero(polygon, 0);

        // --- 1. Scaling Setup ---
        CoordinatesScalerForPGS scaler = new CoordinatesScalerForPGS();
//        int scaleFactor = 1_000_000_000;
        int scaleFactor = Integer.MAX_VALUE;
        CoordinatesScalerForPGS.ScalingData scalingData;
        List<Pair<BigFraction, BigFraction>> scaledBigFractionPoints;
        List<Pair<Double, Double>> allCoordsInput = new ArrayList<>();
        List<Integer> ringSizes = new ArrayList<>(); // Store size of each ring (incl. duplicate end point)
        // This set will be populated by ringToLoop with SCALED coordinates
        final Set<Coordinate> edgeCoordsSet_Scaled = new HashSet<>();


        // --- 2. Extract ORIGINAL Coordinates (for scaling) ---
        // We do NOT populate edgeCoordsSet_Scaled here. Original coords are only needed for scaler input.
        try {
            // Process Exterior Ring for scaling input
            LinearRing exteriorRing = polygon.getExteriorRing();
            Coordinate[] exteriorCoords = exteriorRing.getCoordinates();
            if (exteriorCoords.length < 4) {
                System.err.println("Error: Exterior ring has too few points (<4).");
                return new PShape(PConstants.GROUP);
            }
            ringSizes.add(exteriorCoords.length);
            for (Coordinate c : exteriorCoords) {
                allCoordsInput.add(Pair.of(c.getX(), c.getY()));
            }

            // Process Interior Rings (Holes) for scaling input
            for (int i = 0; i < polygon.getNumInteriorRing(); i++) {
                LinearRing holeRing = polygon.getInteriorRingN(i);
                Coordinate[] holeCoords = holeRing.getCoordinates();
                if (holeCoords.length < 4) {
                    System.out.println("Warning: Skipping hole with too few points (<4).");
                    continue;
                }
                ringSizes.add(holeCoords.length);
                for (Coordinate c : holeCoords) {
                    allCoordsInput.add(Pair.of(c.getX(), c.getY()));
                }
            }

            if (allCoordsInput.isEmpty()) {
                 System.err.println("No valid coordinates found in the polygon after filtering.");
                 return new PShape(PConstants.GROUP);
            }

            // --- 3. Perform Scaling ---
            Pair<CoordinatesScalerForPGS.ScalingData, List<Pair<BigFraction, BigFraction>>> result = scaler.scalePoints(scaleFactor, allCoordsInput);
            scalingData = result.getLeft();
            scaledBigFractionPoints = result.getRight();

        } catch (ParseExceptionForPGS e) { // Adjust exception type if needed
            System.err.println("Error during coordinate scaling: " + e.getMessage());
             return new PShape(PConstants.GROUP);
        } catch (Exception e) {
             System.err.println("Unexpected error during coordinate extraction or scaling: " + e.getMessage());
             e.printStackTrace();
             return new PShape(PConstants.GROUP);
        }

        // --- 4. Reconstruct SCALED Rings and Prepare Input for twak ---
        final LoopL<Edge> loops = new LoopL<>(); // Use LoopL to hold multiple Loop<Edge>
		final Machine speed = new Machine(10);
        int currentPointIndex = 0;

        try {
            // Process Exterior Ring (isHole = false)
            int exteriorRingSize = ringSizes.get(0);
            Coordinate[] scaledExteriorCoords = new Coordinate[exteriorRingSize];
            for (int j = 0; j < exteriorRingSize; j++) {
                 Pair<BigFraction, BigFraction> bfPair = scaledBigFractionPoints.get(currentPointIndex + j);
                 scaledExteriorCoords[j] = new Coordinate(bfPair.getLeft().doubleValue(), bfPair.getRight().doubleValue());
            }
            LinearRing scaledExteriorRing = GEOMETRY_FACTORY.createLinearRing(scaledExteriorCoords);
            // Call ringToLoop for exterior, passing the set to be populated
            Loop<Edge> exteriorLoop = ringToLoop(scaledExteriorRing, false, edgeCoordsSet_Scaled, speed);
             if (exteriorLoop.count() > 0) { // Only add if loop is valid
                loops.add(exteriorLoop);
             }
            currentPointIndex += exteriorRingSize;

            // Process Interior Rings (isHole = true)
            for (int i = 1; i < ringSizes.size(); i++) {
                int holeRingSize = ringSizes.get(i);
                Coordinate[] scaledHoleCoords = new Coordinate[holeRingSize];
                 for (int j = 0; j < holeRingSize; j++) {
                    Pair<BigFraction, BigFraction> bfPair = scaledBigFractionPoints.get(currentPointIndex + j);
                    scaledHoleCoords[j] = new Coordinate(bfPair.getLeft().doubleValue(), bfPair.getRight().doubleValue());
                 }
                 try {
                      LinearRing scaledHoleRing = GEOMETRY_FACTORY.createLinearRing(scaledHoleCoords);
                      // Call ringToLoop for hole, passing the SAME set to add more points
                      Loop<Edge> holeLoop = ringToLoop(scaledHoleRing, true, edgeCoordsSet_Scaled, speed);
                      if (holeLoop.count() > 0) { // Only add if loop is valid
                         loops.add(holeLoop);
                      }
                 } catch (IllegalArgumentException illegalArgEx) {
                      System.err.println("Error creating scaled LinearRing for hole " + (i-1) + ". Skipping hole. Message: " + illegalArgEx.getMessage());
                 }
                 currentPointIndex += holeRingSize;
            }
        } catch(IllegalArgumentException illegalArgEx) {
             System.err.println("Error creating scaled LinearRing for exterior ring. Cannot proceed. Message: " + illegalArgEx.getMessage());
             return new PShape(PConstants.GROUP);
        } catch (Exception e) {
            System.err.println("Error reconstructing scaled rings or preparing twak input: " + e.getMessage());
            e.printStackTrace();
            return new PShape(PConstants.GROUP);
        }


        // --- 5. Run Skeletonization (on SCALED data) ---
        final PShape lines = new PShape(PConstants.GROUP);
		final PShape faces = new PShape(PConstants.GROUP);
        final Skeleton skeleton;
        // These sets store PEdge objects with ORIGINAL coordinates
		final Set<PEdge> branchEdges = new HashSet<>();
		final Set<PEdge> boneEdges = new HashSet<>();

        // --- 6. Get Inverse Transform ---
        Function<Pair<Double, Double>, Pair<BigFraction, BigFraction>> invTransform = CoordinatesScalerForPGS.inverseTransform(scalingData);

		try {
            if (loops.isEmpty()) {
                 System.err.println("No valid loops provided to Skeleton constructor. Skipping skeletonization.");
                 return new PShape(PConstants.GROUP);
            }
			skeleton = new Skeleton(loops, true);
			skeleton.skeleton();

            // --- 7. Process Results and UNSCALE ---
            final double tolerance = 1e-9; // Tolerance for comparing scaled coordinates

			skeleton.output.faces.values().forEach(f -> {
//                org.twak.utils.collections.Loop<?> faceLoopRaw = f.getLoopL().iterator().next();


                // --- MODIFICATION START ---
                LoopL<?> faceLoopList = f.getLoopL(); // Get the LoopL first

                // Check if the LoopL itself is null or empty OR contains no loops
                if (faceLoopList == null || faceLoopList.isEmpty()) {
                    System.err.println("Warning: Skipping face from twak output because it has no associated boundary loops. Face ID (hashCode): " + f.hashCode());
                    // Optionally add more context here if available, like the input polygon WKT if you have it accessible
                    return; // Continue to the next face in the forEach loop
                }

                org.twak.utils.collections.Loop<?> faceLoopRaw;
                try {
                    // Now that we know faceLoopList is not empty, try to get the first loop.
                    // We assume standard faces have exactly one loop in their LoopL.
                    faceLoopRaw = faceLoopList.iterator().next();

                    // Add an extra check: is the obtained loop itself null or empty/degenerate?
                    if (faceLoopRaw == null || faceLoopRaw.count() < 3) {
                         System.err.println("Warning: Skipping face from twak output because its boundary loop is null, empty, or degenerate (less than 3 points). Face ID (hashCode): " + f.hashCode());
                         return; // Continue to the next face
                    }

                } catch (NoSuchElementException e) {
                    // This catch block is now less likely due to the isEmpty() check,
                    // but acts as a safeguard against unexpected iterator behavior.
                    System.err.println("Warning: Skipping face due to unexpected NoSuchElementException while getting loop iterator.next(). Face ID (hashCode): " + f.hashCode());
                    return; // Continue to the next face
                }
                // --- MODIFICATION END ---

                
                List<Point3d> vertices_scaled = new ArrayList<>();
                 for (Object node : faceLoopRaw) { // Extract scaled vertices from twak output
                     if (node instanceof org.twak.camp.Corner) {
                         org.twak.camp.Corner corner = (org.twak.camp.Corner) node;
                         vertices_scaled.add(new Point3d(corner.x, corner.y, 0));
                     } else if (node instanceof Point3d) {
                         vertices_scaled.add((Point3d)node);
                     }
                 }

				List<PVector> faceVertices_original = new ArrayList<>(); // For drawing face shapes

				for (int i = 0; i < vertices_scaled.size(); i++) {
					final Point3d p1_scaled_pt = vertices_scaled.get(i);
					final Point3d p2_scaled_pt = vertices_scaled.get((i + 1) % vertices_scaled.size());

                    // Create Coordinate objects for SCALED points for checking
                    Coordinate p1_scaled_coord = new Coordinate(p1_scaled_pt.x, p1_scaled_pt.y);
                    Coordinate p2_scaled_coord = new Coordinate(p2_scaled_pt.x, p2_scaled_pt.y);

                    // Check if SCALED coordinates are near the SCALED boundary coordinates
                    final boolean a = edgeCoordsSet_Scaled.stream().anyMatch(orig_scaled -> orig_scaled.distance(p1_scaled_coord) < tolerance);
                    final boolean b = edgeCoordsSet_Scaled.stream().anyMatch(orig_scaled -> orig_scaled.distance(p2_scaled_coord) < tolerance);

                    // --- Unscale points back to original coordinate system for drawing ---
                    Pair<BigFraction, BigFraction> p1_orig_bf = invTransform.apply(Pair.of(p1_scaled_pt.x, p1_scaled_pt.y));
                    Pair<BigFraction, BigFraction> p2_orig_bf = invTransform.apply(Pair.of(p2_scaled_pt.x, p2_scaled_pt.y));

                    double p1x_orig = p1_orig_bf.getLeft().doubleValue();
                    double p1y_orig = p1_orig_bf.getRight().doubleValue();
                    double p2x_orig = p2_orig_bf.getLeft().doubleValue();
                    double p2y_orig = p2_orig_bf.getRight().doubleValue();

                    // Add unscaled vertex for PShape face
					faceVertices_original.add(new PVector((float)p1x_orig, (float)p1y_orig)); // I use double for modified PShape core

                    // Add edge to appropriate set IF the points are distinct enough
                    // Use ORIGINAL coordinates for the PEdge objects
                    if (p1_scaled_coord.distance(p2_scaled_coord) > tolerance) { // Check distinctness using scaled coords
                        if (a ^ b) { // branch (one point on scaled boundary, one not)
                            branchEdges.add(new PEdge(p1x_orig, p1y_orig, p2x_orig, p2y_orig));
                        } else {
                            if (!a) { // bone (neither point on scaled boundary)
                                boneEdges.add(new PEdge(p1x_orig, p1y_orig, p2x_orig, p2y_orig));
                            }
                            // else: Both a and b are true (within tolerance) - segment lies on original boundary, ignore.
                        }
                    }
				}

                // Create face PShape using ORIGINAL coordinates
                if (!faceVertices_original.isEmpty()) {
    				PShape face = PGS_Conversion.fromPVector(faceVertices_original);
                    if(face != null) {
        				face.setStroke(true);
        				face.setStrokeWeight(2);
        				face.setStroke(ColorUtils.composeColor(147, 112, 219)); // Purple using micycle
                        faces.addChild(face);
                    }
                }
			});
		} catch (Exception e) {
            System.err.println("Error during skeleton computation or output processing: " + e.getMessage());
			e.printStackTrace();
                      for (StackTraceElement element : e.getStackTrace()) {
        System.err.println("    at " + element.getClassName() + "." +
                           element.getMethodName() + "(" + element.getFileName() +
                           ":" + element.getLineNumber() + ")");
    }
		}

        // --- 8. Assemble Final PShape using ORIGINAL coordinates ---
		final PShape bones = prepareLinesPShape(null, null, 4);
		boneEdges.forEach(e -> {
			bones.vertex(e.a.x, e.a.y);
			bones.vertex(e.b.x, e.b.y);
		});
		bones.endShape();

		final PShape branches = prepareLinesPShape(ColorUtils.composeColor(40, 235, 180), null, 1); // Cyan using micycle
		branchEdges.forEach(e -> {
			branches.vertex(e.a.x, e.a.y);
			branches.vertex(e.b.x, e.b.y);
		});
		branches.endShape();

		lines.addChild(faces);
		lines.addChild(branches);
		lines.addChild(bones);

		return lines;
	}

[micycle1]

What's the code for collide() now? The WIP PGS v2.1 points at a new version of campskeleton/twak jutils that has come changes I've made. I changed collide() from using a matrix solver to a direct solution (though precision more-or-less the same, but much it's much faster): https://github.com/twak/jutils/blob/9b771cb488b29a38a82e1cea5400f3c77c1772fe/src/org/twak/utils/geom/LinearForm3D.java#L179 So I'm curious what a more robust version would look like.

…

________________________________ From: istinnstudio ***@***.***> Sent: 03 April 2025 11:41 To: micycle1/PGS ***@***.***> Cc: Michael Carleton ***@***.***>; Author ***@***.***> Subject: Re: [micycle1/PGS] Rounding errors and side effects with straight skeleton algorithm (Discussion #133) I have managed to enhance precision a bit further by using a combination of an other coordinate scaler (affine) with a specific and "sensitive" small scale factor of 10 (e.g using 100 breaks some faces again) this time and an implementation of Shewchuk robust calculations for org.twak.utils.geom LinearForm3D Tuple3d collide method. I need more tests to ensure this. I still get some symmetrical overlaid faces so I guess this is a different hidden mystery in twak. But now results are pretty much consistent. I cannot consider this as a universal solution, just a reference. Those notes have been focused on my very special case and is a trial and error experimental process. — Reply to this email directly, view it on GitHub<#133 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACG7RKTXC5EI4H5Z4KNFNSL2XUF5NAVCNFSM6AAAAABUL5MC72VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTENZRGIZTMOA>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[istinnstudio]

robust_collide_experimental.zip
I cannot test it universally or if it has issues, it is AI work that I use experimentally. I think it is more stable than the simple one but I cannot ensure it. Probably needs tests and improvements. Luckily combined with an affine transformation scaler and inversion it produces the best results for me even with scale set to 1. I do not know why that happens but in those sensitive areas of unstable numbers unpredictable things can happen.

[istinnstudio]

I wonder if I can expose weights of the twak now, do you think this is possible? Is this the speed parameter in Machine class? If I change this from the default value of 10 it cannot work, in fact values other than 9.5 to 10.5 cannot work. I might look at this as it could introduce interesting results if used right. Maybe using an optional map of weights.

[micycle1]

I do believe it's the speed parameter.

Variable weights do work (this is from the camp org.twajk.camp.debug.CampSkeleton -- yellow points denote edge weights)

[istinnstudio]

it is interesting, especially if used in round symmetry or something.

[istinnstudio]

I also have a simpler version of the robust collide that seems to also work for me, here it is, might be faster :

 package org.twak.utils.geom.robust;


/**
 * A simplified Java translation of key parts of Shewchuk's robust arithmetic,
 * adapted for our 3×3 determinant needs.
 */
public class ShewchukRobustSimpler {

    /**
     * Computes the sum of two doubles with a compensated error.
     * @param a first value
     * @param b second value
     * @return an array where index 0 is the sum and index 1 is the round-off error.
     */
    public static double[] twoSum(double a, double b) {
        double sum = a + b;
        double bVirtual = sum - a;
        double aVirtual = sum - bVirtual;
        double bRoundoff = b - bVirtual;
        double aRoundoff = a - aVirtual;
        double error = aRoundoff + bRoundoff;
        return new double[]{sum, error};
    }

    /**
     * Computes the product of two doubles with a compensated error.
     * @param a first value
     * @param b second value
     * @return an array where index 0 is the product and index 1 is the round-off error.
     */
    public static double[] twoProduct(double a, double b) {
        double product = a * b;
        // Using Math.fma (if available) to compute the error term.
        double error = Math.fma(a, b, -product);
        return new double[]{product, error};
    }
    
    /**
     * Computes the determinant of a 3x3 matrix robustly.
     * The matrix m is expected to be a 3x3 array.
     * 
     * The determinant is computed as:
     *   m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1])
     * - m[0][1]*(m[1][0]*m[2][2] - m[1][2]*m[2][0])
     * + m[0][2]*(m[1][0]*m[2][1] - m[1][1]*m[2][0])
     * 
     * @param m the 3x3 matrix
     * @return the robust determinant value
     */
    public static double robustDeterminant3x3(double[][] m) {
        // First term: m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1])
        double[] prod1 = twoProduct(m[1][1], m[2][2]);
        double[] prod2 = twoProduct(m[1][2], m[2][1]);
        double term1 = m[0][0] * (prod1[0] - prod2[0]);
        
        // Second term: m[0][1]*(m[1][0]*m[2][2] - m[1][2]*m[2][0])
        double[] prod3 = twoProduct(m[1][0], m[2][2]);
        double[] prod4 = twoProduct(m[1][2], m[2][0]);
        double term2 = m[0][1] * (prod3[0] - prod4[0]);
        
        // Third term: m[0][2]*(m[1][0]*m[2][1] - m[1][1]*m[2][0])
        double[] prod5 = twoProduct(m[1][0], m[2][1]);
        double[] prod6 = twoProduct(m[1][1], m[2][0]);
        double term3 = m[0][2] * (prod5[0] - prod6[0]);
        
        // Combine the terms robustly
        double det = term1 - term2 + term3;
        return det;
    }
}


/**
 * Finds the intersection point between this plane and two others in 3D space,
 * using a robust determinant computation based on Shewchuk's arithmetic.
 * 
 * @param b second plane
 * @param c third plane
 * @return the intersection point as a Point3d
 * @throws Error if the planes are nearly degenerate (i.e. do not intersect at a single point)
 */
public Tuple3d collide(final LinearForm3D b, final LinearForm3D c) {
    final LinearForm3D a = this;
    
    if (a.hasNaN() || b.hasNaN() || c.hasNaN()) {
        throw new Error("Invalid plane data.");
    }
    
    // Construct the matrix from the normals (first three coefficients) of the planes
    double[][] detMatrix = {
        {a.A, a.B, a.C},
        {b.A, b.B, b.C},
        {c.A, c.B, c.C}
    };
    double det = ShewchukRobustSimpler.robustDeterminant3x3(detMatrix);
    
    // Use an appropriate tolerance for near-degenerate cases
    if (Math.abs(det) < 1e-10) {
        throw new Error("Planes do not intersect at a single point.");
    }
    
    // Build submatrices for Cramer's rule to solve for x, y, z.
    double[][] Mx = {
        {-a.D, a.B, a.C},
        {-b.D, b.B, b.C},
        {-c.D, c.B, c.C}
    };
    double[][] My = {
        {a.A, -a.D, a.C},
        {b.A, -b.D, b.C},
        {c.A, -c.D, c.C}
    };
    double[][] Mz = {
        {a.A, a.B, -a.D},
        {b.A, b.B, -b.D},
        {c.A, c.B, -c.D}
    };
    
    double x = ShewchukRobustSimpler.robustDeterminant3x3(Mx) / det;
    double y = ShewchukRobustSimpler.robustDeterminant3x3(My) / det;
    double z = ShewchukRobustSimpler.robustDeterminant3x3(Mz) / det;
    
    return new Point3d(x, y, z);
}

[istinnstudio]

I would like also to note that the result is also affected by the usage of clipper2. I mean positively. The coordinates are more unstable without using it, even with very small offsets that do not affect visuals. This is probably due to the fact that clipper2 uses internally big integer conversions, if I remember correctly. So the final result is a combination of factors that is indeed unpredictable and the strategy is when values are getting stable, note the status of those "intermediate" calculations and parameters that have been used. It is indeed trial and error. Anything can interfere. So because my case is on the extreme side, I am not sure this calculation alchemy can be universally accepted.

[micycle1]

I've done some more digging...

I think most of your problems are probably caused by the default tolerance value of 0.01 (1e-2) in HeightCollision.process().

Reducing this (even to 1e-3) fixes the problems in the demo you gave me.

It make sense that increasing the scale would help, as you've experienced, since this tolerance value becomes relatively smaller compared to the coordinates -- however too small, and deformities can occur again, since events are missed.

[istinnstudio]

Thanks for this, I have seen those magic numbers, there are more of those in twak. I have tested 1e-3, there is a change but unfortunately I still get missing faces, it is a point for investigation, that is a useful finding. Maybe relate it directly to coordinate range or normalize it. I have also tested weights, this is an other monster that can give outstanding results when used carefully but hard to manage. There are magic numbers there also related to coords. For example the threshold for my cords is down below to 0.05, with a very narrow range of speed that is useable. So this is also coords dependent. Speed of 10 is on the safe side, so I guess it is ok as default value for unweighted. There are also issues of possibly wrong weighted skeletonization on the same shapes that appear rotated in the same grid. This looks wrong but can introduce "organic" results. Reordering points can change this but achieving full symmetry is challenging especially on odd sided polygons.

[micycle1]

Try other values, but it is indeed related to coordinate range (my polygons are screen units -- in the order of about 1000 wide/high).
Oddly, for smaller, unit-sized polygons I found a larger number is needed.

The best approach would be for the skeleton code to normalise the input (perhaps to unit square), and use suitable magic numbers based off that. But I suppose it's an weakness of this algorithm (in lieu of exact arithmetic) that is has to group co-sited collisions based on proximity, which is inherently indeterminate.