lib/gfxpoly/moments.c

   1 #include "moments.h"
   2
   3 #define MOMENTS
   4
   5 void moments_update(moments_t*moments, actlist_t*actlist, int32_t y1, int32_t y2)
   6 {
   7     segment_t* s = actlist_leftmost(actlist);
   8     segment_t* l = 0;
   9
  10     /* the actual coordinate of grid points is at the bottom right, hence
  11        we have to add 1.0 to both coordinates (or just 1.0 to the sum) */
  12     double mid = (y1+y2)/2.0 + 1.0;
  13
  14     double area = 0.0;
  15
  16     while(s) {
  17         if(l && l->wind.is_filled) {
  18             area += (XPOS(s,mid) - XPOS(l,mid));
  19
  20 #ifdef MOMENTS
  21             double dx1 = (l->b.x - l->a.x) / (double)(l->b.y - l->a.y);
  22             double o1 = l->a.x - l->a.y*dx1;
  23
  24             double dx2 = (s->b.x - s->a.x) / (double)(s->b.y - s->a.y);
  25             double o2 = s->b.x - s->b.y*dx2;
  26
  27 #ifdef DEBUG
  28             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",
  29                     y1*0.05, y2*0.05,
  30                     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,
  31                     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);
  32 #endif
  33
  34 #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))
  35             double m1x = S1(y2)-S1(y1);
  36 #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))
  37             double m2x = S2(y2)-S2(y1);
  38             moments->m[0][0] += (XPOS(s,mid) - XPOS(l,mid))*(y2-y1);
  39             moments->m[1][0] += m1x;
  40             moments->m[2][0] += m2x;
  41 #endif
  42         }
  43         l = s;
  44         s = s->right;
  45     }
  46
  47     area *= (y2-y1);
  48
  49     moments->area += area;
  50 }
  51 void moments_normalize(moments_t*moments, double gridsize)
  52 {
  53     moments->area *= gridsize*gridsize;
  54     moments->m[0][0] *= gridsize*gridsize;
  55     moments->m[1][0] *= gridsize*gridsize*gridsize*gridsize;
  56     moments->m[2][0] *= gridsize*gridsize*gridsize*gridsize*gridsize*gridsize;
  57 }