-
Notifications
You must be signed in to change notification settings - Fork 0
/
utilities.py
169 lines (138 loc) · 6.6 KB
/
utilities.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 3 17:09:56 2014
@author: robinsonm
"""
import csv
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import array
import vtk
def import_columns(filename):
reader = csv.reader(open(filename), delimiter=' ',skipinitialspace=True)
first_line = reader.next()
print first_line
decrement = 0
if (first_line[0][0]=='#'):
decrement = decrement+1
if (first_line[-1]==''):
decrement = decrement+1
n = len(first_line)-decrement
columns = [array.array('f') for _ in xrange(n)]
if first_line[0][0]!='#':
for i in xrange(n):
columns[i].append(float(first_line[i]))
for line in reader:
if line[0][0]=='#':
continue
for i in xrange(n):
columns[i].append(float(line[i]))
for i in xrange(n):
columns[i] = np.array(columns[i])
return columns
def mesh(x,y,z):
nx = 100
ny = 100
xmin = np.min(x)
xmax = np.max(x)
ymin = np.min(y)
ymax = np.max(y)
xi = np.linspace(xmin, xmax, nx)
yi = np.linspace(ymin, ymax, ny)
xi, yi = np.meshgrid(xi, yi)
x_new = (x - xmin) / (xmax - xmin)
xi_new = (xi - xmin) / (xmax - xmin)
y_new = (y - ymin) / (ymax - ymin)
yi_new = (yi - ymin) / (ymax - ymin)
print x_new.shape
print y_new.shape
print z.shape
print xi_new.shape
print yi_new.shape
zi = plt.mlab.griddata(x_new, y_new, z, xi_new, yi_new)
return xi,yi,zi
def volume_plot(input_data_matrix):
nx,ny,nz = input_data_matrix.shape
min_data = input_data_matrix.min()
max_data = input_data_matrix.max()
print 'volume plot of data matrix with shape ',input_data_matrix.shape,' and (min,max) = (',min_data,',',max_data,")"
data_matrix = np.uint8(255.0*(input_data_matrix-min_data)/(max_data-min_data))
# data_matrix[0:35, 0:35, 0:35] = 50
# data_matrix[25:55, 25:55, 25:55] = 100
# data_matrix[45:69, 45:69, 45:69] = 150
nx,ny,nz = data_matrix.shape
min_data = data_matrix.min()
max_data = data_matrix.max()
#print 'volume plot of data matrix with shape ',data_matrix.shape,' and (min,max) = (',min_data,',',max_data,")"
# For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which
# imports raw data and stores it.
dataImporter = vtk.vtkImageImport()
# The preaviusly created array is converted to a string of chars and imported.
data_string = data_matrix.tostring()
dataImporter.CopyImportVoidPointer(data_string, len(data_string))
# The type of the newly imported data is set to unsigned char (uint8)
dataImporter.SetDataScalarTypeToUnsignedChar()
# Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer
# must be told this is the case.
dataImporter.SetNumberOfScalarComponents(1)
# The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this
# simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case.
# I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although
# VTK complains if not both are used.
dataImporter.SetDataExtent(0, nx-1, 0, ny-1, 0, nz-1)
dataImporter.SetWholeExtent(0, nx-1, 0, ny-1, 0, nz-1)
# The following class is used to store transparencyv-values for later retrival. In our case, we want the value 0 to be
# completly opaque whereas the three different cubes are given different transperancy-values to show how it works.
alphaChannelFunc = vtk.vtkPiecewiseFunction()
alphaChannelFunc.AddPoint(min_data, 0.0)
alphaChannelFunc.AddPoint(min_data+1, 0.2)
#alphaChannelFunc.AddPoint(0.5*(max_data+min_data), 0.1)
#alphaChannelFunc.AddPoint(max_data, 1.0)
# alphaChannelFunc.AddPoint(0, 0.0)
# alphaChannelFunc.AddPoint(50, 0.05)
# alphaChannelFunc.AddPoint(100, 0.1)
# alphaChannelFunc.AddPoint(150, 0.2)
# This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes
# to be of the colors red green and blue.
colorFunc = vtk.vtkColorTransferFunction()
print min_data,max_data
colorFunc.AddRGBPoint(min_data, 1.0, 0.0, 0.0)
colorFunc.AddRGBPoint(max_data, 0.0, 0.0, 1.0)
# The preavius two classes stored properties. Because we want to apply these properties to the volume we want to render,
# we have to store them in a class that stores volume prpoperties.
volumeProperty = vtk.vtkVolumeProperty()
volumeProperty.SetColor(colorFunc)
volumeProperty.SetScalarOpacity(alphaChannelFunc)
# This class describes how the volume is rendered (through ray tracing).
compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
# We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered.
volumeMapper = vtk.vtkVolumeRayCastMapper()
volumeMapper.SetVolumeRayCastFunction(compositeFunction)
volumeMapper.SetInputConnection(dataImporter.GetOutputPort())
# The class vtkVolume is used to pair the preaviusly declared volume as well as the properties to be used when rendering that volume.
volume = vtk.vtkVolume()
volume.SetMapper(volumeMapper)
volume.SetProperty(volumeProperty)
# With almost everything else ready, its time to initialize the renderer and window, as well as creating a method for exiting the application
renderer = vtk.vtkRenderer()
renderWin = vtk.vtkRenderWindow()
renderWin.AddRenderer(renderer)
renderInteractor = vtk.vtkRenderWindowInteractor()
renderInteractor.SetRenderWindow(renderWin)
# We add the volume to the renderer ...
renderer.AddVolume(volume)
# ... set background color to white ...
renderer.SetBackground(1, 1, 1)
# ... and set window size.
renderWin.SetSize(400, 400)
# A simple function to be called when the user decides to quit the application.
def exitCheck(obj, event):
if obj.GetEventPending() != 0:
obj.SetAbortRender(1)
# Tell the application to use the function as an exit check.
renderWin.AddObserver("AbortCheckEvent", exitCheck)
renderInteractor.Initialize()
# Because nothing will be rendered without any input, we order the first render manually before control is handed over to the main-loop.
renderWin.Render()
renderInteractor.Start()