From 1f743f502443a6c6e5e16e1ed84de73e75cd3e45 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Thu, 14 Feb 2019 20:49:00 +0100
Subject: [PATCH] 2018 -> 2019

---
 Energy.py                          |   2 +-
 Energy_Regularization.py           |   2 +-
 Energy_WarpedImage.py              |   2 +-
 ImageIterator.py                   |   2 +-
 ImageSeries.py                     |   2 +-
 NonlinearSolver.py                 |   2 +-
 Problem.py                         |   2 +-
 Problem_ImageRegistration.py       |   2 +-
 compute_displacement_error.py      |   2 +-
 compute_displacements.py           |   2 +-
 compute_quadrature_degree.py       |   2 +-
 compute_strains.py                 |   2 +-
 fedic.py                           |   2 +-
 fedic2.py                          |   2 +-
 generate_images.py                 |   2 +-
 generate_undersampled_images.py    |   2 +-
 generated_grad_expressions_cpp.py  | 176 +++++++++++++++++++++++++++++
 generated_image_expressions_cpp.py |   2 +-
 image_expressions_cpp.py           |   2 +-
 image_expressions_py.py            |   2 +-
 plot_binned_strains_vs_radius.py   |   2 +-
 plot_displacement_error.py         |   2 +-
 plot_regional_strains.py           |   2 +-
 plot_strains.py                    |   2 +-
 plot_strains_vs_radius.py          |   2 +-
 write_VTU_file.py                  |   2 +-
 26 files changed, 201 insertions(+), 25 deletions(-)
 create mode 100644 generated_grad_expressions_cpp.py

diff --git a/Energy.py b/Energy.py
index a645ef4..f0502d8 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 651d89d..1596d14 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 454865a..9411303 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 e650c94..ae45597 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 89eea3b..c3c219f 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 86fc7d6..a071751 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 cdbba9b..a402876 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 36fe471..40e31e5 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 8e7b6ac..74afebb 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 17332be..8b539df 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 78755d5..e507fc5 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 a6c2a27..a1f37c4 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 ee505f8..d7bfac8 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 af560c5..cd2e049 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 163cc8d..ffcf501 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 7deda20..bcd2333 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 0000000..fc64124
--- /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 548e326..7a62201 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 a360e32..541e41d 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 7e7585b..dabb8f1 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 25bfc6f..8d0033e 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 727e80d..57a5756 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 4237804..ab3472b 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 9087239..ffa967b 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 73b4e8d..8294398 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 9a71f67..512d1b4 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                                   ###
 ###                                                                          ###
-- 
GitLab