analytics

Monday, December 29, 2014

Shortnet definition and source code

I'm defining the Shortnet of a group of points to be the subset S of the collection of every possible line connecting pairs of points C, such that every line crossing a line shorter than itself has been removed.
Example: You have your points...
 Then you create all the lines connecting every point pairwise...
Then you remove any line that crosses a line shorter than itself...
Here is one with 50 points:

There will never be any lines that cross left as one or the other of each pair of lines that cross is removed.
Source Code using PIL drawing library and Python:
import random
import Image, ImageDraw
def cross(a1, b1, c1, d1, a2, b2, c2, d2):
    if (-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)== 0:
        return False
    if (-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1) == 0:
        return False

    s = 1.0*(a1*d1-d2*a1+d2*c1-c1*b1+c2*b1-c2*d1)/(-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)
    t = 1.0*(-b2*c1+a2*d1+b2*c2-a2*d2+d2*c1-c2*d1)/(-a2*b1+a2*d1+c2*b1-c2*d1+b2*a1-b2*c1-d2*a1+d2*c1)
    if s > 0 and s < 1.0:
        if t > 0 and t < 1.0:
            return True
    return False
def gen(num):
    points = []
    for i in range(0, num):
        point = [random.randint(100, 700), random.randint(100, 700)]
        points.append(point)
    return points
def shorter(a1, b1, c1, d1, a2, b2, c2, d2):
    d1 = ((c1 - a1)**2.0 + (d1 - b1)**2.0)**.5
    d2 = ((c2 - a2)**2.0 + (d2 - b2)**2.0)**.5
    if d1 <= d2:
        return True
    else:
        return False
def makelines(points):
    lines = []
    for i in range(0, len(points)-1):
        for j in range(i+1, len(points)):
            lines.append([points[i][0], points[i][1], points[j][0], points[j][1]])
    return lines
def picklines(lines):
    marked = []
    for i in range(0, len(lines)):
        marked.append(True)
    for i in range(0, len(lines)-1):
        for j in range(i+1, len(lines)):
            if cross(lines[i][0], lines[i][1], lines[i][2], lines[i][3], lines[j][0], lines[j][1], lines[j][2], lines[j][3]):
                if shorter(lines[i][0], lines[i][1], lines[i][2], lines[i][3], lines[j][0], lines[j][1], lines[j][2], lines[j][3]):
                    marked[j] = False
                else:
                    marked[i] = False
    for i in range(0, len(lines)):
        if marked[i] == False:
            lines[i] = False
    return lines
def drawlines(d, lines):
    for line in lines:
        if line != False:
            d.line(line)
def drawpoints(d, points):
    for point in points:
        d.ellipse([point[0]-10, point[1]-10, point[0]+10, point[1]+10], fill=(255,0,0))
def main():
    points = gen(10)
    plot = Image.new("RGB", [800, 800])
    d = ImageDraw.Draw(plot)
    drawpoints(d, points)
    plot.save("dots.png")
    lines = makelines(points)
    drawlines(d, lines)
    plot.save("all.png")
    plot = Image.new("RGB", [800, 800])
    d = ImageDraw.Draw(plot)
    drawpoints(d, points)
    lines = picklines(lines)
    drawlines(d, lines)
    plot.save("shortnet.png")
main()

Friday, December 19, 2014

Parallelizable convex hull algorithm

Suppose you have a number of points in any number of dimensions, here I will just show 2.
First find the centroid, and the point farthest from that centroid Pc and construct the circle or n-dimensional sphere to it from the centroid... (Centroid shown in blue)
Now consider all the straight lines C to P(i) to B(i) from the Centroid to each Point and continuing on to a point on the Bounding region. Make a list of Pc and all points P(i) such that P(i) to B(i) is shorter than any other line from a different point P(j) to B(i).
This last step is the part that can be done in parallel, passing the entire set to as many processors as there are points and each one finding the solution for its given point. The points that are closest to their respective bounding point than any other point is to that bounding point are the hull. 

Tuesday, December 9, 2014

recursive formula for points on the unit circle of regular 2^n gon's

I found that formulas of this type:


describe the points on a circle corresponding to regular 2^(n+1)-gons where n is the number of square root symbols in the formulas and a,b,c... are each 1 or 2. For example a 16-gon
I haven't quite figured out where each point will be on the unit circle by the signs of the terms though. 

Monday, December 8, 2014

Fast parallelizable way to generate approximate points for a circle without using trig functions

This method starts with the four points of the bounding square...
Then replaces each point with two points (1- root 2/2) away from the sides... This is basically like truncating every corner to the right proportion...
Then iterates this process but using .25 instead of (1-root 2/2) that I found by experiment works very well. I suppose for perfect accuracy there is a recursive series involving root 2 but it approaches .25 so quickly that I didn't bother with it.


I found 4 iterations is enough for fairly large circles...
And in the end you have an list with all the points in counter clockwise order around the circle. 
Python source code (Using PIL for drawing)
import Image, ImageDraw
def plot(points, draw):
    for i in range(0, len(points)-1):
        draw.line([points[i], points[i+1]])
    draw.line([points[len(points)-1], points[0]])
def truncate(points):
    newpoints = []
    param = .29
    if len(points)>4:
        param = .25
    for i in range(0, len(points)-1):
        x1 = points[i+1][0]*(param)+points[i][0]*(1.0-param)
        y1 = points[i+1][1]*(param)+points[i][1]*(1.0-param)
        x2 = points[i+1][0]*(1.0-param)+points[i][0]*(param)
        y2 = points[i+1][1]*(1.0-param)+points[i][1]*(param)
        newpoints.append((x1, y1))
        newpoints.append((x2, y2))
    x1 = points[0][0]*(param)+points[len(points)-1][0]*(1.0-param)
    y1 = points[0][1]*(param)+points[len(points)-1][1]*(1.0-param)
    x2 = points[0][0]*(1.0-param)+points[len(points)-1][0]*(param)
    y2 = points[0][1]*(1.0-param)+points[len(points)-1][1]*(param)
    newpoints.append((x1, y1))
    newpoints.append((x2, y2))
    return newpoints
def main():

    points = [(100,100),(100,700),(700,700),(700,100)]
    
    for i in range(1, 5):
        p = Image.new("RGB", [800,800])
        draw = ImageDraw.Draw(p)
        points = truncate(points)  
        plot(points, draw)
        p.save("plots"+str(i)+".png")
main()

Thursday, December 4, 2014

Using physics to try and find a Hamilton Circuit of a graph with odd number of vertices

Supposing you have a simply connected planar graph of even number of vertices with no loops edges labeled like so:
I found one can set one edge leading to one vertex, say A arbitrarily to 1, then write equations from Kirtchoff's laws for the graph...
http://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws

Kirchoff's laws say that the sum of the edges connected to each vertex is 0 and the sum around each face is 0.

It's kind of interesting that after you set one edge to 1, this will always be the right number of equations because of the graph equation V+F=E+2, and we have V+F equations to solve for E unknowns, The plus two in this formula is canceled by the fact that usually the formula counts the region outside of the graph as a face, but we don't use write an equation for that, and we've set one edge value to 1 so there is one less unknown edge.

Note there will be no solution if c,g,l, or m in this example is the edge chosen to be 1. But for any the other edges it will either be this solution or -1*all the values.

Then we can plot these solutions...
Now we can follow the flow of the graph around through every vertex exactly once by starting with A which is 1, then following the path from that vertex that is -1, and from there to the vertex following the path of 1, and so on alternating values. Since there are an odd number of vertices this works out evenly.
If you had an odd number of vertices to start out with, you could break the edge that you set as 1 into two edges one -1 and 1...

**Reasoning**
I don't know how to prove it, it just seems to work in every case. It makes sense to me electrically that if there is a nonzero flow around each face's edges, that some faces will combine and form one flow around the outside. I don't know how to prove that that will always go through every vertex though.

Friday, November 28, 2014

Branching poisson disk algorithm

My program outputs a distribution like this:
Every dot is at least a distance of 10 pixels from any other.
The way it works is to start at the middle, and branch out in 3 evenly spaced directions at a random angle from the center dot, then from those 3 points branch in a similar way generating more points and continue in that fashion for a certain number of iterations, above is 25 generations. I think this picture helps explain it a lot.
The speed is fairly quick, generating these graphs in python took only about 1.5 seconds for around 1300 points, I think it could be speeded up considerally using parallel processing for example on a gpu with a faster language. 
**Python source code**
import math
import random
import Image, ImageDraw
def check(d, pointlist, point, gridsize):
    tooclose = False
    if point[0] < 10 or point[0] > gridsize[0]-10 or point[1] < 10 or point[1] > gridsize[1]-10:
        return True
    for p in range(0, len(pointlist)):
        if ((point[0]-pointlist[p][0])**2.0 + (point[1]-pointlist[p][1])**2.0)**.5 < d-1:
            return True
    return tooclose
def findspot(pointlist, added, d, gridsize, draw):
    newlist = []
    newadded =[]
    for point in pointlist:
        newlist.append(point)
    for point in added:
        tx = point[0]
        ty = point[1]
        theta = random.random()*2*3.14159
        for i in range(0, 3):
            px = tx + d*math.cos(theta+i*(2*3.14159)/3)
            py = ty + d*math.sin(theta+i*(2*3.14159)/3)
            if check(d, newlist, [px,py], gridsize) == False:
                newlist.append([px,py])
                newadded.append([px,py])
    return newlist, newadded
class graph():
    def __init__(this, d, gridsize, draw):
        this.pointlist = [[400,400]]
        this.added = [[400,400]]
        this.d = d
        this.gridsize = gridsize
        this.draw = draw
    def tree(this, limit):
        newlist, newadded = findspot(this.pointlist,this.added,this.d, this.gridsize, this.draw)
        if limit < 25:
            this.pointlist = newlist
            this.added = newadded
            this.tree(limit+1)
def main():
    gridsize = [790,790]
    d = 10.0
    point = [400,400]
    plot = Image.new("RGB", [800,800])
    draw = ImageDraw.Draw(plot)
    G = graph(d, gridsize, draw)
    G.tree(0)
    for point in G.pointlist:
        for i in range(-1, 2):
            for j in range(-1, 2):
                plot.putpixel((int(point[0]+i), int(point[1])+j), (255,255,255))
    plot.save("plotdot.png")
main()

Sunday, November 23, 2014

recursively moving around the unit circle with constant angle to the origin

Starting with solving this general formula for two vectors [x,y] and [a,b]:
The above says the dot product between them equals to c and the magnitude of the cross product equals to d (counterclockwise). And then it is solved for a,b.
So we begin by thinking about if [x,y] starts as the vector [1,0]... then we can find possible values for c and d like so, picking any value as close or far away from 1 as we'd like for x and solving...
So plug those in for d and c respectively...Or alternatively use a rational value from the rational parameterization of the circle. These are the values for the dot and cross product from the vector to the generated point to the vector [1, 0] because the dot product will just be the x coordinate of the point we generated and the cross product will just be the y.   Let's call this Formula 1.
Now c and d is also our first vector after [1,0] so we can plug that in for x and y as well and get...
This vector [a,b] is the same angle to [.9801,.141] as [.9801,.141] was to [1,0], because we made the dot and cross products stay the same.

Now we can plug this [a,b] in for x and y in formula 1 and get...
And in this way work our way around the unit circle by that small angle (which we actually don't know) without ever having to figure a cosine or sine! An approximation of the angle can be found by seeing how long it takes for this process to get close to a known value like [0,1]