Subdividing square problems due to floating point precision errors

Question

I have a problem (in JAVA) which I unable to solve due to floating point precision errors.

I have an axis-aligned square class which is defined by specifying 2 points. In the constructor it is determined what the largest distance is (either in x- or y-direction) and this is used to create the square, with the points equal distance form the center. This means the points are on the boundary.

As an example, if I define a square with points (0, 2) and (3, 3) the largest distance is the x-distance (3) and the square will be defines as follows:

• lower left point = (0, 1)
• upper right point = (3, 4)

As can be seen the points are on the edge and the midpoint of the points is exactly the center of the square (1.5, 2.5).

I created a method to check if a square contains a certain point, which looks like something like this:

see added code sample below

Which means that if a point is on the boundary it is 'contained'.

Now I want to subdivide the square into 4 'equally' sized sections (NorthEast, NorthWest, SouthWest and SouthEast) and logically the initial points that defined the original square must be contained in at least 1 of the sections. But when unit testing this which random points it fails and it looks like it is because of double precision floating point errors.

I tried different methods of solving it and my last iteration was as follows:

• use the initial points to define the midpoint (p0) of the square
• use the initial points + midpoint to define the 4 corner points (p1, p2, p3 p4 in clockwise direction, starting bottom-left) of the square

This guarantees that the original points are contained and I can easily create the subdivisions and I thought it guarantees that the original points are contained it at least 1 of the sections if I ensure no math is performed on the characteristics of the points defining the square (since they are on the boundary). My subdivision routine is as follows:

see added code sample below

But when running a unit test of a large amount of iterations with random points almost 1 in 16 still fails and I have no idea why the edge points are changed. In all those tests the initial containment check (whether or not the parent square contains the points although they are on the edge) passes 100%.

EDIT Some actual code showing my implementation:

public class Point implements IPoint {
    double x, y;

    public Point(double x, double y) {
        this.x = x;
        this.y = y;
    }

    @Override
    public double x() {
        return x;
    }

    @Override
    public double y() {
        return y;
    }

    @Override
    public IPoint midPoint(IPoint other) {
        return new Point((x() + other.x()) / 2, (y() + other.y()) / 2);
    }
}

public class Rectangle implements IRectangle {
    IPoint p0, p1, p2, p3, p4;

    public Rectangle(IPoint v1, IPoint v2) {
        double dx, dy, dl;
        IPoint v0;

        // calculate dominant length
        dx = Math.abs(v1.x() - v2.x());
        dy = Math.abs(v1.y() - v2.y());
        dl = dx >= dy ? dx : dy;

        if (dx >= dy) {
            // make sure v0 = left-most
            if (v1.x() <= v2.x()) {
                v0 = v1;
                v1 = v2;
            } else {
               v0 = v2;
            }
        } else {
            // make sure v0 = bottom-most
            if (v1.y() <= v2.y()) {
                v0 = v1;
                v1 = v2;
            } else {
                v0 = v2;
            }
        }

        this.p0 = v0.midPoint(v1);
        if (dx >= dy) {
            // this way v0 and v1 are always on the vertical boundaries
            this.p1 = new Point(v0.x(), this.p0.y() - dl / 2);
            this.p2 = new Point(v0.x(), this.p0.y() + dl / 2);
            this.p3 = new Point(v1.x(), this.p0.y() + dl / 2);
            this.p4 = new Point(v1.x(), this.p0.y() - dl / 2);
        } else {
            // this way v0 and v1 are always on the horizontal boundaries
            this.p1 = new Point(this.p0.x() - dl / 2, v0.y());
            this.p2 = new Point(this.p0.x() - dl / 2, v1.y());
            this.p3 = new Point(this.p0.x() + dl / 2, v1.y());
            this.p4 = new Point(this.p0.x() + dl / 2, v0.y());
        }
    }

    @Override
    public boolean contains(IPoint p) {
        if (p.x() < p1.x() || p.x() > p4.x()) return false;
        if (p.y() < p1.y() || p.y() > p2.y()) return false;
        return true;
    }

    @Override
    public IRectangle[] subdivide() {
        return new Rectangle[] {
            new Rectangle(p0, p2),
            new Rectangle(p0, p3),
            new Rectangle(p0, p4),
            new Rectangle(p0, p1)
        };
    }
}

And here is the test case:

@Test
public void testMassedSubdivide() throws Exception {
    Random r = new Random();
    IPoint p1, p2;
    IRectangle[] rects;
    boolean check1, check2;
    for (int i = 0; i < 100000; i++) {
        p1 = new Point(r.nextDouble(), r.nextDouble());
        p2 = new Point(r.nextDouble(), r.nextDouble());
        q = new Rectangle(p1, p2);
        assertTrue(q.contains(p1));
        assertTrue(q.contains(p2));

        rects = q.subdivide();
        check1 = rects[0].contains(p1) || rects[1].contains(p1) || rects[2].contains(p1) || rects[3].contains(p1);
        check2 = rects[0].contains(p2) || rects[1].contains(p2) || rects[2].contains(p2) || rects[3].contains(p2);
        assertTrue(check1);
        assertTrue(check2);
    }
}

One of the failure cases resulting from my randomized test:

p1 = (0.31587198758298796,  0.12796964677511913)
p2 = (0.04837609765424089,  0.6711236142940149)

This one fails, because p1 should be in the SouthEast sector, but that one is defined as:

p0=(0.31791253449833834,    0.2637581386548431),    
p1=(0.18212404261861442,    0.12796964677511916), <- wrong, last 6 should be 3  
p2=(0.18212404261861442,    0.39954663053456707),   
p3=(0.4537010263780623,     0.39954663053456707),   
p4=(0.4537010263780623,     0.12796964677511916)  <- wrong, last 6 should be 3

Answer 1

After looking at your code, it doesn't make any sense that it would fail, given that as far as I can tell, you are not applying any operations to the Y values in that failure case--you're just passing it through without any manipulation so floating-point precision loss is irrelevant. It was a bit hard for me to follow though, when p1-p4 could represent any corner, so I rewrote the Rectangle class as follows in the hopes of being a bit clearer:

public class Rectangle implements IRectangle {
    IPoint centroid, bottomLeft, topLeft, bottomRight, topRight;

    public Rectangle(IPoint v0, IPoint v1) {
        IPoint bottomLeft = new Point(Math.min(v0.x, v1.x), Math.min(v0.y, v1.y));
        IPoint topRight   = new Point(Math.max(v0.x, v1.x), Math.max(v0.y, v1.y));

        // calculate dominant length
        double dx = topRight.x - bottomLeft.x;
        double dy = topRight.y - bottomLeft.y;  // Assumes (0, 0) is in the bottom-left.
        double dl = dx >= dy ? dx : dy;

        this.centroid = bottomLeft.midPoint(topRight);
        if (dx >= dy) {
            // this way bottomLeft and topRight are always on the vertical boundaries
            this.bottomLeft  = new Point(bottomLeft.x(), this.centroid.y() - dl / 2);
            this.topLeft     = new Point(bottomLeft.x(), this.centroid.y() + dl / 2);
            this.bottomRight = new Point(topRight.x(), this.centroid.y() - dl / 2);
            this.topRight    = new Point(topRight.x(), this.centroid.y() + dl / 2);
        } else {
            // this way bottomLeft and topRight are always on the horizontal boundaries
            this.bottomLeft  = new Point(this.centroid.x() - dl / 2, bottomLeft.y());
            this.topLeft     = new Point(this.centroid.x() - dl / 2, topLeft.y());
            this.bottomRight = new Point(this.centroid.x() + dl / 2, bottomLeft.y());
            this.topRight    = new Point(this.centroid.x() + dl / 2, topLeft.y());
        }
    }

    @Override
    public boolean contains(IPoint p) {
        if (p.x() < bottomLeft.x() || p.x() > topRight.x()) return false;
        if (p.y() < bottomLeft.y() || p.y() > topRight.y()) return false;
        return true;
    }

    @Override
    public IRectangle[] subdivide() {
        return new Rectangle[] {
            new Rectangle(centroid, bottomLeft),
            new Rectangle(centroid, topLeft),
            new Rectangle(centroid, bottomRight),
            new Rectangle(centroid, topRight)
        };
    }
}

If that doesn't fix the issue, perhaps some logging will suss out at point the 0.12796964677511913 changes to 0.12796964677511916