Sum area between groups of overlapping forms

Question

I am having trouble finding a clean way to calculate the remaining area between groups of curves plotted in Python in a given zone.

Here is an image to illustrate:

Each form is made of 2 half ellipses with different parameters for the height with parametric equation:

X     =Xc+A *cos(Theta)
Y_down=Yc+B1*sin(Theta)
Y_up  =Yc+B2*sin(Theta)

Along 1 line (X direction), the parametric equations are the same except for Xc. Along the Y axis (vertical direction), A varies along with Xc and Yc.

All the forms are made by iteration on the X axis and iteration on the Y axis. I used Zorder in the plot so that they overlap in the order of their creation.

The problem is that even if I can calculate the area of every form, I don't see how I can find the red area since those forms are overlapping in every way possible. For the moment, I can access the red area by plotting all the curves and binarizing the output figure and summing. But I would like to find a more analytical/elegant solution that does not depend on the DPI of the output figure. Or is there something else I could do?

Thank you! I hope I was clear.

Answer 1

This sounds like a line scan problem.

Imagine a line that moves from top to bottom, left to right, whatever direction is easier for you to handle. Let's say we move it top to bottom.

At every point the shapes that overlap with the line will generate intervals (open where the line encounters the shape and closes where the line leaves the shape. Given the intervals you can easily compute how much of the line is read (outside of any interval). This is O(N) in complexity.

Ok so now we need to move it top to bottom and sum the areas. But you don't want to move it pixel by pixel since this makes it DPI dependent. So note that the intervals move slightly as the line moves but they don't change shape/unless the shapes intersect at that point, there's a new shape that the line encounters and where you leave a shape above.

Take the min&max y for every shape and sort them. You also need intersections between every pair of overlapping shapes. But instead of computing them you can only compute the ones that are near each other on the scan line (keep them in a heap since you need to find the next one at every step).

So one way to do it is to move the line a bit (to the next point of interest), compute the area you've just swept with the line and add it to your total. This will be O(N^3), O(N) per line move and you might have O(N^2) line moves (assuming every shape overlaps with every other shape).

You can do it a bit faster O(N^2 log N) if you speed up the area computation. The only part of it you need to compute is the around the intervals that swap. This means that for every free interval you need to remember were you last updated the are.

I leave the math to you, but it's simple. If you consider an interval it's basically a rectangle padded by two ellipses. Area of the rectangle is simple and you just need to add the are outside of the ellipse slice.

Answer 2

This is geometry problem. You need to remove overlapping area from the total area sum. That is not an easy task. You can try to solve it by integration between the elliptic arc curves. For 2 ellipses is that easy but for more you would need first to decide which arc to use and if it is interior or exterior. That could devolve to tedious code.

Instead you can try to divide your scene to vertical slices with small enough width (accuracy of intergration) and just using hit test determine/sum all the non filled areas by ray casting.

1. process scene by some dx step

   the dx step will determine the result accuracy.

2. for each x collect all intersection points

   this is simple O(N) search where N is number of ellipses. Just compute two y coordinates of each ellipse per the current x and add it to list (if valid).

3. remove all intersection points that are inside of another ellipse

   This is O(n^2) where n is number of intersections. simply test if any intersection is inside any ellipse intersection points. If yes flag it for removal ... and when all done remove them. Do not delete them during the search as that would invalidate the results. This way you r list will contain only outer intersection points. Beware there might be duplicates present !!!

4. sort intersections by y

   This will speed up / simplify the next process.

5. cast ray from y=y_start

6. find first intersection in path

   it is simply first intersection point in the list. Just check if it is inside integrated range and skip/stop if not. If valid then add the distance to last point (y_start or y from last iteration). Also do not forget to clamp the y value to your y_end otherwise you could add areas not visible on screen to the area...

7. find end of this ellipse

   simply increment index pointing to the intersections list until y coordinate jumps (to handle duplicate values). If end of list use y_end as a y value ... add to area and stop.

8. loop #6 until y_end coordinate is hit or crossed

Here my C++/VCL implementation for this:

//---------------------------------------------------------------------------
struct _ellipse
    {
    double x,y,rx,ry0,ry1;
    _ellipse()  {}
    _ellipse(_ellipse& a)   { *this=a; }
    ~_ellipse() {}
    _ellipse* operator = (const _ellipse *a) { *this=*a; return this; }
    //_ellipse* operator = (const _ellipse &a) { ...copy... return this; }
    };
struct _hit
    {
    double y;   // hit y coordinate
    int ix;     // ellipse index
    int f;      // 0 means ry0, 1 means ry1 edge, >=2 means deleted
    _hit()  {}
    _hit(_hit& a) { *this=a; }
    ~_hit() {}
    _hit* operator = (const _hit *a) { *this=*a; return this; }
    //_hit* operator = (const _hit &a) { ...copy... return this; }
    };
const int N=50;
_ellipse ell[N];
//---------------------------------------------------------------------------
void sort_asc_bubble(_hit *a,int n)
    {
    int i,e; _hit q;
    for (e=1;e;n--)
     for (e=0,i=1;i<n;i++)
      if (a[i-1].y>a[i].y)
       { q=a[i-1]; a[i-1]=a[i]; a[i]=q; e=1; }
    }
//---------------------------------------------------------------------------
void init()
    {
    int i; double xs=512,ys=512;
    _ellipse a;
    RandSeed=587654321;
    for (i=0;i<N;i++)
        {
        a.x=xs*Random();
        a.y=ys*Random();
        a.rx =xs*(0.02+(0.25*Random()));
        a.ry0=ys*(0.02+(0.09*Random()));
        a.ry1=ys*(0.02+(0.09*Random()));
        ell[i]=a;
        }
    }
//---------------------------------------------------------------------------
double area()
    {
    double area,dx=1.0;
    double x,y,xx,y0,y1,rxx,ryy0,ryy1;
    _ellipse *p;
    _hit  hit[N+N];     // intersection points
    int  i,j,n;         // n = number of intersections
    int _in;
    for (area=0.0,x=0.0;x<scr.xs;x+=dx)     // all vertical slices
        {
        // [prepare data]
        // O(N) precompute intersection points for ray/ellipses
        for (n=0,p=ell,i=0;i<N;i++,p++)
         if ((x>=p->x-p->rx)&&(x<=p->x+p->rx))
            {
            xx=x-p->x; xx*=xx;
            rxx =p->rx *p->rx ;
            ryy0=p->ry0*p->ry0;
            ryy1=p->ry1*p->ry1;
            hit[n].ix=i; hit[n].f=0; hit[n].y=p->y-sqrt(ryy0*(1.0-(xx/rxx))); n++;
            hit[n].ix=i; hit[n].f=1; hit[n].y=p->y+sqrt(ryy1*(1.0-(xx/rxx))); n++;
            }
        // O(n^2) flag inside edges
        for (i=0;i+1<n;i+=2)
            {
            y0=hit[i+0].y;
            y1=hit[i+1].y;
            for (j=0;j<n;j++)
             if ((i!=j)&&(i+1!=j))
              if ((hit[j].y>y0)&&(hit[j].y<y1))
               hit[j].f|=2;
            }
        // O(n) delete flagged edges
        for (i=0,j=0;i<n;i++)
         if (hit[i].f<2)
          { hit[j]=hit[i]; j++; } n=j;
        // O(n^2) sort y asc and original indexes are in i0,i1
        sort_asc_bubble(hit,n);

        // [integrate area]
        for (i=-1,y0=0.0,y1=0.0;;)
            {
            i++; if (i<n) y=hit[i].y; else y=scr.ys;
            if (y>scr.ys) y=y>scr.ys;
            if (y>y1)
                {
                y0=y1; y1=y;
                area+=y1-y0;
                // debug draw for visual check (render red pixel at x,y0 ... x,y1)
                //for (y=y0;y<=y1;y++) scr.pyx[int(y)][int(x)]=0x00FF0000;
                }
            if (i>=n) break;
            y=hit[i].y;
            for (;i<n;i++) if (hit[i].y>y) { y1=hit[i].y; break; }
            }
        } area*=dx; // rectangular rule
    return area;
    }
//---------------------------------------------------------------------------

Here result for dx=1.0 and resolution 512x512:

The number in the Window Caption is the computed area in [pixels^2] roughly 31% of the window surface. You can use any integration rule or dx step. Also you can change the bubble sort for something different but often for such small n the bubble beats others ... and is simple enough to code.

The code itself is not yet optimized ...