Python API: batch_hull fails bizarrely

[briansturgill]

batch_hull from Python is failing very strangely.

Further experimentation suggests that the number of segments being a multiple
of the y length, is heavily related to making it occur.

import manifold3d as m

def rounded_rectangle(x, y, radius, segments):
    circ = m.CrossSection.circle(radius, segments)
    o = m.CrossSection.batch_hull([
        circ.translate((radius, radius)),
        circ.translate((x - radius, radius)),
        circ.translate((x - radius, y - radius)),
        circ.translate((radius, y - radius)),
    ])
    return o

# You might find this helpful to view the problem
def _save_svg(filename, obj):
    txt = []
    bb = obj.bounds()
    width = round(bb[2] - bb[0], 5)
    height = round(bb[3] - bb[1], 5)
    color = (0, 128, 255)
    units = "mm"

    txt.append("<?xml version=\"1.0\" encoding=\"UTF-8\"?>")
    txt.append("<!-- Created by Piecad. -->")
    txt.append("<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1 Tiny//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11-tiny.dtd\">")
    txt.append(f"<svg width=\"{width}{units}\" height=\"{height}{units}\" viewBox=\"0 0 {width} {height}\" fill=\"rgb({color[0]},{color[1]},{color[2]})\" xmlns=\"http://www.w3.org/2000/svg\" fill-rule=\"evenodd\">")

    txt.append(f"<g><path d=\"")

    off_x = 0 - bb[0]
    off_y = 0 - bb[1]
    paths = obj.to_polygons()
    for path in paths:
        for i in range(len(path)):
            vert = path[i]
            x = round(vert[0] + off_x, 5)
            y = round(vert[1] + off_y, 5)
            if i == 0:
                txt.append(f"M{x} {y}")
            else:
                txt.append(f"L{x} {y}")

    txt.append("\"/></g>\n")

    txt.append("</svg>")

    with open(filename, "w") as f:
        f.write("\n".join(txt))


radius = 9.0
x = 51
y = 36
segments = 36

# These don't work
o = rounded_rectangle(x, y, radius, segments)
print(o.num_vert())
_save_svg("test.svg", o)
o = rounded_rectangle(x, y, radius+3, segments)
print(o.num_vert())
o = rounded_rectangle(x+1, y, radius, segments)
print(o.num_vert())
o = rounded_rectangle(x-1, y, radius, segments)
print(o.num_vert())

# These work fine
o = rounded_rectangle(x, y, radius-6, segments)
print(o.num_vert())
o = rounded_rectangle(x, y, radius-3, segments)
print(o.num_vert())
o = rounded_rectangle(x, y, radius-1, segments)
print(o.num_vert())
o = rounded_rectangle(x, y, radius+1, segments)
print(o.num_vert())
o = rounded_rectangle(x, y+1, radius, segments)
print(o.num_vert())
o = rounded_rectangle(x, y+2, radius, segments)
print(o.num_vert())
o = rounded_rectangle(x, y-1, radius, segments)
print(o.num_vert())
o = rounded_rectangle(x+1, y+1, radius, segments)
print(o.num_vert())

[elalish]

Thanks - when you say "don't work", do you mean it crashes or you get the wrong output? Any chance you can toss a few screenshots in?

[briansturgill]

@elalish My apologies, I've been looking at it too long.
By fail, I mean incorrect output... if you run the script, it will produce test.svg which will show you an incorrect result.
With segment size fixed, the number of verts should be the same no matter what the x, y or radius.
The working ones have 40 (proper for 36 segments). The broken ones have a number of verts other than 40.

[briansturgill]

Here's a screen shot from the _save_svg in the code above.

[elalish]

I wouldn't bet on having exactly the same number of verts - floating point error, etc. But those big gouges look like a serious problem. We use quickhull for this, so maybe an upstream bug? There's also been talk of switching to a better dependency for more numerically-stable hulls - this may be good impetus.

[briansturgill]

You use quick hull for a 3D hull. The example I am giving is 2D as is written in cross_section.cpp.
It is not a part of Clipper2.
Here's an implementation I made of a Graham scan for CSharpCAD.
I use a fake atan2 which greatly speeds the algorithm.

namespace CSharpCAD;

public static partial class CSCAD
{
    private class pt
    {
        internal readonly Vec2 point;
        internal readonly double angle;
        internal readonly double distSq;
        internal pt(Vec2 point, double angle, double distSq)
        {
            this.point = point;
            this.angle = angle;
            this.distSq = distSq;
        }
    };
    /*
     * Create a convex hull of the given geom2 geometries.
     * Uses https://en.wikipedia.org/wiki/Graham_scan
     */
    internal static Geom2 HullGeom2(params Geom2[] geometries)
    {
        // Extract unique points from the geometries.
        // To avoid a second pass, also determine the minimum point.
        var uniquePoints = new HashSet<Vec2>();
        var min = new Vec2(double.PositiveInfinity, double.PositiveInfinity);
        foreach (var g in geometries)
        {
            var outlines = g.ToOutlines();
            foreach (var outline in outlines)
            {
                foreach (var point in outline)
                {
                    uniquePoints.Add(point);
                    if (point.Y < min.Y || (point.Y == min.Y && point.X < min.X))
                    {
                        min = point;
                    }
                }
            }
        }

        // Returned "angle" is really 1/tan (inverse of slope) made negative to increase with angle.
        // This function is strictly for sorting in this algorithm.
        double fakeAtan2(double y, double x)
        {
            // The "if" is a special case for when the minimum vector found in loop above is present.
            // We need to insure that it sorts as the minimum point. Otherwise this becomes NaN.
            if (y == 0 && x == 0) { return double.NegativeInfinity; }
            return -(x / y);
        }

        // Gather information for sorting by polar coordinates.
        var points = new List<pt>(uniquePoints.Count);
        foreach (var v in uniquePoints)
        {
            // Use of fakeAtan2 avoids use of Math.Atan2 which slows things down.
            // A simple timing test suggests it saves about 10% of total time.
            var angle = fakeAtan2((v.Y - min.Y), (v.X - min.X));
            //var angle = Math.Atan2((v.Y - min.Y), (v.X - min.X));
            var distSq = min.SquaredDistance(v);
            points.Add(new pt(v, angle, distSq));
        }
        points.Sort((pt pt1, pt pt2) => pt1.angle < pt2.angle ? -1 : pt1.angle > pt2.angle ? 1 :
            pt1.distSq < pt2.distSq ? -1 : pt1.distSq > pt2.distSq ? 1 : 0);

        // ccw returns:  < 0 clockwise, 0 colinear, > 0 counter-clockwise.
        double ccw(Vec2 v1, Vec2 v2, Vec2 v3) => (v2.X - v1.X) * (v3.Y - v1.Y) - (v2.Y - v1.Y) * (v3.X - v1.X);
        var stack = new List<Vec2>(); // Start with empty stack
        foreach (var point in points)
        {
            var cnt = stack.Count;
            while (cnt > 1 && ccw(stack[cnt - 2], stack[cnt - 1], point.point) <= C.EPSILON)
            {
                stack.RemoveAt(stack.Count - 1); // Pop - gets rid of colinear and interior (clockwise) points.
                cnt = stack.Count;
            }
            stack.Add(point.point); // Push
        }

        if (stack.Count < 3) return new Geom2();

        return new Geom2(stack);
    }
}

[geoffder]

Actually, that would be the fault of the simple 2D hull algorithm (or my porting of it) used in CrossSection, we don't use quickhull there. I don't have time to contribute to a fix for it right now though unfortunately.

[elalish]

Ah, good to know. Well, now that we have quickhull as a dependency, perhaps we should switch to that (assuming it does 2D as well as 3D?). Or whatever other dependency we decide to switch to. I'd rather not maintain our own 2D hulling code if I can avoid it.

[briansturgill]

@elalish QuickHull is 3D only. The 2D algorithm is simple as you can see above. If you like I can take a stab at a PR.

[elalish]

Yeah, ideally see if you can fix our existing code, rather than replace it wholesale, but whatever works. Thanks!