RandomComplex.py

## Generate Random Complex K (for circle packing algorithm)
import random
import bpy
## 1rst Generating Random polygon 
##   
## Starting by generating point maxima and minima.
##  In this method, for n>=3, iteratively start by generating 3 points.

## Iterated point procedures randomly chooses an edge, then with minima,
## and maxima constraints picks a point randomly between edge minima and
## maxima, forming a new edge.
## A limit can be placed on subdivisions such that no further subdivision
## of any newly formed edge can take place without the subdivision
## of any older higher Queue priority edges in the randomization process.

MaxSize = 10
PolygonSize = 20  ## must be 3 or higher

def clockwisewalktest(walk):
    ## works with nonconvex polygons should be safe
    ## I believe for the primitive polygon type (3 vertices)
    ## constructed in this algorithm.
    prev = None
    samt = 0.0
    newwalk = walk[0:len(walk)]
    newwalk.append(walk[0])
    for vert in newwalk:
        if prev == None:
            prev = vert
            continue
##        samt += (vert[0]-prev[0])*(vert[1]+prev[1])
        samt += (prev[0]*vert[1]-prev[1]*vert[0])
##        if walk.index(vert) == len(walk)-1:
##            samt += (walk[0][0] - vert[0])*(walk[0][1]+vert[1])
        prev = vert
    print('samt: ', samt)
    if samt < 0:
        print('original walk is clockwise')
        return True
    else:
        print('original walk is counter clockwise')
        return False

def polygonwalk(vert,last,target,vedges, walk):
    for edge in vedges[vert]:
        va, vb = edge
        prev = vert
        nextv = None
        if va == vert:
            if vb != last:
                walk.append(vb)
                if vb != target:
                    nextv = vb
                    polygonwalk(nextv,prev,target,vedges,walk)

        else:
            if va != last:
                walk.append(va)
                if va != target:
                    nextv = va
                    polygonwalk(nextv,prev,target,vedges,walk)
def convextest(v1,v2,v3,v4):
    def crossproduct(p1,p2):
        return (p1[0]*p2[1] - p1[1]*p2[0])
    ## assumed v1,v2,v3,v4 are sequentially ordered on the
    ##polygon walk
    ## this is cross product comparison
    t1 = crossproduct(v1,v2) == crossproduct(v2,v3)
    t2 = crossproduct(v3,v4) == crossproduct(v1,v2)
    if t1 and t2:
        return True
    else:
        return False
    

def generateRandomVertex():
    return (random.random()*MaxSize, random.random()*MaxSize)

def generateRandomVertexMM(minX,maxX,minY,maxY):
    return (random.uniform(minX,maxX),random.uniform(minY,maxY))

def cubicInterpolate (p, x):
    return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])))

def distance(a,b):
    ax,ay = a
    bx,by = b
    return ((ax-bx)**2+(ay-by)**2)**.5

def getMinMax(edge):
    edgex = edge[0:len(edge)]
    edgey = edge[0:len(edge)]
    edgex = list(edgex)
    edgey = list(edgey)
    edgex.sort(key = lambda tup:tup[0])
    edgey.sort(key = lambda tup:tup[1])
    minx = edgex[0][0]
    maxx = edgex[1][0]
    miny = edgey[0][1]
    maxy = edgey[1][1]
    return (minx,maxx,miny,maxy)

def getXScale(minx,maxx):
    return 1.0/abs(maxx-minx)

def scale(scale, points):
    rpoints = []
    for point in points:
        rpoints.append((scale*point[0], scale*point[1]))
    return rpoints

def translateX(tr, points):
    rpoints = []
    for point in points:
        rpoints.append((tr + point[0], point[1]))
    return rpoints

def slope(edge):
    a, b = edge
    ax,ay = a
    bx,by = b
    return (by - ay)/(bx - ax)

def midpoint(edge):
    a,b = edge
    return ((a[0]+b[0])/2,(a[1]+b[1])/2)

def slopenormal(edgeslope):
    return - 1/edgeslope

def testdirection(edge1, edge2):
    #assumed edge1 = (a,b) and edge2 = (b,c)
    # where b intersect edge 1 and 2
    a,b = edge1
    b,c = edge2
    if c > b:
        if a > b:
            return False
        else:
            return True
    else:
        if a > b:
            return True
        else:
            return False

def setRotation(edge, rotheir):
    ## closest distance to walk pair root determines
    ## direction of the vector
    root = rotheir[edge]
    ra, rb = root
    rbx, rby = rb
    rax, ray = ra
    vec = [rbx-rax, rby-ray]
##    a,b = edge
##    ax,ay = a
##    bx,by = b
##    ## find which vertex is closest to root a
##    dara = distance(a,ra)
##    dbra = distance(b,ra)
##    if dara < dbra:
##    vec = [bx-ax, by-ay]
##    else:
##    vec = [ax-bx, ay-by]
    ## 90 degree rotation
    print('ab vector: ', vec)
    vec = [-vec[1], vec[0]]
    print('rotation edge: ', edge)
    print('rotation vector: ', vec)
    return vec

def getY(point, slope, x):
    return slope*(x - point[0]) + point[1]

def getneighborverts(edge,vedges):
    a,b = edge
    n1 = None
    n2 = None
    ne1 = None
    ne2 = None
    for nedge in vedges[a]:
        if edge != nedge:
            na,nb = nedge
            ne1 = nedge
            if na == a:
                n1 = nb
            else:
                n1 = na
    for nedge in vedges[b]:
        if edge != nedge:
            na,nb = nedge
            ne2 = nedge
            if na == b:
                n2 = nb
            else:
                n2 = na
    return n1,n2,ne1,ne2
    

def updateEdges(a,b,edges,dedge,vedges):
    edges.append((a,b))
    d = distance(a,b)
    if d in dedge:
        dedge[d].append((a,b))
    else:
        dedge[d] = [(a,b)]
    if a in vedges:
        vedges[a].append((a,b))
    else:
        vedges[a] = [(a,b)]
    if b in vedges:
        vedges[b].append((a,b))
    else:
        vedges[b] = [(a,b)]

def updateRotatheir(edge, parent, rotheir):
    root = rotheir[parent]
    ra, rb = root
    a,b = edge
    ax,ay = a
    bx,by = b
    ## find which vertex is closest to root a
    dara = distance(a,ra)
    dbra = distance(b,ra)
    if dara < dbra:
        rotheir[edge] = (a,b)
    else:
        rotheir[edge] = (b,a)

def deleteEdge(edge,edges,dedge,vedges):
    a,b = edge
    d = distance(a,b)
    edges.remove(edge)
    dedge[d].remove(edge)
    if len(dedge[d]) == 0:
        del dedge[d]
##    print(edge)
##    print(vedges)
    vedges[a].remove(edge)
    vedges[b].remove(edge)

edgecount = 0
vertices = []
edges = []
dedge = {}
vedges = {}
edged = {}
rotheir = {}
for i in range(0,3):
    vertices.append(generateRandomVertex())
    if len(vertices) > 0 and len(vertices) != 1:
        a = vertices[0]
        b = vertices[len(vertices)-1]
        updateEdges(a,b,edges,dedge,vedges)
        rotheir[(a,b)] = (a,b)

    edgecount += 1
a = vertices[2]
b = vertices[1]
updateEdges(a,b,edges,dedge,vedges)
rotheir[(a,b)] = (a,b)
edgecount += 1

## create walk order used in determining rotations
verts = list(vedges.keys())
print(vedges)
a = verts[0]
print('a',a)
tedge = vedges[a][0]
last = None
target = None
walk = [a]
print('walk:', walk)
for vert in tedge:
    if vert != a:
        target = vert
        last = vert
polygonwalk(a,last,target,vedges,walk)

print('original walk: ', walk)
newwalk = []
if not clockwisewalktest(walk):
    walk0 = walk[0]
    walk1 = walk[1:len(walk)]
    walk1 = walk1[::-1]  ## reverse the order for clockwise
    newwalk.append(walk0)
    newwalk += walk1
    walk = newwalk
    
print('rotation order walk: ', walk)

def setrotatheirorder(walk, rotheir):
    prev = None
    for vert in walk:
        if prev == None:
            prev = vert
            continue
        print('vert,prev pair: ', (vert,prev))
        if (vert, prev) in rotheir:
            print('found rev order rotheir key')
            print('vert,prev pair: ', (vert,prev))
            rotheir[(vert,prev)] = (prev,vert)
        if walk.index(vert) == len(walk)-1:
            if (walk[0], vert) in rotheir:
                rotheir[(walk[0],vert)] = (vert,walk[0])
        prev = vert

setrotatheirorder(walk,rotheir)
print('rotheir: ', rotheir)

Q = []
edgec = 0
print(dedge)
qedges = None
pedge = None
parents = []
i = 0
while (edgecount < PolygonSize+1):
    if len(Q) == 0:
       ##fill Q
       ##edgescopy = edges[0:len(edges)]
       ##random.shuffle(edgescopy)
       edgekeys = list(dedge.keys())
       edgekeys.sort(reverse = True)
       Q = edgekeys
       i += 1
    if edgec == 0:
        qedges = dedge[Q[0]][0:len(dedge[Q[0]])]   
        pedge = dedge[Q[0]][edgec]
    else:
        pedge = qedges[edgec]
    minx,maxx,miny,maxy = getMinMax(pedge)
    ##nvert = generateRandomVertexMM(minx,maxx,miny,maxy)
    x = (maxx+minx)/2.0
    n1,n2,ne1,ne2 = getneighborverts(pedge,vedges)
##    if testdirection(ne1,pedge) and testdirection(pedge,ne2):
    if 1 == 0:
        a,b = pedge
        xscale = getXScale(minx,maxx)
    ##    print('n1:',n1)
    ##    print('n2:', n2)
    ##    print('ne1:', ne1)
    ##    print('ne2:', ne2)

        p = [n1,a,b,n2]
        ## we need to set up interpolation which means scaling
        ## and translating positions to end up on interval [0,1]
        ## for p1 and p2, will also need to compute position x =-1
        ## for p0 at y.
        p = scale(xscale, p)
        tr = -minx*xscale
        sxt = xscale*x+tr
        p = translateX(tr, p)
        if p[1][0] != 0.0:
            p = [p[3],p[2],p[1],p[0]]
        ne1 = (p[0],p[1])
        ne2 = (p[2],p[3])
        sne1 = slope(ne1)
        sne2 = slope(ne2)
        ny1 = getY(p[0], sne1, -1)
        ny2 = getY(p[3], sne2, 2)
        print('ny1:', ny1)
        print('ny2:', ny2)
        print('p:',p)
        py = [ny1,p[1][1],p[2][1],ny2]
        syt = cubicInterpolate (py, sxt)
        ##rescale y back to original coordinate
        ## note: we don't worry about retranslating since this isn't an xcoordinate
        y = syt*1/xscale
    else:
        mpoint = midpoint(pedge)
##        pslope = slope(pedge)
##        nslope = slopenormal(pslope)
        rvec = setRotation(pedge, rotheir)
##        y = getY(pedge[0], pslope, x)
        x, y = mpoint
        print('y at midpoint: ', y)
        x = x + rvec[0]/(3.5) ##+ i)
        y = y + rvec[1]/(3.5 )##+ i)
    vertices.append((x,y))
    nvert = (x,y)
    updateEdges(pedge[0],nvert,edges,dedge,vedges)
    updateEdges(pedge[1],nvert,edges,dedge,vedges)

    deleteEdge(pedge,edges,dedge,vedges)
    nedge1 = (pedge[0],nvert)
    nedge2 = (pedge[1],nvert)
    updateRotatheir(nedge1, pedge, rotheir)
    updateRotatheir(nedge2, pedge, rotheir)
    del rotheir[pedge]
    edgecount += 1
    print('Q[0]', Q[0])
    print('qedges', qedges)
    print('length qedges - 1: ',len(qedges)-1)
    if edgec == len(qedges)-1:
        del Q[0]
        edgec = 0
    else:    
        edgec += 1




verts = list(vedges.keys())
print(vedges)
a = verts[0]
print('a',a)
tedge = vedges[a][0]
last = None
target = None
walk = [a]
print('walk:', walk)
for vert in tedge:
    if vert != a:
        target = vert
        last = vert
polygonwalk(a,last,target,vedges,walk)
face = []
for vert in walk:
    face.append(vertices.index(vert))
bvertices = []
for vert in vertices:
    x,y = vert
    bvertices.append((x,y,0.0))
faces = []
faces.append(tuple(face))
meshName = "Polygon"
obName = "PolygonObj"
me = bpy.data.meshes.new(meshName)
ob = bpy.data.objects.new(obName, me)
ob.location = bpy.context.scene.cursor_location
bpy.context.scene.objects.link(ob)
me.from_pydata(bvertices,[],faces)      
me.update(calc_edges=True)