3

I have a 2D convex polygon in 3D space and a function to measure the area of the polygon.

public double area() {
    if (vertices.size() >= 3) {
        double area = 0;
        Vector3 origin = vertices.get(0);
        Vector3 prev = vertices.get(1).clone();
        prev.sub(origin);
        for (int i = 2; i < vertices.size(); i++) {
            Vector3 current = vertices.get(i).clone();
            current.sub(origin);
            Vector3 cross = prev.cross(current);
            area += cross.magnitude();
            prev = current;
        }
        area /= 2;
        return area;
    } else {
        return 0;
    }
}

To test that this method works at all orientations of the polygon I had my program rotate it a little bit each iteration and calculate the area. Like so...

Face f = poly.getFaces().get(0);
        for (int i = 0; i < f.size(); i++) {
            Vector3 v = f.getVertex(i);
            v.rotate(0.1f, 0.2f, 0.3f);
        }
        if (blah % 1000 == 0)
            System.out.println(blah + ":\t" + f.area());

My method seems correct when testing with a 20x20 square. However the rotate method (a method in the Vector3 class) seems to introduce some error into the position of each vertex in the polygon, which affects the area calculation. Here is the Vector3.rotate() method

public void rotate(double xAngle, double yAngle, double zAngle) {
    double oldY = y;
    double oldZ = z;
    y = oldY * Math.cos(xAngle) - oldZ * Math.sin(xAngle);
    z = oldY * Math.sin(xAngle) + oldZ * Math.cos(xAngle);

    oldZ = z;
    double oldX = x;
    z = oldZ * Math.cos(yAngle) - oldX * Math.sin(yAngle);
    x = oldZ * Math.sin(yAngle) + oldX * Math.cos(yAngle);

    oldX = x;
    oldY = y;
    x = oldX * Math.cos(zAngle) - oldY * Math.sin(zAngle);
    y = oldX * Math.sin(zAngle) + oldY * Math.cos(zAngle);
}

Here is the output for my program in the format "iteration: area":

0:  400.0
1000:   399.9999999999981
2000:   399.99999999999744
3000:   399.9999999999959
4000:   399.9999999999924
5000:   399.9999999999912
6000:   399.99999999999187
7000:   399.9999999999892
8000:   399.9999999999868
9000:   399.99999999998664
10000:  399.99999999998386
11000:  399.99999999998283
12000:  399.99999999998215
13000:  399.9999999999805
14000:  399.99999999998016
15000:  399.99999999997897
16000:  399.9999999999782
17000:  399.99999999997715
18000:  399.99999999997726
19000:  399.9999999999769
20000:  399.99999999997584

Since this is intended to eventually be for a physics engine I would like to know how I can minimise the cumulative error since the Vector3.rotate() method will be used on a very regular basis.

Thanks!

A couple of odd notes:

  • The error is proportional to the amount rotated. ie. bigger rotation per iteration -> bigger error per iteration.

  • There is more error when passing doubles to the rotate function than when passing it floats.

1
  • The advantage of using floats maybe due to some truncation/rounding. This would hint that some kind of rounding scheme would give you less error. Also these are numbers are not particularly anything to worry about.
    – Mikhail
    Apr 8, 2012 at 10:22

2 Answers 2

9

You'll always have some cumulative error with repeated floating point trig operations — that's just how they work. To deal with it, you basically have two options:

  1. Just ignore it. Note that, in your example, after 20,000 iterations(!) the area is still accurate down to 13 decimal places. That's not bad, considering that doubles can only store about 16 decimal places to begin with.

    Indeed, plotting your graph, the area of your square seems to be going down more or less linearly:
    PLOT
    This makes sense, assuming that the effective determinant of your approximate rotation matrix is about 1 − 3.417825 × 10-18, which is well within normal double precision floating point error range of one. If that's the case, the area of your square would continue a very slow exponential decay towards zero, such that you'd need about two billion (2 × 109) 7.3 × 1014 iterations to get the area down to 399. Assuming 100 iterations per second, that's about seven and a half months 230 thousand years.

    Edit: When I first calculated how long it would take for the area to reach 399, it seems I made a mistake and somehow managed to overestimate the decay rate by a factor of about 400,000(!). I've corrected the mistake above.

  2. If you still feel you don't want any cumulative error, the answer is simple: don't iterate floating point rotations. Instead, have your object store its current orientation in a member variable, and use that information to always rotate the object from its original orientation to its current one.

    This is simple in 2D, since you just have to store an angle. In 3D, I'd suggest storing either a quaternion or a matrix, and occasionally rescaling it so that its norm / determinant stays approximately one (and, if you're using a matrix to represent the orientation of a rigid body, that it remains approximately orthogonal).

    Of course, this approach won't eliminate cumulative error in the orientation of the object, but the rescaling does ensure that the volume, area and/or shape of the object won't be affected.

2
  • Excellent and informed answer - far better than my 'laymans' understanding. Apr 8, 2012 at 11:31
  • Wow, thanks for putting the time for the area decay into perspective, I hadn't even thought of how long it would take. The rotation calls will never be as fast as 100 per second per body once the full simulation is happening. And besides, I already store a polyhedron as a set of base coords for the vertices and the orientation around each axis in a Vector3. I didn't make the logical connection that that would mitigate any accumulation. Excellent answer. Sorry @PaulSullivan, I'm going to have to change the correct answer to this one. Apr 8, 2012 at 11:43
1

You say there is cumulative error but I don't believe there is (note how your output desn't always go down) and the rest of the error is just due to rounding and loss of precision in a float.

I did work on a 2d physics engine in university (in java) and found double to be more precise (of course it is see oracles datatype sizes

In short you will never get rid of this behaviour you just have to accept the limitations of precision

EDIT:

Now I look at your .area function there is possibly some cumulative due to

+= cross.magnitude

but I have to say that whole function looks a bit odd. Why does it need to know the previous vertices to calculate the current area?

4
  • Thanks, I think you're right that I'll just have to live with these errors. I think the error is still cumulative since it is rotating an evermore inaccurate position vector each iteration. The area function takes the vector from the first vertex to vertices i and i+1. It then takes the area of the parallelogram made up of those 2 vectors (which works out to be the magnitude of the cross product) Then is just sums that for every i and then halves at the end since we wanted the sum of the triangles, not the parallelograms. Storing the previous vertex just speeds it up ever so slightly :P Apr 8, 2012 at 10:42
  • Ahh I think I see - actually what you are seeing is what led to chaos theory and the term 'butterfly effect' - A computer scientist used this imprecision to show that complex systems with very small variations quickly become chaotic @null0pointer Apr 8, 2012 at 11:09
  • That's pretty cool, I didn't know chaos theory arose out of computer science :) Apr 8, 2012 at 11:44
  • @PaulSullivan - Note that the error will be essentially a random walk, which will diverge in multiple dimensions. And, yes, the error is cumulative.
    – user85109
    Apr 9, 2012 at 0:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.