fixed moments calculcation

4 files changed, 27 insertions, 11 deletions

diff --git a/config.h.in b/config.h.in
index 60bf3227..9f122a03 100644
--- a/config.h.in
+++ b/config.h.in
@@ -277,10 +277,6 @@
 #endif
 #endif
 
-#ifndef WIN32
-#define CHECKS
-#endif
-
 // supply a substitute calloc function if necessary
 #ifndef HAVE_CALLOC
 #define calloc rfx_calloc_replacement
diff --git a/lib/gfxpoly/moments.c b/lib/gfxpoly/moments.c
index 19be93a3..f8db1a13 100644
--- a/lib/gfxpoly/moments.c
+++ b/lib/gfxpoly/moments.c
@@ -1,7 +1,8 @@
 #include "moments.h"
 
 #define MOMENTS
-void update_moments(moments_t*moments, actlist_t*actlist, int32_t y1, int32_t y2)
+
+void moments_update(moments_t*moments, actlist_t*actlist, int32_t y1, int32_t y2)
 {
     segment_t* s = actlist_leftmost(actlist);
     segment_t* l = 0;
@@ -18,16 +19,23 @@ void update_moments(moments_t*moments, actlist_t*actlist, int32_t y1, int32_t y2
 
 #ifdef MOMENTS
             double dx1 = (l->b.x - l->a.x) / (double)(l->b.y - l->a.y);
-            double o1 = - y1 * dx1;
+            double o1 = l->a.x - l->a.y*dx1;
 
             double dx2 = (s->b.x - s->a.x) / (double)(s->b.y - s->a.y);
-            double o2 = - y2 * dx2;
+            double o2 = s->b.x - s->b.y*dx2;
+
+#ifdef DEBUG
+            printf("y=%f-%f\n\tl1=([%f,%f,%f,%f] dx=%f o=%f)\n\tl2=([%f,%f,%f,%f] dx=%f o=%f)\n",
+                    y1*0.05, y2*0.05,
+                    l->a.x*0.05, l->a.y*0.05, l->b.x*0.05, l->b.y*0.05, dx1*0.05, o1*0.05,
+                    s->a.x*0.05, s->a.y*0.05, s->b.x*0.05, s->b.y*0.05, dx2*0.05, o2*0.05);
+#endif
 
 #define S1(y) 0.5*(1/3.0*(dx2*dx2-dx1*dx1)*(y)*(y)*(y)+1/2.0*(2*dx2*o2-2*dx1*o1)*(y)*(y)+(o2*o2-o1*o1)*(y))
             double m1x = S1(y2)-S1(y1);
 #define S2(y) (1/3.0)*(1/4.0*(dx2*dx2*dx2-dx1*dx1*dx1)*(y)*(y)*(y)*(y)+1/3.0*(3*dx2*dx2*o2-3*dx1*dx1*o1)*(y)*(y)*(y)+1/2.0*(3*dx2*o2*o2-3*dx1*o1*o1)*(y)*(y)+(o2*o2*o2-o1*o1*o1)*(y))
             double m2x = S2(y2)-S2(y1);
-            moments->m[0][0] += (XPOS(s,mid) - XPOS(l,mid));
+            moments->m[0][0] += (XPOS(s,mid) - XPOS(l,mid))*(y2-y1);
             moments->m[1][0] += m1x;
             moments->m[2][0] += m2x;
 #endif
@@ -40,3 +48,10 @@ void update_moments(moments_t*moments, actlist_t*actlist, int32_t y1, int32_t y2
 
     moments->area += area;
 }
+void moments_normalize(moments_t*moments, double gridsize)
+{
+    moments->area *= gridsize*gridsize;
+    moments->m[0][0] *= gridsize*gridsize;
+    moments->m[1][0] *= gridsize*gridsize*gridsize*gridsize;
+    moments->m[2][0] *= gridsize*gridsize*gridsize*gridsize*gridsize*gridsize;
+}
diff --git a/lib/gfxpoly/poly.c b/lib/gfxpoly/poly.c
index d911af16..d1c069a5 100644
--- a/lib/gfxpoly/poly.c
+++ b/lib/gfxpoly/poly.c
@@ -1546,7 +1546,7 @@ gfxpoly_t* gfxpoly_process(gfxpoly_t*poly1, gfxpoly_t*poly2, windrule_t*windrule
         actlist_verify(status.actlist, status.y-1);
 #endif
         if(moments && lasty > INT_MIN) {
-            update_moments(moments, status.actlist, lasty, status.y);
+            moments_update(moments, status.actlist, lasty, status.y);
         }
 
         xrow_reset(status.xrow);
@@ -1620,16 +1620,21 @@ double gfxpoly_area(gfxpoly_t*p)
     moments_t moments;
     gfxpoly_t*p2 = gfxpoly_process(p, 0, &windrule_evenodd, &onepolygon, &moments);
     gfxpoly_destroy(p2);
-    return moments.area * p->gridsize * p->gridsize;
+    moments_normalize(&moments, p->gridsize);
+    return moments.area;
 }
 double gfxpoly_intersection_area(gfxpoly_t*p1, gfxpoly_t*p2)
 {
     moments_t moments;
     gfxpoly_t*p3 = gfxpoly_process(p1, p2, &windrule_intersect, &twopolygons, &moments);
     gfxpoly_destroy(p3);
+
+    moments_normalize(&moments, p1->gridsize);
+
     printf("%f %f %f\n",
             moments.m[0][0],
             moments.m[1][0],
             moments.m[2][0]);
+
     return moments.area * p1->gridsize * p2->gridsize;
 }
diff --git a/lib/gfxpoly/test.c b/lib/gfxpoly/test.c
index 7238225c..7aeb3470 100644
--- a/lib/gfxpoly/test.c
+++ b/lib/gfxpoly/test.c
@@ -323,7 +323,7 @@ int test_area()
         gfxpoly_t*poly1 = gfxpoly_from_fill(line1, 0.05);
         gfxpoly_t*poly2 = gfxpoly_from_fill(line2, 0.05);
         double area = gfxpoly_intersection_area(poly1, poly2);
-        printf("%f\n", area);
+        //printf("%f\n", area);
 
         gfxpoly_destroy(poly1);
         gfxpoly_destroy(poly2);