converter.py
Go to the documentation of this file.
1 '''
2 A module to produce paraview / vtk objects
3 '''
4 
5 import math
6 import numpy
7 from collections import defaultdict
8 
9 
10 def extrude(pts, dx):
11  '''
12  make a 3d set of cells based on ring of pts extruded along X axis by dx
13 
14  Return points and "relative cells"
15 
16 
17  '''
18  pts2 = [ [pt[0]+dx,pt[1],pt[2]] for pt in pts] # the other face
19  all_pts = pts + pts2
20 
21  n = len(pts)
22  top_cell = range(n)
23  bot_cell = range(n, 2*n)
24  cells = [top_cell, bot_cell]
25 
26  # enumerate the sides
27  for ind in range(n):
28  ind2 = (ind+1)%n
29  cell = [top_cell[ind], top_cell[ind2], bot_cell[ind2], bot_cell[ind]]
30  cells.append(cell)
31 
32  return all_pts, cells
33 
34 
35 
36 def orderpoints(pointset):
37  c = [0.0,0.0,0.0]
38  for p in pointset:
39  for i in range(3):
40  c[i] += p[i]
41  n = len(pointset)
42  for i in range(3):
43  c[i] /= n
44 
45  byang = list()
46  for p in pointset:
47  ang = math.atan2(p[2]-c[2], p[1]-c[1]);
48  byang.append((ang, p))
49  byang.sort()
50  return [p for a,p in byang]
51 
52 
53 def depos2pts(arr):
54  '''
55  Convert numpy array like which comes from 'depo_data_0' key of npz
56  file from NumpyDepoSaver to tvtk unstructured grid.
57  '''
58  from tvtk.api import tvtk, write_data
59 
60  npts = arr.shape[1]
61  # t,q,x,y,z,dl,dt
62  q = arr[1,:].reshape(npts)
63  pts = arr[2:5,:].T
64 
65  indices = list(range(npts))
66 
67  ret = tvtk.PolyData(points=pts)
68  verts = numpy.arange(0, npts, 1)
69  verts.shape = (npts,1)
70  ret.verts = verts
71  ret.point_data.scalars = indices[:npts]
72  ret.point_data.scalars.name = 'indices'
73 
74  ret.point_data.add_array(q)
75  ret.point_data.get_array(1).name = 'charge'
76 
77  return ret
78 
79 
80 
82  '''
83  Given a graph object return a tvtk data object with blbos.
84  '''
85  from tvtk.api import tvtk, write_data
86 
87  all_points = list()
88  blob_cells = list()
89  datasetnames = set()
90  values = list()
91  for node, ndata in gr.nodes.data():
92  if ndata['code'] != 'b':
93  continue;
94  vals = dict();
95  thickness = 1.0
96  for key,val in ndata.items():
97  if key == 'corners':
98  pts = orderpoints(val)
99  continue
100  if key == 'span':
101  thickness = val
102  continue
103  if key == 'code':
104  continue
105  datasetnames.add(key)
106  vals[key] = val;
107  pts,cells = extrude(pts, thickness)
108  all_points += pts
109  blob_cells.append((len(pts), cells))
110  values.append(vals)
111 
112  ugrid = tvtk.UnstructuredGrid(points = all_points);
113  ptype = tvtk.Polyhedron().cell_type
114  offset = 0
115  for npts,cells in blob_cells:
116  cell_ids = [len(cells)]
117  for cell in cells:
118  cell_ids.append(len(cell))
119  cell_ids += [offset+cid for cid in cell]
120  ugrid.insert_next_cell(ptype, cell_ids)
121  offset += npts
122 
123  ugrid.cell_data.scalars = list(range(len(values)))
124  ugrid.cell_data.scalars.name = "indices"
125 
126  narrays = 1
127  for datasetname in sorted(datasetnames):
128  arr = numpy.asarray([vals.get(datasetname, 0.0) for vals in values], dtype=float)
129  ugrid.cell_data.add_array(arr)
130  ugrid.cell_data.get_array(narrays).name = datasetname
131  narrays += 1
132 
133  return ugrid
134 
135 
136 
137 def get_blob(gr,node):
138  for other in gr[node]:
139  odat = gr.nodes[other]
140  if odat['code'] == 'b':
141  return other
142  return None
143 
144 def get_slice(gr, bnode):
145  for other in gr[bnode]:
146  odat = gr.nodes[other]
147  if odat['code'] == 's':
148  return other
149  return None
150 
151 
153  from tvtk.api import tvtk, write_data
154 
155  class Perwpid:
156  def __init__(self):
157  self.allind = list()
158  self.allchs = list()
159  self.values = list()
160  perwpid = defaultdict(Perwpid)
161  for node, ndata in gr.nodes.data():
162  if ndata['code'] != 'm':
163  continue;
164  bnode = get_blob(gr, node)
165  if bnode is None:
166  raise ValueError("bad graph structure")
167  #x = gr.nodes[bnode]['corners'][0][0] # "t"
168  #eckses.append(x)
169  snode = get_slice(gr, bnode)
170  if snode is None:
171  raise ValueError("bad graph structure")
172  wpid = ndata['wpid']
173 
174 
175  sdat = gr.nodes[snode]
176  snum = sdat['ident']
177  perwpid[wpid].allind.append(snum)
178  sact = sdat['activity']
179  val = 0.0
180  chids = ndata['chids']
181  perwpid[wpid].allchs += chids
182  for chid in chids:
183  val += float(sact[str(chid)])
184  perwpid[wpid].values.append((snum,chids,val))
185  all_imgdat=dict()
186  for wpid, dat in perwpid.items():
187  smin = min(dat.allind)
188  smax = max(dat.allind)
189  cmin = min(dat.allchs)
190  cmax = max(dat.allchs)
191  arr = numpy.zeros((smax-smin+1, cmax-cmin+1))
192  for ind,chids,val in dat.values:
193  for ch in chids:
194  arr[ind - smin, ch - cmin] += val
195  imgdat = tvtk.ImageData(spacing=(1,1,1), origin=(0,0,0))
196  imgdat.point_data.scalars = arr.T.flatten()
197  imgdat.point_data.scalars.name = 'activity'
198  imgdat.dimensions = list(arr.shape)+[1]
199  all_imgdat[wpid] = imgdat
200  return all_imgdat
201 
202 def blob_center(bdat):
203  '''
204  Return an array of one point at the center of the blob
205  '''
206  thickness = bdat['span']
207  value = bdat['value']
208  arr = numpy.asarray(bdat['corners'])
209 
210  npts = arr.shape[0]
211  center = numpy.array([0.0]*4, dtype=float)
212  center[:3] = numpy.sum(arr, 0) / npts
213  center[0] += 0.5*thickness
214  center[3] = value;
215  return center
216 
217 
218 def blob_uniform_sample(bdat, density):
219  '''
220  Return an array of points uniformly sampled in the blob
221  '''
222  import random
223  from shapely.geometry import Polygon, Point
224  thickness = bdat['span']
225  value = bdat['value']
226 
227  # z is x, y is y
228  xstart = bdat['corners'][0][0]
229  corners = [(cp[2],cp[1]) for cp in orderpoints(bdat['corners'])]
230  npts = len(corners);
231  pgon = Polygon(corners)
232  nwant = max(1, int(pgon.area * thickness * density))
233  pts = list()
234  min_x, min_y, max_x, max_y = pgon.bounds
235 
236  while len(pts) != nwant:
237  p = Point([random.uniform(min_x, max_x), random.uniform(min_y, max_y)])
238  if (p.within(pgon)):
239  pts.append([random.uniform(xstart, xstart+thickness), p.y, p.x, value/nwant]);
240  return numpy.asarray(pts)
241 
242 
243 
244 def blobpoints(gr, sample_method=blob_center):
245  '''
246  return Nx4 array with rows made of x,y,z,q.
247  '''
248  arr = None
249 
250  for node, ndata in gr.nodes.data():
251  if ndata['code'] != 'b':
252  continue;
253  one = sample_method(ndata)
254  if arr is None:
255  arr = one
256  else:
257  arr = numpy.vstack((arr, one))
258 
259  return arr
def blob_uniform_sample(bdat, density)
Definition: converter.py:218
def get_blob(gr, node)
Definition: converter.py:137
def blobpoints(gr, sample_method=blob_center)
Definition: converter.py:244
def get_slice(gr, bnode)
Definition: converter.py:144
def extrude(pts, dx)
Definition: converter.py:10
static int max(int a, int b)
def orderpoints(pointset)
Definition: converter.py:36
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
static QCString str