dolfin_warp/generate_images.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import glob
+import numpy
+import os
+import shutil
+import vtk
+
+import myPythonLibrary    as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_warp as dwarp
+
+from .generate_images_Image   import Image
+from .generate_images_Mapping import Mapping
+
+################################################################################
+
+def set_I_woGrad(
+        image,
+        X,
+        I,
+        scalars,
+        k_point,
+        G=None,
+        Finv=None,
+        vectors=None):
+
+    image.I0(X, I)
+    scalars.SetTuple(k_point, I)
+
+def set_I_wGrad(
+        image,
+        X,
+        I,
+        scalars,
+        k_point,
+        G,
+        Finv,
+        vectors):
+
+    image.I0_wGrad(X, I, G)
+    scalars.SetTuple(k_point, I)
+    G = numpy.dot(G, Finv)
+    vectors.SetTuple(k_point, G)
+
+################################################################################
+
+def generate_images(
+        images,
+        structure,
+        texture,
+        noise,
+        deformation,
+        evolution,
+        generate_gradient=0,
+        keep_temp_images=0,
+        verbose=0):
+
+    mypy.my_print(verbose, "*** generate_images ***")
+
+    assert ("n_integration" not in images),\
+        "\"n_integration\" has been deprecated. Use \"upsampling\" instead. Aborting."
+
+    if ("upsampling_factors" not in images):
+        images["upsampling_factors"] = [1]*images["n_dim"]
+    if ("downsampling_factors" not in images):
+        images["downsampling_factors"] = [1]*images["n_dim"]
+    if ("zfill" not in images):
+        images["zfill"] = len(str(images["n_frames"]))
+    if ("ext" not in images):
+        images["ext"] = "vti"
+
+    if not os.path.exists(images["folder"]):
+        os.mkdir(images["folder"])
+
+    for filename in glob.glob(images["folder"]+"/"+images["basename"]+"_*.*"):
+        os.remove(filename)
+
+    image = Image(
+        images,
+        structure,
+        texture,
+        noise,
+        generate_gradient)
+    mapping = Mapping(
+        images,
+        structure,
+        deformation,
+        evolution,
+        generate_gradient)
+
+    vtk_image = myvtk.createImageFromSizeAndRes(
+        dim  = images["n_dim"],
+        size = images["L"],
+        res  = images["n_voxels"],
+        up   = images["upsampling_factors"])
+    n_points_upsampled = vtk_image.GetNumberOfPoints()
+    vtk_scalars = vtk_image.GetPointData().GetScalars()
+
+    if (generate_gradient):
+        vtk_gradient = vtk.vtkImageData()
+        vtk_gradient.DeepCopy(vtk_image)
+
+        vtk_vectors = myvtk.createFloatArray(
+            name="ImageScalarsGradient",
+            n_components=3,
+            n_tuples=n_points_upsampled,
+            verbose=verbose-1)
+        vtk_gradient.GetPointData().SetScalars(vtk_vectors)
+    else:
+        vtk_gradient = None
+        vtk_vectors  = None
+
+    x = numpy.empty(3)
+    X = numpy.empty(3)
+    I = numpy.empty(1)
+    global_min = float("+Inf")
+    global_max = float("-Inf")
+    if (generate_gradient):
+        G     = numpy.empty(3)
+        Finv  = numpy.empty((3,3))
+        set_I = set_I_wGrad
+    else:
+        G     = None
+        Finv  = None
+        set_I = set_I_woGrad
+
+    for k_frame in range(images["n_frames"]):
+        t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
+        mypy.my_print(verbose, "t = "+str(t))
+        mapping.init_t(t)
+        for k_point in range(n_points_upsampled):
+            vtk_image.GetPoint(k_point, x)
+            #print("x0 = "+str(x))
+            mapping.X(x, X, Finv)
+            #print("X = "+str(X))
+            set_I(image, X, I, vtk_scalars, k_point, G, Finv, vtk_vectors)
+            global_min = min(global_min, I[0])
+            global_max = max(global_max, I[0])
+        #print(vtk_image)
+        myvtk.writeImage(
+            image=vtk_image,
+            filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
+            verbose=verbose-1)
+        if (generate_gradient):
+            #print(vtk_gradient)
+            myvtk.writeImage(
+                image=vtk_gradient,
+                filename=images["folder"]+"/"+images["basename"]+"-grad"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
+                verbose=verbose-1)
+    # mypy.my_print(verbose, "global_min = "+str(global_min))
+    # mypy.my_print(verbose, "global_max = "+str(global_max))
+
+    if (images["upsampling_factors"] == [1]*images["n_dim"]):
+        downsampling = False
+    else:
+        downsampling = True
+
+        if (keep_temp_images):
+            for k_frame in range(images["n_frames"]):
+                shutil.copy(
+                    src=images["folder"]+"/"+images["basename"]             +"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
+                    dst=images["folder"]+"/"+images["basename"]+"_upsampled"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"])
+
+        dwarp.compute_downsampled_images(
+            images_folder=images["folder"],
+            images_basename=images["basename"],
+            downsampling_factors=images["upsampling_factors"],
+            keep_resolution=0,
+            verbose=verbose)
+
+    if (images["downsampling_factors"] == [1]*images["n_dim"]):
+        downsampling = False
+    else:
+        downsampling = True
+
+        if (keep_temp_images):
+            for k_frame in range(images["n_frames"]):
+                shutil.copy(
+                    src=images["folder"]+"/"+images["basename"]                  +"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
+                    dst=images["folder"]+"/"+images["basename"]+"_predownsampled"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"])
+
+        dwarp.compute_downsampled_images(
+            images_folder=images["folder"],
+            images_basename=images["basename"],
+            downsampling_factors=images["downsampling_factors"],
+            keep_resolution=1,
+            verbose=verbose)
+
+    if (images["data_type"] in ("float")):
+        normalizing = False
+    elif (images["data_type"] in ("unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float", "uint8", "uint16", "uint32", "uint64", "ufloat")):
+        normalizing = True
+
+        if (keep_temp_images):
+            for k_frame in range(images["n_frames"]):
+                shutil.copy(
+                    src=images["folder"]+"/"+images["basename"]                 +"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"],
+                    dst=images["folder"]+"/"+images["basename"]+"_prenormalized"+"_"+str(k_frame).zfill(images["zfill"])+"."+images["ext"])
+
+        dwarp.compute_normalized_images(
+            images_folder=images["folder"],
+            images_basename=images["basename"],
+            images_datatype=images["data_type"],
+            verbose=verbose)

dolfin_warp/generate_images_Image.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import math
+import numpy
+import random
+
+################################################################################
+
+class Image():
+
+    def __init__(
+            self,
+            images,
+            structure,
+            texture,
+            noise,
+            generate_image_gradient=0):
+
+        self.L = images["L"]
+
+        # structure
+        if (structure["type"] == "no"):
+            self.I0_structure = self.I0_structure_no_wGrad if (generate_image_gradient) else self.I0_structure_no
+        elif (structure["type"] == "box"):
+            self.I0_structure = self.I0_structure_box_wGrad if (generate_image_gradient) else self.I0_structure_box
+            self.Xmin = structure["Xmin"]+[float("-Inf")]*(3-images["n_dim"])
+            self.Xmax = structure["Xmax"]+[float("+Inf")]*(3-images["n_dim"])
+        elif (structure["type"] in ("ring", "heart")):
+            self.R = float()
+            self.Ri = structure["Ri"]
+            self.Re = structure["Re"]
+            self.X0 = structure["X0"] if ("X0" in structure) else [images["L"][0]/2, images["L"][1]/2]
+            if (images["n_dim"] == 2):
+                self.I0_structure = self.I0_structure_heart_2_wGrad if (generate_image_gradient) else self.I0_structure_heart_2
+            elif (images["n_dim"] == 3):
+                self.I0_structure = self.I0_structure_heart_3_wGrad if (generate_image_gradient) else self.I0_structure_heart_3
+                self.Zmin = structure.Zmin if ("Zmin" in structure) else 0.
+                self.Zmax = structure.Zmax if ("Zmax" in structure) else images["L"][2]
+            else:
+                assert (0), "n_dim must be \"2\" or \"3 for \"ring\"/\"heart\" type structure. Aborting."
+        else:
+            assert (0), "structure type must be \"no\", \"box\", \"ring\" or \"heart\". Aborting."
+
+        # texture
+        if (texture["type"] == "no"):
+            self.I0_texture = self.I0_texture_no_wGrad if (generate_image_gradient) else self.I0_texture_no
+        elif (texture["type"].startswith("tagging")):
+            if   (images["n_dim"] == 1):
+                if ("-signed" in texture["type"]):
+                    self.I0_texture = self.I0_texture_tagging_signed_X_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_X
+                else:
+                    self.I0_texture = self.I0_texture_tagging_X_wGrad if (generate_image_gradient) else self.I0_texture_tagging_X
+            elif (images["n_dim"] == 2):
+                if ("-signed" in texture["type"]):
+                    if   ("-addComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_signed_XY_wAdditiveCombination_wGrad       if (generate_image_gradient) else self.I0_texture_tagging_signed_XY_wAdditiveCombination
+                    elif ("-diffComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_signed_XY_wDifferentiableCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_XY_wDifferentiableCombination
+                    else:
+                        self.I0_texture = self.I0_texture_tagging_signed_XY_wMultiplicativeCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_XY_wMultiplicativeCombination
+                else:
+                    if   ("-addComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_XY_wAdditiveCombination_wGrad       if (generate_image_gradient) else self.I0_texture_tagging_XY_wAdditiveCombination
+                    elif ("-diffComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_XY_wDifferentiableCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_XY_wDifferentiableCombination
+                    else:
+                        self.I0_texture = self.I0_texture_tagging_XY_wMultiplicativeCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_XY_wMultiplicativeCombination
+            elif (images["n_dim"] == 3):
+                if ("-signed" in texture["type"]):
+                    if   ("-addComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_signed_XYZ_wAdditiveCombination_wGrad       if (generate_image_gradient) else self.I0_texture_tagging_signed_XYZ_wAdditiveCombination
+                    elif ("-diffComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_signed_XYZ_wDifferentiableCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_XYZ_wDifferentiableCombination
+                    else:
+                        self.I0_texture = self.I0_texture_tagging_signed_XYZ_wMultiplicativeCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_XYZ_wMultiplicativeCombination
+                else:
+                    if   ("-addComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_XYZ_wAdditiveCombination_wGrad       if (generate_image_gradient) else self.I0_texture_tagging_XYZ_wAdditiveCombination
+                    elif ("-diffComb" in texture["type"]):
+                        self.I0_texture = self.I0_texture_tagging_XYZ_wDifferentiableCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_XYZ_wDifferentiableCombination
+                    else:
+                        self.I0_texture = self.I0_texture_tagging_XYZ_wMultiplicativeCombination_wGrad if (generate_image_gradient) else self.I0_texture_tagging_XYZ_wMultiplicativeCombination
+            else:
+                assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
+            self.s = texture["s"]
+        elif (texture["type"].startswith("taggX")):
+            if ("-signed" in texture["type"]):
+                self.I0_texture = self.I0_texture_tagging_signed_X_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_X
+            else:
+                self.I0_texture = self.I0_texture_tagging_X_wGrad if (generate_image_gradient) else self.I0_texture_tagging_X
+            self.s = texture["s"]
+        elif (texture["type"].startswith("taggY")):
+            if ("-signed" in texture["type"]):
+                self.I0_texture = self.I0_texture_tagging_signed_Y_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_Y
+            else:
+                self.I0_texture = self.I0_texture_tagging_Y_wGrad if (generate_image_gradient) else self.I0_texture_tagging_Y
+            self.s = texture["s"]
+        elif (texture["type"].startswith("taggZ")):
+            if ("-signed" in texture["type"]):
+                self.I0_texture = self.I0_texture_tagging_signed_Z_wGrad if (generate_image_gradient) else self.I0_texture_tagging_signed_Z
+            else:
+                self.I0_texture = self.I0_texture_tagging_Z_wGrad if (generate_image_gradient) else self.I0_texture_tagging_Z
+            self.s = texture["s"]
+        else:
+            assert (0), "texture type must be \"no\", \"tagging\", \"taggX\", \"taggY\" or \"taggZ\". Aborting."
+
+        # noise (MG20220818: Should use dwarp.Noise)
+        if (noise["type"] == "no"):
+            self.I0_noise = self.I0_noise_no_wGrad if (generate_image_gradient) else self.I0_noise_no
+        elif (noise["type"] == "normal"):
+            self.I0_noise = self.I0_noise_normal_wGrad if (generate_image_gradient) else self.I0_noise_normal
+            self.avg = noise["avg"] if ("avg" in noise) else 0.
+            self.std = noise["stdev"]
+        else:
+            assert (0), "noise type must be \"no\" or \"normal\". Aborting."
+
+    def I0(self, X, I):
+        self.I0_structure(X, I)
+        self.I0_texture(X, I)
+        self.I0_noise(I)
+
+    def I0_wGrad(self, X, I, G):
+        self.I0_structure_wGrad(X, I, G)
+        self.I0_texture_wGrad(X, I, G)
+        self.I0_noise_wGrad(I, G)
+
+    def I0_structure_no(self, X, I):
+        I[0] = 1.
+
+    def I0_structure_no_wGrad(self, X, I, G):
+        self.I0_structure_no(X, I)
+        G[:] = 1. # MG 20180806: gradient is given by texture; here it is just indicator function
+
+    def I0_structure_box(self, X, I):
+        if all(numpy.greater_equal(X, self.Xmin)) and all(numpy.less_equal(X, self.Xmax)):
+            I[0] = 1.
+        else:
+            I[0] = 0.
+
+    def I0_structure_box_wGrad(self, X, I, G):
+        if all(numpy.greater_equal(X, self.Xmin)) and all(numpy.less_equal(X, self.Xmax)):
+            I[0] = 1.
+            G[:] = 1. # MG 20180806: gradient is given by texture; here it is just indicator function
+        else:
+            I[0] = 0.
+            G[:] = 0. # MG 20180806: gradient is given by texture; here it is just indicator function
+
+    def I0_structure_heart_2(self, X, I):
+        self.R = ((X[0]-self.X0[0])**2 + (X[1]-self.X0[1])**2)**(1./2)
+        if (self.R >= self.Ri) and (self.R <= self.Re):
+            I[0] = 1.
+        else:
+            I[0] = 0.
+
+    def I0_structure_heart_2_wGrad(self, X, I, G):
+        self.R = ((X[0]-self.X0[0])**2 + (X[1]-self.X0[1])**2)**(1./2)
+        if (self.R >= self.Ri) and (self.R <= self.Re):
+            I[0] = 1.
+            G[:] = 1. # MG 20180806: gradient is given by texture; here it is just indicator function
+        else:
+            I[0] = 0.
+            G[:] = 0. # MG 20180806: gradient is given by texture; here it is just indicator function
+
+    def I0_structure_heart_3(self, X, I):
+        self.R = ((X[0]-self.X0[0])**2 + (X[1]-self.X0[1])**2)**(1./2)
+        if (self.R >= self.Ri) and (self.R <= self.Re) and (X[2] >= self.Zmin) and (X[2] <= self.Zmax):
+            I[0] = 1.
+        else:
+            I[0] = 0.
+
+    def I0_structure_heart_3_wGrad(self, X, I, G):
+        self.R = ((X[0]-self.X0[0])**2 + (X[1]-self.X0[1])**2)**(1./2)
+        if (self.R >= self.Ri) and (self.R <= self.Re) and (X[2] >= self.Zmin) and (X[2] <= self.Zmax):
+            I[0] = 1.
+            G[:] = 1. # MG 20180806: gradient is given by texture; here it is just indicator function
+        else:
+            I[0] = 0.
+            G[:] = 0. # MG 20180806: gradient is given by texture; here it is just indicator function
+
+    def I0_texture_no(self, X, I):
+        I[0] *= 1.
+
+    def I0_texture_no_wGrad(self, X, I, G):
+        self.I0_texture_no(X, I)
+        G[:] *= 0.
+
+    def I0_texture_tagging_X(self, X, I):
+        I[0] *= abs(math.sin(math.pi*X[0]/self.s))
+
+    def I0_texture_tagging_X_wGrad(self, X, I, G):
+        self.I0_texture_tagging_X(X, I)
+        G[0] *= math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s)
+        G[1] *= 0.
+        G[2] *= 0.
+
+    def I0_texture_tagging_Y(self, X, I):
+        I[0] *= abs(math.sin(math.pi*X[1]/self.s))
+
+    def I0_texture_tagging_Y_wGrad(self, X, I, G):
+        self.I0_texture_tagging_Y(X, I)
+        G[0] *= 0.
+        G[1] *= math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s)
+        G[2] *= 0.
+
+    def I0_texture_tagging_Z(self, X, I):
+        I[0] *= abs(math.sin(math.pi*X[2]/self.s))
+
+    def I0_texture_tagging_Z_wGrad(self, X, I, G):
+        self.I0_texture_tagging_Z(X, I)
+        G[0] *= 0.
+        G[1] *= 0.
+        G[2] *= math.copysign(1, math.sin(math.pi*X[2]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[2]/self.s)
+
+    def I0_texture_tagging_XY_wAdditiveCombination(self, X, I):
+        I[0] *= (abs(math.sin(math.pi*X[0]/self.s))
+               + abs(math.sin(math.pi*X[1]/self.s)))/2
+
+    def I0_texture_tagging_XY_wAdditiveCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XY_wAdditiveCombination(X, I)
+        G[0] *= math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) / 2
+        G[1] *= math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) / 2
+        G[2] *= 0.
+
+    def I0_texture_tagging_XY_wMultiplicativeCombination(self, X, I):
+        I[0] *= (abs(math.sin(math.pi*X[0]/self.s))
+             *   abs(math.sin(math.pi*X[1]/self.s)))**(1./2)
+
+    def I0_texture_tagging_XY_wMultiplicativeCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XY_wMultiplicativeCombination(X, I)
+        G[0] *= math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) * abs(math.sin(math.pi*X[1]/self.s)) / 2 / I[0]
+        G[1] *= math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) / 2 / I[0]
+        G[2] *= 0.
+
+    def I0_texture_tagging_XY_wDifferentiableCombination(self, X, I):
+        I[0] *= (1 + 3 * abs(math.sin(math.pi*X[0]/self.s))
+                       * abs(math.sin(math.pi*X[1]/self.s)))**(1./2) - 1
+
+    def I0_texture_tagging_XY_wDifferentiableCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XY_wDifferentiableCombination(X, I)
+        G[0] *= 3 * math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) * abs(math.sin(math.pi*X[1]/self.s)) / 2 / (I[0] + 1)
+        G[1] *= 3 * math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) / 2 / (I[0] + 1)
+        G[2] *= 0.
+
+    def I0_texture_tagging_XYZ_wAdditiveCombination(self, X, I):
+        I[0] *= (abs(math.sin(math.pi*X[0]/self.s))
+               + abs(math.sin(math.pi*X[1]/self.s))
+               + abs(math.sin(math.pi*X[2]/self.s)))/3
+
+    def I0_texture_tagging_XYZ_wAdditiveCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XYZ_wAdditiveCombination(X, I)
+        G[0] *= math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) / 3
+        G[1] *= math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) / 3
+        G[2] *= math.copysign(1, math.sin(math.pi*X[2]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[2]/self.s) / 3
+
+    def I0_texture_tagging_XYZ_wMultiplicativeCombination(self, X, I):
+        I[0] *= (abs(math.sin(math.pi*X[0]/self.s))
+             *   abs(math.sin(math.pi*X[1]/self.s))
+             *   abs(math.sin(math.pi*X[2]/self.s)))**(1./3)
+
+    def I0_texture_tagging_XYZ_wMultiplicativeCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XYZ_wMultiplicativeCombination(X, I)
+        G[0] *= math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) * abs(math.sin(math.pi*X[1]/self.s)) * abs(math.sin(math.pi*X[2]/self.s)) / 3 / I[0]**2
+        G[1] *= math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) * abs(math.sin(math.pi*X[2]/self.s)) / 3 / I[0]**2
+        G[2] *= math.copysign(1, math.sin(math.pi*X[2]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[2]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) * abs(math.sin(math.pi*X[1]/self.s)) / 3 / I[0]**2
+
+    def I0_texture_tagging_XYZ_wDifferentiableCombination(self, X, I):
+        I[0] *= (1 + 7 * abs(math.sin(math.pi*X[0]/self.s))
+                       * abs(math.sin(math.pi*X[1]/self.s))
+                       * abs(math.sin(math.pi*X[2]/self.s)))**(1./3) - 1
+
+    def I0_texture_tagging_XYZ_wDifferentiableCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XYZ_wDifferentiableCombination(X, I)
+        G[0] *= 7 * math.copysign(1, math.sin(math.pi*X[0]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s) * abs(math.sin(math.pi*X[1]/self.s)) * abs(math.sin(math.pi*X[2]/self.s)) / 3 / (I[0] + 1)
+        G[1] *= 7 * math.copysign(1, math.sin(math.pi*X[1]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) * abs(math.sin(math.pi*X[2]/self.s)) / 3 / (I[0] + 1)
+        G[2] *= 7 * math.copysign(1, math.sin(math.pi*X[2]/self.s)) * (math.pi/self.s) * math.cos(math.pi*X[2]/self.s) * abs(math.sin(math.pi*X[0]/self.s)) * abs(math.sin(math.pi*X[1]/self.s)) / 3 / (I[0] + 1)
+
+    def I0_texture_tagging_signed_X(self, X, I):
+        I[0] *= (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+
+    def I0_texture_tagging_signed_X_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_X(X, I)
+        G[0] *= (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2) / 2
+        G[1] *= 0.
+        G[2] *= 0.
+
+    def I0_texture_tagging_signed_Y(self, X, I):
+        I[0] *= (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2
+
+    def I0_texture_tagging_signed_Y_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_Y(X, I)
+        G[0] *= 0.
+        G[1] *= (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2) / 2
+        G[2] *= 0.
+
+    def I0_texture_tagging_signed_Z(self, X, I):
+        I[0] *= (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2
+
+    def I0_texture_tagging_signed_Z_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_Z(X, I)
+        G[0] *= 0.
+        G[1] *= 0.
+        G[2] *= (math.pi/self.s) * math.cos(math.pi*X[2]/self.s-math.pi/2) / 2
+
+    def I0_texture_tagging_signed_XY_wAdditiveCombination(self, X, I):
+        I[0] *= ((1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+              +  (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2) / 2
+
+    def I0_texture_tagging_signed_XY_wAdditiveCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_XY_wAdditiveCombination(X, I)
+        G[0] *= (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 / 2
+        G[1] *= (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 / 2
+        G[2] *= 0.
+
+    def I0_texture_tagging_signed_XY_wMultiplicativeCombination(self, X, I):
+        I[0] *= ((1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+             *   (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2)**(1./2)
+
+    def I0_texture_tagging_signed_XY_wMultiplicativeCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XY_wMultiplicativeCombination(X, I)
+        G[0] *= (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 / 2 / I[0]
+        G[1] *= (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 / 2 / I[0]
+        G[2] *= 0.
+
+    def I0_texture_tagging_signed_XY_wDifferentiableCombination(self, X, I):
+        I[0] *= (1 + 3 * (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+                       * (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2)**(1./2) - 1
+
+    def I0_texture_tagging_signed_XY_wDifferentiableCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_XY_wDifferentiableCombination(X, I)
+        G[0] *= 3 * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 * (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2 / 2 / (I[0] + 1)
+        G[1] *= 3 * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 * (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2 / 2 / (I[0] + 1)
+        G[2] *= 0.
+
+    def I0_texture_tagging_signed_XYZ_wAdditiveCombination(self, X, I):
+        I[0] *= ((1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+              +  (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2
+              +  (1+math.sin(math.pi*X[2]/self.s-math.pi/2))/2) / 3
+
+    def I0_texture_tagging_signed_XYZ_wAdditiveCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_XYZ_wAdditiveCombination(X, I)
+        G[0] *= (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 / 3
+        G[1] *= (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 / 3
+        G[2] *= (math.pi/self.s) * math.cos(math.pi*X[2]/self.s-math.pi/2)/2 / 3
+
+    def I0_texture_tagging_signed_XYZ_wMultiplicativeCombination(self, X, I):
+        I[0] *= ((1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+             *   (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2
+             *   (1+math.sin(math.pi*X[2]/self.s-math.pi/2))/2)**(1./3)
+
+    def I0_texture_tagging_signed_XYZ_wMultiplicativeCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_XYZ_wMultiplicativeCombination(X, I)
+        G[0] *= (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 / 3 / I[0]
+        G[1] *= (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 / 3 / I[0]
+        G[2] *= (math.pi/self.s) * math.cos(math.pi*X[2]/self.s-math.pi/2)/2 / 3 / I[0]
+
+    def I0_texture_tagging_signed_XYZ_wDifferentiableCombination(self, X, I):
+        I[0] *= (1 + 7 * (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2
+                       * (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2
+                       * (1+math.sin(math.pi*X[2]/self.s-math.pi/2))/2)**(1./3) - 1
+
+    def I0_texture_tagging_signed_XYZ_wDifferentiableCombination_wGrad(self, X, I, G):
+        self.I0_texture_tagging_signed_XYZ_wDifferentiableCombination(X, I)
+        G[0] *= 7 * (math.pi/self.s) * math.cos(math.pi*X[0]/self.s-math.pi/2)/2 * (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2 * (1+math.sin(math.pi*X[2]/self.s-math.pi/2))/2 / 3 / (I[0] + 1)
+        G[1] *= 7 * (math.pi/self.s) * math.cos(math.pi*X[1]/self.s-math.pi/2)/2 * (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2 * (1+math.sin(math.pi*X[2]/self.s-math.pi/2))/2 / 3 / (I[0] + 1)
+        G[2] *= 7 * (math.pi/self.s) * math.cos(math.pi*X[2]/self.s-math.pi/2)/2 * (1+math.sin(math.pi*X[0]/self.s-math.pi/2))/2 * (1+math.sin(math.pi*X[1]/self.s-math.pi/2))/2 / 3 / (I[0] + 1)
+
+    def I0_noise_no(self, I):
+        pass
+
+    def I0_noise_no_wGrad(self, I, G):
+        pass
+
+    def I0_noise_normal(self, I):
+        I[0] += random.normalvariate(self.avg, self.std)
+
+    def I0_noise_normal_wGrad(self, I, G):
+        self.I0_noise_normal(I)
+        G[:] += [2*random.normalvariate(self.avg, self.std) for k in range(len(G))]

dolfin_warp/generate_images_Mapping.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import math
+import numpy
+
+################################################################################
+
+class Mapping():
+
+    def __init__(
+            self,
+            images,
+            structure,
+            deformation,
+            evolution,
+            generate_image_gradient=0):
+
+        self.deformation = deformation
+        if (self.deformation["type"] == "no"):
+            self.init_t = self.init_t_no
+            self.X = self.X_no
+            self.x = self.x_no
+        elif (self.deformation["type"] == "translation"):
+            self.init_t = self.init_t_translation
+            self.X = self.X_translation
+            self.x = self.x_translation
+            self.D = numpy.empty(3)
+        elif (self.deformation["type"] == "rotation"):
+            self.init_t = self.init_t_rotation
+            self.X = self.X_rotation
+            self.x = self.x_rotation
+            self.C = numpy.empty(3)
+            self.R = numpy.empty((3,3))
+            self.Rinv = numpy.empty((3,3))
+        elif (self.deformation["type"] == "homogeneous"):
+            self.init_t = self.init_t_homogeneous
+            self.X = self.X_homogeneous
+            self.x = self.x_homogeneous
+            self.X0 = numpy.empty(3)
+        elif (self.deformation["type"] == "heart"):
+            assert (structure["type"] in ("ring", "heart")), "structure type must be \"ring\" or \"heart\" for \"heart\" type deformation, not \""+str(structure["type"])+"\". Aborting."
+            self.init_t = self.init_t_heart
+            self.X = self.X_heart
+            self.x = self.x_heart
+            self.x_inplane = numpy.empty(2)
+            self.X_inplane = numpy.empty(2)
+            self.rt = numpy.empty(2)
+            self.RT = numpy.empty(2)
+            self.L = images["L"]
+            self.Ri = structure["Ri"]
+            self.Re = structure["Re"]
+            self.R = numpy.empty((3,3))
+        else:
+            assert (0), "deformation type must be \"no\", \"translation\", \"rotation\", \"homogeneous\" or \"heart\", not \""+str(self.deformation["type"])+"\". Aborting."
+
+        if (evolution["type"] == "linear"):
+            self.phi = self.phi_linear
+        elif (evolution["type"] == "sine"):
+            self.phi = self.phi_sine
+            self.T = evolution["T"]
+        else:
+            assert (0), "evolution type must be \"linear\" or \"sine\", not \""+str(evolution["type"])+"\". Aborting."
+
+    def phi_linear(self, t):
+        return t
+
+    def phi_sine(self, t):
+        return math.sin(math.pi*t/self.T)**2
+
+    def init_t_no(self, t):
+        pass
+
+    def init_t_translation(self, t):
+        self.D[0] = self.deformation["Dx"] if ("Dx" in self.deformation) else 0.
+        self.D[1] = self.deformation["Dy"] if ("Dy" in self.deformation) else 0.
+        self.D[2] = self.deformation["Dz"] if ("Dz" in self.deformation) else 0.
+        self.D *= self.phi(t)
+
+    def init_t_rotation(self, t):
+        self.C[0] = self.deformation["Cx"] if ("Cx" in self.deformation) else 0.
+        self.C[1] = self.deformation["Cy"] if ("Cy" in self.deformation) else 0.
+        self.C[2] = self.deformation["Cz"] if ("Cz" in self.deformation) else 0.
+        Rx = self.deformation["Rx"]*math.pi/180*self.phi(t) if ("Rx" in self.deformation) else 0.
+        Ry = self.deformation["Ry"]*math.pi/180*self.phi(t) if ("Ry" in self.deformation) else 0.
+        Rz = self.deformation["Rz"]*math.pi/180*self.phi(t) if ("Rz" in self.deformation) else 0.
+        RRx = numpy.array([[          1. ,           0. ,           0. ],
+                           [          0. , +math.cos(Rx), -math.sin(Rx)],
+                           [          0. , +math.sin(Rx), +math.cos(Rx)]])
+        RRy = numpy.array([[+math.cos(Ry),           0. , +math.sin(Ry)],
+                           [          0. ,           1. ,           0. ],
+                           [-math.sin(Ry),           0. , +math.cos(Ry)]])
+        RRz = numpy.array([[+math.cos(Rz), -math.sin(Rz),           0. ],
+                           [+math.sin(Rz), +math.cos(Rz),           0. ],
+                           [          0. ,           0. ,           1. ]])
+        self.R[:,:] = numpy.dot(numpy.dot(RRx, RRy), RRz)
+        self.Rinv[:,:] = numpy.linalg.inv(self.R)
+
+    def init_t_homogeneous(self, t):
+        self.X0[0] = self.deformation["X0"] if ("X0" in self.deformation) else 0.
+        self.X0[1] = self.deformation["Y0"] if ("Y0" in self.deformation) else 0.
+        self.X0[2] = self.deformation["Z0"] if ("Z0" in self.deformation) else 0.
+        if (any(E in self.deformation for E in ("Exx", "Eyy", "Ezz"))): # build F from E
+            Exx = self.deformation["Exx"] if ("Exx" in self.deformation) else 0.
+            Eyy = self.deformation["Eyy"] if ("Eyy" in self.deformation) else 0.
+            Ezz = self.deformation["Ezz"] if ("Ezz" in self.deformation) else 0.
+            self.F = numpy.array([[Exx,  0.,  0.],
+                                  [ 0., Eyy,  0.],
+                                  [ 0.,  0., Ezz]])*self.phi(t)
+            self.F *= 2
+            self.F += numpy.eye(3)
+            w, v = numpy.linalg.eig(self.F)
+            # assert (numpy.diag(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v)) == w).all(), str(numpy.dot(numpy.dot(numpy.transpose(v), self.F), v))+" ≠ "+str(numpy.diag(w))+". Aborting."
+            self.F = numpy.dot(numpy.dot(v, numpy.diag(numpy.sqrt(w))), numpy.transpose(v))
+        else:
+            Fxx = self.deformation["Fxx"] if ("Fxx" in self.deformation) else 1.
+            Fyy = self.deformation["Fyy"] if ("Fyy" in self.deformation) else 1.
+            Fzz = self.deformation["Fzz"] if ("Fzz" in self.deformation) else 1.
+            Fxy = self.deformation["Fxy"] if ("Fxy" in self.deformation) else 0.
+            Fyx = self.deformation["Fyx"] if ("Fyx" in self.deformation) else 0.
+            Fyz = self.deformation["Fyz"] if ("Fyz" in self.deformation) else 0.
+            Fzy = self.deformation["Fzy"] if ("Fzy" in self.deformation) else 0.
+            Fzx = self.deformation["Fzx"] if ("Fzx" in self.deformation) else 0.
+            Fxz = self.deformation["Fxz"] if ("Fxz" in self.deformation) else 0.
+            self.F = numpy.eye(3) + (numpy.array([[Fxx, Fxy, Fxz],
+                                                  [Fyx, Fyy, Fyz],
+                                                  [Fzx, Fzy, Fzz]]) - numpy.eye(3))*self.phi(t)
+        self.Finv = numpy.linalg.inv(self.F)
+
+    def init_t_heart(self, t):
+        self.dRi = self.deformation["dRi"]*self.phi(t) if ("dRi" in self.deformation) else 0.
+        self.dRe = self.deformation["dRi"]*self.phi(t) if ("dRi" in self.deformation) else 0.
+        self.dTi = self.deformation["dTi"]*self.phi(t) if ("dTi" in self.deformation) else 0.
+        self.dTe = self.deformation["dTe"]*self.phi(t) if ("dTe" in self.deformation) else 0.
+        self.A = numpy.array([[1.-(self.dRi-self.dRe)/(self.Re-self.Ri), 0.],
+                              [  -(self.dTi-self.dTe)/(self.Re-self.Ri), 1.]])
+        self.Ainv = numpy.linalg.inv(self.A)
+        self.B = numpy.array([(1.+self.Ri/(self.Re-self.Ri))*self.dRi-self.Ri/(self.Re-self.Ri)*self.dRe,
+                              (1.+self.Ri/(self.Re-self.Ri))*self.dTi-self.Ri/(self.Re-self.Ri)*self.dTe])
+
+    def X_no(self, x, X, Finv=None):
+        X[:] = x
+        if (Finv is not None): Finv[:,:] = numpy.identity(numpy.sqrt(numpy.size(Finv)))
+
+    def X_translation(self, x, X, Finv=None):
+        X[:] = x - self.D
+        if (Finv is not None): Finv[:,:] = numpy.identity(numpy.sqrt(numpy.size(Finv)))
+
+    def X_rotation(self, x, X, Finv=None):
+        X[:] = numpy.dot(self.Rinv, x - self.C) + self.C
+        if (Finv is not None): Finv[:,:] = self.Rinv
+
+    def X_homogeneous(self, x, X, Finv=None):
+        X[:] = numpy.dot(self.Finv, x - self.X0) + self.X0
+        if (Finv is not None): Finv[:,:] = self.Finv
+
+    def X_heart(self, x, X, Finv=None):
+        #print("x = "+str(x))
+        self.x_inplane[0] = x[0] - self.L[0]/2
+        self.x_inplane[1] = x[1] - self.L[1]/2
+        #print("x_inplane = "+str(self.x_inplane))
+        self.rt[0] = numpy.linalg.norm(self.x_inplane)
+        self.rt[1] = math.atan2(self.x_inplane[1], self.x_inplane[0])
+        #print("rt = "+str(self.rt))
+        self.RT[:] = numpy.dot(self.Ainv, self.rt-self.B)
+        #print("RT = "+str(self.RT))
+        X[0] = self.RT[0] * math.cos(self.RT[1]) + self.L[0]/2
+        X[1] = self.RT[0] * math.sin(self.RT[1]) + self.L[1]/2
+        X[2] = x[2]
+        #print("X = "+str(X))
+        if (Finv is not None):
+            Finv[0,0] = 1.+(self.dRe-self.dRi)/(self.Re-self.Ri)
+            Finv[0,1] = 0.
+            Finv[0,2] = 0.
+            Finv[1,0] = (self.dTe-self.dTi)/(self.Re-self.Ri)*self.rt[0]
+            Finv[1,1] = self.rt[0]/self.RT[0]
+            Finv[1,2] = 0.
+            Finv[2,0] = 0.
+            Finv[2,1] = 0.
+            Finv[2,2] = 1.
+            #print("F = "+str(Finv))
+            Finv[:,:] = numpy.linalg.inv(Finv)
+            #print("Finv = "+str(Finv))
+            self.R[0,0] = +math.cos(self.RT[1])
+            self.R[0,1] = +math.sin(self.RT[1])
+            self.R[0,2] = 0.
+            self.R[1,0] = -math.sin(self.RT[1])
+            self.R[1,1] = +math.cos(self.RT[1])
+            self.R[1,2] = 0.
+            self.R[2,0] = 0.
+            self.R[2,1] = 0.
+            self.R[2,2] = 1.
+            #print("R = "+str(self.R))
+            Finv[:] = numpy.dot(numpy.transpose(self.R), numpy.dot(Finv, self.R))
+            #print("Finv = "+str(Finv))
+
+    def x_no(self, X, x, F=None):
+        x[:] = X
+        if (F is not None): F[:,:] = numpy.identity(numpy.sqrt(numpy.size(F)))
+
+    def x_translation(self, X, x, F=None):
+        x[:] = X + self.D
+        if (F is not None): F[:,:] = numpy.identity(numpy.sqrt(numpy.size(F)))
+
+    def x_rotation(self, X, x, F=None):
+        x[:] = numpy.dot(self.R, X - self.C) + self.C
+        if (F is not None): F[:,:] = self.R
+
+    def x_homogeneous(self, X, x, F=None):
+        x[:] = numpy.dot(self.F, X - self.X0) + self.X0
+        if (F is not None): F[:,:] = self.F
+
+    def x_heart(self, X, x, F=None):
+        #print("X = "+str(X))
+        self.X_inplane[0] = X[0] - self.L[0]/2
+        self.X_inplane[1] = X[1] - self.L[1]/2
+        #print("X_inplane = "+str(self.X_inplane))
+        self.RT[0] = numpy.linalg.norm(self.X_inplane)
+        self.RT[1] = math.atan2(self.X_inplane[1], self.X_inplane[0])
+        #print("RT = "+str(self.RT))
+        self.rt[:] = numpy.dot(self.A, self.RT) + self.B
+        #print("rt = "+str(self.rt))
+        x[0] = self.rt[0] * math.cos(self.rt[1]) + self.L[0]/2
+        x[1] = self.rt[0] * math.sin(self.rt[1]) + self.L[1]/2
+        x[2] = X[2]
+        #print("x = "+str(x))
+        if (F is not None):
+            F[0,0] = 1.+(self.dRe-self.dRi)/(self.Re-self.Ri)
+            F[0,1] = 0.
+            F[0,2] = 0.
+            F[1,0] = (self.dTe-self.dTi)/(self.Re-self.Ri)*self.rt[0]
+            F[1,1] = self.rt[0]/self.RT[0]
+            F[1,2] = 0.
+            F[2,0] = 0.
+            F[2,1] = 0.
+            F[2,2] = 1.
+            #print("F = "+str(F))
+            self.R[0,0] = +math.cos(self.RT[1])
+            self.R[0,1] = +math.sin(self.RT[1])
+            self.R[0,2] = 0.
+            self.R[1,0] = -math.sin(self.RT[1])
+            self.R[1,1] = +math.cos(self.RT[1])
+            self.R[1,2] = 0.
+            self.R[2,0] = 0.
+            self.R[2,1] = 0.
+            self.R[2,2] = 1.
+            F[:] = numpy.dot(numpy.transpose(self.R), numpy.dot(F, self.R))
+            #print("F = "+str(F))

dolfin_warp/generate_images_Noise.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+from __future__ import annotations # MG20220819: Necessary list[float] type hints in python < 3.10
+
+import random
+
+################################################################################
+
+class Noise():
+    def __init__(self,
+            params: dict = {}) -> None:
+
+        noise_type = params.get("type", "no")
+        if (noise_type == "no"):
+            self.add_noise = self.add_noise_no
+        elif (noise_type == "normal"):
+            self.add_noise = self.add_noise_normal
+            self.avg = params.get("avg", 0.)
+            self.std = params.get("stdev", 0.)
+        else:
+            assert (0), "noise type must be \"no\" or \"normal\". Aborting."
+
+    def add_noise_no(self,
+            I: list[float]) -> None:
+        pass
+
+    def add_noise_normal(self,
+            I: list[float]) -> None:
+        I[0] += random.normalvariate(self.avg, self.std)

dolfin_warp/generate_undersampled_images.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import numpy
+import os
+import vtk
+
+import myPythonLibrary    as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_warp as dwarp
+
+################################################################################
+
+def generateUndersampledImages(
+        images,
+        structure,
+        texture,
+        noise,
+        deformation,
+        evolution,
+        undersampling_level,
+        keep_temporary_images=0,
+        verbose=0):
+
+    mypy.my_print(verbose, "*** generateUndersampledImages ***")
+
+    images_basename = images["basename"]
+    images_n_voxels = images["n_voxels"][:]
+
+    images["n_voxels"][1] /= undersampling_level
+    if (images["n_dim"] == 3):
+        images["n_voxels"][2] /= undersampling_level
+    images["basename"] = images_basename+"-X"
+    texture["type"] = "taggX"
+    dwarp.generate_images(
+        images=images,
+        structure=structure,
+        texture=texture,
+        noise=noise,
+        deformation=deformation,
+        evolution=evolution,
+        verbose=verbose-1)
+    images["n_voxels"][:] = images_n_voxels[:]
+    images["basename"] = images_basename
+
+    images["n_voxels"][0] /= undersampling_level
+    if (images["n_dim"] == 3):
+        images["n_voxels"][2] /= undersampling_level
+    images["basename"] = images_basename+"-Y"
+    texture["type"] = "taggY"
+    dwarp.generate_images(
+        images=images,
+        structure=structure,
+        texture=texture,
+        noise=noise,
+        deformation=deformation,
+        evolution=evolution,
+        verbose=verbose-1)
+    images["n_voxels"][:] = images_n_voxels[:]
+    images["basename"] = images_basename
+
+    if (images["n_dim"] == 3):
+        images["n_voxels"][0] /= undersampling_level
+        images["n_voxels"][1] /= undersampling_level
+        images["basename"] = images_basename+"-Z"
+        texture["type"] = "taggZ"
+        dwarp.generate_images(
+            images=images,
+            structure=structure,
+            texture=texture,
+            noise=noise,
+            deformation=deformation,
+            evolution=evolution,
+            verbose=verbose-1)
+        images["n_voxels"][:] = images_n_voxels[:]
+        images["basename"] = images_basename
+
+    if ("zfill" not in images.keys()):
+        images["zfill"] = len(str(images["n_frames"]))
+    for k_frame in range(images["n_frames"]):
+        imageX = myvtk.readImage(
+            filename=images["folder"]+"/"+images["basename"]+"-X_"+str(k_frame).zfill(images["zfill"])+".vti",
+            verbose=verbose-1)
+        interpolatorX = myvtk.getImageInterpolator(
+            image=imageX,
+            verbose=verbose-1)
+        iX = numpy.empty(1)
+
+        imageY = myvtk.readImage(
+            filename=images["folder"]+"/"+images["basename"]+"-Y_"+str(k_frame).zfill(images["zfill"])+".vti",
+            verbose=verbose-1)
+        interpolatorY = myvtk.getImageInterpolator(
+            image=imageY,
+            verbose=verbose-1)
+        iY = numpy.empty(1)
+
+        if (images["n_dim"] == 2):
+            imageXY = myvtk.createImage(
+                origin=[images["L"][0]/images["n_voxels"][0]/2, images["L"][1]/images["n_voxels"][1]/2, 0.],
+                spacing=[images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], 1.],
+                extent=[0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, 0])
+            scalars = imageXY.GetPointData().GetScalars()
+            x = numpy.empty(3)
+            for k_point in range(imageXY.GetNumberOfPoints()):
+                imageXY.GetPoint(k_point, x)
+                interpolatorX.Interpolate(x, iX)
+                interpolatorY.Interpolate(x, iY)
+                scalars.SetTuple(k_point, (iX*iY)**(1./2))
+            myvtk.writeImage(
+                image=imageXY,
+                filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+                verbose=verbose-1)
+            if not (keep_temporary_images):
+                os.system("rm "+images["folder"]+"/"+images["basename"]+"-X_"+str(k_frame).zfill(images["zfill"])+".vti")
+                os.system("rm "+images["folder"]+"/"+images["basename"]+"-Y_"+str(k_frame).zfill(images["zfill"])+".vti")
+
+        elif (images["n_dim"] == 3):
+            imageZ = myvtk.readImage(
+                filename=images["folder"]+"/"+images["basename"]+"-Z_"+str(k_frame).zfill(images["zfill"])+".vti",
+                verbose=verbose-1)
+            interpolatorZ = myvtk.getImageInterpolator(
+                image=imageZ,
+                verbose=verbose-1)
+            iZ = numpy.empty(1)
+
+            imageXYZ = myvtk.createImage(
+                origin=[images["L"][0]/images["n_voxels"][0]/2, images["L"][1]/images["n_voxels"][1]/2, images["L"][2]/images["n_voxels"][2]/2],
+                spacing=[images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], images["L"][2]/images["n_voxels"][2]],
+                extent=[0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, images["n_voxels"][2]-1])
+            scalars = imageXYZ.GetPointData().GetScalars()
+            x = numpy.empty(3)
+            for k_point in range(imageXYZ.GetNumberOfPoints()):
+                imageXYZ.GetPoint(k_point, x)
+                interpolatorX.Interpolate(x, iX)
+                interpolatorY.Interpolate(x, iY)
+                interpolatorZ.Interpolate(x, iZ)
+                scalars.SetTuple(k_point, (iX*iY*iZ)**(1./3))
+            myvtk.writeImage(
+                image=imageXYZ,
+                filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+                verbose=verbose-1)
+            if not (keep_temporary_images):
+                os.system("rm "+images["folder"]+"/"+images["basename"]+"-X_"+str(k_frame).zfill(images["zfill"])+".vti")
+                os.system("rm "+images["folder"]+"/"+images["basename"]+"-Y_"+str(k_frame).zfill(images["zfill"])+".vti")
+                os.system("rm "+images["folder"]+"/"+images["basename"]+"-Z_"+str(k_frame).zfill(images["zfill"])+".vti")

dolfin_warp/get_centroids.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Ezgi Berberoğlu, 2017-2021                                    ###
+###                                                                          ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland         ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+### And Martin Genet, 2016-2024                                              ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import numpy
+
+import myVTKPythonLibrary as myvtk
+
+################################################################################
+
+def get_centroids(mesh):
+
+    assert mesh.GetCellData().HasArray("sector_id"), \
+        "There is no field named sector_id. Aborting."
+    sector_id = mesh.GetCellData().GetArray("sector_id")
+    n_sector_ids = int(sector_id.GetRange()[1]+1)
+    n_cells = mesh.GetNumberOfCells()
+
+    sector_centroids = numpy.zeros([n_sector_ids,3])
+    n_cells_per_sector = numpy.zeros(n_sector_ids)
+    cell_centers = myvtk.getCellCenters(
+        mesh=mesh)
+
+    for k_cell in range(n_cells):
+        k_cell_sector_id = int(sector_id.GetTuple(k_cell)[0])
+        sector_centroids[k_cell_sector_id,:] = numpy.add(
+            sector_centroids[k_cell_sector_id,:],
+            cell_centers.GetPoint(k_cell))
+        n_cells_per_sector[k_cell_sector_id] += 1
+
+    for k_sector in range(n_sector_ids):
+        sector_centroids[k_sector,:] = sector_centroids[k_sector,:]/n_cells_per_sector[k_sector]
+
+    return sector_centroids

plot_strains.py → dolfin_warp/plot_binned_strains_vs_radius.py

 #coding=utf8
 
-########################################################################
-###                                                                  ###
-### Created by Martin Genet, 2016                                    ###
-###                                                                  ###
-### École Polytechnique, Palaiseau, France                           ###
-###                                                                  ###
-########################################################################
-
-import math
-import numpy
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
 import os
 import sys
 
-########################################################################
+################################################################################
 
-def plot_strains(
+def plot_binned_strains_vs_radius(
         working_folder,
         working_basenames,
         ref_folder=None,
@@ -24,22 +22,30 @@ def plot_strains(
         suffix="",
         verbose=1):
 
-    assert (components in ("all", "circ-long", "rad-circ"))
-
     if (ref_folder is not None) and (ref_basename is not None):
-        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()[1:]
     else:
-        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()
-    n_frames = len(lines)-1
-    n_sectors = (len(lines[-1].split(" "))-1)/12
-    #print "n_frames = " + str(n_frames)
-    #print "n_sectors = " + str(n_sectors)
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
 
-    plotfile = open("plot_strains"+("-"+suffix)*(suffix!="")+".plt", "w")
+    assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basenames[0]+"-binned_strains_vs_radius"
+    else:
+        plotfile_basename = "plot_binned_strains_vs_radius"+("-"+suffix)*(suffix!="")
+
+    plotfile = open(plotfile_basename+".plt", "w")
     plotfile.write('''\
 set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
 
-set output "plot_strains'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+set output "'''+plotfile_basename+'''.pdf"
 
 load "Set1.plt"
 set linestyle 1 pointtype 1
@@ -57,20 +63,14 @@ set key box textcolor variable width +0
 set grid
 
 ''')
-
-    for k_sector in range(n_sectors):
+    for k_frame in range(n_frames):
         plotfile.write('''\
 set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
 
+set xrange [0.:1.]
+set xtics("endocardium" 0, "epicardium" 1)
+
 ''')
-        if (components == "all"):
-            comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
-        elif (components == "circ-long"):
-            comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
-        elif (components == "rad-circ"):
-            comp_names = ["radial", "circumferential", "radial-circumferential"]
-        else:
-            assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
         for comp_name in comp_names:
             if   (comp_name == "radial"                      ): k_comp = 0
             elif (comp_name == "circumferential"             ): k_comp = 1
@@ -79,29 +79,27 @@ set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in
             elif (comp_name == "radial-longitudinal"         ): k_comp = 4
             elif (comp_name == "circumferential-longitudinal"): k_comp = 5
             plotfile.write('''\
-set xrange [0:'''+str(n_frames)+''']
-
 set ylabel "'''+comp_name+''' strain (%)"
-set yrange [-30:30]
+set yrange [-50:50]
 
 plot 0 linecolor rgb "black" notitle,\\
 ''')
             if (ref_folder is not None) and (ref_basename is not None):
                 plotfile.write('''\
-     "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linecolor "black" linewidth 5 notitle,\
-     "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linecolor "black" pointtype 1 linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
+     "'''+ref_folder+'''/'''+ref_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+''') index '''+str(k_frame)+''' with lines linecolor "black" linewidth 3 notitle,\\
+     "'''+ref_folder+'''/'''+ref_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+'''):(100*$'''+str(3+2*k_comp+1)+''') index '''+str(k_frame)+''' with errorbars linecolor "black" linewidth 1 notitle'''+(len(working_basenames)>0)*(''',\\
 ''')+(len(working_basenames)==0)*('''
 
 '''))
             for k_basename in range(len(working_basenames)):
                 working_basename = working_basenames[k_basename]
                 plotfile.write('''\
-     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ($1+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\
-     "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ($1+'''+str(k_basename)+'''./100):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+     "'''+working_folder+'''/'''+working_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+''') index '''+str(k_frame)+''' with lines linestyle '''+str(k_basename+1)+''' linewidth 3 title "'''+working_basename+'''",\\
+     "'''+working_folder+'''/'''+working_basename+'''-binned_strains_vs_radius.dat" using ($2):(100*$'''+str(3+2*k_comp)+'''):(100*$'''+str(3+2*k_comp+1)+''') index '''+str(k_frame)+''' with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
 ''')+(k_basename==len(working_basenames)-1)*('''
 
 '''))
 
     plotfile.close()
 
-    os.system("gnuplot plot_strains"+("-"+suffix)*(suffix!="")+".plt")
+    os.system("gnuplot "+plotfile_basename+".plt")

plot_displacement_error.py → dolfin_warp/plot_displacement_error.py

 #coding=utf8
 
-########################################################################
-###                                                                  ###
-### Created by Martin Genet, 2016                                    ###
-###                                                                  ###
-### École Polytechnique, Palaiseau, France                           ###
-###                                                                  ###
-########################################################################
-
-import math
-import numpy
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
 import os
-import sys
 
-########################################################################
+################################################################################
 
 def plot_displacement_error(
         working_folder,
@@ -22,7 +19,7 @@ def plot_displacement_error(
         verbose=1):
 
     n_frames = len(open(working_folder+"/"+working_basenames[0]+"-error.dat").readlines())-1
-    #print "n_frames = " + str(n_frames)
+    #print("n_frames = " + str(n_frames))
 
     plotfile = open("plot_displacement_error"+("-"+suffix)*(suffix!="")+".plt", "w")
     plotfile.write('''\

dolfin_warp/plot_regional_strains.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import math
+import matplotlib.pyplot
+import numpy
+
+################################################################################
+
+def plot_regional_strains(
+        working_folder,
+        working_basename,
+        k_frame=None,
+        components="all", # all, circ-long, or rad-circ
+        yranges=[0]*6,
+        n_sectors_c=6,
+        n_sectors_l=3,
+        suffix="",
+        verbose=1):
+
+    assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+        n_cols = 3
+        n_rows = 2
+    elif (components == "circ-long"):
+        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
+        n_cols = 3
+        n_rows = 1
+    elif (components == "rad-circ"):
+        comp_names = ["radial", "circumferential", "radial-circumferential"]
+        n_cols = 3
+        n_rows = 1
+    n_comp = len(comp_names)
+
+    strains_all = 100*numpy.loadtxt(working_folder+"/"+working_basename+"-strains.dat")[:,1:]
+
+    n_frames = len(strains_all)
+    if (k_frame is None):
+        #k_frame = n_frames//2
+        k_frame = numpy.argmin(strains_all[:,2])
+
+    strains_es = strains_all[k_frame]
+    #print("len(strains_es) = "+str(len(strains_es)))
+
+    n_sectors = n_sectors_l * n_sectors_c
+    assert (len(strains_es)//2 == (1+n_sectors)*n_comp), "Number of strain components ("+str(len(strains_es)//2)+") inconsistent with number of sectors (n_sectors_c="+str(n_sectors_c)+", n_sectors_l="+str(n_sectors_l)+"). Aborting."
+
+    strains_es_avg = [[strains_es[(k_sector+1)*2*n_comp+2*k_comp] for k_sector in range(n_sectors)] for k_comp in range(n_comp)]
+
+    size = 4
+    fig = matplotlib.pyplot.figure(figsize=(n_cols*size,n_rows*size))
+
+    for k_comp in range(n_comp):
+
+        subplot = fig.add_subplot(n_rows, n_cols, k_comp+1, projection='polar')
+        subplot.set_title(comp_names[k_comp]+" strain (%)")
+
+        subplot.set_xticks([])
+        subplot.set_yticks([])
+
+        strains_comp = strains_es_avg[k_comp]
+
+        yrange = yranges[k_comp]
+        if (yrange == 0):
+            strains_comp_min = min(strains_comp)
+            strains_comp_max = max(strains_comp)
+            strains_comp_min = min(strains_comp_min, -strains_comp_max)
+            strains_comp_max = -strains_comp_min
+        else:
+            strains_comp_min = -yrange
+            strains_comp_max = +yrange
+        assert (strains_comp_max>strains_comp_min), "strains_comp_max ("+str(strains_comp_max)+") <= strains_comp_min (strains_comp_min). Aborting."
+
+        cmap = matplotlib.pyplot.cm.get_cmap('coolwarm')
+        smap = matplotlib.pyplot.cm.ScalarMappable(
+            cmap=cmap,
+            norm=matplotlib.pyplot.Normalize(vmin=strains_comp_min, vmax=strains_comp_max))
+        smap._A = []
+
+        cbar = matplotlib.pyplot.colorbar(
+            mappable=smap,
+            #orientation="horizontal",
+            format="%+g",
+            shrink=2./3)
+        #cbar.set_label(comp_names[k_comp]+" (%)")
+        cbar.solids.set_edgecolor("face")
+
+        k_sector = 0
+        for k_l in range(n_sectors_l):
+            for k_c in range(n_sectors_c):
+                subplot.bar(
+                    left = k_c * 2*math.pi/n_sectors_c,
+                    height = 1./n_sectors_l,
+                    width = 2*math.pi/n_sectors_c,
+                    bottom = 1.-float(k_l+1)/n_sectors_l,
+                    color = cmap(float(strains_comp[k_sector]-strains_comp_min)/(strains_comp_max - strains_comp_min)),
+                    linewidth=0)
+                #subplot.annotate(
+                    #"{:+2.1f}".format(strains_comp[k_sector]),
+                    #xy=    [(k_c+0.5) * 2*math.pi/n_sectors_c, 1.-float(k_l+0.5)/n_sectors_l],
+                    #xytext=[(k_c+0.5) * 2*math.pi/n_sectors_c, 1.-float(k_l+0.5)/n_sectors_l],
+                    #xycoords='polar',
+                    #textcoords='data',
+                    #horizontalalignment='center',
+                    #verticalalignment='center')
+                subplot.annotate(
+                    str(k_sector+1),
+                    xy=    [(k_c+0.5) * 2*math.pi/n_sectors_c, 1.-float(k_l+0.5)/n_sectors_l],
+                    xytext=[(k_c+0.5) * 2*math.pi/n_sectors_c, 1.-float(k_l+0.5)/n_sectors_l],
+                    xycoords='polar',
+                    textcoords='data',
+                    horizontalalignment='center',
+                    verticalalignment='center')
+                k_sector += 1
+
+        subplot.annotate(
+            "anterior",
+            textcoords='data',
+            xycoords='polar',
+            xy=    [1*math.pi/4,1.1],
+            xytext=[1*math.pi/4,1.1],
+            horizontalalignment='center',
+            verticalalignment='center',
+            rotation=-45)
+        subplot.annotate(
+            "septal",
+            textcoords='data',
+            xycoords='polar',
+            xy=    [3*math.pi/4,1.1],
+            xytext=[3*math.pi/4,1.1],
+            horizontalalignment='center',
+            verticalalignment='center',
+            rotation=45)
+        subplot.annotate(
+            "inferior",
+            textcoords='data',
+            xycoords='polar',
+            xy=    [5*math.pi/4,1.1],
+            xytext=[5*math.pi/4,1.1],
+            horizontalalignment='center',
+            verticalalignment='center',
+            rotation=-45)
+        subplot.annotate(
+            "lateral",
+            textcoords='data',
+            xycoords='polar',
+            xy=    [7*math.pi/4,1.1],
+            xytext=[7*math.pi/4,1.1],
+            horizontalalignment='center',
+            verticalalignment='center',
+            rotation=45)
+
+    matplotlib.pyplot.tight_layout()
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basename+"-regional_strains"
+    else:
+        plotfile_basename = "plot_regional_strains"+("-"+suffix)*(suffix!="")
+    #matplotlib.pyplot.savefig(plotfile_basename+".pdf", format='pdf')
+    matplotlib.pyplot.savefig(plotfile_basename+".pdf", format='pdf', bbox_inches="tight")

dolfin_warp/plot_strains.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import os
+
+################################################################################
+
+def plot_strains(
+        working_folder,
+        working_basenames,
+        working_scalings=None,
+        ref_folder=None,
+        ref_basename=None,
+        components="all", # all, circ-long, or rad-circ
+        yranges=[0]*6,
+        suffix="",
+        verbose=1):
+
+    if (working_scalings is not None):
+        assert (len(working_scalings) == len(working_basenames))
+    else:
+        working_scalings = [1.] * len(working_basenames)
+
+    if (ref_folder is not None) and (ref_basename is not None):
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()[1:]
+    else:
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
+    n_sectors = (len(lines[0].split(" "))-1)//12
+    #print("n_frames = " + str(n_frames))
+    #print("n_sectors = " + str(n_sectors))
+
+    comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
+    assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_indx = [0,1,2,3,4,5]
+        n_cols = 3
+        n_rows = 2
+    elif (components == "circ-long"):
+        comp_indx = [1,2,5]
+        n_cols = 3
+        n_rows = 1
+    elif (components == "rad-circ"):
+        comp_indx = [0,1,3]
+        n_cols = 3
+        n_rows = 1
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basenames[0]+"-strains"
+    else:
+        plotfile_basename = "plot_strains"+("-"+suffix)*(suffix!="")
+    plotfile = open(plotfile_basename+".plt", "w")
+
+    size_x = 5
+    size_y = 3
+    plotfile.write('''\
+set terminal pdf enhanced size '''+str(n_cols*size_x)+''','''+str(n_rows*size_y)+'''
+
+set output "'''+plotfile_basename+'''.pdf"
+
+load "Set1.plt"
+
+set key box textcolor variable width +0
+
+set grid
+
+''')
+
+    for k_sector in range(n_sectors):
+        plotfile.write('''\
+set multiplot layout '''+str(n_rows)+''','''+str(n_cols)+(''' title "sector '''+str(k_sector)+'''"''')*(k_sector > 0)+'''
+
+set xlabel "frame ()"
+set xrange [0:'''+str(n_frames)+''']
+
+''')
+        for k_comp in comp_indx:
+            plotfile.write('''\
+set ylabel "'''+comp_names[k_comp]+''' strain (%)"
+''')
+            yrange = yranges[k_comp]
+            if (yrange > 0):
+                plotfile.write('''\
+set yrange [-'''+str(yrange)+''':'''+str(yrange)+''']
+''')
+            plotfile.write('''\
+
+plot 0 linecolor rgb "black" notitle,\\
+''')
+            if (ref_folder is not None) and (ref_basename is not None):
+                plotfile.write('''\
+    "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linecolor "black" linewidth 5 notitle,\\
+    "'''+ref_folder+'''/'''+ref_basename+'''-strains.dat" using ($1):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linecolor "black" linewidth 1 pointtype 1 notitle'''+(len(working_basenames)>0)*(''',\\
+''')+(len(working_basenames)==0)*('''
+
+'''))
+            for k_basename in range(len(working_basenames)):
+                working_basename = working_basenames[k_basename]
+                working_scaling = working_scalings[k_basename]
+                plotfile.write('''\
+    "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with lines linestyle '''+str(k_basename+1)+''' linewidth 5 title "'''+working_basename+'''",\\
+    "'''+working_folder+'''/'''+working_basename+'''-strains.dat" using ('''+str(working_scaling)+'''*($1)+'''+str(k_basename)+'''./10):(100*$'''+str(2+12*k_sector+2*k_comp)+'''):(100*$'''+str(2+12*k_sector+2*k_comp+1)+''') with errorbars linestyle '''+str(k_basename+1)+''' linewidth 1 pointtype 1 notitle'''+(k_basename<len(working_basenames)-1)*(''',\\
+''')+(k_basename==len(working_basenames)-1)*('''
+
+'''))
+
+    plotfile.close()
+
+    os.system("gnuplot \""+plotfile_basename+".plt\"")

plot_strains_vs_radius.py → dolfin_warp/plot_strains_vs_radius.py

 #coding=utf8
 
-########################################################################
-###                                                                  ###
-### Created by Martin Genet, 2016                                    ###
-###                                                                  ###
-### École Polytechnique, Palaiseau, France                           ###
-###                                                                  ###
-########################################################################
-
-import math
-import numpy
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
 import os
-import sys
 
-########################################################################
+################################################################################
 
 def plot_strains_vs_radius(
         working_folder,
         working_basenames,
-        index=10,
         ref_folder=None,
         ref_basename=None,
-        components="all",
+        components="all", # all, circ-long, or rad-circ
+        yranges=[0]*6,
         suffix="",
         verbose=1):
 
+    if (ref_folder is not None) and (ref_basename is not None):
+        lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()[1:]
+    else:
+        lines = open(working_folder+"/"+working_basenames[0]+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
+
+    comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
     assert (components in ("all", "circ-long", "rad-circ"))
+    if (components == "all"):
+        comp_indx = [0,1,2,3,4,5]
+        n_cols = 3
+        n_rows = 2
+    elif (components == "circ-long"):
+        comp_indx = [1,2,5]
+        n_cols = 3
+        n_rows = 1
+    elif (components == "rad-circ"):
+        comp_indx = [0,1,3]
+        n_cols = 3
+        n_rows = 1
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basenames[0]+"-strains_vs_radius"
+    else:
+        plotfile_basename = "plot_strains_vs_radius"+("-"+suffix)*(suffix!="")
+    plotfile = open(plotfile_basename+".plt", "w")
 
-    plotfile = open("plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt", "w")
+    size_x = 5
+    size_y = 3
     plotfile.write('''\
-set terminal pdf enhanced size 15,'''+('''6''')*(components=="all")+('''3''')*(components in ("circ-long", "rad-circ"))+'''
+set terminal pdf enhanced size '''+str(n_cols*size_x)+''','''+str(n_rows*size_y)+'''
 
-set output "plot_strains_vs_radius'''+('''-'''+suffix)*(suffix!="")+'''.pdf"
+set output "'''+plotfile_basename+'''.pdf"
 
 load "Set1.plt"
-set linestyle 1 pointtype 1
-set linestyle 2 pointtype 1
-set linestyle 3 pointtype 1
-set linestyle 4 pointtype 1
-set linestyle 5 pointtype 1
-set linestyle 6 pointtype 1
-set linestyle 7 pointtype 1
-set linestyle 8 pointtype 1
-set linestyle 9 pointtype 1
 
 set key box textcolor variable width +0
 
 set grid
 
-set multiplot layout '''+('''2''')*(components=="all")+('''1''')*(components in ("circ-long", "rad-circ"))+''',3
-
 ''')
-    if (components == "all"):
-        comp_names = ["radial", "circumferential", "longitudinal", "radial-circumferential", "radial-longitudinal", "circumferential-longitudinal"]
-    elif (components == "circ-long"):
-        comp_names = ["circumferential","longitudinal","circumferential-longitudinal"]
-    elif (components == "rad-circ"):
-        comp_names = ["radial", "circumferential", "radial-circumferential"]
-    else:
-        assert (0), "components must be \"all\", \"circ-long\" or \"rad-circ\". Aborting."
-    for comp_name in comp_names:
-        if   (comp_name == "radial"                      ): k_comp = 0
-        elif (comp_name == "circumferential"             ): k_comp = 1
-        elif (comp_name == "longitudinal"                ): k_comp = 2
-        elif (comp_name == "radial-circumferential"      ): k_comp = 3
-        elif (comp_name == "radial-longitudinal"         ): k_comp = 4
-        elif (comp_name == "circumferential-longitudinal"): k_comp = 5
+    for k_frame in range(n_frames):
         plotfile.write('''\
+set multiplot layout '''+str(n_rows)+''','''+str(n_cols)+'''
 
-set xlabel "r (mm)"
-set xrange [0.2:0.4]
+set xrange [0.:1.]
+set xtics("endocardium" 0, "epicardium" 1)
 
-set ylabel "'''+comp_name+''' strain (%)"
-set yrange [-40:40]
+''')
+        for k_comp in comp_indx:
+            plotfile.write('''\
+set ylabel "'''+comp_names[k_comp]+''' strain (%)"
+''')
+            yrange = yranges[k_comp]
+            if (yrange > 0):
+                plotfile.write('''\
+set yrange [-'''+str(yrange)+''':'''+str(yrange)+''']
 
 plot 0 linecolor rgb "black" notitle,\\
 ''')
-        if (ref_folder is not None) and (ref_basename is not None):
-            plotfile.write('''\
-     "'''+ref_folder+'''/'''+ref_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linecolor "black" pointtype 1 pointsize 1 linewidth 3 notitle'''+(len(working_basenames)>0)*(''',\\
+            if (ref_folder is not None) and (ref_basename is not None):
+                plotfile.write('''\
+     "'''+ref_folder+'''/'''+ref_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(k_frame)+''' with points linecolor "black" linewidth 3 pointtype 1 pointsize 1 notitle'''+(len(working_basenames)>0)*(''',\\
 ''')+(len(working_basenames)==0)*('''
 
 '''))
-        for k_basename in range(len(working_basenames)):
-            working_basename = working_basenames[k_basename]
-            plotfile.write('''\
-     "'''+working_folder+'''/'''+working_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(index)+''' with points linestyle '''+str(k_basename+1)+''' pointtype 1 pointsize 1 linewidth 3 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
+            for k_basename in range(len(working_basenames)):
+                working_basename = working_basenames[k_basename]
+                plotfile.write('''\
+     "'''+working_folder+'''/'''+working_basename+'''-strains_vs_radius.dat" using ($2):(100*$'''+str(3+k_comp)+''') index '''+str(k_frame)+''' with points linestyle '''+str(k_basename+1)+''' linewidth 3 pointtype 1 pointsize 1 title "'''+working_basename+'''"'''+(k_basename<len(working_basenames)-1)*(''',\\
 ''')+(k_basename==len(working_basenames)-1)*('''
 
 '''))
 
     plotfile.close()
 
-    os.system("gnuplot plot_strains_vs_radius"+("-"+suffix)*(suffix!="")+".plt")
+    os.system("gnuplot "+plotfile_basename+".plt")

dolfin_warp/plot_twist_vs_height.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import os
+
+################################################################################
+
+def plot_twist_vs_height(
+        working_folder,
+        working_basename,
+        beta_range=15.,
+        suffix="",
+        verbose=1):
+
+    lines = open(working_folder+"/"+working_basename+"-strains.dat").readlines()[1:]
+    n_frames = len(lines)
+
+    if (suffix is None):
+        plotfile_basename = working_folder+"/"+working_basename+"-twist_vs_height"
+    else:
+        plotfile_basename = "plot_twist_vs_height"+("-"+suffix)*(suffix!="")
+    plotfile = open(plotfile_basename+".plt", "w")
+
+    plotfile.write('''\
+set terminal pdf enhanced
+
+load "Set1.plt"
+
+set output "'''+plotfile_basename+'''.pdf"
+
+set key off
+# set key box textcolor variable width +0
+
+set grid
+
+set xlabel "beta (deg)"
+xrange = '''+str(beta_range)+'''
+set xrange [-xrange:+xrange]
+
+# set ylabel "ll ()"
+set yrange [0.:1.]
+set ytics("apex" 0, "base" 1)
+
+# f(x) = a*x + b
+
+# g(x) = x/c - d/c
+# h(x) = c*x + d
+
+g(x) = c*x + d
+h(x) = x/c - d/c
+
+FIT_MAXITER = 10
+
+datafile = "'''+working_folder+'''/'''+working_basename+'''-twist_vs_height.dat"
+
+''')
+    for k_frame in range(n_frames):
+        plotfile.write('''\
+# a = 1.
+# b = 0.5
+
+# c = 1.
+# d = 0.5
+
+c = 1.
+d = -0.5*c
+
+set title "index '''+str(k_frame)+'''"
+# fit f(x) datafile using 3:2 index '''+str(k_frame)+''' via a,b
+'''+('''# ''')*(k_frame==0)+'''fit g(x) datafile using 2:3 index '''+str(k_frame)+''' via c,d
+plot datafile using 3:2 index '''+str(k_frame)+''' notitle, h(x) notitle
+
+''')
+
+    plotfile.close()
+
+    os.system("gnuplot "+plotfile_basename+".plt")

dolfin_warp/warp.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import dolfin
+
+import dolfin_warp as dwarp
+
+################################################################################
+
+def warp(
+        working_folder                              : str,
+        working_basename                            : str,
+        images_folder                               : str,
+        images_basename                             : str,
+        images_grad_basename                        : str         = None                                ,
+        images_ext                                  : str         = "vti"                               , # vti, vtk
+        images_n_frames                             : int         = None                                ,
+        images_ref_frame                            : int         = 0                                   ,
+        images_quadrature                           : int         = None                                ,
+        images_quadrature_from                      : str         = "points_count"                      , # points_count, integral
+        images_static_scaling                       : bool        = False                               ,
+        images_dynamic_scaling                      : bool        = False                               ,
+        images_char_func                            : bool        = True                                ,
+        images_is_cone                              : bool        = False                               ,
+        mesh                                        : dolfin.Mesh = None                                ,
+        mesh_folder                                 : str         = None                                ,
+        mesh_basename                               : str         = None                                ,
+        mesh_degree                                 : int         = 1                                   ,
+        kinematics_type                             : str         = "full"                              , # full, reduced
+        reduced_kinematics_model                    : str         = "translation+rotation+scaling+shear", # translation, rotation, scaling, shear, translation+rotation+scaling+shear, etc.
+        regul_type                                  : str         = "continuous-equilibrated"           , # continuous-linear-equilibrated, continuous-linear-elastic, continuous-equilibrated, continuous-elastic, continuous-hyperelastic, discrete-simple-equilibrated, discrete-simple-elastic, discrete-linear-equilibrated, discrete-linear-tractions, discrete-linear-tractions-normal, discrete-linear-tractions-tangential, discrete-linear-tractions-normal-tangential, discrete-equilibrated, discrete-tractions, discrete-tractions-normal, discrete-tractions-tangential, discrete-tractions-normal-tangential
+        regul_types                                 : list        = None                                ,
+        regul_model                                 : str         = "ogdenciarletgeymonatneohookean"    , # hooke, kirchhoff, ogdenciarletgeymonatneohookean, ogdenciarletgeymonatneohookeanmooneyrivlin
+        regul_models                                : list        = None                                ,
+        regul_quadrature                            : int         = None                                ,
+        regul_level                                 : float       = 0.                                  ,
+        regul_levels                                : list        = None                                ,
+        regul_poisson                               : float       = 0.                                  ,
+        regul_b                                     : float       = None                                ,
+        regul_volume_subdomain_data                               = None                                ,
+        regul_volume_subdomain_id                                 = None                                ,
+        regul_surface_subdomain_data                              = None                                ,
+        regul_surface_subdomain_id                                = None                                ,
+        relax_type                                  : str         = None                                , # None, constant, aitken, backtracking, gss
+        relax_backtracking_factor                   : float       = None                                ,
+        relax_tol                                   : float       = None                                ,
+        relax_n_iter_max                            : int         = None                                ,
+        relax_must_advance                          : bool        = None                                ,
+        normalize_energies                          : bool        = False                               ,
+        initialize_reduced_U_from_file              : bool        = False                               ,
+        initialize_reduced_U_filename               : str         = None                                ,
+        initialize_U_from_file                      : bool        = False                               ,
+        initialize_U_folder                         : str         = None                                ,
+        initialize_U_basename                       : str         = None                                ,
+        initialize_U_ext                            : str         = "vtu"                               ,
+        initialize_U_array_name                     : str         = "displacement"                      ,
+        initialize_U_method                         : str         = "dofs_transfer"                     , # dofs_transfer, interpolation, projection
+        register_ref_frame                          : bool        = False                               ,
+        iteration_mode                              : str         = "normal"                            , # normal, loop
+        gimic                                       : bool        = False                               ,
+        gimic_texture                               : str         = "no"                                ,
+        gimic_resample                              : int         = 1                                   ,
+        nonlinearsolver                             : str         = "newton"                            , # None, newton, CMA
+        tol_res_rel                                 : float       = None                                ,
+        tol_dU                                      : float       = None                                ,
+        tol_dU_rel                                  : float       = None                                ,
+        n_iter_max                                  : int         = 100                                 ,
+        continue_after_fail                         : bool        = False                               ,
+        write_qois_limited_precision                : bool        = False                               ,
+        write_VTU_files                             : bool        = True                                ,
+        write_VTU_files_with_preserved_connectivity : bool        = False                               ,
+        write_XML_files                             : bool        = False                               ,
+        print_iterations                            : bool        = False                               ,
+        silent                                      : bool        = False                               ):
+
+    if (regul_types is not None):
+        if (regul_models is not None):
+            assert (len(regul_models) == len(regul_types))
+        else:
+            regul_models = [regul_model]*len(regul_types)
+        if (regul_levels is not None):
+            assert (len(regul_levels) == len(regul_types))
+        else:
+            regul_levels = [regul_level]*len(regul_types)
+    else:
+        assert (regul_type is not None)
+        if ("tractions" in regul_type):
+            if (regul_type.startswith("discrete-linear-equilibrated-")):
+                regul_types = ["discrete-linear-equilibrated"]
+                if (regul_type == "discrete-linear-equilibrated-tractions"):
+                    regul_types += ["discrete-linear-tractions"]
+                elif (regul_type == "discrete-linear-equilibrated-tractions-normal"):
+                    regul_types += ["discrete-linear-tractions-normal"]
+                elif (regul_type == "discrete-linear-equilibrated-tractions-tangential"):
+                    regul_types += ["discrete-linear-tractions-tangential"]
+                elif (regul_type == "discrete-linear-equilibrated-tractions-normal-tangential"):
+                    regul_types += ["discrete-linear-tractions-normal-tangential"]
+            elif (regul_type.startswith("discrete-equilibrated-")):
+                regul_types = ["discrete-equilibrated"]
+                if (regul_type == "discrete-equilibrated-tractions"):
+                    regul_types += ["discrete-tractions"]
+                elif (regul_type == "discrete-equilibrated-tractions-normal"):
+                    regul_types += ["discrete-tractions-normal"]
+                elif (regul_type == "discrete-equilibrated-tractions-tangential"):
+                    regul_types += ["discrete-tractions-tangential"]
+                elif (regul_type == "discrete-equilibrated-tractions-normal-tangential"):
+                    regul_types += ["discrete-tractions-normal-tangential"]
+            else: assert (0), "Unknown regul_type ("+str(regul_type)+"). Aborting."
+            regul_models = [regul_model  ]*2
+            regul_levels = [regul_level/2]*2
+        else:
+            regul_types  = [regul_type ]
+            regul_models = [regul_model]
+            regul_levels = [regul_level]
+    # print (regul_types)
+    # print (regul_models)
+    # print (regul_levels)
+
+    if (kinematics_type == "full"):
+        assert (initialize_reduced_U_from_file is False),\
+            "Should use \"initialize_U_from_file\", not \"initialize_reduced_U_from_file\". Aborting."
+        problem = dwarp.FullKinematicsWarpingProblem(
+            mesh=mesh,
+            mesh_folder=mesh_folder,
+            mesh_basename=mesh_basename,
+            U_degree=mesh_degree,
+            silent=silent)
+    elif (kinematics_type == "reduced"):
+        assert (initialize_U_from_file is False),\
+            "Not implemented. Aborting."
+        problem = dwarp.ReducedKinematicsWarpingProblem(
+            mesh=mesh,
+            mesh_folder=mesh_folder,
+            mesh_basename=mesh_basename,
+            model=reduced_kinematics_model,
+            silent=silent)
+    else:
+        assert (0), "\"kinematics_type\" (="+str(kinematics_type)+") must be \"full\" or \"reduced\". Aborting."
+
+    images_series = dwarp.ImagesSeries(
+        folder=images_folder,
+        basename=images_basename,
+        grad_basename=images_grad_basename,
+        n_frames=images_n_frames,
+        ext=images_ext,
+        printer=problem.printer)
+
+    if (images_quadrature is None):
+        problem.printer.print_str("Computing quadrature degree…")
+        problem.printer.inc()
+        if (images_quadrature_from == "points_count"):
+            images_quadrature = dwarp.compute_quadrature_degree_from_points_count(
+                image_filename=images_series.get_image_filename(k_frame=images_ref_frame),
+                mesh=problem.mesh,
+                verbose=1)
+        elif (images_quadrature_from == "integral"):
+            images_quadrature = dwarp.compute_quadrature_degree_from_integral(
+                image_filename=images_series.get_image_filename(k_frame=images_ref_frame),
+                mesh=problem.mesh,
+                verbose=1)
+        else:
+            assert (0), "\"images_quadrature_from\" (="+str(images_quadrature_from)+") must be \"points_count\" or \"integral\". Aborting."
+        problem.printer.print_var("images_quadrature",images_quadrature)
+        problem.printer.dec()
+
+    image_w = 1.-sum(regul_levels)
+    assert (image_w > 0.),\
+        "1.-sum(regul_levels) must be positive. Aborting."
+
+    if (gimic):
+        generated_image_energy = dwarp.GeneratedImageContinuousEnergy(
+            problem=problem,
+            images_series=images_series,
+            quadrature_degree=images_quadrature,
+            texture=gimic_texture,
+            w=image_w,
+            ref_frame=images_ref_frame,
+            resample=gimic_resample)
+        problem.add_image_energy(generated_image_energy)
+    else:
+        warped_image_energy = dwarp.WarpedImageContinuousEnergy(
+            problem=problem,
+            images_series=images_series,
+            quadrature_degree=images_quadrature,
+            w=image_w,
+            ref_frame=images_ref_frame,
+            w_char_func=images_char_func,
+            im_is_cone=images_is_cone,
+            static_scaling=images_static_scaling,
+            dynamic_scaling=images_dynamic_scaling)
+        problem.add_image_energy(warped_image_energy)
+
+    for regul_type, regul_model, regul_level in zip(regul_types, regul_models, regul_levels):
+        if (regul_level>0):
+            name_suffix  = ""
+            name_suffix += ("_"+    regul_type  )*(len(regul_types )>1)
+            name_suffix += ("_"+    regul_model )*(len(regul_models)>1)
+            name_suffix += ("_"+str(regul_level))*(len(regul_levels)>1)
+            regul_b_ = None
+            if regul_type.startswith("continuous"):
+                regularization_energy_type = dwarp.RegularizationContinuousEnergy
+                if regul_type.startswith("continuous-linear"):
+                    regul_type_ = regul_type.split("-",2)[2]
+                else:
+                    regul_type_ = regul_type.split("-",1)[1]
+            elif regul_type.startswith("discrete-simple"):
+                regularization_energy_type = dwarp.SimpleRegularizationDiscreteEnergy
+                regul_type_ = regul_type.split("-",2)[2]
+            elif regul_type.startswith("discrete"):
+                if ("equilibrated" in regul_type):
+                    regularization_energy_type = dwarp.VolumeRegularizationDiscreteEnergy
+                    regul_b_ = regul_b
+                elif ("tractions" in regul_type):
+                    regularization_energy_type = dwarp.SurfaceRegularizationDiscreteEnergy
+                else: assert (0), "regul_type (= "+str(regul_type)+") unknown. Aborting."
+                if regul_type.startswith("discrete-linear"):
+                    regul_type_ = regul_type.split("-",2)[2]
+                else:
+                    regul_type_ = regul_type.split("-",1)[1]
+            else: assert (0), "regul_type (= "+str(regul_type)+") unknown. Aborting."
+            regularization_energy = regularization_energy_type(
+                name="reg"+name_suffix,
+                problem=problem,
+                w=regul_level,
+                type=regul_type_,
+                model=regul_model,
+                poisson=regul_poisson,
+                b_fin=regul_b_,
+                volume_subdomain_data=regul_volume_subdomain_data,
+                volume_subdomain_id=regul_volume_subdomain_id,
+                surface_subdomain_data=regul_surface_subdomain_data,
+                surface_subdomain_id=regul_surface_subdomain_id,
+                quadrature_degree=regul_quadrature)
+            problem.add_regul_energy(regularization_energy)
+
+    if (normalize_energies):
+        dwarp.compute_energies_normalization(
+            problem=problem,
+            verbose=1)
+
+    if (nonlinearsolver == "newton"):
+        solver = dwarp.NewtonNonlinearSolver(
+            problem=problem,
+            parameters={
+                "working_folder":working_folder,
+                "working_basename":working_basename,
+                "relax_type":relax_type,
+                "relax_backtracking_factor":relax_backtracking_factor,
+                "relax_tol":relax_tol,
+                "relax_n_iter_max":relax_n_iter_max,
+                "relax_must_advance":relax_must_advance,
+                "tol_res_rel":tol_res_rel,
+                "tol_dU":tol_dU,
+                "tol_dU_rel":tol_dU_rel,
+                "n_iter_max":n_iter_max,
+                "write_iterations":print_iterations})
+    elif (nonlinearsolver == "cma"):
+        assert (relax_type is None),\
+            "Not implemented. Aborting."
+        solver = dwarp.CMANonlinearSolver(
+            problem=problem,
+            parameters={
+                "working_folder":working_folder,
+                "working_basename":working_basename,
+                "write_iterations":print_iterations})
+
+    image_iterator = dwarp.ImageIterator(
+        problem=problem,
+        solver=solver,
+        parameters={
+            "working_folder":working_folder,
+            "working_basename":working_basename,
+            "register_ref_frame":register_ref_frame,
+            "initialize_reduced_U_from_file":initialize_reduced_U_from_file,
+            "initialize_reduced_U_filename":initialize_reduced_U_filename,
+            "initialize_U_from_file":initialize_U_from_file,
+            "initialize_U_folder":initialize_U_folder,
+            "initialize_U_basename":initialize_U_basename,
+            "initialize_U_ext":initialize_U_ext,
+            "initialize_U_array_name":initialize_U_array_name,
+            "initialize_U_method":initialize_U_method,
+            "write_qois_limited_precision":write_qois_limited_precision,
+            "write_VTU_files":write_VTU_files,
+            "write_VTU_files_with_preserved_connectivity":write_VTU_files_with_preserved_connectivity,
+            "write_XML_files":write_XML_files,
+            "iteration_mode":iteration_mode,
+            "continue_after_fail":continue_after_fail})
+
+    success = image_iterator.iterate()
+
+    problem.close()
+
+    return success
+
+fedic2 = warp
+
+########################################################################
+
+if (__name__ == "__main__"):
+    import fire
+    fire.Fire(warp)

dolfin_warp/warp_and_refine.py

+#coding=utf8
+
+################################################################################
+###                                                                          ###
+### Created by Martin Genet, 2016-2024                                       ###
+###                                                                          ###
+### École Polytechnique, Palaiseau, France                                   ###
+###                                                                          ###
+################################################################################
+
+import dolfin
+
+import dolfin_warp as dwarp
+
+################################################################################
+
+def warp_and_refine(
+        working_folder               : str,
+        working_basename             : str,
+        images_folder                : str,
+        images_basename              : str,
+        images_quadrature            : int         = None                            ,
+        images_quadrature_from       : str         = "points_count"                  , # points_count, integral
+        mesh                         : dolfin.Mesh = None                            ,
+        refinement_levels            : list        = [0]                             ,
+        meshes                       : list        = None                            ,
+        mesh_folder                  : str         = None                            ,
+        mesh_basenames               : list        = None                            ,
+        regul_type                   : str         = "continuous-equilibrated"       , # continuous-equilibrated, continuous-elastic, continuous-hyperelastic, discrete-linear-equilibrated, discrete-linear-elastic, discrete-equilibrated, discrete-tractions, discrete-tractions-normal, discrete-tractions-tangential, discrete-tractions-normal-tangential
+        regul_types                  : list        = None                            ,
+        regul_model                  : str         = "ogdenciarletgeymonatneohookean", # hooke, kirchhoff, ogdenciarletgeymonatneohookean, ogdenciarletgeymonatneohookeanmooneyrivlin
+        regul_models                 : list        = None                            ,
+        regul_level                  : float       = 0.                              ,
+        regul_levels                 : list        = None                            ,
+        regul_poisson                : float       = 0.                              ,
+        regul_b                      : float       = None                            ,
+        regul_volume_subdomain_data                = None                            ,
+        regul_volume_subdomain_id                  = None                            ,
+        regul_surface_subdomain_data               = None                            ,
+        regul_surface_subdomain_id                 = None                            ,
+        relax_type                   : str         = None                            , # constant, aitken, backtracking, gss
+        relax_tol                    : float       = None                            ,
+        relax_n_iter_max             : int         = None                            ,
+        normalize_energies           : bool        = False                           ,
+        tol_dU                       : float       = None                            ,
+        n_iter_max                   : int         = 100                             ,
+        continue_after_fail          : bool        = False                           ,
+        write_qois_limited_precision : bool        = False                           ,
+        print_iterations             : bool        = False                           ,
+        silent                       : bool        = False                           ):
+
+    if (meshes is None):
+        meshes = []
+        if (mesh is not None) and (refinement_levels is not None):
+            for refinement_level in refinement_levels:
+                mesh_for_warp = dolfin.Mesh(mesh)
+                for _ in range(refinement_level):
+                    mesh_for_warp = dolfin.refine(mesh_for_warp)
+                meshes += [mesh_for_warp]
+        elif (mesh_folder is not None) and (mesh_basenames is not None):
+            for mesh_basename in mesh_basenames:
+                mesh_filename = mesh_folder+"/"+mesh_basename+".xml"
+                meshes += [dolfin.Mesh(mesh_filename)]
+    
+    regul_volume_subdomain_data_lst = regul_volume_subdomain_data
+    regul_volume_subdomain_id_lst = regul_volume_subdomain_id
+    regul_surface_subdomain_data_lst = regul_surface_subdomain_data
+    regul_surface_subdomain_id_lst = regul_surface_subdomain_id
+    
+    for k_mesh, mesh_for_warp in enumerate(meshes):
+        working_basename_for_warp  = working_basename
+        working_basename_for_warp += "-refine="+str(k_mesh)
+
+        if (k_mesh == 0):
+            initialize_U_from_file = False
+
+            working_basename_for_init = None
+ 
+        else:
+            initialize_U_from_file = True
+
+            working_basename_for_init  = working_basename
+            working_basename_for_init += "-refine="+str(k_mesh-1)
+
+
+        if regul_volume_subdomain_data:
+            regul_volume_subdomain_data = regul_volume_subdomain_data_lst[k_mesh]   
+            regul_volume_subdomain_id = regul_volume_subdomain_id_lst[k_mesh] 
+        if regul_surface_subdomain_data:
+            regul_surface_subdomain_data = regul_surface_subdomain_data_lst[k_mesh]   
+            regul_surface_subdomain_id = regul_surface_subdomain_id_lst[k_mesh]         
+
+        success = dwarp.warp(
+            working_folder                              = working_folder,
+            working_basename                            = working_basename_for_warp,
+            images_folder                               = images_folder,
+            images_basename                             = images_basename,
+            images_quadrature                           = images_quadrature,
+            images_quadrature_from                      = images_quadrature_from,
+            mesh                                        = mesh_for_warp,
+            regul_type                                  = regul_type,
+            regul_types                                 = regul_types,
+            regul_model                                 = regul_model,
+            regul_models                                = regul_models,
+            regul_level                                 = regul_level,
+            regul_levels                                = regul_levels,
+            regul_poisson                               = regul_poisson,
+            regul_b                                     = regul_b,
+            regul_volume_subdomain_data                 = regul_volume_subdomain_data,
+            regul_volume_subdomain_id                   = regul_volume_subdomain_id,
+            regul_surface_subdomain_data                = regul_surface_subdomain_data,
+            regul_surface_subdomain_id                  = regul_surface_subdomain_id,
+            relax_type                                  = relax_type,
+            relax_tol                                   = relax_tol,
+            relax_n_iter_max                            = relax_n_iter_max,
+            normalize_energies                          = normalize_energies,
+            tol_dU                                      = tol_dU,
+            n_iter_max                                  = n_iter_max,
+            continue_after_fail                         = continue_after_fail,
+            initialize_U_from_file                      = initialize_U_from_file,
+            initialize_U_folder                         = working_folder,
+            initialize_U_basename                       = working_basename_for_init,
+            initialize_U_ext                            = "vtu",
+            initialize_U_array_name                     = "displacement",
+            initialize_U_method                         = "projection", # dofs_transfer, interpolation, projection
+            write_qois_limited_precision                = write_qois_limited_precision,
+            write_VTU_files                             = True,
+            write_VTU_files_with_preserved_connectivity = True,
+            write_XML_files                             = True,
+            print_iterations                            = print_iterations,
+            silent                                      = silent)
+
+        if not (success) and not (continue_after_fail):
+            break
+
+    return success
+
+########################################################################
+
+if (__name__ == "__main__"):
+    import fire
+    fire.Fire(warp_and_refine)

fedic.py

-#coding=utf8
-
-########################################################################
-###                                                                  ###
-### Created by Martin Genet, 2016                                    ###
-###                                                                  ###
-### École Polytechnique, Palaiseau, France                           ###
-###                                                                  ###
-########################################################################
-
-import dolfin
-import glob
-import numpy
-import os
-
-import myFEniCSPythonLibrary as myFEniCS
-import myVTKPythonLibrary    as myVTK
-
-########################################################################
-
-def fedic(
-        working_folder,
-        working_basename,
-        images_folder,
-        images_basename,
-        mesh_folder,
-        mesh_basename,
-        images_dimension=3,
-        k_ref=0,
-        penalty=0.9,
-        use_I0_tangent=0,
-        relax_type="const",
-        relax_init=0.9,
-        tol_disp=1e-3,
-        n_iter_max=100,
-        print_iterations=0,
-        continue_after_fail=0,
-        verbose=1):
-
-    if not os.path.exists(working_folder):
-        os.mkdir(working_folder)
-
-    print "Loading mesh…"
-    mesh = dolfin.Mesh(mesh_folder+"/"+mesh_basename+".xml")
-
-    fs = dolfin.VectorFunctionSpace(
-        mesh=mesh,
-        family="Lagrange",
-        degree=1)
-    dX = dolfin.dx(mesh)
-
-    print "Computing quadrature degree…"
-    degree = myFEniCS.compute_quadrature_degree(
-        image_filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-        dX=dX,
-        image_dimension=images_dimension)
-    print "degree = " + str(degree)
-
-    fe = dolfin.FiniteElement(
-        family="Quadrature",
-        degree=degree)
-
-    print "Loading reference image…"
-    if (images_dimension == 2):
-        I0 = myFEniCS.ExprIm2(
-            filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-            element=fe)
-        if (use_I0_tangent):
-            DXI0 = myFEniCS.ExprGradXIm2(
-                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-                element=fe)
-            DYI0 = myFEniCS.ExprGradYIm2(
-                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-                element=fe)
-    elif (images_dimension == 3):
-        I0 = myFEniCS.ExprIm3(
-            filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-            element=fe)
-        if (use_I0_tangent):
-            DXI0 = myFEniCS.ExprGradXIm3(
-                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-                element=fe)
-            DYI0 = myFEniCS.ExprGradYIm3(
-                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-                element=fe)
-            DZI0 = myFEniCS.ExprGradZIm3(
-                filename=images_folder+"/"+images_basename+"_"+str(k_ref).zfill(2)+".vti",
-                element=fe)
-    else:
-        assert (0), "images_dimension must be 2 or 3. Aborting."
-    I0_norm = dolfin.assemble(I0**2 * dX)**(1./2)
-
-    print "Defining functions…"
-    penalty = dolfin.Constant(penalty)
-    U = dolfin.Function(
-        fs,
-        name="displacement")
-    dU = dolfin.Function(
-        fs,
-        name="ddisplacement")
-    ddU = dolfin.TrialFunction(fs)
-    V = dolfin.TestFunction(fs)
-
-    print "Printing initial solution…"
-    file_pvd = dolfin.File(working_folder+"/"+working_basename+"_.pvd")
-    for vtu in glob.glob(working_folder+"/"+working_basename+"_*.vtu"):
-        os.remove(vtu)
-    file_pvd << (U, float(k_ref))
-
-    print "Defining variational forms…"
-    if (images_dimension == 2):
-        I1 = myFEniCS.ExprDefIm2(
-            U=U,
-            element=fe)
-        DXI1 = myFEniCS.ExprGradXDefIm2(
-            U=U,
-            element=fe)
-        DYI1 = myFEniCS.ExprGradYDefIm2(
-            U=U,
-            element=fe)
-        if (use_I0_tangent):
-            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0)), ddU),
-                                            dolfin.dot(dolfin.as_vector((DXI0, DYI0)),   V))*dX\
-              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
-                                            dolfin.grad(  V))*dX
-        else:
-            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1)), ddU),
-                                            dolfin.dot(dolfin.as_vector((DXI1, DYI1)),   V))*dX\
-              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
-                                            dolfin.grad(  V))*dX
-        b =     penalty  * dolfin.inner(I0-I1,
-                                        dolfin.dot(dolfin.as_vector((DXI1, DYI1)), V))*dX\
-          - (1.-penalty) * dolfin.inner(dolfin.grad(U),
-                                        dolfin.grad(V))*dX
-    elif (images_dimension == 3):
-        I1 = myFEniCS.ExprDefIm3(
-            U=U,
-            element=fe)
-        DXI1 = myFEniCS.ExprGradXDefIm3(
-            U=U,
-            element=fe)
-        DYI1 = myFEniCS.ExprGradYDefIm3(
-            U=U,
-            element=fe)
-        DZI1 = myFEniCS.ExprGradZDefIm3(
-            U=U,
-            element=fe)
-        if (use_I0_tangent):
-            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)), ddU),
-                                            dolfin.dot(dolfin.as_vector((DXI0, DYI0, DZI0)),   V))*dX\
-              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
-                                            dolfin.grad(  V))*dX
-        else:
-            a =     penalty  * dolfin.inner(dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), ddU),
-                                            dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)),   V))*dX\
-              + (1.-penalty) * dolfin.inner(dolfin.grad(ddU),
-                                            dolfin.grad(  V))*dX
-        b =     penalty  * dolfin.inner(I0-I1,
-                                        dolfin.dot(dolfin.as_vector((DXI1, DYI1, DZI1)), V))*dX\
-          - (1.-penalty) * dolfin.inner(dolfin.grad(U),
-                                        dolfin.grad(V))*dX
-    else:
-        assert (0), "images_dimension must be 2 or 3. Aborting."
-
-    # linear system
-    A = None
-    B = None
-    if (use_I0_tangent):
-        A = dolfin.assemble(a, tensor=A)
-
-    print "Checking number of frames…"
-    n_frames = len(glob.glob(images_folder+"/"+images_basename+"_*.vti"))
-    print "n_frames = " + str(n_frames)
-
-    assert (abs(k_ref) < n_frames)
-    k_ref = k_ref%n_frames
-
-    n_iter_tot = 0
-    for forward_or_backward in ["forward", "backward"]: # not sure if this still works…
-        print "forward_or_backward = " + forward_or_backward
-
-        if (forward_or_backward == "forward"):
-            k_next_frames = range(k_ref+1, n_frames, +1)
-        elif (forward_or_backward == "backward"):
-            k_next_frames = range(k_ref-1, -1, -1)
-        print "k_next_frames = " + str(k_next_frames)
-
-        for k_next_frame in k_next_frames:
-            print "k_next_frame = " + str(k_next_frame)
-
-            if (print_iterations):
-                if (os.path.exists(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf")):
-                    os.remove(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf")
-                file_dat_iter = open(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".dat", "w")
-
-                file_pvd_iter = dolfin.File(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_.pvd")
-                for vtu in glob.glob(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_*.vtu"):
-                    os.remove(vtu)
-                file_pvd_iter << (U, 0.)
-
-            print "Loading image and image gradient…"
-            if (images_dimension == 2):
-                I1.init_image(  filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-                DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-                DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-            elif (images_dimension == 3):
-                I1.init_image(  filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-                DXI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-                DYI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-                DZI1.init_image(filename=images_folder+"/"+images_basename+"_"+str(k_next_frame).zfill(2)+".vti")
-            else:
-                assert (0), "images_dimension must be 2 or 3. Aborting."
-
-            print "Running registration…"
-            k_iter = 0
-            relax = relax_init
-            while (True):
-                print "k_iter = " + str(k_iter)
-                n_iter_tot += 1
-
-                # linear system
-                if not (use_I0_tangent):
-                    A = dolfin.assemble(a, tensor=A)
-                #print "A = " + str(A.array())
-                if (k_iter > 0):
-                    if (k_iter == 1):
-                        B_old = B.copy()
-                    elif (k_iter > 1):
-                        B_old[:] = B[:]
-                    #print "B_old = " + str(B_old[0])
-                B = dolfin.assemble(b, tensor=B)
-                #print "B = " + str(B[0])
-                dolfin.solve(A, dU.vector(), B)
-                #print "dU = " + str(dU.vector().array())
-
-                # relaxation
-                if (relax_type == "aitken") and (k_iter > 0):
-                    if (k_iter == 1):
-                        dB = B - B_old
-                    elif (k_iter > 1):
-                        dB[:] = B[:] - B_old[:]
-                    relax *= (-1.) * B_old.inner(dB) / dB.inner(dB)
-                print "relax = " + str(relax)
-
-                # solution update
-                U.vector()[:] += relax * dU.vector()[:]
-                #U.vector().axpy(relax, dU.vector())
-
-                if (print_iterations):
-                    #print "U = " + str(U.vector().array())
-                    file_pvd_iter << (U, float(k_iter+1))
-
-                # displacement error
-                dU_max = numpy.linalg.norm(dU.vector().array(), ord=numpy.Inf)
-                U_max = numpy.linalg.norm(U.vector().array(), ord=numpy.Inf)
-                err_disp_abs = max(dU_max, U_max)
-                U_norm = numpy.linalg.norm(U.vector().array())
-                err_disp_rel = dU_max/U_norm
-                print "err_disp_abs = " + str(err_disp_abs)
-                print "err_disp_rel = " + str(err_disp_rel)
-
-                # image error
-                I1I0_norm = dolfin.assemble((I1-I0)**2 * dX)**(1./2)
-                err_im = I1I0_norm/I0_norm
-                print "err_im = " + str(err_im)
-
-                if (print_iterations):
-                    file_dat_iter.write(" ".join([str(val) for val in [k_iter, err_disp_abs, err_disp_rel, err_im]])+"\n")
-
-                # exit test
-                if (err_disp_rel < tol_disp) or (err_disp_abs < tol_disp):
-                    print "Nonlinear solver converged…"
-                    success = True
-                    break
-
-                if (k_iter >= n_iter_max-1):
-                    print "Warning! Nonlinear solver failed to converge…"
-                    success = False
-                    break
-
-                # increment counter
-                k_iter += 1
-
-            if (print_iterations):
-                os.remove(working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+"_.pvd")
-                file_dat_iter.close()
-                os.system("gnuplot -e \"set terminal pdf; set output '"+working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".pdf'; set logscale y; plot '"+working_folder+"/"+working_basename+"-frame="+str(k_next_frame).zfill(2)+".dat' using 1:3\"")
-
-            if not (success) and not (continue_after_fail):
-                break
-
-            print "Printing solution…"
-            file_pvd << (U, float(k_next_frame))
-
-        if not (success) and not (continue_after_fail):
-            break
-
-    print "n_iter_tot = " + str(n_iter_tot)
-
-    os.remove(working_folder+"/"+working_basename+"_.pvd")
-
-    return success

image_expressions.py

-#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])
-        #self.evalcount = 0
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        self.s = getScalingFactor(self.image.GetScalarTypeAsString())
-        self.interpolator = myVTK.createImageInterpolator(
-            image=self.image,
-            verbose=0)
-
-    def eval(self, Im, X):
-        #print "Im"
-        self.X[0] = X[0]
-        self.X[1] = X[1]
-        #print "    X = " + str(X)
-        self.interpolator.Interpolate(self.X, Im)
-        #print "    Im = " + str(Im)
-        Im /= self.s
-        #print "    Im = " + str(Im)
-        #self.evalcount += 1
-
-class ExprIm3(dolfin.Expression):
-    def __init__(self, filename, **kwargs):
-        if filename is not None:
-            self.init_image(filename=filename)
-        #self.evalcount = 0
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        self.s = getScalingFactor(self.image.GetScalarTypeAsString())
-        self.interpolator = myVTK.createImageInterpolator(
-            image=self.image,
-            verbose=0)
-
-    def eval(self, Im, X):
-        #print "Im"
-        #print "    X = " + str(X)
-        self.interpolator.Interpolate(X, Im)
-        #print "    Im = " + str(Im)
-        Im /= self.s
-        #print "    Im = " + str(Im)
-        #self.evalcount += 1
-
-class ExprGradIm2(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
-    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)
-        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,)
-
-    def eval(self, GradIm, X):
-        self.X[0] = X[0]
-        self.X[1] = X[1]
-        self.interpolator.Interpolate(self.X, GradIm)
-        GradIm /= self.s
-
-class ExprGradXIm2(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)
-        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 eval(self, GradXIm, X):
-        self.X[0] = X[0]
-        self.X[1] = X[1]
-        GradXIm[0] = self.interpolator.Interpolate(self.X[0], self.X[1], self.X[2], 0)
-        GradXIm /= self.s
-
-class ExprGradYIm2(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)
-        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 eval(self, GradYIm, X):
-        self.X[0] = X[0]
-        self.X[1] = X[1]
-        GradYIm[0] = self.interpolator.Interpolate(self.X[0], self.X[1], self.X[2], 1)
-        GradYIm /= self.s
-
-class ExprGradIm3(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
-    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)
-        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,)
-
-    def eval(self, GradIm, X):
-        self.interpolator.Interpolate(X, GradIm)
-        GradIm /= self.s
-
-class ExprGradXIm3(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)
-        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 eval(self, GradXIm, X):
-        GradXIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 0)
-        GradXIm /= self.s
-
-class ExprGradYIm3(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)
-        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 eval(self, GradYIm, X):
-        GradYIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 1)
-        GradYIm /= self.s
-
-class ExprGradZIm3(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)
-        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 eval(self, GradZIm, X):
-        GradZIm[0] = self.interpolator.Interpolate(X[0], X[1], X[2], 2)
-        GradZIm /= 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.array([float()]*2)
-        self.x = numpy.array([float()]*2+[Z])
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        self.s = getScalingFactor(self.image.GetScalarTypeAsString())
-        self.interpolator = myVTK.createImageInterpolator(
-            image=self.image,
-            verbose=0)
-
-    def eval(self, DefIm, X):
-        #print "DefIm"
-        #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] = X[0] + self.UX[0]
-        self.x[1] = X[1] + self.UX[1]
-        #print "    x = " + str(self.x)
-        #print "    DefIm = " + str(DefIm)
-        self.interpolator.Interpolate(self.x, DefIm)
-        #print "    DefIm = " + str(DefIm)
-        DefIm /= self.s
-        #print "    DefIm = " + str(DefIm)
-
-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.array([float()]*3)
-        self.x = numpy.array([float()]*3)
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        self.s = getScalingFactor(self.image.GetScalarTypeAsString())
-        self.interpolator = myVTK.createImageInterpolator(
-            image=self.image,
-            verbose=0)
-
-    def eval(self, DefIm, X):
-        #print "DefIm"
-        #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)
-        #print "    DefIm = " + str(DefIm)
-        self.interpolator.Interpolate(self.x, DefIm)
-        #print "    DefIm = " + str(DefIm)
-        DefIm /= self.s
-        #print "    DefIm = " + str(DefIm)
-
-class ExprGradDefIm2(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
-    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.array([float()]*2)
-        self.x = numpy.array([float()]*2+[Z])
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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,)
-
-    def eval(self, GradDefIm, X):
-        #print "GradDefIm"
-        #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] = X[0] + self.UX[0]
-        self.x[1] = X[1] + self.UX[1]
-        #print "    x = " + str(self.x)
-        #print "    GradDefIm = " + str(GradDefIm)
-        self.interpolator.Interpolate(self.x, GradDefIm)
-        #print "    GradDefIm = " + str(GradDefIm)
-        GradDefIm /= self.s
-        #print "    GradDefIm = " + str(GradDefIm)
-
-class ExprGradXDefIm2(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.array([float()]*2)
-        self.x = numpy.array([float()]*2+[Z])
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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 eval(self, GradXDefIm, X):
-        #print "GradXDefIm"
-        #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] = X[0] + self.UX[0]
-        self.x[1] = X[1] + self.UX[1]
-        #print "    x = " + str(self.x)
-        #print "    GradXDefIm = " + str(GradXDefIm)
-        GradXDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 0)
-        #print "    GradXDefIm = " + str(GradXDefIm)
-        GradXDefIm /= self.s
-        #print "    GradXDefIm = " + str(GradXDefIm)
-
-class ExprGradYDefIm2(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.array([float()]*2)
-        self.x = numpy.array([float()]*2+[Z])
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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 eval(self, GradYDefIm, X):
-        #print "GradYDefIm"
-        #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] = X[0] + self.UX[0]
-        self.x[1] = X[1] + self.UX[1]
-        #print "    x = " + str(self.x)
-        #print "    GradYDefIm = " + str(GradYDefIm)
-        GradYDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 1)
-        #print "    GradYDefIm = " + str(GradYDefIm)
-        GradYDefIm /= self.s
-        #print "    GradYDefIm = " + str(GradYDefIm)
-
-class ExprGradDefIm3(dolfin.Expression): # for some reason this does not work with quadrature finite elements, hence the following scalar gradients definition
-    def __init__(self, U, filename=None, **kwargs):
-        if filename is not None:
-            self.init_image(filename=filename)
-        self.U = U
-        self.UX = numpy.array([float()]*3)
-        self.x = numpy.array([float()]*3)
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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,)
-
-    def eval(self, GradDefIm, X):
-        #print "GradDefIm"
-        #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)
-        #print "    GradDefIm = " + str(GradDefIm)
-        self.interpolator.Interpolate(self.x, GradDefIm)
-        #print "    GradDefIm = " + str(GradDefIm)
-        GradDefIm /= self.s
-        #print "    GradDefIm = " + str(GradDefIm)
-
-class ExprGradXDefIm3(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.array([float()]*3)
-        self.x = numpy.array([float()]*3)
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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 eval(self, GradXDefIm, X):
-        #print "GradXDefIm"
-        #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)
-        #print "    GradXDefIm = " + str(GradXDefIm)
-        GradXDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 0)
-        #print "    GradXDefIm = " + str(GradXDefIm)
-        GradXDefIm /= self.s
-        #print "    GradXDefIm = " + str(GradXDefIm)
-
-class ExprGradYDefIm3(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.array([float()]*3)
-        self.x = numpy.array([float()]*3)
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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 eval(self, GradYDefIm, X):
-        #print "GradYDefIm"
-        #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)
-        #print "    GradYDefIm = " + str(GradYDefIm)
-        GradYDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 1)
-        #print "    GradYDefIm = " + str(GradYDefIm)
-        GradYDefIm /= self.s
-        #print "    GradYDefIm = " + str(GradYDefIm)
-
-class ExprGradZDefIm3(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.array([float()]*3)
-        self.x = numpy.array([float()]*3)
-
-    def init_image(self, filename):
-        self.image = myVTK.readImage(
-            filename=filename,
-            verbose=0)
-        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 eval(self, GradZDefIm, X):
-        #print "GradZDefIm"
-        #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)
-        #print "    GradZDefIm = " + str(GradZDefIm)
-        GradZDefIm[0] = self.interpolator.Interpolate(self.x[0], self.x[1], self.x[2], 2)
-        #print "    GradZDefIm = " + str(GradZDefIm)
-        GradZDefIm /= self.s
-        #print "    GradZDefIm = " + str(GradZDefIm)
-
-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();
-};
-
-}'''

setup.py

+import datetime
+import os
+import setuptools
+
+# version = os.environ['CI_COMMIT_TAG']
+version = datetime.date.today().strftime("%Y.%m.%d")
+
+setuptools.setup(
+    name="dolfin_warp",
+    version=version,
+    author="Martin Genet",
+    author_email="martin.genet@polytechnique.edu",
+    description=open("README.md", "r").readlines()[1][:-1],
+    long_description = open("README.md", "r").read(),
+    long_description_content_type="text/markdown",
+    url="https://gitlab.inria.fr/mgenet/dolfin_warp",
+    packages=["dolfin_warp"],
+    license="GPLv3",
+    classifiers=[
+        "Programming Language :: Python :: 3",
+        "License :: OSI Approved :: GNU General Public License v3 (GPLv3)",
+        "Operating System :: OS Independent",
+    ],
+    install_requires=["dolfin_mech", "matplotlib", "meshio", "myPythonLibrary", "myVTKPythonLibrary", "numpy", "pandas", "scipy", "vtk", "vtkpython_cbl"],
+)