diff --git a/Energy.py b/Energy.py
index a645ef42bb0639043bb6ac83b2253ad86360a0c6..f0502d8cec80554eb0e0717feeafb13945992c2d 100644
--- a/Energy.py
+++ b/Energy.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Energy_Regularization.py b/Energy_Regularization.py
index 651d89d813210528acfebed45ec15370954f6b26..1596d14bba92707c0faa4cca6839ff4c6295b62d 100644
--- a/Energy_Regularization.py
+++ b/Energy_Regularization.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Energy_WarpedImage.py b/Energy_WarpedImage.py
index 454865a43e676fe073247e1a9de1c554daeead63..94113039c4a2d02de32a012b73895701e88df701 100644
--- a/Energy_WarpedImage.py
+++ b/Energy_WarpedImage.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/ImageIterator.py b/ImageIterator.py
index e650c94a79d970e8f14af73361e1ed66a9e0ff24..ae45597f9583554fd2805b8aff645da4f233010a 100644
--- a/ImageIterator.py
+++ b/ImageIterator.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/ImageSeries.py b/ImageSeries.py
index 89eea3b992c6fd2b86b054d77f3fe6051043f2ef..c3c219f2d6ea2720acbdbccc64258493c04e50a7 100644
--- a/ImageSeries.py
+++ b/ImageSeries.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/NonlinearSolver.py b/NonlinearSolver.py
index 86fc7d62c1b758a2e402db7034be537e47fe665c..a07175117b3af7c9f4b8df37bd41f4f19ba7bc06 100644
--- a/NonlinearSolver.py
+++ b/NonlinearSolver.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Problem.py b/Problem.py
index cdbba9bc1076127ebeb4b5fd8467c4274c38d03a..a4028765883fc4571e4de301aaf64dd3830ac8bb 100644
--- a/Problem.py
+++ b/Problem.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/Problem_ImageRegistration.py b/Problem_ImageRegistration.py
index 36fe4712a5c2cdd4f1c682a36f3b6e46814829d3..40e31e5da8f8dc31a45287d5ca598bc95fd39c87 100644
--- a/Problem_ImageRegistration.py
+++ b/Problem_ImageRegistration.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_displacement_error.py b/compute_displacement_error.py
index 8e7b6acef8b31dd91f13539e33ee576436b1f34d..74afebb4d0aececa71b7e612315fc99fb533ea2b 100644
--- a/compute_displacement_error.py
+++ b/compute_displacement_error.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_displacements.py b/compute_displacements.py
index 17332be5a0a4d7f1a45f20ffac4f3aa3705687cc..8b539df001f97dbcfcf6d08b6fdcac08998d24d1 100644
--- a/compute_displacements.py
+++ b/compute_displacements.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_quadrature_degree.py b/compute_quadrature_degree.py
index 78755d55ebbca9ecfc145024a437b94374b02cc4..e507fc5e5b3e9dd93471ec1c476ebf4da6564c94 100644
--- a/compute_quadrature_degree.py
+++ b/compute_quadrature_degree.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/compute_strains.py b/compute_strains.py
index a6c2a27749a7dcffec21da78f3a21d63e1a97f92..a1f37c4c2e7d3290662e14e7378cc12414c41c98 100644
--- a/compute_strains.py
+++ b/compute_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/fedic.py b/fedic.py
index ee505f8193732d74952ae7c3d9f1207472a357af..d7bfac8586cb43c885905babdc8a649bbbb0c1bf 100644
--- a/fedic.py
+++ b/fedic.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/fedic2.py b/fedic2.py
index af560c58da46f45725ea79cafa251276ba141983..cd2e049990f8040a8bf3230f0cb8b033f5da8b6d 100644
--- a/fedic2.py
+++ b/fedic2.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/generate_images.py b/generate_images.py
index 163cc8db67e6600344dad211b3318eedf7739123..ffcf5017505b02462e51095b3336ea617e00016c 100644
--- a/generate_images.py
+++ b/generate_images.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/generate_undersampled_images.py b/generate_undersampled_images.py
index 7deda2065669f79fd0344e841c06dfa30800d3c5..bcd23330727d1434c3800bd410a834aa49fa30b8 100644
--- a/generate_undersampled_images.py
+++ b/generate_undersampled_images.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/generated_grad_expressions_cpp.py b/generated_grad_expressions_cpp.py
new file mode 100644
index 0000000000000000000000000000000000000000..fc64124bbb8141e2c149c1d2feb565aa8c94a9f0
--- /dev/null
+++ b/generated_grad_expressions_cpp.py
@@ -0,0 +1,176 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Martin Genet, 2016-2019 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+################################################################################
+
+def get_ExprGenGrad_cpp(
+ im_dim,
+ im_type="grad",
+ im_is_def=0,
+ disp_type="fenics", # "vtk"
+ verbose=0):
+
+ assert (im_dim in (2,3))
+ assert (im_type in ("grad"))
+
+ ExprGenGrad_cpp = '''\
+#include <string.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkStructuredPointsReader.h>
+#include <vtkXMLImageDataReader.h>
+#include <vtkImageData.h>
+#include <vtkImageGradient.h>
+#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);
+}
+
+namespace dolfin
+{
+
+class MyExpr : public Expression
+{
+ vtkSmartPointer<vtkImageData> image;
+ vtkSmartPointer<vtkImageGradient> gradient;
+ vtkSmartPointer<vtkImageInterpolator> interpolator;
+ double static_scaling;
+ std::shared_ptr<dolfin::Mesh> mesh;
+ std::shared_ptr<dolfin::Function> U;
+ mutable Array<double> UX;
+ mutable Array<double> x;
+
+public:
+
+ MyExpr():
+ Expression('''+str(im_dim)*(im_type in ("grad"))+'''),
+ UX('''+str(im_dim)+'''),
+ x(3),
+ interpolator(vtkSmartPointer<vtkImageInterpolator>::New())
+
+ {
+ }
+
+ void init_mesh(
+ const std::shared_ptr<dolfin::Mesh> mesh_)
+ {
+ mesh = mesh_;
+
+ std::cout << mesh->num_vertices() << std::endl;
+ std::cout << mesh->num_cells() << std::endl;
+ }
+
+ void init_disp(
+ const std::shared_ptr<dolfin::Function> U_)
+ {
+ U = U_;
+ }
+
+ void init_image()
+ {
+ image = vtkSmartPointer<vtkImageData>::New(); // MG20180913: another way to instantiate object
+ gradient = vtkSmartPointer<vtkImageGradient>::New();
+ }
+
+
+
+ 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());'''+('''
+ 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)+'''
+ }
+
+
+ void prepareImage( int n1, int n2, int n3, float origin1, float origin2, float origin3, float spac1, float spac2, float spac3)
+ {
+ image->SetDimensions(n1,n2,n3);
+ image->AllocateScalars(VTK_FLOAT,1);
+
+ image->SetSpacing(spac1, spac2, spac3);
+ image->SetOrigin(origin1, origin2, origin3);
+
+ static_scaling = getStaticScalingFactor("double");
+ }
+
+
+ void assignValues( int i, int j, int k, float value)
+ {
+ image->SetScalarComponentFromFloat(i,j,k,0,value);
+ }
+
+ void setGradient_Interpol()
+ {
+ interpolator->SetInterpolationModeToLinear();
+ interpolator->SetOutValue(.0);
+
+ gradient->SetInputData(image);
+ gradient->SetDimensionality('''+str(im_dim)+''');
+ gradient->Update();
+ interpolator->Initialize(gradient->GetOutput());
+
+ interpolator->Update();
+ }
+
+
+ void prepareImage_2D( int n1, int n2, float origin1, float origin2, float spac1, float spac2)
+ {
+ const double &Z=0.;
+ image->SetExtent(0, n1-1, 0, n2-1, 0, 0);
+
+ image->AllocateScalars(VTK_FLOAT,1);
+
+ image->SetSpacing(spac1, spac2, 1);
+ image->SetOrigin(origin1, origin2, 0);
+
+ static_scaling = getStaticScalingFactor("double");
+ x[2] = Z;
+
+ }
+
+
+ void assignValues_2D( int i, int j, float value)
+ {
+ image->SetScalarComponentFromFloat(i,j,0,0,value);
+ }
+};
+
+}'''
+ #print ExprIm_cpp
+ return ExprGenGrad_cpp
diff --git a/generated_image_expressions_cpp.py b/generated_image_expressions_cpp.py
index 548e3265a8f3ad2429758168613aa11dec7deea7..7a62201923049f3d97cbf369e05ed4c328ffdb37 100644
--- a/generated_image_expressions_cpp.py
+++ b/generated_image_expressions_cpp.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index a360e3283bb9ec7dd0d718145e2cf1df0cad5857..541e41d8ae2f766a5bd96100abd4e9976ce70200 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/image_expressions_py.py b/image_expressions_py.py
index 7e7585bc81d09ab4128da02d85643cfc37fee449..dabb8f138c3cc5a6ba4b0f1cd839cfaa54c4d89b 100644
--- a/image_expressions_py.py
+++ b/image_expressions_py.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_binned_strains_vs_radius.py b/plot_binned_strains_vs_radius.py
index 25bfc6fce9bc2c39451dd575ef8f0e7315fc5edd..8d0033e959bb4783d2e2e2216af647a8e770917b 100644
--- a/plot_binned_strains_vs_radius.py
+++ b/plot_binned_strains_vs_radius.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_displacement_error.py b/plot_displacement_error.py
index 727e80dc5f294f19f764cf04604f0af4a9d4d017..57a575678a2f47592b05baed7e082325e8bed488 100644
--- a/plot_displacement_error.py
+++ b/plot_displacement_error.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_regional_strains.py b/plot_regional_strains.py
index 4237804bb5dc51375a0de51e06b3db3b86824d54..ab3472bfeb05f4dc27ad81d0e58c0075dffa9368 100644
--- a/plot_regional_strains.py
+++ b/plot_regional_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_strains.py b/plot_strains.py
index 90872397ca9841c7a74ca58e27bd9f777001b90c..ffa967bed80b6d20056eaa742c987f587529a6df 100644
--- a/plot_strains.py
+++ b/plot_strains.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/plot_strains_vs_radius.py b/plot_strains_vs_radius.py
index 73b4e8d8f6d3fc07d5ca45a1ad24d510546b5f6d..8294398de7e47abaf53d647c77817fa6cf015cc6 100644
--- a/plot_strains_vs_radius.py
+++ b/plot_strains_vs_radius.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###
diff --git a/write_VTU_file.py b/write_VTU_file.py
index 9a71f67414a5b5c34e6861c7171a0fc721998709..512d1b46e966665464213e2d098a789eb5e30c1e 100644
--- a/write_VTU_file.py
+++ b/write_VTU_file.py
@@ -2,7 +2,7 @@
################################################################################
### ###
-### Created by Martin Genet, 2016-2018 ###
+### Created by Martin Genet, 2016-2019 ###
### ###
### École Polytechnique, Palaiseau, France ###
### ###