- ScatterPlotItem disables render cache during export
- Fixes for SVG exporter - functions.isosurface() is a bazillion times faster (API change: return value format has changed)
This commit is contained in:
parent
3de5719011
commit
ecca8855df
@ -45,9 +45,9 @@ data = np.abs(np.fromfunction(psi, (50,50,100)))
|
||||
|
||||
|
||||
print("Generating isosurface..")
|
||||
verts = pg.isosurface(data, data.max()/4.)
|
||||
verts, faces = pg.isosurface(data, data.max()/4.)
|
||||
|
||||
md = gl.MeshData(vertexes=verts)
|
||||
md = gl.MeshData(vertexes=verts, faces=faces)
|
||||
|
||||
colors = np.ones((md.faceCount(), 4), dtype=float)
|
||||
colors[:,3] = 0.2
|
||||
|
@ -73,7 +73,8 @@ class Exporter(object):
|
||||
|
||||
def getSourceRect(self):
|
||||
if isinstance(self.item, pg.GraphicsScene):
|
||||
return self.item.getViewWidget().viewRect()
|
||||
w = self.item.getViewWidget()
|
||||
return w.viewportTransform().inverted()[0].mapRect(w.rect())
|
||||
else:
|
||||
return self.item.sceneBoundingRect()
|
||||
|
||||
|
@ -36,11 +36,14 @@ class SVGExporter(Exporter):
|
||||
return
|
||||
self.svg = QtSvg.QSvgGenerator()
|
||||
self.svg.setFileName(fileName)
|
||||
self.svg.setSize(QtCore.QSize(100,100))
|
||||
#self.svg.setResolution(600)
|
||||
dpi = QtGui.QDesktopWidget().physicalDpiX()
|
||||
## not really sure why this works, but it seems to be important:
|
||||
self.svg.setSize(QtCore.QSize(self.params['width']*dpi/90., self.params['height']*dpi/90.))
|
||||
self.svg.setResolution(dpi)
|
||||
#self.svg.setViewBox()
|
||||
targetRect = QtCore.QRect(0, 0, self.params['width'], self.params['height'])
|
||||
sourceRect = self.getSourceRect()
|
||||
|
||||
painter = QtGui.QPainter(self.svg)
|
||||
try:
|
||||
self.setExportMode(True)
|
||||
|
181
functions.py
181
functions.py
@ -1145,26 +1145,33 @@ def isocurve(data, level):
|
||||
|
||||
return lines ## a list of pairs of points
|
||||
|
||||
|
||||
IsosurfaceDataCache = None
|
||||
def isosurface(data, level):
|
||||
"""
|
||||
Generate isosurface from volumetric data using marching cubes algorithm.
|
||||
See Paul Bourke, "Polygonising a Scalar Field"
|
||||
(http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/)
|
||||
(http://paulbourke.net/geometry/polygonise/)
|
||||
|
||||
*data* 3D numpy array of scalar values
|
||||
*level* The level at which to generate an isosurface
|
||||
|
||||
Returns an array of vertex coordinates (N, 3, 3);
|
||||
|
||||
This function is SLOW; plenty of room for optimization here.
|
||||
Returns an array of vertex coordinates (Nv, 3) and an array of
|
||||
per-face vertex indexes (Nf, 3)
|
||||
"""
|
||||
## For improvement, see:
|
||||
##
|
||||
## Efficient implementation of Marching Cubes' cases with topological guarantees.
|
||||
## Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan Tavares.
|
||||
## Journal of Graphics Tools 8(2): pp. 1-15 (december 2003)
|
||||
|
||||
## Precompute lookup tables on the first run
|
||||
global IsosurfaceDataCache
|
||||
if IsosurfaceDataCache is None:
|
||||
## map from grid cell index to edge index.
|
||||
## grid cell index tells us which corners are below the isosurface,
|
||||
## edge index tells us which edges are cut by the isosurface.
|
||||
## (Data stolen from Bourk; see above.)
|
||||
edgeTable = [
|
||||
edgeTable = np.array([
|
||||
0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
|
||||
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
|
||||
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
|
||||
@ -1196,7 +1203,7 @@ def isosurface(data, level):
|
||||
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
|
||||
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
|
||||
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
|
||||
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 ]
|
||||
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 ], dtype=np.uint16)
|
||||
|
||||
## Table of triangles to use for filling each grid cell.
|
||||
## Each set of three integers tells us which three edges to
|
||||
@ -1460,27 +1467,43 @@ def isosurface(data, level):
|
||||
[0, 3, 8],
|
||||
[]
|
||||
]
|
||||
edgeShifts = np.array([ ## maps edge ID (0-11) to (x,y,z) cell offset and edge ID (0-2)
|
||||
[0, 0, 0, 0],
|
||||
[1, 0, 0, 1],
|
||||
[0, 1, 0, 0],
|
||||
[0, 0, 0, 1],
|
||||
[0, 0, 1, 0],
|
||||
[1, 0, 1, 1],
|
||||
[0, 1, 1, 0],
|
||||
[0, 0, 1, 1],
|
||||
[0, 0, 0, 2],
|
||||
[1, 0, 0, 2],
|
||||
[1, 1, 0, 2],
|
||||
[0, 1, 0, 2],
|
||||
#[9, 9, 9, 9] ## fake
|
||||
], dtype=np.ubyte)
|
||||
nTableFaces = np.array([len(f)/3 for f in triTable], dtype=np.ubyte)
|
||||
faceShiftTables = [None]
|
||||
for i in range(1,6):
|
||||
## compute lookup table of index: vertexes mapping
|
||||
faceTableI = np.zeros((len(triTable), i*3), dtype=np.ubyte)
|
||||
faceTableInds = np.argwhere(nTableFaces == i)
|
||||
faceTableI[faceTableInds[:,0]] = np.array([triTable[j] for j in faceTableInds])
|
||||
faceTableI = faceTableI.reshape((len(triTable), i, 3))
|
||||
faceShiftTables.append(edgeShifts[faceTableI])
|
||||
|
||||
## translation between edge index and
|
||||
## the vertex indexes that bound the edge
|
||||
edgeKey = [
|
||||
[(0,0,0), (1,0,0)],
|
||||
[(1,0,0), (1,1,0)],
|
||||
[(1,1,0), (0,1,0)],
|
||||
[(0,1,0), (0,0,0)],
|
||||
[(0,0,1), (1,0,1)],
|
||||
[(1,0,1), (1,1,1)],
|
||||
[(1,1,1), (0,1,1)],
|
||||
[(0,1,1), (0,0,1)],
|
||||
[(0,0,0), (0,0,1)],
|
||||
[(1,0,0), (1,0,1)],
|
||||
[(1,1,0), (1,1,1)],
|
||||
[(0,1,0), (0,1,1)],
|
||||
]
|
||||
## Let's try something different:
|
||||
#faceTable = np.empty((256, 5, 3, 4), dtype=np.ubyte) # (grid cell index, faces, vertexes, edge lookup)
|
||||
#for i,f in enumerate(triTable):
|
||||
#f = np.array(f + [12] * (15-len(f))).reshape(5,3)
|
||||
#faceTable[i] = edgeShifts[f]
|
||||
|
||||
|
||||
IsosurfaceDataCache = (faceShiftTables, edgeShifts, edgeTable, nTableFaces)
|
||||
else:
|
||||
faceShiftTables, edgeShifts, edgeTable, nTableFaces = IsosurfaceDataCache
|
||||
|
||||
|
||||
facets = []
|
||||
|
||||
## mark everything below the isosurface level
|
||||
mask = data < level
|
||||
@ -1494,35 +1517,93 @@ def isosurface(data, level):
|
||||
for k in [0,1]:
|
||||
fields[i,j,k] = mask[slices[i], slices[j], slices[k]]
|
||||
vertIndex = i - 2*j*i + 3*j + 4*k ## this is just to match Bourk's vertex numbering scheme
|
||||
#print i,j,k," : ", fields[i,j,k], 2**vertIndex
|
||||
index += fields[i,j,k] * 2**vertIndex
|
||||
#print index
|
||||
#print index
|
||||
|
||||
## add facets
|
||||
for i in range(index.shape[0]): # data x-axis
|
||||
for j in range(index.shape[1]): # data y-axis
|
||||
for k in range(index.shape[2]): # data z-axis
|
||||
tris = triTable[index[i,j,k]]
|
||||
for l in range(0, len(tris), 3): ## faces for this grid cell
|
||||
edges = tris[l:l+3]
|
||||
pts = []
|
||||
for m in [0,1,2]: # points in this face
|
||||
p1 = edgeKey[edges[m]][0]
|
||||
p2 = edgeKey[edges[m]][1]
|
||||
v1 = data[i+p1[0], j+p1[1], k+p1[2]]
|
||||
v2 = data[i+p2[0], j+p2[1], k+p2[2]]
|
||||
f = (level-v1) / (v2-v1)
|
||||
fi = 1.0 - f
|
||||
p = ( ## interpolate between corners
|
||||
p1[0]*fi + p2[0]*f + i + 0.5,
|
||||
p1[1]*fi + p2[1]*f + j + 0.5,
|
||||
p1[2]*fi + p2[2]*f + k + 0.5
|
||||
)
|
||||
pts.append(p)
|
||||
facets.append(pts)
|
||||
### Generate table of edges that have been cut
|
||||
cutEdges = np.zeros([x+1 for x in index.shape]+[3], dtype=np.uint32)
|
||||
edges = edgeTable[index]
|
||||
for i, shift in enumerate(edgeShifts[:12]):
|
||||
slices = [slice(shift[j],cutEdges.shape[j]+(shift[j]-1)) for j in range(3)]
|
||||
cutEdges[slices[0], slices[1], slices[2], shift[3]] += edges & 2**i
|
||||
|
||||
return np.array(facets)
|
||||
## for each cut edge, interpolate to see where exactly the edge is cut and generate vertex positions
|
||||
m = cutEdges > 0
|
||||
vertexInds = np.argwhere(m) ## argwhere is slow!
|
||||
vertexes = vertexInds[:,:3].astype(np.float32)
|
||||
dataFlat = data.reshape(data.shape[0]*data.shape[1]*data.shape[2])
|
||||
|
||||
## re-use the cutEdges array as a lookup table for vertex IDs
|
||||
cutEdges[vertexInds[:,0], vertexInds[:,1], vertexInds[:,2], vertexInds[:,3]] = np.arange(vertexInds.shape[0])
|
||||
|
||||
for i in [0,1,2]:
|
||||
vim = vertexInds[:,3] == i
|
||||
vi = vertexInds[vim, :3]
|
||||
viFlat = (vi * (np.array(data.strides[:3]) / data.itemsize)[np.newaxis,:]).sum(axis=1)
|
||||
v1 = dataFlat[viFlat]
|
||||
v2 = dataFlat[viFlat + data.strides[i]/data.itemsize]
|
||||
vertexes[vim,i] += (level-v1) / (v2-v1)
|
||||
|
||||
### compute the set of vertex indexes for each face.
|
||||
|
||||
## This works, but runs a bit slower.
|
||||
#cells = np.argwhere((index != 0) & (index != 255)) ## all cells with at least one face
|
||||
#cellInds = index[cells[:,0], cells[:,1], cells[:,2]]
|
||||
#verts = faceTable[cellInds]
|
||||
#mask = verts[...,0,0] != 9
|
||||
#verts[...,:3] += cells[:,np.newaxis,np.newaxis,:] ## we now have indexes into cutEdges
|
||||
#verts = verts[mask]
|
||||
#faces = cutEdges[verts[...,0], verts[...,1], verts[...,2], verts[...,3]] ## and these are the vertex indexes we want.
|
||||
|
||||
|
||||
## To allow this to be vectorized efficiently, we count the number of faces in each
|
||||
## grid cell and handle each group of cells with the same number together.
|
||||
## determine how many faces to assign to each grid cell
|
||||
nFaces = nTableFaces[index]
|
||||
totFaces = nFaces.sum()
|
||||
faces = np.empty((totFaces, 3), dtype=np.uint32)
|
||||
ptr = 0
|
||||
#import debug
|
||||
#p = debug.Profiler('isosurface', disabled=False)
|
||||
|
||||
## this helps speed up an indexing operation later on
|
||||
cs = np.array(cutEdges.strides)/cutEdges.itemsize
|
||||
cutEdges = cutEdges.flatten()
|
||||
|
||||
## this, strangely, does not seem to help.
|
||||
#ins = np.array(index.strides)/index.itemsize
|
||||
#index = index.flatten()
|
||||
|
||||
for i in range(1,6):
|
||||
### expensive:
|
||||
#p.mark('1')
|
||||
cells = np.argwhere(nFaces == i) ## all cells which require i faces (argwhere is expensive)
|
||||
#p.mark('2')
|
||||
if cells.shape[0] == 0:
|
||||
continue
|
||||
#cellInds = index[(cells*ins[np.newaxis,:]).sum(axis=1)]
|
||||
cellInds = index[cells[:,0], cells[:,1], cells[:,2]] ## index values of cells to process for this round
|
||||
#p.mark('3')
|
||||
|
||||
### expensive:
|
||||
verts = faceShiftTables[i][cellInds]
|
||||
#p.mark('4')
|
||||
verts[...,:3] += cells[:,np.newaxis,np.newaxis,:] ## we now have indexes into cutEdges
|
||||
verts = verts.reshape((verts.shape[0]*i,)+verts.shape[2:])
|
||||
#p.mark('5')
|
||||
|
||||
### expensive:
|
||||
#print verts.shape
|
||||
verts = (verts * cs[np.newaxis, np.newaxis, :]).sum(axis=2)
|
||||
#vertInds = cutEdges[verts[...,0], verts[...,1], verts[...,2], verts[...,3]] ## and these are the vertex indexes we want.
|
||||
vertInds = cutEdges[verts]
|
||||
#p.mark('6')
|
||||
nv = vertInds.shape[0]
|
||||
#p.mark('7')
|
||||
faces[ptr:ptr+nv] = vertInds #.reshape((nv, 3))
|
||||
#p.mark('8')
|
||||
ptr += nv
|
||||
|
||||
return vertexes, faces
|
||||
|
||||
|
||||
|
||||
|
@ -233,7 +233,7 @@ class ScatterPlotItem(GraphicsObject):
|
||||
self.bounds = [None, None] ## caches data bounds
|
||||
self._maxSpotWidth = 0 ## maximum size of the scale-variant portion of all spots
|
||||
self._maxSpotPxWidth = 0 ## maximum size of the scale-invariant portion of all spots
|
||||
self.opts = {'pxMode': True, 'useCache': True} ## If useCache is False, symbols are re-drawn on every paint.
|
||||
self.opts = {'pxMode': True, 'useCache': True, 'exportMode': False} ## If useCache is False, symbols are re-drawn on every paint.
|
||||
|
||||
self.setPen(200,200,200, update=False)
|
||||
self.setBrush(100,100,150, update=False)
|
||||
@ -664,10 +664,14 @@ class ScatterPlotItem(GraphicsObject):
|
||||
rect = QtCore.QRectF(y, x, h, w)
|
||||
self.fragments.append(QtGui.QPainter.PixmapFragment.create(pos, rect))
|
||||
|
||||
def setExportMode(self, enabled, opts):
|
||||
self.opts['exportMode'] = enabled
|
||||
|
||||
|
||||
def paint(self, p, *args):
|
||||
#p.setPen(fn.mkPen('r'))
|
||||
#p.drawRect(self.boundingRect())
|
||||
if self.opts['pxMode']:
|
||||
if self.opts['pxMode'] is True:
|
||||
atlas = self.fragmentAtlas.getAtlas()
|
||||
#arr = fn.imageToArray(atlas.toImage(), copy=True)
|
||||
#if hasattr(self, 'lastAtlas'):
|
||||
@ -681,7 +685,7 @@ class ScatterPlotItem(GraphicsObject):
|
||||
|
||||
p.resetTransform()
|
||||
|
||||
if not USE_PYSIDE and self.opts['useCache']:
|
||||
if not USE_PYSIDE and self.opts['useCache'] and self.opts['exportMode'] is False:
|
||||
p.drawPixmapFragments(self.fragments, atlas)
|
||||
else:
|
||||
for i in range(len(self.data)):
|
||||
|
Loading…
Reference in New Issue
Block a user