finding the area of a closed 2d uniform cubic B-sp

问题:

I have a list of 2d points which are the control vertices (Dx) for a closed uniform cubic B-spline. I am assuming a simple curve (non-self-intersecting, all control points are distinct).

I am trying to find the area enclosed by the curve:

If I calculate the knot points (Px), I can treat the curve as a polygon; then I "just" need to find the remaining delta areas, for each segment, between the actual curve and the straight line joining the knot points.

I understand that the shape (and therefore area) of a Bspline is invariant under rotation and translation - so for each segment, I can find a translation to put the t=0 knot at the origin and a rotation to put the t=1 knot on the +x axis:

I can find the equation for the curve by plugging the points in and re-grouping:

P(t) = (
    (t**3)*(-Dm1 + 3*D0 - 3*D1 + D2)
    + (t**2)*(3*Dm1 - 6*D0 + 3*D1)
    + t*(-3*Dm1 + 3*D1)
    + (Dm1 + 4*D0 + D1)
) / 6.

but I'm tearing my hair out trying to integrate it - I can do

 1
/
| Py(t) dt
/
t=0

but that doesn't give me area. I think what I need is

 Px(t=1)
/
| Py(t) (dPx(t) / dt) dt
/
x = Px(t=0)

but before I go further, I'd really like to know:

1. Is this the correct calculation for area? Ideally, an analytic solution would make my day!

2. Once I find this area, how can I tell whether I need to add or subtract it from the base poly (red vs green areas in the first diagram)?

3. Are there any Python modules which would do this calculation for me? Numpy has some methods for evaluating cubic Bsplines, but none which seem to deal with area.

4. Is there an easier way to do this? I am thinking about maybe evaluating P(t) at a bunch of points - something like t = numpy.arange(0.0, 1.0, 0.05) - and treating the whole thing as a polygon. Any idea how many subdivisions are needed to guarantee a given level of accuracy (I'd like error < 1%)?

回答1:

1. Pick some arbitrary point as pivot p₀ (e.g. the origin (0,0))
2. Pick some point along the curve p₁ = (x,y)
3. Differentiate the curve at that point, to get a velocity v = <vx,vy>
4. Form a triangle from the three points, and calculate the area. Easiest done with a cross product between the vectors p₀p₁ and v, and dividing by two.
5. Integrate this area over t, from 0 to 1.

The result you get is the area for one segment of the whole curve. Some will be negative, as they face the opposite direction. If you sum up the areas for all segments, you get the area of the whole curve.

The result is:

Area_i = (31 x_i-1 y_i + 28 x_i-1 y_i+1 + x_i-1 y_i+2 - 31 x_i y_i-1 + 183 x_i y_i+1 + 28 x_i y_i+2 - 28 x_i+1 y_i-1 - 183 x_i+1 y_i + 31 x_i+1 y_i+2 - x_i+2 y_i-1 - 28 x_i+2 y_i - 31 x_i+2 y_i+1) / 720

If you convert it into matrix form, you get:

Area_i = <x_i-1   x_i   x_i+1   x_i+2> · P · <y_i-1   y_i   y_i+1   y_i+2>^T

where P is

[    0   31   28    1]
[  -31    0  183   28] / 720
[  -28 -183    0   31]
[   -1  -28  -31    0]

---

If the control points are [(0,0) (1,0) (1,1) (0,1)], the resulting areas become:

[(0,0), (1,0), (1,1), (0,1)] -> 242/720
[(1,0), (1,1), (0,1), (0,0)] -> 242/720
[(1,1), (0,1), (0,0), (1,0)] ->   2/720
[(0,1), (0,0), (1,0), (1,1)] ->   2/720

The sum is 488/720 = 61/90.

回答2:

Personally, I would use the splines to their best advantage and rewrite the area integral as a contour integral using Green's theorem.

Since you already know the curve, it'll be an easy matter to do the integration using Gaussian quadrature of sufficient order. No shenanigans to estimate the extra areas needed. I'll bet it'll be computationally efficient as well, because Gaussian quadrature over polynomials is so well behaved. Cubic B-splines will integrate nicely.

I would write your code in such a way that the first and last points must be coincident. That's an invariant on the problem.

This approach will even work for areas with holes in them. You integrate along the outer curve, make an imaginary straight line that connects the outer to the inner curve, integrate along the inner curve, then pass back to the outer along the straight line that got you there.

回答3:

Well... it looks like you already know what to do.

1. Yes, this is correct way of calculating area under the segment, indeed

   dS=Y*dX=Y*(dX/dt)*dt
   

2. You don't need to care about adding or subtracting. You need to add all the time, but some of the integrals (red segments) will be negative (if you always set put your P_n in the beginning of coordinates and P_n+1 along X axis positive direction). So for every segment you need to calculate translation and rotation and then integral (all done analytically).

3. I don't know about Python module, but it seems that making such a module will take a day at most.
4. I think analytic solution will be better and is not so hard.

Anyway the graphics is absolutely stunning