Estimating area of irregular region

Suppose you have an irregular region with continuous boundary like this:

I found that you can estimate the area pretty well by finding points roughly evenly spread around the perimeter xi,yi...

Then using this formula over all the points xi,yi:

Using the computer to count black pixels in the image gave an area of:

185604

and using the formula over the points gave:

182466

Using more points on the perimeter may give closer approximations.

My reasoning in coming up with the formula was first I find xc,yc which is the centroid of the sample points around the perimeter of the curve. Then find the average distance from that centroid to each of the sample points. Then I reason that the area of the irregular region should be roughly equal to a circle of radius that average distance to the centroid. It's as if you somehow reshape the irregular region into a circle in a way that preserves center of mass. Then find the area of the resulting circle. 

I don't know if it always works but it seems likely given how random the shape I started out with was and that it ended up within 1% of the right answer. 

**Three dimensions**

where the list of points could be in some 3 dimensional configuration like so:

The constraint is that it can only calculate volumes of solids topologically equivalent to a sphere...

plotting differential equations with a recurrence relationship

So imagine we have the differential equation:
y = -dy_2
saying the y value is related to -1 times the second derivative.
Suppose y, y1, and y2 are 3 points separated by a certain delta x. Then we can solve for y in terms of y1 and y2 like so:

Here I've substituted an expression for the slope of change in slopes as you go from y2 to y1 and then to y, and placed that instead of -dy_2 in the original differential equation. 

Then if you know two starting points, like y1=0.8572989892, y2=0.8521080219 for a delta x of 1/100. Then you can generate a new y value, then make y2=y1, y1=y and repeat generating the graph....

find the next prime after a list of known primes

Ok suppose you have a list of known primes like:
{2,3,5,7}
You'd like to know if the next odd number is prime:
So assume you're list of primes is:
{2,3,5,7,9}
Of the 9 numbers from 1 to 9 we know that roughly 9/2 -1 will be composite with a factor of 2.
So 9-3 = 6
Of these 6, 6/3 - 1 will divide 3 but not 2.
6-1=5 possible remaining primes
of these 5  5/5 -1 will divide 5 but not any previous prime
5-0=5
And then the same kind of thing for 9 so:
5-0 = 5.... saying that there are 5 primes up to 9, but one more thing is if you think of the 5 numbers this process leaves in the list, the number 1 hasn't been accounted for. So 5-1 = 4, saying there are only 4 primes in our list contradicting the fact that there are 5. Therefore 9 is not prime.  

Now suppose we skip 9 as a possible prime and move on to 11...
{2,3,5,7,11}
11/2 -1 = 4
11-4 = 7
7/3-1 = 1
7-1 = 6
6/5 - 1 is 0
6-0 is still 6 and
6/7 -1 is 0
and 6/11 is 0,
so then we do the final step as we did above of removing the number 1 from the list of possibilities, giving that there are 5 primes in the list, confirming our assumption.

So I don't know if this always works, I'll write a program to try and generate the list of the first 1000 primes or something in this way so expect an update...

new drawings

Smooth curve generator implemented in python

Given a list of ordered points as in the bottom of this image, my program produces the top part of the image:

For more about the theory behind it see my last blog post.
The run time for the program was on the order of a hundredth of a second but I was plotting the pixels one by one I think there are faster ways to plot the output.

The code;s not commented extensively but short and straightforward. if there are any questions let me know... it requires the Image module from PIL to make the image shown above (the only difference is I made the dots bigger in the image above with GIMP to make them stand out....

import Image
import math
import time

def plot(curve, blank):
    for point in curve:
        blank.putpixel(point, (255,255,255))
def bplot(curve, blank):
    for point in curve:
        blank.putpixel((point[0], point[1]+250), (255, 0, 0))
def circ(point1, point2, point3):
    ax = point1[0]
    ay = point1[1]
    ax = point1[0]
    ay = point1[1]
    bx = point2[0]
    by = point2[1]
    cx = point3[0]
    cy = point3[1]
    '''find the x,y and radius for the circle through the 3 points'''
    if (ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay) != 0:
        x=.5*(-pow(ay, 2)*cy+pow(ay, 2)*by-ay*pow(bx, 2)\
        -ay*pow(by, 2)+ay*pow(cy, 2)+ay*pow(cx, 2)-\
        pow(cx, 2)*by+pow(ax, 2)*by+pow(bx, 2)*\
        cy-pow(ax, 2)*cy-pow(cy, 2)*by+cy*pow(by, 2))\
        /(ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay)
        y=-.5*(-pow(ax, 2)*cx+pow(ax, 2)*bx-ax*pow(by, 2)\
        -ax*pow(bx, 2)+ax*pow(cx, 2)+ax*pow(cy, 2)-\
        pow(cy, 2)*bx+pow(ay, 2)*bx+pow(by, 2)*cx\
        -pow(ay, 2)*cx-pow(cx, 2)*bx+cx*pow(bx, 2))\
        /(ax*by-ax*cy-cx*by+cy*bx-bx*ay+cx*ay)
    else:
        return False
    
    r=pow(pow(x-ax, 2)+pow(y-ay, 2), .5)
    
    return x, y, r
def findpoint(eq1, eq2, point1, point2):
    '''find the centroid of the overlapping part of two circles
    from their equations'''
    thetabeg = math.acos((point1[0]-eq1[0])/eq1[2])
    thetaend = math.acos((point2[0]-eq1[0])/eq1[2])
    mid1x = eq1[2]*math.cos((thetabeg+thetaend)/2)+eq1[0]
    thetaybeg = math.asin((point1[1]-eq1[1])/eq1[2])
    thetayend = math.asin((point2[1]-eq1[1])/eq1[2])
    mid1y = eq1[2]*math.sin((thetaybeg+thetayend)/2)+eq1[1]

    thetabeg2 = math.acos((point1[0]-eq2[0])/eq2[2])
    thetaend2 = math.acos((point2[0]-eq2[0])/eq2[2])
    mid2x = eq2[2]*math.cos((thetabeg2+thetaend2)/2)+eq2[0]
    thetaybeg2 = math.asin((point1[1]-eq2[1])/eq2[2])
    thetayend2 = math.asin((point2[1]-eq2[1])/eq2[2])
    mid2y = eq2[2]*math.sin((thetaybeg2+thetayend2)/2)+eq2[1]
    return [(mid2x+mid1x)/2, (mid2y+mid1y)/2]
def main():
    starttime = time.time()
    ''' Ordered list of points, the first and last affect the shape
    of the curve but are not connected though drawn'''
    curve = [[23,150],[43,127],[100,49],[179,83],[253,155],[362,215],[445,122],[551,57],[560,23], [540, 10], [520,20],[500,50],[450,70]]
    blank = Image.new("RGB", (600, 600))
    bplot(curve, blank)
    
    for j in range(0, 6):
        newpoints = []
        for i in range(0, len(curve)-3):
            eq = circ(curve[i], curve[i+1], curve[i+2])
            eq2 = circ(curve[i+1], curve[i+2], curve[i+3])
            if eq == False or eq2 == False:
                newpoints.append([(curve[i+1][0]+curve[i+2][0])/2, (curve[i+1][1]+curve[i+2][1])/2])
            else:    
                newpoints.append(findpoint(eq, eq2, curve[i+1], curve[i+2]))
        for point in newpoints:
            point[0] = int(round(point[0]))
            point[1] = int(round(point[1]))
        for m in range(0, len(newpoints)):
            curve.insert(2*m+2, newpoints[m])
        
    plot(curve, blank)
    print time.time()-starttime
    blank.save("smoothcurve.png")
    
main()

**The above code runs on Windows 7 Home, Python 2.7.3, PIL 32 bit

I think it's nice to be able to do like the following, usually selection in image programs follows polygonal lines but with this algorithm you can do this instead:

And it's fast enough that you could make adjustments to add new nodes or move individual nodes and it recalculate in real time, which I think would be handy....
***LOL I know it looks like I just hand drew that outline around the cloud but I really gave it the coordinates of the red dots and it did all the rest
**Comment from Ofnuts from the gimp mailing list:
Ofnuts: The curve doesn't include the two end points and it seems it doesn't even get through one of the points (next before last at top right):

Me: Right I should have made it more clear, I left it as a comment in the code, but the first and last point figure into the calculation and do have some effect but are not directly connected to the rest of the curve. They could be left undrawn or automatically be selected as the top left of the bounding rectangle and the bottom right or as I was saying put anywhere to give the beginning and end of the curve the "feel" that it was coming from a certain direction before and after  the points that actually are all connected.

The second issue is that I was placing the red dots in gimp and reading off their coordinates also with gimp to put in the program, so it wasn't exact, but the coordinates actually put in the program are all plotted as pixels so in the end are drawn as part of the curve... I just checked and made sure 520,20 given in the program is where it should be, the red dots were wrong...

Constructing a smooth curve from a set of points geometrically

Suppose you start with a set of points and you know the order they should be connected...

The first step of this method is to construct arcs between each set of 3 consecutive points in the geometrical sense.  Here I've done the first 3 points...

Doing each set of 3 consecutive points gives...

Now the next step is to place new points in the centroids of any enclosed areas like so...

Now all the arcs currently drawn can be deleted...

And the process reiterates... more arcs are drawn for every 3 consecutive points...

And then once again a point is placed for each centroid, and more arcs are drawn, here is where we are after 3 iterations...

As you can see the overall path becomes smoother and smoother. 

It can go on forever or a certain number of iterations. 

**EDIT**

Ofnuts from the gimp forum asked me a good question about how many iterations would be necessary... On a computer screen with pixels you would need to go log_2(pixel arc length) iterations.  The pixel arc length is the number of pixels in the final curve. An upper bound for the pixel arc length is the pixel length of the circular arcs placed in the first iteration.This is because the number of pixels you have placed doubles each iteration and the final curve fits inside of the first iterations arcs.