Newer
Older
#coding=utf8
########################################################################
### ###
### Created by Martin Genet, 2016 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
########################################################################
import dolfin
import numpy
import myVTKPythonLibrary as myVTK
########################################################################
def getScalingFactor(scalar_type_as_string):
if (scalar_type_as_string == 'unsigned char' ): return float(2**8 -1)
elif (scalar_type_as_string == 'unsigned short'): return float(2**16-1)
elif (scalar_type_as_string == 'unsigned int' ): return float(2**32-1)
elif (scalar_type_as_string == 'unsigned long' ): return float(2**64-1)
elif (scalar_type_as_string == 'float' ): return 1.
elif (scalar_type_as_string == 'double' ): return 1.
else: assert (0), "Wrong image scalar type. Aborting."
class ExprIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
self.X[0:2] = X[0:2]
Martin Genet
committed
self.interpolator.Interpolate(self.X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprIm3(dolfin.Expression):
def __init__(self, filename, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
Martin Genet
committed
self.interpolator.Interpolate(X, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (2,)
Martin Genet
committed
def eval(self, Expr, X):
self.X[0:2] = X[0:2]
Martin Genet
committed
self.interpolator.Interpolate(self.X, Expr)
Expr /= self.s
class ExprGradIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (3,)
Martin Genet
committed
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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def eval(self, Expr, X):
self.interpolator.Interpolate(X, Expr)
Expr /= self.s
class ExprHessIm2(dolfin.Expression):
def __init__(self, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.X = numpy.array([float()]*2+[Z])
self.H_vec = numpy.empty(4)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.X[0:2] = X[0:2]
self.interpolator.Interpolate(self.X, self.H_vec)
cvec4_to_mat22(self.H_vec, Expr)
Expr /= self.s
class ExprHessIm3(dolfin.Expression):
def __init__(self, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.H_vec = numpy.empty(9)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.interpolator.Interpolate(X, self.H_vec)
cvec9_to_mat33(self.H_vec, Expr)
Expr /= self.s
class ExprDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
Martin Genet
committed
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
Martin Genet
committed
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (2,)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[0:2] = X[0:2] + self.UX[0:2]
Martin Genet
committed
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprGradDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
Martin Genet
committed
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageGradient(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (3,)
Martin Genet
committed
def eval(self, Expr, X):
#print "Expr"
#print " U = " + str(U.vector().array())
#print " X = " + str(X)
#print " UX = " + str(self.UX)
self.U.eval(self.UX, X)
#print " UX = " + str(self.UX)
self.x[:] = X + self.UX
#print " x = " + str(self.x)
Martin Genet
committed
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
#print " Expr = " + str(Expr)
self.interpolator.Interpolate(self.x, Expr)
#print " Expr = " + str(Expr)
Expr /= self.s
#print " Expr = " + str(Expr)
class ExprHessDefIm2(dolfin.Expression):
def __init__(self, U, filename=None, Z=0., **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(2)
self.x = numpy.array([float()]*2+[Z])
self.H_vec = numpy.empty(4)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (2,2)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[0:2] = X[0:2] + self.UX[0:2]
self.interpolator.Interpolate(self.x, self.H_vec)
cvec4_to_mat22(self.H_vec, Expr)
Expr /= self.s
class ExprHessDefIm3(dolfin.Expression):
def __init__(self, U, filename=None, **kwargs):
if filename is not None:
self.init_image(filename=filename)
self.U = U
self.UX = numpy.empty(3)
self.x = numpy.empty(3)
self.H_vec = numpy.empty(9)
def init_image(self, filename):
self.image = myVTK.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
self.image.GetScalarTypeAsString())
self.image = myVTK.computeImageHessian(
image=self.image,
verbose=0)
self.interpolator = myVTK.createImageInterpolator(
image=self.image,
verbose=0)
def value_shape(self):
return (3,3)
def eval(self, Expr, X):
self.U.eval(self.UX, X)
self.x[:] = X + self.UX
self.interpolator.Interpolate(self.x, self.H_vec)
cvec9_to_mat33(self.H_vec, Expr)
Expr /= self.s
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
cppcode_ExprIm='''
#include <vtkSmartPointer.h>
#include <vtkXMLImageDataReader.h>
#include <vtkImageData.h>
#include <vtkImageInterpolator.h>
namespace dolfin
{
class MyFun : public Expression
{
public:
MyFun(): Expression()
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
reader->SetFileName("<<filename>>");
reader->Update();
interpolator->Initialize(reader->GetOutput());
interpolator->Update();
};
void eval(Array<double>& values, const Array<double>& x) const
{
interpolator->Interpolate(x.data(), values.data());
}
private:
vtkSmartPointer<vtkImageInterpolator> interpolator = vtkSmartPointer<vtkImageInterpolator>::New();
};
}'''