CGAL users discussion list

[Cody Rose <>]

Hello,

I've stumbled into some CGAL behavior I can't explain and I was hoping 
that someone could help me understand it. In my application, I'm using 
Newell's method to calculate the approximate normal for a polygon-like 
series of points (http://cs.haifa.ac.il/~gordon/plane.pdf). (The idea is 
that if all the points are actually coplanar, it will yield the true 
normal, but if they aren't quite, it will give a "best fit.") The 
problem I'm having is essentially this, given a particular affine 
transformation T and the kernel Cartesian<leda_real>:

newell(T(points)) != T(newell(points))

And I can't figure out why not, since I was under the impression that 
this kernel is exact (and that this should work in an exact kernel). I 
can't chalk it up to variations in how Newell's method approximates a 
plane fit before and after the transformation, because in the example I 
attached, all the points are actually coplanar after all. (Honestly that 
explanation wouldn't make sense to me anyway.)

So I feel like I'm missing something about the way exactness (or this 
affine transformation) works. Can anyone clear this up?

Cody Rose
#define CGAL_USE_LEDA
#define CGAL_LEDA_VERSION 630
#define LEDA_DLL

#include <vector>
#include <cstdio>

#include <CGAL/leda_real.h>
#include <CGAL/Cartesian.h>

typedef leda_real                                               NT;
typedef CGAL::Cartesian<NT>                             K;
typedef CGAL::Direction_3<K>                    direction_3;
typedef CGAL::Point_3<K>                                point_3;
typedef CGAL::Aff_transformation_3<K>   transformation_3;

direction_3 newell(const std::vector<point_3> & loop) {
        NT a(0.0), b(0.0), c(0.0);
        for (size_t i = 0; i < loop.size(); ++i) {
                auto & curr = loop[i];
                auto & next = loop[(i+1) % loop.size()];
                a += (curr.y() - next.y()) * (curr.z() + next.z());
                b += (curr.z() - next.z()) * (curr.x() + next.x());
                c += (curr.x() - next.x()) * (curr.y() + next.y());
        }
        return direction_3(a, b, c);
}

int main() {
        transformation_3 t(0.00000000000000000, -0.97175545666443730, 
0.23599011090062880, -6363.1356980000000,
                                           1.0000000000000000,  
0.00000000000000000,  0.00000000000000000, 871.99942999999996,
                                           0.00000000000000000, 
0.23599011090062882,  0.97175545666443719, 216.00000000000000);

        // These points are all coplanar...
        std::vector<point_3> v;
        v.push_back(point_3(0, 0, 0));
        v.push_back(point_3(1, 0, 0));
        v.push_back(point_3(1, 1, 0));
        v.push_back(point_3(0, 1, 0));

        /// so normal_before = <0, 0, 1>.
        direction_3 normal_before = newell(v);
        for (auto p = v.begin(); p != v.end(); ++p) { *p = t(*p); }
        direction_3 normal_after = newell(v);

        direction_3 transformed = t(normal_before);

        if (t(normal_before) == normal_after) {
                // This is what I expected to happen.
                printf("equal!\n");
        }
        else {
                // This is what actually happens!
                printf("not equal!\n");
                NT ddx = transformed.dx() - normal_after.dx();
                NT ddy = transformed.dy() - normal_after.dy();
                NT ddz = transformed.dz() - normal_after.dz();
                if (CGAL::is_zero(ddx)) { printf("dx difference is zero\n"); }
                else { printf("dx difference is %+2.25f\n", 
CGAL::to_double(ddx)); }
                if (CGAL::is_zero(ddy)) { printf("dy difference is zero\n"); }
                else { printf("dy difference is %+2.25f\n", 
CGAL::to_double(ddy)); }
                if (CGAL::is_zero(ddz)) { printf("dz difference is zero\n"); }
                else { printf("dz difference is %+2.25f\n", 
CGAL::to_double(ddz)); }
        }

        return 0;
}