- Author:
- rjag008 <rjag008@auckland.ac.nz>
- Date:
- 2018-08-18 17:01:42+12:00
- Desc:
- Final Release
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/51d/rawfile/9ea33dda840d259caf68e26d21300bece4eeff3f/digitiser/mapping.py
'''
Created on 12/06/2018
@author: rjag008
'''
import pickle
import numpy as np
from scipy.interpolate import splprep, splev
from opencmiss.zinc.context import Context
from collections import OrderedDict
from opencmiss.zinc.node import Node
from opencmiss.zinc.element import Element, Elementbasis
def pointOnLine(points, line, rtol=1.0e-5, atol=1.0e-8):
"""Determine if a point is on a line segment
Input
points: Coordinates of either
* one point given by sequence [x, y]
* multiple points given by list of points or Nx2 array
line: Endpoint coordinates [[x0, y0], [x1, y1]] or
the equivalent 2x2 numeric array with each row corresponding
to a point.
rtol: Relative error for how close a point must be to be accepted
atol: Absolute error for how close a point must be to be accepted
Output
True or False
Notes
Line can be degenerate and function still works to discern coinciding
points from non-coinciding.
Tolerances rtol and atol are used with numpy.allclose()
"""
if len(points.shape) == 1:
# One point only - make into 1 x 2 array
points = points[np.newaxis, :]
one_point = True
else:
one_point = False
msg = 'Argument points must be either [x, y] or an Nx[>2] array of points'
if len(points.shape) != 2:
raise Exception(msg)
if not points.shape[0] > 0:
raise Exception(msg)
if points.shape[1] < 2:
raise Exception(msg)
N = points.shape[0] # Number of points
x = points[:, 0]
y = points[:, 1]
x0, y0 = line[0]
x1, y1 = line[1]
# Vector from beginning of line to point
a0 = x - x0
a1 = y - y0
# It's normal vector
a_normal0 = a1
a_normal1 = -a0
# Vector parallel to line
b0 = x1 - x0
b1 = y1 - y0
# Dot product
nominator = abs(a_normal0 * b0 + a_normal1 * b1)
denominator = b0 * b0 + b1 * b1
# Determine if point vector is parallel to line up to a tolerance
is_parallel = np.zeros(N, dtype=np.bool) # All False
is_parallel[nominator <= atol + rtol * denominator] = True
# Determine for points parallel to line if they are within end points
a0p = a0[is_parallel]
a1p = a1[is_parallel]
len_a = np.sqrt(a0p * a0p + a1p * a1p)
len_b = np.sqrt(b0 * b0 + b1 * b1)
cross = a0p * b0 + a1p * b1
# Initialise result to all False
result = np.zeros(N, dtype=np.bool)
# Result is True only if a0 * b0 + a1 * b1 >= 0 and len_a <= len_b
result[is_parallel] = (cross >= 0) * (len_a <= len_b)
# Return either boolean scalar or boolean vector
if one_point:
return result[0]
else:
return result
class MeshMapper(object):
'''
Map material coordinates to anatomical coordinates to region type
and vice versa
'''
def __init__(self,discret=50,boundaryDiscret=50):
'''
Constructor
'''
self.discret = discret
self.boundaryDiscret = boundaryDiscret
self.colors = [[147,231,255],[168,209,141],[162,198,255],[255,0,0],[255,218,104],[255,255,255]]
def setupByFile(self,soureMeshFile,circumferentialElements=8,axialElements=11):
self.context = Context('MeshMapper')
region = self.context.getDefaultRegion()
region.readFile(soureMeshFile)
self.region = region
with open(soureMeshFile,'r') as mesh:
self.mesh = mesh.read()
self.generateMaps(circumferentialElements,axialElements)
def setupByRegion(self,source,circumferentialElements=8,axialElements=11):
self.region = source
#As of 13 June 2018, no python api for writing to memory from zinc region
import tempfile,os
ofile = tempfile.NamedTemporaryFile(delete=False)
self.region.writeFile(ofile.name)
with open(ofile.name,'r') as mesh:
self.mesh = mesh.read()
os.remove(ofile.name)
self.generateMaps(circumferentialElements,axialElements)
def setupByPickle(self,source):
self.context = Context('MeshMapper')
region = self.context.getDefaultRegion()
self.region = region
with open(source,'rb') as ser:
discret,boundaryDiscret,lb,rb,ladder,surfaceCoordinates,anatomicalCoordinates,regionMappings,mesh=pickle.load(ser)
self.discret = discret
self.boundaryDiscret = boundaryDiscret
self.leftWallCoordinates = lb
self.rightWallCoordinates = rb
self.ladderCoordinates = ladder
self.surfaceCoordinates = surfaceCoordinates
self.anatomicalCoordinates = anatomicalCoordinates
self.regionMappings = regionMappings
sir = region.createStreaminformationRegion()
sir.createStreamresourceMemoryBuffer(str(mesh))
region.read(sir)
self.mesh = mesh
def pointsInPolygon(self,surfaceCoordinates,polygon,closed=True):
'''
Based on https://github.com/inasafe/python-safe/blob/master/safe/engine/polygon.py
'''
# Suppress numpy warnings (as we'll be dividing by zero)
original_numpy_settings = np.seterr(invalid='ignore', divide='ignore')
rtol = 0.0
atol = 0.0
scShape = surfaceCoordinates.shape
points = surfaceCoordinates.reshape((-1,scShape[-1]))
# Get polygon extents to quickly rule out points that
# are outside its bounding box
maxpx,maxpy = np.max(polygon,axis=0)
minpx,minpy = np.min(polygon,axis=0)
M = points.shape[0]
N = polygon.shape[0]
x = points[:, 0]
y = points[:, 1]
# Vector keeping track of which points are inside
inside = np.zeros(M, dtype=np.int) # All assumed outside initially
# Mask for points can be considered for inclusion
candidates = np.ones(M, dtype=np.bool) # All True initially
# Only work on those that are inside polygon bounding box
outsideBox = (x > maxpx) + (x < minpx) + (y > maxpy) + (y < minpy)
insideBox = np.logical_not(outsideBox)
candidates *= insideBox
# Don't continue if all points are outside bounding box
if not np.sometrue(candidates):
return inside>0
# Restrict computations to candidates only
ix = np.where(candidates==True)[0]
inCandidates = candidates[ix]
cPoints = points[ix]
cInside = inside[ix]
cy = y[ix]
cx = x[ix]
# Find points on polygon boundary
if closed:
for i in range(N):
# Loop through polygon edges
j = (i + 1) % N
edge = [polygon[i, :], polygon[j, :]]
# Select those that are on the boundary
boundaryPoints = pointOnLine(cPoints, edge, rtol, atol)
cInside[boundaryPoints] = True
# Remove boundary point from further analysis
inCandidates[boundaryPoints] = False
else:
for i in range(N):
# Loop through polygon edges
j = (i + 1) % N
edge = [polygon[i, :], polygon[j, :]]
# Select those that are on the boundary
boundaryPoints = pointOnLine(cPoints, edge, rtol, atol)
cInside[boundaryPoints] = False
# Remove boundary point from further analysis
inCandidates[boundaryPoints] = False
# Algorithm for finding points inside polygon
for i in range(N):
# Loop through polygon edges
j = (i + 1) % N
px_i, py_i = polygon[i, :]
px_j, py_j = polygon[j, :]
# Intersection formula
sigma = (cy - py_i) / (py_j - py_i) * (px_j - px_i)
seg_i = (py_i < cy) * (py_j >= cy)
seg_j = (py_j < cy) * (py_i >= cy)
mask = (px_i + sigma < cx) * (seg_i + seg_j) * inCandidates
cInside[mask] = 1 - cInside[mask]
# Restore numpy warnings
np.seterr(**original_numpy_settings)
inside[ix] = cInside
if len(scShape)==3:
return inside.reshape((scShape[0],scShape[1]))>0
return inside>0
def generateMaps(self,circumferentialElements,axialElements):
fieldModule = self.region.getFieldmodule()
coordinatesField = fieldModule.findFieldByName('coordinates')
#Get the ladder coordinates
nodeset = fieldModule.findNodesetByFieldDomainType(coordinatesField.DOMAIN_TYPE_NODES)
nodeIterator = nodeset.createNodeiterator ()
nodes = dict()
node = nodeIterator.next()
while node.isValid():
nodes[int(node.getIdentifier())-1] = node
node = nodeIterator.next()
fieldCache = fieldModule.createFieldcache()
coordinates = np.zeros((len(nodes),3))
for nd,node in nodes.items():
fieldCache.setNode(node)
_,coord = coordinatesField.evaluateReal(fieldCache,3)
coordinates[nd-1] = coord
ladderCoordinates = np.zeros((axialElements+1,6))
#Fit splines through each axial ring and determine the coordinates for which z = 0
cds = np.zeros((circumferentialElements,3))
for a in range(axialElements+1):
axoff = a*circumferentialElements
#Find the coordinates of a cross-section nodes number go from start-start+numCircumferentialNodes
for nd in range(circumferentialElements):
nid = nd + axoff
cds[nd] = coordinates[nid]
#Spline interpolate
tck, u = splprep(cds.T, u=None, s=0.0, per=1)
#Get the midpoint
xs, ys, zs = splev([0.5], tck, der=0)
ladderCoordinates[a,:3]=cds[0]
ladderCoordinates[a,3:]=[xs[0],ys[0],zs[0]]
#Compute the mesh walls using ladder Coordinates
#Interpolate of left
tck, u = splprep(ladderCoordinates[:,:2].T, u=None, s=0.0)
u_new = np.linspace(u.min(), u.max(), self.boundaryDiscret)
xs, ys = splev(u_new, tck, der=0)
self.leftWallCoordinates = np.c_[xs,ys]
#Interpolate of right
tck, u = splprep(ladderCoordinates[:,3:5].T, u=None, s=0.0)
u_new = np.linspace(u.min(), u.max(), self.boundaryDiscret)[::-1]
xs, ys = splev(u_new, tck, der=0)
self.rightWallCoordinates = np.c_[xs,ys]
self.ladderCoordinates = ladderCoordinates
mesh = fieldModule.findMeshByDimension(3)
ei = mesh.createElementiterator()
elemDict = dict()
elem = ei.next()
while elem.isValid():
elemDict[int(elem.getIdentifier())] = elem
elem = ei.next()
threeDMesh = True
if len(elemDict)==0: #Surface meshes
threeDMesh = False
mesh = fieldModule.findMeshByDimension(2)
ei = mesh.createElementiterator()
elemDict = dict()
elem = ei.next()
while elem.isValid():
elemDict[int(elem.getIdentifier())] = elem
elem = ei.next()
fieldCache = fieldModule.createFieldcache()
# following method uses the mesh, but this causes issues with ladder intersection continuity
# x2 = np.linspace(0.00001, 1.0-0.00001, self.boundaryDiscret)
# x1 = x2*0.0
# xmid = x1+0.5
# lbcoordXI = np.c_[x1,x2]
# rbcoordXI = np.c_[xmid,x2]
#
# for i,coordXI in enumerate([lbcoordXI,rbcoordXI]):
# coords = []
# for val in coordXI:
# xpos = int(val[0]*circumferentialElements)
# xi1 = val[0]*circumferentialElements - xpos
# ypos = int(val[1]*axialElements)
# xi2 = val[1]*axialElements - ypos
# eid = xpos + ypos*circumferentialElements + 1
# fieldCache.setMeshLocation(elemDict[eid],[xi1,xi2,0.0])
# _,v = coordinatesField.evaluateReal(fieldCache,3)
# coords.append(v)
# if i==0:
# lbcoords = coords
# else:
# rbcoords = coords
#
# self.leftWallCoordinates = np.np.array(lbcoords)
# self.rightWallCoordinates = np.np.array(rbcoords)
# self.ladderCoordinates = ladderCoordinates
x2 = np.linspace(0.00001, 1.0-0.00001, self.discret)
surfaceCoordinates = np.zeros((self.discret,self.discret,3))
xi1,xi2 = np.meshgrid(x2,x2)
if threeDMesh:
x = [0.0,0.0,0.0]
else:
x = [0.0,0.0]
for i in range(self.discret):
for j in range(self.discret):
xpos = int(xi1[i,j]*circumferentialElements)
xi1s = xi1[i,j]*circumferentialElements - xpos
ypos = int(xi2[i,j]*axialElements)
xi2s = xi2[i,j]*axialElements - ypos
eid = xpos + ypos*circumferentialElements + 1
x[0] = xi1s
x[1] = xi2s
fieldCache.setMeshLocation(elemDict[eid],x)
_,v = coordinatesField.evaluateReal(fieldCache,3)
surfaceCoordinates[i,j] = v
anatomicalCoordinates = np.zeros((self.discret,self.discret,2))
#Columns - x
for j in range(1,self.discret):
anatomicalCoordinates[:,j,0] = anatomicalCoordinates[:,j-1,0] + np.linalg.norm(surfaceCoordinates[:,j-1]-surfaceCoordinates[:,j],axis=1)
#Rows - y
for j in range(1,self.discret):
anatomicalCoordinates[j,:,1] = anatomicalCoordinates[j-1,:,1] + np.linalg.norm(surfaceCoordinates[j-1,:]-surfaceCoordinates[j,:],axis=1)
xmax = anatomicalCoordinates[:,:,0].max()
ymax = anatomicalCoordinates[:,:,1].max()
#Scaling
anatomicalCoordinates[:,:,0] = anatomicalCoordinates[:,:,0]/xmax
anatomicalCoordinates[:,:,1] = anatomicalCoordinates[:,:,1]/ymax
'''
#Centering
for j in range(self.discret):
xm = anatomicalCoordinates[j,-1,0]/2.0
anatomicalCoordinates[j,:,0] = anatomicalCoordinates[j,:,0] - xm + 0.5
for j in range(self.discret):
ym = anatomicalCoordinates[-1,j,1]/2.0
anatomicalCoordinates[:,j,1] = anatomicalCoordinates[:,j,1] - ym + 0.5
'''
self.surfaceCoordinates = surfaceCoordinates
self.anatomicalCoordinates = anatomicalCoordinates
def generateSegmentLengthsAndWeights(self,circumferentialElements,axialElements):
#Generate node to element xi mapping
numberOfCircumfrentialNodes = circumferentialElements
# create node element map
ev = 0
nodes = dict()
#Compute the material coordinates of nodes
for lengthElementIdx in range(1,axialElements+1):
for circumfrentialElementIdx in range(1,circumferentialElements+1):
#Element coord with respect to elements along each dim
localNode1 = circumfrentialElementIdx + (lengthElementIdx - 1)*circumferentialElements
if circumfrentialElementIdx == circumferentialElements:
localNode2 = 1 + (lengthElementIdx-1)*numberOfCircumfrentialNodes
else:
localNode2 = localNode1 + 1
localNode3 = localNode1 + numberOfCircumfrentialNodes
localNode4 = localNode2 + numberOfCircumfrentialNodes
l1xi = [0.0,0.0]
l2xi = [1.0,0.0]
l3xi = [0.0,1.0]
l4xi = [1.0,1.0]
ev +=1
nds = [localNode1,localNode2,localNode3,localNode4]
xis = [l1xi,l2xi,l3xi,l4xi]
for i,nd in enumerate(nds):
if nd in nodes:
if xis[i] == [0.0,0.0]:
nodes[nd] = [ev,xis[i]]
if nodes[nd][1] != [0.0,0.0] and xis[i]==[0.0,1.0]:
nodes[nd] = [ev,xis[i]]
else:
nodes[nd] = [ev,xis[i]]
def arcLength(x):
npts = x.shape[0]
arc = 0.0
for k in range(1, npts):
arc = arc + np.linalg.norm(x[k]-x[k-1])
return arc
segmentLengths = np.zeros((len(nodes),2))
ndis = 20
xdiscret = np.linspace(0.0,1.0,ndis)
fieldModule = self.region.getFieldmodule()
coordinatesField = fieldModule.findFieldByName('coordinates')
mesh = fieldModule.findMeshByDimension(3)
ei = mesh.createElementiterator()
elemDict = dict()
elem = ei.next()
while elem.isValid():
elemDict[int(elem.getIdentifier())] = elem
elem = ei.next()
fieldCache = fieldModule.createFieldcache()
coords = np.zeros((ndis,3))
for nd, xi in nodes.items():
ev = xi[0]
xiv = xi[1]
elem = elemDict[ev]
if xiv[0] == 0.0:
if xiv[1] == 0.0: #bottom left node
#Get the x1 length
for i in range(ndis):
fieldCache.setMeshLocation(elem,[xdiscret[i],0.0,0.0])
_,v = coordinatesField.evaluateReal(fieldCache,3)
coords[i] = v
segmentLengths[nd-1,0] = arcLength(coords)
#Get the x2 length
for i in range(ndis):
fieldCache.setMeshLocation(elem,[0.0,xdiscret[i],0.0])
_,v = coordinatesField.evaluateReal(fieldCache,3)
coords[i] = v
segmentLengths[nd-1,1] = arcLength(coords)
elif xiv[1] == 1.0: #top left node
#Get the x1 length
for i in range(ndis):
fieldCache.setMeshLocation(elem,[xdiscret[i],0.0,0.0])
_,v = coordinatesField.evaluateReal(fieldCache,3)
coords[i] = v
segmentLengths[nd-1,0] = arcLength(coords)
circumferentialLengths = np.zeros((axialElements+1,1))
axialLengths = np.zeros((circumferentialElements,1))
for a in range(axialElements+1):
snode = a*numberOfCircumfrentialNodes
for i in range(numberOfCircumfrentialNodes):
circumferentialLengths[a] += segmentLengths[snode+i,0]
axialLengths[i] +=segmentLengths[snode+i,1]
segmentWeights = np.zeros(segmentLengths.shape)
for a in range(axialElements+1):
snode = a*numberOfCircumfrentialNodes
for i in range(numberOfCircumfrentialNodes):
segmentWeights[snode+i,0] = segmentLengths[snode+i,0]/circumferentialLengths[a]
segmentWeights[snode+i,1] = segmentLengths[snode+i,1]/axialLengths[i]
self.segmentWeights = segmentWeights
self.segmentLenths = segmentLengths
self.circumferentialLengths = circumferentialLengths
self.axialLengths = axialLengths
def generateFlatMesh(self,region,circumferentialElements,axialElements):
elementsCount1 = circumferentialElements
elementsCount2 = axialElements
coordinateDimensions = 3
#context = Context('pm')
#region = context.getDefaultRegion()
fm = region.getFieldmodule()
fm.beginChange()
coordinates = fm.createFieldFiniteElement(coordinateDimensions)
coordinates.setName('coordinates')
coordinates.setManaged(True)
coordinates.setTypeCoordinate(True)
coordinates.setCoordinateSystemType(coordinates.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN)
coordinates.setComponentName(1, 'x')
coordinates.setComponentName(2, 'y')
coordinates.setComponentName(3, 'z')
texture = fm.createFieldFiniteElement(2)
texture.setName('texture')
texture.setManaged(True)
texture.setCoordinateSystemType(coordinates.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN)
texture.setComponentName(1, 'x')
texture.setComponentName(2, 'y')
nodes = fm.findNodesetByFieldDomainType(coordinates.DOMAIN_TYPE_NODES)
nodetemplate = nodes.createNodetemplate()
nodetemplate.defineField(coordinates)
nodetemplate.defineField(texture)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
nodetemplate.setValueNumberOfVersions(texture, -1, Node.VALUE_LABEL_VALUE, 1)
nodetemplate.setValueNumberOfVersions(texture, -1, Node.VALUE_LABEL_D_DS1, 1)
nodetemplate.setValueNumberOfVersions(texture, -1, Node.VALUE_LABEL_D_DS2, 1)
mesh = fm.findMeshByDimension(2)
bicubicHermiteBasis = fm.createElementbasis(2, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE)
eft = mesh.createElementfieldtemplate(bicubicHermiteBasis)
for n in range(4):
eft.setFunctionNumberOfTerms(n*4 + 4, 0)
elementtemplate = mesh.createElementtemplate()
elementtemplate.setElementShapeType(Element.SHAPE_TYPE_SQUARE)
result = elementtemplate.defineField(coordinates, -1, eft)
result = elementtemplate.defineField(texture, -1, eft)
cache = fm.createFieldcache()
# create nodes
nodeIdentifier = 1
x = [ 0.0, 0.0, 0.0 ]
dx_ds1 = [ 1.0 / elementsCount1, 0.0, 0.0 ]
dx_ds2 = [ 0.0, 1.0 / elementsCount2, 0.0 ]
numberOfCircumferentialNodes = elementsCount1 +1
allNodes = dict()
for n2 in range(elementsCount2 + 1):
x[1] = float(n2) / elementsCount2 - 0.5
yv = 2*x[1]
for n1 in range(numberOfCircumferentialNodes):
x[0] = float(n1) / (elementsCount1) - 0.5
#Use square to circle map
xv = 2*x[0]
#Flip along x to match the element orientated opening of the mesh
xx = xv*np.sqrt(1.0-0.5*yv**2)
yy = yv*np.sqrt(1.0-0.5*xv**2)
node = nodes.createNode(nodeIdentifier, nodetemplate)
allNodes[nodeIdentifier] = node
cache.setNode(node)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [xx,yy,0.0])
texture.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, [x[0]+0.5,x[1]+0.5])
texture.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
texture.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2)
nodeIdentifier = nodeIdentifier + 1
#Introduce some kink at the top and bottom
'''
#Does not produce expected mesh
#Normalize x1-distribution of nodes
for a in range(elementsCount2+1):
snode = a*(elementsCount1 +1 ) + 1
cylinderNode = a*(elementsCount1) + 1
enode = snode + elementsCount1
cache.setNode(allNodes[snode])
_,scoord = coordinates.evaluateReal(cache,3)
cache.setNode(allNodes[enode])
_,ecoord = coordinates.evaluateReal(cache,3)
scoord = np.array(scoord)
vec = np.array(ecoord) - scoord
length = np.linalg.norm(vec)
vec = vec/length
#slength = length/elementsCount1
for n1 in range(1,elementsCount1):
cnode = snode + n1
cache.setNode(allNodes[cnode])
nc = scoord + vec*length*n1*self.segmentWeights[cylinderNode-1,0]
coordinates.assignReal(cache,list(nc))
#Normalize x2-distribution of nodes
for i in range(elementsCount1+1):
for a in range(1,elementsCount2+1):
snode = (a-1)*(elementsCount1+1) +1+i
cnode = snode + numberOfCircumferentialNodes
cylinderNode = (a-1)*(elementsCount1) + 1 + i
print snode,cylinderNode
cache.setNode(allNodes[snode])
_,scoord = coordinates.evaluateReal(cache,3)
cache.setNode(allNodes[cnode])
_,ecoord = coordinates.evaluateReal(cache,3)
scoord = np.array(scoord)
ecoord = np.array(ecoord)
vec = scoord-ecoord
ecoord = scoord +vec*self.segmentWeights[cylinderNode-1,1]
coordinates.assignReal(cache,list(ecoord))
'''
#Manually distort
NumberOfElements = elementsCount1*elementsCount2
# create elements
elementIdentifier = 1
no2 = (elementsCount1 + 1)
for e2 in range(elementsCount2):
for e1 in range(elementsCount1):
#Note that the xi ordering is different with the top having the lowest element number to the bottom
element = mesh.createElement(NumberOfElements-elementIdentifier+1, elementtemplate)
bni = e2*no2 + e1 + 1
nodeIdentifiers = [ bni, bni + 1, bni + no2, bni + no2 + 1 ]
result = element.setNodesByIdentifier(eft, nodeIdentifiers)
elementIdentifier = elementIdentifier + 1
fm.defineAllFaces()
smooth = fm.createFieldsmoothing()
coordinates.smooth(smooth)
fm.endChange()
def assignRegions(self,regionMappings,colors):
'''
Region assignment is in the descending order of regionvalue
'''
rnm = list(reversed(regionMappings.keys()))
shape =self.surfaceCoordinates.shape
rshape =(shape[0],shape[1])
if len(shape)==2:
rshape=(shape[0],1)
self.regionMappings = np.ones(rshape,dtype=np.int)*len(rnm)
self.colors = []
for i,rn in enumerate(rnm):
self.colors.append(colors[rn])
poly = regionMappings[rn]
if isinstance(poly, list):
poly= np.array(poly)
spts = self.pointsInPolygon(self.surfaceCoordinates, poly)
self.regionMappings[spts] = i
#Add a color for unmapped regions
self.colors.append([255,255,255])
def getBoundaryPoints(self):
'''
Return the boundary coordinates for the left, right wall and ladder
'''
#Senfd copies do that modification are not reflected
return np.array(self.leftWallCoordinates),np.array(self.rightWallCoordinates),np.array(self.ladderCoordinates)
def serialize(self,filename):
with open(filename,'wb') as ser:
if hasattr(self,'regionMappings'):
sobj = [self.discret,self.boundaryDiscret,self.leftWallCoordinates,self.rightWallCoordinates, self.ladderCoordinates,\
self.surfaceCoordinates,self.anatomicalCoordinates,self.regionMappings,self.mesh]
else:
sobj = [self.discret,self.boundaryDiscret,self.leftWallCoordinates,self.rightWallCoordinates, self.ladderCoordinates,\
self.surfaceCoordinates,self.anatomicalCoordinates,None,self.mesh]
pickle.dump(sobj,ser,2)
def createImage(self,filename):
pic = np.zeros((self.discret,self.discret,3),dtype=np.uint8)
if hasattr(self, 'regionMappings'):
for i in range(self.discret):
for j in range(self.discret):
pic[i,j] = self.colors[self.regionMappings[i,j]] #Note that the coordinate system for numpy is different from qimage
else:
for i in range(self.discret):
for j in range(self.discret):
x,y = self.anatomicalCoordinates[i,j]
x = int(x*self.discret)+1
y = int(y*self.discret)+1
pic[x,y] = [255, 255, 255] #Note that the coordinate system for numpy is different from qimage
from PIL import Image
img = Image.fromarray(pic, 'RGB')
img = img.resize((400,400),Image.BICUBIC)
img.save(filename,format="png")
#img.show()
