''' 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()