"""
PySCeS - Python Simulator for Cellular Systems (http://pysces.sourceforge.net)
Copyright (C) 2004-2009 B.G. Olivier, J.M. Rohwer, J.-H.S Hofmeyr all rights reserved,
Brett G. Olivier (bgoli@users.sourceforge.net)
Triple-J Group for Molecular Cell Physiology
Stellenbosch University, South Africa.
Permission to use, modify, and distribute this software is given under the
terms of the PySceS (BSD style) license. See LICENSE.txt that came with
this distribution for specifics.
NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
Brett G. Olivier
"""
import numpy, scipy, scipy.io
class DataTools:
"""
Data manipulation tools, requires scipy, scipy.io
"""
def data_filter_func(row):
print 'You should overwrite this (data_filter_func) method with your own or supply a function to data_filter'
return False
def DataFilter(self, data, filter_func=None, report=True):
"""Takes in an array and a function which returns a boolean for filter==True"""
if filter_func == None:
filter_func = self.data_filter_func
data_out = numpy.zeros(data.shape,'d')
rr = 0
filtered = []
for r in range(data_out.shape[0]):
if filter_func(data[r,:]):
filtered.append('Filtered: %s' % data[r,:])
else:
data_out[rr,:] = data[r,:]
rr += 1
# 'remove filtered rows'
data_out = data_out[:rr,:]
if report:
for a in filtered:
print(a)
return data_out
def GridSortLR(self, arr):
"""
Sort irregular unstructured grid data into monotonically rising grid data
(where possible) Work from left column to right (uses ultra cool NumPy functions :-)
"""
print 'arr.shape:', arr.shape
sorted = numpy.rot90(arr)
print 'rotated.shape:', sorted.shape
sorted_idx = numpy.lexsort(sorted)
print 'sorted_idx.shape', sorted_idx.shape
sorted = sorted.take(sorted_idx, axis=-1)
sorted = numpy.flipud(sorted)
arr = numpy.transpose(sorted)
print 'arr.shape:', arr.shape
return arr
def writeGnuPlot_3D(self, arr, fname):
if arr.shape[0] > 5e6:
print "warning: writeGnuPlot_3D does not work for large arrays"
try:
scipy.io.write_array(fname+'_result_gplt.dat', arr, precision=9, separator=' ')
except Exception, ex:
print 'writeGnuPlot_3D exception raised'
print ex
def vtk_scalar_func(self, val):
"""
Manipulate the scalar (or z-axis) values
Accepts and returns a double
"""
return val
def writeVTK_UnstructuredGrid(self, arr, fname, scalar_func=None):
"""
Write arr as a VTK unstructured grid to file "fname.vtk"
arr must either be in:
x,y,z format (in which case z is used for scalar data) or
z,y,z,s format (s being point scalar value)
Scalar_func(val) function to manipulate the scalar (z-axis) values
Original version http://mayavi.sourceforge.net/mwiki/PointClouds
Hacked quite a bit by brett 20070221
"""
assert arr.shape[1] == 3 or arr.shape[1] == 4, '\nneed 3 or 4 columns for this'
if scalar_func == None:
scalar_func = self.vtk_scalar_func
if arr.shape[1] == 4:
HAVE_SCALARS = 1
else:
HAVE_SCALARS = 0
print 'No scalar values supplied Z axis values will be used'
n=arr.shape[0]
print "n:",n
# write data to vtk polydata file
# write header
out = file(fname+'.vtk', 'w')
h1 = "# vtk DataFile Version 2.0\n"
h1 += "%s\n" % fname
h1 += "ASCII\n"
h1 += "DATASET UNSTRUCTURED_GRID\n"
h1 += "POINTS " + str(n) + " double\n"
out.write(h1)
# write xyz data
for r in range(n):
#s = '%15.2f %15.2f %15.2f' % (x[i], y[i], z[i])
out.write(str(arr[r,0])+" "+str(arr[r,1])+" "+str(arr[r,2])+'\n')
# write cell data
out.write("CELLS "+ str(n)+ " "+ str(2*n)+'\n')
for r in range(n):
#s = '1 %d \n' % (i)
out.write("1 "+str(r)+"\n")
# write cell types
out.write("CELL_TYPES " + str(n)+'\n')
for r in range(n):
out.write("1 \n")
# write z scalar values
h2 = '\n' + """POINT_DATA """ + str(n) + "\n"
h3 = "SCALARS %s double 1\n" % fname
h3 += "LOOKUP_TABLE default\n"
out.write(h2 + h3)
for r in range(n):
if HAVE_SCALARS:
sc=(scalar_func(arr[r,3]))
else:
sc=(scalar_func(arr[r,2]))
out.write(str(sc)+ "\n")
out.write('\n')
out.close()
class IntersectionAnalysis(DataTools):
zmin = 0.01
intersect_tol = 0.04
coords = None
vtkcoords = None
intersection = None
vtkintersection = None
def setCoords(self, coords):
self.coords = coords
def setVTKCoords(self):
self.vtkcoords = numpy.zeros((self.coords.shape[0], self.coords.shape[1]+1),'d')
self.vtkcoords[:,:self.coords.shape[1]] = self.coords.copy()
def data_filter_func(self, row):
val = False
## if (numpy.abs(row) < 1.0e-10).any() or (row < 0.0).any():
if (row < self.zmin).any():
val = True
return val
def CalculateIntersect(self, tol=None):
"""calc and add scalar value for intersection to vtkcoords"""
if tol != None:
self.intersect_tol = tol
intersect = []
cntr = 0
for r in range(self.coords.shape[0]):
if numpy.abs(self.coords[r, 2] - self.coords[r, 3]) < self.intersect_tol:
intersect.append(self.coords[r,:])
cntr += 1
print '\n*****\n%s intersections found.\n*****\n' % cntr
self.intersection = numpy.array(intersect)
self.vtkintersection = self.intersection.copy()
del intersect
def LogVTKdata(self):
print "Logging VTK coords"
# log for vtk
for r in range(self.vtkcoords.shape[0]):
self.vtkcoords[r,:self.coords.shape[1]] = scipy.log10(self.vtkcoords[r,:self.coords.shape[1]])
def LogVTKIntersectiondata(self):
print "Logging VTK intersection"
# log for vtk intersect
for r in range(self.vtkintersection.shape[0]):
self.vtkintersection[r,:] = scipy.log10(self.intersection[r,:])
|