Newer
Older
################################################################################
### ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
################################################################################
################################################################################
def get_ExprIm_cpp(
im_dim,
im_type="im",
im_is_def=0,
verbose=0):
assert (im_dim in (2,3))
assert (im_type in ("im","grad"))
#include <string.h>
#include <vtkSmartPointer.h>
#include <vtkStructuredPointsReader.h>
#include <vtkXMLImageDataReader.h>
#include <vtkImageData.h>'''+('''
#include <vtkImageGradient.h>''')*(im_type=="grad")+'''
#include <vtkImageInterpolator.h>
double getStaticScalingFactor(const char* scalar_type_as_string)
{
if (strcmp(scalar_type_as_string, "unsigned char" ) == 0) return pow(2, 8)-1;
if (strcmp(scalar_type_as_string, "unsigned short") == 0) return pow(2, 16)-1;
if (strcmp(scalar_type_as_string, "unsigned int" ) == 0) return pow(2, 32)-1;
if (strcmp(scalar_type_as_string, "unsigned long" ) == 0) return pow(2, 64)-1;
if (strcmp(scalar_type_as_string, "float" ) == 0) return 1.;
if (strcmp(scalar_type_as_string, "double" ) == 0) return 1.;
assert (0);
}
vtkSmartPointer<vtkImageInterpolator> interpolator;
double static_scaling;'''+('''
Array<double>* dynamic_scaling; // does not work
double dynamic_scaling_a; // should not be needed
double dynamic_scaling_b; // should not be needed
mutable Array<double> UX;
mutable Array<double> x;''')*(im_is_def)+('''
mutable Array<double> X3D;''')*(not im_is_def)*(im_dim==2)+'''
Expression('''+str(im_dim)*(im_type=="grad")+''')'''+(''',
Martin Genet
committed
dynamic_scaling_a(1.), // should not be needed
dynamic_scaling_b(0.), // should not be needed
UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''',
X3D(3)''')*(not im_is_def)*(im_dim==2)+'''
void init_dynamic_scaling(
const Array<double> &scaling)
{
//dynamic_scaling = scaling; // does not work
// std::cout << "dynamic_scaling = " << dynamic_scaling->str(1) << std::endl; // does not work
dynamic_scaling_a = scaling[0]; // should not be needed
dynamic_scaling_b = scaling[1]; // should not be needed
void update_dynamic_scaling( // should not be needed
const Array<double> &scaling) // should not be needed
{ // should not be needed
dynamic_scaling_a = scaling[0]; // should not be needed
dynamic_scaling_b = scaling[1]; // should not be needed
} // should not be needed
void init_disp(
Function* UU)
{
U = UU;
}''')*(im_is_def)+'''
void init_image(
const char* filename,
const char* interpol_mode="'''+('''nearest''')*(im_type=="im")+('''nearest''')*(im_type=="grad")+'''",
Martin Genet
committed
const double &interpol_out_value='''+('''0.''')*(im_type=="im")+('''0.''')*(im_type=="grad")+(''',
const double &Z=0.''')*(im_dim==2)+''')
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
reader->SetFileName(filename);
reader->Update();
static_scaling = getStaticScalingFactor(reader->GetOutput()->GetScalarTypeAsString());'''+('''
vtkSmartPointer<vtkImageGradient> gradient = vtkSmartPointer<vtkImageGradient>::New();
#if VTK_MAJOR_VERSION >= 6
gradient->SetInputData(reader->GetOutput());
#else
gradient->SetInput(reader->GetOutput());
#endif
gradient->SetDimensionality('''+str(im_dim)+''');
gradient->Update();''')*(im_type=="grad")+'''
interpolator = vtkSmartPointer<vtkImageInterpolator>::New();
if (strcmp(interpol_mode, "nearest") == 0)
interpolator->SetInterpolationModeToNearest();
}
else if (strcmp(interpol_mode, "linear") == 0)
interpolator->SetInterpolationModeToLinear();
else if (strcmp(interpol_mode, "cubic") == 0)
{
interpolator->SetInterpolationModeToCubic();
}
else
{
std::cout << "Interpolator interpol_mode (" << interpol_mode << ") must be \\"nearest\\", \\"linear\\" or \\"cubic\\". Aborting." << std::endl;
assert(0);
}
interpolator->SetOutValue(interpol_out_value);
interpolator->Initialize('''+('''reader->GetOutput()''')*(im_type=="im")+('''gradient->GetOutput()''')*(im_type=="grad")+''');
interpolator->Update();'''+('''
x[2] = Z;''')*(im_is_def)*(im_dim==2)+('''
X3D[2] = Z;''')*(not im_is_def)*(im_dim==2)+'''
void eval(Array<double>& expr, const Array<double>& X) const
std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+('''
U->eval(UX, X);'''+('''
std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+('''
x[0] = X[0] + UX[0];
x[1] = X[1] + UX[1];''')*(im_dim==2)+('''
x[0] = X[0] + UX[0];
x[1] = X[1] + UX[1];
x[2] = X[2] + UX[2];''')*(im_dim==3)+('''
std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+'''
interpolator->Interpolate(x.data(), expr.data());''')*(im_is_def)+('''
X3D[0] = X[0];
X3D[1] = X[1];'''+('''
std::cout << "X3D = " << X3D.str(1) << std::endl;''')*(verbose)+'''
interpolator->Interpolate(X3D.data(), expr.data());''')*(not im_is_def)*(im_dim==2)+('''
interpolator->Interpolate(X.data(), expr.data());''')*(not im_is_def)*(im_dim==3)+('''
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+('''
expr[0] /= static_scaling;''')*(im_type=="im")+('''
expr[0] /= static_scaling;
expr[1] /= static_scaling;''')*(im_type=="grad")*(im_dim==2)+('''
expr[0] /= static_scaling;
expr[1] /= static_scaling;
expr[2] /= static_scaling;''')*(im_type=="grad")*(im_dim==3)+('''
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+(('''
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
// std::cout << "in (im)" << std::endl;
// std::cout << "expr = " << expr.str(1) << std::endl;
// expr[0] *= (*dynamic_scaling)[0]; // does not work
// expr[0] += (*dynamic_scaling)[1]; // does not work
expr[0] *= dynamic_scaling_a; // should not be needed
expr[0] += dynamic_scaling_b; // should not be needed
// std::cout << "expr = " << expr.str(1) << std::endl;
// std::cout << "out (im)" << std::endl;''')*(im_type=="im")+('''
// std::cout << "in (grad)" << std::endl;
// std::cout << "expr = " << expr.str(1) << std::endl;
// expr[0] *= (*dynamic_scaling)[0]; // does not work
// expr[1] *= (*dynamic_scaling)[0]; // does not work
expr[0] *= dynamic_scaling_a; // should not be needed
expr[1] *= dynamic_scaling_a; // should not be needed
// std::cout << "expr = " << expr.str(1) << std::endl;
// std::cout << "out (grad)" << std::endl;''')*(im_type=="grad")*(im_dim==2)+('''
// std::cout << "in (grad)" << std::endl;
// std::cout << "expr = " << expr.str(1) << std::endl;
// expr[0] *= (*dynamic_scaling)[0]; // does not work
// expr[1] *= (*dynamic_scaling)[0]; // does not work
// expr[2] *= (*dynamic_scaling)[0]; // does not work
expr[0] *= dynamic_scaling_a; // should not be needed
expr[1] *= dynamic_scaling_a; // should not be needed
expr[2] *= dynamic_scaling_a; // should not be needed
// std::cout << "expr = " << expr.str(1) << std::endl;
// std::cout << "out (grad)" << std::endl;''')*(im_type=="grad")*(im_dim==3)+('''
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose))*(im_is_def)+'''
Martin Genet
committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
def get_ExprCharFuncIm_cpp(
im_dim,
verbose=0):
assert (im_dim in (2,3))
ExprCharFuncIm_cpp = '''\
#include <string.h>
#include <vtkSmartPointer.h>
#include <vtkXMLImageDataReader.h>
#include <vtkImageData.h>
namespace dolfin
{
class MyExpr : public Expression
{
double xmin, xmax, ymin, ymax, zmin, zmax;
public:
MyExpr():
Expression()
{
}
void init_image(
const char* filename)
{
vtkSmartPointer<vtkXMLImageDataReader> reader = vtkSmartPointer<vtkXMLImageDataReader>::New();
reader->SetFileName(filename);
reader->Update();
vtkSmartPointer<vtkImageData> image = reader->GetOutput();
double* bounds = image->GetBounds();
xmin = bounds[0];
xmax = bounds[1];
ymin = bounds[2];
ymax = bounds[3];
zmin = bounds[4];
zmax = bounds[5];'''+('''
std::cout << "bounds = " << bounds[0] << " " << bounds[1] << " " << bounds[2] << " " << bounds[3] << " " << bounds[4] << " " << bounds[5] << std::endl;
std::cout << "xmin = " << xmin << std::endl;
std::cout << "xmax = " << xmax << std::endl;
std::cout << "ymin = " << ymin << std::endl;
std::cout << "ymax = " << ymax << std::endl;
std::cout << "zmin = " << zmin << std::endl;
std::cout << "zmax = " << zmax << std::endl;''')*(verbose)+'''
}
void eval(Array<double>& expr, const Array<double>& X) const
{'''+('''
std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+('''
if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax))''')*(im_dim==2)+('''
if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax) && (X[2] >= zmin) && (X[2] <= zmax))''')*(im_dim==3)+'''
{
expr[0] = 1.;
}
else
{
expr[0] = 0.;
}'''+('''
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+'''
}
};
}'''
#print ExprCharFuncIm_cpp
return ExprCharFuncIm_cpp