ProfileGeneration.py

About

This module primarily houses the ProfileGeneration() class which manages everything for the generation of a single geometry profile. Geometry profiles can be expressed in the following forms,

  • Cartesian coordinate list
  • polar coordinate list
  • Binary encoded string
  • Image array

Radial Parameterization

In the context of the downstream genetic algorithm these geometry profiles describe different realization of parameter space. An example profile geometry is shown below:

drawing

each red dot indicates an individual element of a profile in polar coordinates, represented as an ordered pair, \((\theta_n,r_n)\) where \(r_n\) is the point along the line cast from the center to the edge of the domain that intersects with the polygon, this realization of the polygon is refered to as the "radial parameterization" later on

Profile Generators

there are five profile generator functions implemented, with plans to extend merge rectangleGenerator() and triangleGenerator() into a generalized generator for n-sided polygons.

Generation of circles is by far the simplest process every radial term is set equal resulting in a uniform circle

def circleGenerator(self):
    self.t = np.linspace(0,2*np.pi, self.nVertices)
    self.fourierGenerator(1)
    self.r = np.random.randint(self.rMin, self.rMax)*(self.r / np.max(self.r))

as detailed in the previous section radial parameterization of arbitrary polygons is accomplished by finding the point where a central ray intersects with polygon edges. For rectangles this process has a special simplified case seen in the below code segment.

def rectangleGenerator(self, spanMin = 40, spanMax = 135):
    rectW = np.random.randint(spanMin, spanMax)
    rectH = np.random.randint(spanMin, spanMax)
    self.t = np.linspace(0,2*np.pi, self.nVertices)
    self.r = np.zeros(self.nVertices)
    #calculate rays cast far from the origin
    x, y = 240*np.cos(self.t), 240*np.sin(self.t)
    for v in range(self.nVertices):
        if np.abs(y[v]/x[v]) > rectH/rectW:
            #The far ray intersects with the top or bottom of the rectangle
            self.r[v] = (rectH/2)/np.abs(np.sin(self.t[v]))
        else:
            #The far ray intersects with the side of the rectangle
            self.r[v] = (rectW/2)/np.abs(np.cos(self.t[v]))

In the beginning of this section it was mentioned that a n-sided polygon generator is planned to replace rectangleGenerator() and triangleGenerator(), the blueprint of this would follow the triangleGenerator() code below, so it's a good idea to break down this code

  def triangleGenerator(self):
    self.t = np.linspace(0,2*np.pi, self.nVertices)
    self.r = np.zeros(self.nVertices)

    #places three random vertices within the domain
    v0 = np.random.randint(-70,70, 2); x1, y1 = v0
    v1 = np.random.randint(-70,70, 2); x2, y2 = v1
    v2 = np.random.randint(-70,70, 2); x3, y3 = v2
    area = 0.5 * np.abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
    while area <= 4*np.pi*self.rMin*self.rMin:
        v0 = np.random.randint(-70,70, 2); x1, y1 = v0
        v1 = np.random.randint(-70,70, 2); x2, y2 = v1
        v2 = np.random.randint(-70,70, 2); x3, y3 = v2
        area = 0.5 * np.abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
    #find the center of the triangle
    xC = np.mean([v0[0], v1[0], v2[0]]); yC = np.mean([v0[1], v1[1], v2[1]])
    #shift each vertice so that the triangle is centered at the origin
    vC0 = (v0[0]-xC, v0[1]-yC); vC1 = (v1[0]-xC, v1[1]-yC); vC2 = (v2[0]-xC, v2[1]-yC)

    for v in range(self.nVertices):
        ray = LineString([(0,0), (180*np.cos(self.t[v]), 180*np.sin(self.t[v]))])
        edge01 = LineString([(vC0[0], vC0[1]),(vC1[0], vC1[1])])
        edge12 = LineString([(vC1[0], vC1[1]),(vC2[0], vC2[1])])
        edge20 = LineString([(vC2[0], vC2[1]),(vC0[0], vC0[1])])

        edge01Int = ray.intersection(edge01)
        edge12Int = ray.intersection(edge12)
        edge20Int = ray.intersection(edge20)

        if np.shape(edge01Int.coords)[0] != 0:
            x, y = edge01Int.coords[0][0], edge01Int.coords[0][1]
        elif np.shape(edge12Int.coords)[0] != 0:
            x, y = edge12Int.coords[0][0], edge12Int.coords[0][1]
        elif np.shape(edge20Int.coords)[0] != 0:
            x, y = edge20Int.coords[0][0], edge20Int.coords[0][1]
        else:
            x, y = 0, 0

        self.r[v] = np.sqrt(x**2 + y**2)

Three random vertices are chosen and the area of that triangle is calculated and if the area is below a set limit (four times the area of the minimum circle seems to produce results that are reasonable), then the vertices are re-selected until that condition is passed and then the vertices are shifted so the center of mass of the triangle is centered at the origin

v0 = np.random.randint(-70,70, 2); x1, y1 = v0
v1 = np.random.randint(-70,70, 2); x2, y2 = v1
v2 = np.random.randint(-70,70, 2); x3, y3 = v2

area = 0.5 * np.abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
while area <= 4*np.pi*self.rMin*self.rMin:
    v0 = np.random.randint(-70,70, 2); x1, y1 = v0
    v1 = np.random.randint(-70,70, 2); x2, y2 = v1
    v2 = np.random.randint(-70,70, 2); x3, y3 = v2
    area = 0.5 * np.abs(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))

xC = np.mean([v0[0], v1[0], v2[0]]); yC = np.mean([v0[1], v1[1], v2[1]])
#shift each vertice so that the triangle is centered at the origin
vC0 = (v0[0]-xC, v0[1]-yC); vC1 = (v1[0]-xC, v1[1]-yC); vC2 = (v2[0]-xC, v2[1]-yC)

for each polar angle associated with a radial parameter (generally refered to a vertex but that language is a bit confusing in this context) a line segment from 0 along the polar angle to the edge of the domain is generated. 

ray = LineString([(0,0), (180*np.cos(self.t[v]), 180*np.sin(self.t[v]))])

an the edeges of the desired triangle are drawn from it's vertices

edge01 = LineString([(vC0[0], vC0[1]),(vC1[0], vC1[1])])
edge12 = LineString([(vC1[0], vC1[1]),(vC2[0], vC2[1])])
edge20 = LineString([(vC2[0], vC2[1]),(vC0[0], vC0[1])])

using LineString() from the Shapely python library the each edge is checked for an intersection with the test line segment and the distance from the origin to that coordinate is saved as the radial value. 

edge01Int = ray.intersection(edge01)
edge12Int = ray.intersection(edge12)
edge20Int = ray.intersection(edge20)

if np.shape(edge01Int.coords)[0] != 0:
    x, y = edge01Int.coords[0][0], edge01Int.coords[0][1]
elif np.shape(edge12Int.coords)[0] != 0:
    x, y = edge12Int.coords[0][0], edge12Int.coords[0][1]
elif np.shape(edge20Int.coords)[0] != 0:
    x, y = edge20Int.coords[0][0], edge20Int.coords[0][1]
else:
    x, y = 0, 0

self.r[v] = np.sqrt(x**2 + y**2)

This geometric analysis approach works very well and is reliable in generating an assortment of triangles, but the following generators take a less controlled approach

def polygonGenerator(self):
    self.t = np.linspace(0,2*np.pi, self.nVertices)
    self.r = np.random.randint(self.rMin, self.rMax, self.nVertices)
    self.r = self.r*1.0
    self.r[-1] = self.r[0]

in this method each radial value is assigned a random number in the allowed range, generating widely varied polygons, which is ideal for initializing a chromosome population where high diversity is diserable.

def fourierGenerator(self, nTerms):
    self.r = np.zeros(self.nVertices)
    self.t = np.linspace(0,2*np.pi, self.nVertices)
    #setup random fourier series terms
    c = np.random.rand(nTerms)
    d = np.random.rand(nTerms)
    #fourier series terms should sum to rMax for scaling purposes
    c = self.rMax*(c / sum(c))
    d = self.rMax*(d / sum(d))
    for v in range(self.nVertices):
        for n in range(nTerms):
            self.r[v] += c[n]*np.sin(n*self.t[v])+d[n]*np.cos(n*self.t[v])

    #restrict r to the range [rMin, rMax]
    self.r = self.r + (self.r <= self.rMin)*(self.rMin-self.r)
    self.r = self.r + (self.r >= self.rMax)*(self.rMax-self.r)

an alternative method for polygon generation is to build a fourier series with random weights such that \(R(0) = R(2\pi)\) the structure of these shapes is highly dependent on the number of terms, if only one term is used the same circle will always be drawn and if too many terms are used (more than 5) the shape generated isn't quite useable

drawing

In the above image the affect of increasing terms is demonstrated more than five terms will generally results in profiles that are quite short and increasingly wide along the right hand side.

Binary Representation

in the genetic operations expect profiles to be in a binary representation, that conversion is handled in this module with encodePolygon() and the binary representation can be decoded back to the radial form with decodePolygon() 

def encodePolygon(self):
self.binaryPolygon = ''
for v in range(self.nVertices):
    self.binaryPolygon += self.__encodeRadius(self.r[v])

def decodePolygon(self):
n = int(len(self.binaryPolygon) / self.p)
self.decodedPolygon = np.zeros(n)
for i in range(n):
    temp = self.binaryPolygon[self.p*i: (self.p*i+self.p)]
    self.decodedPolygon[i] = self.__decodeRadius(temp)

these functions are supported by the following additional private functions

def __getCoding(self):
self.c = np.linspace(self.rMin, self.rMax, 2**self.p)

the coding is a key for associating radial values with a binary index value defined by the binary precision, this approach reduces hidden precision truncation when dealing with the floating point radial values and just overal is very simple to implement across many different optimization problems

def __binary2Integer(self,b):
b = b[::-1]
temp = 0
for i in range(len(b)):
    temp += int(b[i])*(2**i)
return temp

converting from binary numbers to integers is a well documented process, each binary element is multiplied by 2 raised to the power of it's digit placement, the sum of all of these digits in this process gives an integer index that can be fed into the coding

def __decodeRadius(self, b):
self.__getCoding()
return self.c[self.__binary2Integer(b)]

decoding a binary number to a radial value follows the previous code and just inputs the index to the coding to recover the allowed radial value

  def __encodeRadius(self, v):
#create binary encoding and decoding key
self.__getCoding()
#find the index of the radial value from the binary key
idx = list(self.c).index(v)
#format that index to a binary string
B = format(idx, "b")
#append zeros to the front of the binary string to preserve precision
for i in range(self.p - len(B)):
    B = '0' + B
return B

encoding a radial value to a binary number is similarly simple, the index of the coding corresponding to the known radial value is found, that index is then formatted to a binary string. When formatting to a binary string the precision is limited to whatever the smallest value that can fit that integer but doesn't necessarily match the defined binary precision so additional '0' are added to the left hand side for precision matching.

for v in range(self.nVertices):
    self.binaryPolygon += self.__encodeRadius(self.r[v])

During encoding this process is repeated in a for-loop across every radial value in the list to give a single large binary string

n = int(len(self.binaryPolygon) / self.p)
self.decodedPolygon = np.zeros(n)

decoding is a bit more complicated, first from the length of the full binary polygon string and the binary precision the number of vertices is recovered.

for i in range(n):
    temp = self.binaryPolygon[self.p*i: (self.p*i+self.p)]
    self.decodedPolygon[i] = self.__decodeRadius(temp)

then each sub-string of the binary polygon string which are single binary radii are individual decoded and stored back inot a decoded polygon list.

Image Rendering

For interoperability with trained neural networks later on radial paramterization needs to be transformed into an image array, this is handled by arrayConversion() and polar2Cartesian()

def arrayConversion(self):
    self.polar2Cartesian()
    self.cartPolygon = np.stack((self.y,self.x), axis=1)
    self.arrayImage = grid_points_in_poly(self.dom, self.cartPolygon).astype(int)
    self.smoothedImage = filters.gaussian(self.arrayImage, self.s)
    self.smoothedImage = np.round(self.smoothedImage / np.max(self.smoothedImage))

def polar2Cartesian(self):
    self.__getCoding()
    for v in range(self.nVertices):
        self.r[v] = Support.closestRadius(self.r[v], self.rMin, self.rMax, self.p)
    self.x = self.r*np.cos(self.t) + self.dom[0]/2
    self.y = self.r*np.sin(self.t) + self.dom[1]/2
    xCenter, yCenter = np.mean(self.x), np.mean(self.y)
    a = np.arctan2(self.y - yCenter, self.x - xCenter)
    sortIdx = np.argsort(a)
    self.x, self.y = self.x[sortIdx], self.y[sortIdx]
    self.x = np.append(self.x, self.x[0]) - (xCenter - self.dom[0]/2)
    self.y = np.append(self.y, self.y[0]) - (yCenter - self.dom[1]/2)

Let's break down what's happening between these two functions, before the generated polar coordiantes can be converted to cartesian coordinates the individual radii need to comply with the coding discussed in the binary represenation section. So a radial value is replaced with the closest allowed radius

self.__getCoding()
for v in range(self.nVertices):
    self.r[v] = Support.closestRadius(self.r[v], self.rMin, self.rMax, self.p)

from there the polar coordinate conversion follows directly with transformation from an origin at \((0,0)\) to an origin in the center of the domain

self.x = self.r*np.cos(self.t) + self.dom[0]/2
self.y = self.r*np.sin(self.t) + self.dom[1]/2

this process could end at this point but an sorting of the vertices by angle is performed to guarantee that the polygon is always 'untangled' which is to say that no two polygon edges should intersect.

xCenter, yCenter = np.mean(self.x), np.mean(self.y)
a = np.arctan2(self.y - yCenter, self.x - xCenter)
sortIdx = np.argsort(a)
self.x, self.y = self.x[sortIdx], self.y[sortIdx]

the polygon should be closed so the first element of the list is appended to the end and all of center of mass of the polygon is shifted to the center of the domain

self.x = np.append(self.x, self.x[0]) - (xCenter - self.dom[0]/2)
self.y = np.append(self.y, self.y[0]) - (yCenter - self.dom[1]/2)

from the cartesian coordinates an image array is drawn with grid_points_in_poly() from scikit-image which takes a look at each polygon coordinate and builds an array with ones within the polygon perimeter and zeros outside the polygon perimeter

self.cartPolygon = np.stack((self.y,self.x), axis=1)
self.arrayImage = grid_points_in_poly(self.dom, self.cartPolygon).astype(int)

the image array is then smoothed by a gaussin filter to remove any small features

self.smoothedImage = filters.gaussian(self.arrayImage, self.s)
self.smoothedImage = np.round(self.smoothedImage / np.max(self.smoothedImage))

this smoothing results in the polygon edges acting more like splines to guide the nanostructure profile shape but still act as a reliable parameterization of the nanostructure.

drawing

Example Testing Code

#Example of using ProfileGeneration.py to generate a triangular polygon
#/testing/generationTriangle.py
import sys
import matplotlib.pyplot as plt

sys.path.append("../modules")
import ProfileGeneration

#initialize profile generation class

sigma = 5       #Images will be smoothed by a gaussian blur with sigma=5
precision = 12  #Encoding will occur with a binary precision of 12
#Polygons generated will have 24 Vertices with radial values between 20 and 80
pg = ProfileGeneration.ProfileGeneration(nVertices=24,rMin=20,rMax=80,s=5,p=12)

#generate a triangular polygon
pg.triangleGenerator()
#convert radial paramters into a smoothed image
pg.arrayConversion()

#plot the different realizations of the triangle
fig = plt.figure(figsize=(10, 10))
ax1 = plt.subplot(221,polar=True)
ax2 = plt.subplot(222)
ax3 = plt.subplot(223)
ax4 = plt.subplot(224)
ax1.set_rlim(0,90)
ax2.set_xlim(0,180); ax2.set_xticks([])
ax2.set_ylim(0,180); ax2.set_yticks([])
ax3.set_xticks([]); ax4.set_xticks([])
ax3.set_yticks([]); ax4.set_yticks([])
ax1.set_title("Radial Parameterization", fontsize=18)
ax2.set_title("Cartesian Coordinates", fontsize=18)
ax3.set_title("Array Image", fontsize=18)
ax4.set_title("Smoothed Image", fontsize=18)
ax1.plot(pg.t, pg.r, c='tab:red', zorder= 1)
ax1.scatter(pg.t, pg.r, c='tab:blue', s=10, zorder= 2)
ax2.plot(pg.x, pg.y, c='tab:red', zorder= 1)
ax2.scatter(pg.x, pg.y, c='tab:blue', s=10, zorder= 2)
ax3.imshow(pg.arrayImage, cmap='binary')
ax3.plot(pg.x, pg.y, c='tab:red', zorder= 1)
ax3.scatter(pg.x, pg.y, c='tab:blue', s=10, zorder= 2)
ax4.imshow(pg.smoothedImage, cmap='binary')
ax4.plot(pg.x, pg.y, c='tab:red', zorder= 1)
ax4.scatter(pg.x, pg.y, c='tab:blue', s=10, zorder= 2)
plt.savefig("temp.png")
plt.close()

#verify that shape encoding/decoding is working properly
#encode polygon to chromosome
pg.encodePolygon()
originalTri = pg.r
encodedTri = pg.binaryPolygon
#decode chromosome back to polygon
pg.decodeChromosome(encodedTri)
decodedTri = pg.decodedPolygon
print(f"Encoded Chromosome:     {encodedTri}")
print(f"Original Radial values: {originalTri}")
print(f"Decoded Radial values:  {decodedTri}")
#check all of the triangle radii and and return false if any radii doesn't match
print(f"Encoding Successful?: {any(decodedTri == originalTri)}")

Which should produce an image similar to the following

drawing