diff options
author | Matthias Kramm <kramm@quiss.org> | 2011-02-14 10:55:59 -0800 |
---|---|---|
committer | Matthias Kramm <kramm@quiss.org> | 2011-02-14 10:55:59 -0800 |
commit | 02f4644e27283162e68d56bde6a4b216fb642d01 (patch) | |
tree | 8e85dfdc2d9991c695fc060c56b7ff58c1f57894 | |
parent | bb2ad19a7b6f340d87faf231b085f0f22bc4c4bc (diff) |
fixed moments calculcation
-rw-r--r-- | config.h.in | 4 | ||||
-rw-r--r-- | lib/gfxpoly/moments.c | 23 | ||||
-rw-r--r-- | lib/gfxpoly/poly.c | 9 | ||||
-rw-r--r-- | lib/gfxpoly/test.c | 2 |
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); |