From b29c99d5b0ddd14e1193229b79822f2338bb3ca3 Mon Sep 17 00:00:00 2001
From: Martin Genet <martin.genet@polytechnique.edu>
Date: Tue, 6 Sep 2016 22:04:50 +0200
Subject: [PATCH] First commit after splitting the main library into multiple
libraries
---
__init__.py | 7 +-
__init__.py.py | 26 --
compute_displacement_error.py | 23 +-
compute_displacements.py | 62 +++
compute_quadrature_degree.py | 43 ++-
compute_strains.py | 106 +++--
compute_unwarped_images.py | 95 +++++
compute_warped_images.py | 110 ++++++
compute_warped_mesh.py | 81 ++++
fedic.py | 305 ++++++++-------
generate_images.py | 665 ++++++++++++++++++++++++++++++++
generate_undersampled_images.py | 164 ++++++++
image_expressions_cpp.py | 79 ++--
image_expressions_py.py | 70 ++--
plot_strains.py | 26 +-
print_tools.py | 42 --
16 files changed, 1576 insertions(+), 328 deletions(-)
delete mode 100755 __init__.py.py
create mode 100644 compute_displacements.py
create mode 100644 compute_unwarped_images.py
create mode 100644 compute_warped_images.py
create mode 100644 compute_warped_mesh.py
create mode 100644 generate_images.py
create mode 100644 generate_undersampled_images.py
delete mode 100644 print_tools.py
diff --git a/__init__.py b/__init__.py
index eb6a001..a4a799d 100644
--- a/__init__.py
+++ b/__init__.py
@@ -2,11 +2,16 @@
from fedic import *
from plot_strains_vs_radius import *
+from generate_undersampled_images import *
from plot_strains import *
from image_expressions_cpp import *
from compute_strains import *
-from print_tools import *
+from compute_warped_mesh import *
+from generate_images import *
from plot_displacement_error import *
+from compute_displacements import *
+from compute_unwarped_images import *
from image_expressions_py import *
from compute_displacement_error import *
+from compute_warped_images import *
from compute_quadrature_degree import *
diff --git a/__init__.py.py b/__init__.py.py
deleted file mode 100755
index fded04a..0000000
--- a/__init__.py.py
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/usr/bin/python
-#coding=utf8
-
-########################################################################
-### ###
-### Created by Martin Genet, 2016 ###
-### ###
-### École Polytechnique, Palaiseau, France ###
-### ###
-########################################################################
-
-import os
-
-########################################################################
-
-os.system("rm -rf *.pyc")
-
-init_file = open("__init__.py", "w")
-init_file.write("#this file was generated by __init__.py.py\n\n")
-
-for (dirpath, dirnames, filenames) in os.walk("."):
- for filename in filenames:
- if filename.endswith(".py") and not filename.startswith("__init__"):
- init_file.write("from " + filename[:-3] + " import *\n")
-
-init_file.close()
diff --git a/compute_displacement_error.py b/compute_displacement_error.py
index e84d6e2..38eaebc 100644
--- a/compute_displacement_error.py
+++ b/compute_displacement_error.py
@@ -8,10 +8,14 @@
### ###
########################################################################
+import dolfin
import glob
import numpy
-import myVTKPythonLibrary as myVTK
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
########################################################################
@@ -26,27 +30,30 @@ def compute_displacement_error(
ref_disp_array_name="displacement",
verbose=1):
- n_frames = len(glob.glob(ref_folder+"/"+ref_basename+"_*."+ref_ext))
- assert (len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)) == n_frames)
- if (verbose): print "n_frames = " + str(n_frames)
+ sol_filenames = glob.glob(sol_folder+"/"+sol_basename+"_[0-9]*."+sol_ext)
+ ref_filenames = glob.glob(ref_folder+"/"+ref_basename+"_[0-9]*."+ref_ext)
- ref_zfill = len(glob.glob(ref_folder+"/"+ref_basename+"_*."+ref_ext)[0].rsplit("_")[-1].split(".")[0])
- sol_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
+ sol_zfill = len(sol_filenames[0].rsplit("_",1)[-1].split(".",1)[0])
+ ref_zfill = len(ref_filenames[0].rsplit("_",1)[-1].split(".",1)[0])
if (verbose): print "ref_zfill = " + str(ref_zfill)
if (verbose): print "sol_zfill = " + str(sol_zfill)
+ n_frames = len(sol_filenames)
+ assert (len(ref_filenames) == n_frames)
+ if (verbose): print "n_frames = " + str(n_frames)
+
error_file = open(sol_folder+"/"+sol_basename+"-displacement_error.dat", "w")
error_file.write("#t e\n")
err_int = numpy.empty(n_frames)
ref_int = numpy.empty(n_frames)
for k_frame in xrange(n_frames):
- ref = myVTK.readUGrid(
+ ref = myvtk.readUGrid(
filename=ref_folder+"/"+ref_basename+"_"+str(k_frame).zfill(ref_zfill)+"."+ref_ext,
verbose=0)
n_points = ref.GetNumberOfPoints()
n_cells = ref.GetNumberOfCells()
- sol = myVTK.readUGrid(
+ sol = myvtk.readUGrid(
filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
verbose=0)
assert (sol.GetNumberOfPoints() == n_points)
diff --git a/compute_displacements.py b/compute_displacements.py
new file mode 100644
index 0000000..f3c1348
--- /dev/null
+++ b/compute_displacements.py
@@ -0,0 +1,62 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2016 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import dolfin
+import glob
+import numpy
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+def compute_displacements(
+ sol_folder,
+ sol_basename,
+ ref_frame,
+ sol_ext="vtu",
+ verbose=1):
+
+ sol_filenames = glob.glob(sol_folder+"/"+sol_basename+"_[0-9]*."+sol_ext)
+ sol_zfill = len(sol_filenames[0].rsplit("_",1)[-1].split(".")[0])
+ n_frames = len(sol_filenames)
+ if (verbose): print "n_frames = "+str(n_frames)
+
+ ref_mesh = myvtk.readUGrid(
+ filename=sol_folder+"/"+sol_basename+"_"+str(ref_frame).zfill(sol_zfill)+"."+sol_ext,
+ verbose=verbose)
+ n_points = ref_mesh.GetNumberOfPoints()
+ n_cells = ref_mesh.GetNumberOfCells()
+
+ ref_disp_farray = myvtk.createFloatArray(name="ref_disp")
+ ref_disp_farray.DeepCopy(ref_mesh.GetPointData().GetVectors())
+
+ warper = vtk.vtkWarpVector()
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ warper.SetInputData(ref_mesh)
+ else:
+ warper.SetInput(ref_mesh)
+ warper.Update()
+ warped_mesh = warper.GetOutput()
+ warped_disp_farray = warped_mesh.GetPointData().GetVectors()
+
+ for k_frame in xrange(n_frames):
+ cur_mesh = myvtk.readUGrid(
+ filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
+ verbose=verbose)
+ cur_disp_farray = cur_mesh.GetPointData().GetVectors()
+ [warped_disp_farray.SetTuple(k_point, numpy.array(cur_disp_farray.GetTuple(k_point)) - numpy.array(ref_disp_farray.GetTuple(k_point))) for k_point in xrange(n_points)]
+ myvtk.writeUGrid(
+ ugrid=warped_mesh,
+ filename=sol_folder+"/"+sol_basename+"__"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
+ verbose=verbose)
diff --git a/compute_quadrature_degree.py b/compute_quadrature_degree.py
index 3c89284..aac60a2 100644
--- a/compute_quadrature_degree.py
+++ b/compute_quadrature_degree.py
@@ -11,8 +11,10 @@
import dolfin
import numpy
-import myVTKPythonLibrary as myVTK
-import myFEniCSPythonLibrary as myFEniCS
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
########################################################################
@@ -24,12 +26,12 @@ def compute_quadrature_degree_from_points_count(
deg_max=20,
verbose=1):
- image = myVTK.readImage(
+ image = myvtk.readImage(
filename=image_filename,
verbose=verbose)
n_points = image.GetNumberOfPoints()
- mesh = myVTK.readUGrid(
+ mesh = myvtk.readUGrid(
filename=mesh_filebasename+"."+mesh_ext,
verbose=verbose)
n_cells = mesh.GetNumberOfCells()
@@ -39,26 +41,26 @@ def compute_quadrature_degree_from_points_count(
generic_cell,
k_cell,
subId,
- dist) = myVTK.getCellLocator(
+ dist) = myvtk.getCellLocator(
mesh=mesh,
verbose=verbose)
point = numpy.empty(3)
- n_pixels = numpy.zeros(n_cells)
+ n_pixels_per_cell = numpy.zeros(n_cells)
for k_point in xrange(n_points):
image.GetPoint(k_point, point)
k_cell = cell_locator.FindCell(point)
if (k_cell == -1): continue
- else: n_pixels[k_cell] += 1
- n_pixels_max = int(max(n_pixels))
- n_pixels_avg = int(sum(n_pixels)/n_cells)
+ else: n_pixels_per_cell[k_cell] += 1
+ n_pixels_per_cell_max = int(max(n_pixels_per_cell))
+ n_pixels_per_cell_avg = int(sum(n_pixels_per_cell)/n_cells)
if (verbose):
- #print "n_pixels = "+str(n_pixels)
- #print "sum(n_pixels) = "+str(sum(n_pixels))
- print "n_pixels_max = "+str(n_pixels_max)
- print "n_pixels_avg = "+str(n_pixels_avg)
+ #print "n_pixels_per_cell = "+str(n_pixels_per_cell)
+ #print "sum(n_pixels_per_cell) = "+str(sum(n_pixels_per_cell))
+ print "n_pixels_per_cell_max = "+str(n_pixels_per_cell_max)
+ print "n_pixels_per_cell_avg = "+str(n_pixels_per_cell_avg)
mesh = dolfin.Mesh(mesh_filebasename+"."+"xml")
@@ -72,8 +74,8 @@ def compute_quadrature_degree_from_points_count(
degree=degree,
quad_scheme="default")).dofmap().dofs())/len(mesh.cells())
if (verbose): print "n_quad = "+str(n_quad)
- #if (n_quad > n_pixels_max): break
- if (n_quad > n_pixels_avg): break
+ #if (n_quad > n_pixels_per_cell_max): break
+ if (n_quad > n_pixels_per_cell_avg): break
return degree
@@ -88,8 +90,11 @@ def compute_quadrature_degree_from_integral(
n_under_tol=1,
verbose=1):
- image_dimension = myVTK.computeImageDimensionality(
- image_filename=image_filename,
+ image = myvtk.readImage(
+ filename=image_filename,
+ verbose=verbose)
+ image_dimension = myvtk.getImageDimensionality(
+ image=image,
verbose=verbose)
if (verbose): print "image_dimension = " + str(image_dimension)
@@ -105,11 +110,11 @@ def compute_quadrature_degree_from_integral(
degree=degree,
quad_scheme="default")
if (image_dimension == 2):
- I0 = myFEniCS.ExprIm2(
+ I0 = ddic.ExprIm2(
filename=image_filename,
element=fe)
elif (image_dimension == 3):
- I0 = myFEniCS.ExprIm3(
+ I0 = ddic.ExprIm3(
filename=image_filename,
element=fe)
else:
diff --git a/compute_strains.py b/compute_strains.py
index 4c4445e..2d25eae 100644
--- a/compute_strains.py
+++ b/compute_strains.py
@@ -8,37 +8,45 @@
### ###
########################################################################
+import dolfin
import glob
import numpy
-import myVTKPythonLibrary as myVTK
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
########################################################################
def compute_strains(
sol_folder,
sol_basename,
- sol_zfill=6,
sol_ext="vtu",
+ ref_frame=None,
disp_array_name="displacement",
- ref_folder=None,
- ref_basename=None,
+ defo_grad_array_name="DeformationGradient",
+ strain_array_name="Strain",
+ mesh_w_local_basis_folder=None,
+ mesh_w_local_basis_basename=None,
CYL_or_PPS="PPS",
write_strains=1,
+ plot_strains=1,
write_strain_vs_radius=0,
verbose=1):
- if (ref_folder is not None) and (ref_basename is not None):
- ref_mesh = myVTK.readUGrid(
- filename=ref_folder+"/"+ref_basename+"-WithLocalBasis.vtk",
+ if (mesh_w_local_basis_folder is not None) and (mesh_w_local_basis_basename is not None):
+ mesh_w_local_basis_filename = mesh_w_local_basis_folder+"/"+mesh_w_local_basis_basename+"-WithLocalBasis.vtk"
+ mesh_w_local_basis = myvtk.readUGrid(
+ filename=mesh_w_local_basis_filename,
verbose=verbose)
- ref_n_cells = ref_mesh.GetNumberOfCells()
- if (verbose): print "ref_n_cells = " + str(ref_n_cells)
+ mesh_w_local_basis_n_cells = mesh_w_local_basis.GetNumberOfCells()
+ if (verbose): print "mesh_w_local_basis_n_cells = " + str(mesh_w_local_basis_n_cells)
- if (ref_mesh.GetCellData().HasArray("sector_id")):
- iarray_sector_id = ref_mesh.GetCellData().GetArray("sector_id")
+ if (mesh_w_local_basis.GetCellData().HasArray("sector_id")):
+ iarray_sector_id = mesh_w_local_basis.GetCellData().GetArray("sector_id")
n_sector_ids = 0
- for k_cell in xrange(ref_n_cells):
+ for k_cell in xrange(mesh_w_local_basis_n_cells):
sector_id = int(iarray_sector_id.GetTuple1(k_cell))
if (sector_id < 0): continue
if (sector_id > n_sector_ids-1):
@@ -47,13 +55,15 @@ def compute_strains(
else:
n_sector_ids = 0
else:
- ref_mesh = None
+ mesh_w_local_basis = None
n_sector_ids = 0
- n_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
- if (verbose): print "n_zfill = " + str(n_zfill)
+ sol_filenames = glob.glob(sol_folder+"/"+sol_basename+"_[0-9]*."+sol_ext)
+
+ sol_zfill = len(sol_filenames[0].rsplit("_",1)[-1].split(".")[0])
+ if (verbose): print "sol_zfill = " + str(sol_zfill)
- n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_"+"[0-9]"*sol_zfill+"."+sol_ext))
+ n_frames = len(sol_filenames)
if (verbose): print "n_frames = " + str(n_frames)
if (write_strains):
@@ -64,31 +74,56 @@ def compute_strains(
strain_vs_radius_file = open(sol_folder+"/"+sol_basename+"-strains_vs_radius.dat", "w")
strain_vs_radius_file.write("#t r Err Ecc Ell Erc Erl Ecl\n")
+ if (ref_frame is not None):
+ ref_mesh_filename = sol_folder+"/"+sol_basename+"_"+str(ref_frame).zfill(sol_zfill)+"."+sol_ext
+ ref_mesh = myvtk.readUGrid(
+ filename=ref_mesh_filename,
+ verbose=verbose)
+ myvtk.addDeformationGradients(
+ mesh=ref_mesh,
+ disp_array_name=disp_array_name,
+ verbose=verbose)
+ farray_F0 = ref_mesh.GetCellData().GetArray(defo_grad_array_name)
+
for k_frame in xrange(n_frames):
- mesh = myVTK.readUGrid(
- filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(n_zfill)+"."+sol_ext,
+ mesh_filename = sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext
+ mesh = myvtk.readUGrid(
+ filename=mesh_filename,
verbose=verbose)
n_cells = mesh.GetNumberOfCells()
- if (ref_mesh is not None):
- assert (n_cells == ref_n_cells)
- mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("part_id"))
- mesh.GetCellData().AddArray(ref_mesh.GetCellData().GetArray("sector_id"))
- myVTK.computeStrainsFromDisplacements(
+ if (mesh_w_local_basis is not None):
+ assert (n_cells == mesh_w_local_basis_n_cells)
+ mesh.GetCellData().AddArray(mesh_w_local_basis.GetCellData().GetArray("part_id"))
+ mesh.GetCellData().AddArray(mesh_w_local_basis.GetCellData().GetArray("sector_id"))
+ myvtk.addDeformationGradients(
mesh=mesh,
disp_array_name=disp_array_name,
- ref_mesh=ref_mesh,
+ defo_grad_array_name=defo_grad_array_name,
verbose=verbose)
- myVTK.writeUGrid(
+ if (ref_frame is not None):
+ farray_F = mesh.GetCellData().GetArray(defo_grad_array_name)
+ for k_cell in xrange(n_cells):
+ F = numpy.reshape(farray_F.GetTuple(k_cell) , (3,3), order='C')
+ F0 = numpy.reshape(farray_F0.GetTuple(k_cell), (3,3), order='C')
+ F = numpy.dot(F, numpy.linalg.inv(F0))
+ farray_F.SetTuple(k_cell, numpy.reshape(F, 9, order='C'))
+ myvtk.addStrainsFromDeformationGradients(
+ mesh=mesh,
+ defo_grad_array_name=defo_grad_array_name,
+ strain_array_name=strain_array_name,
+ mesh_w_local_basis=mesh_w_local_basis,
+ verbose=verbose)
+ myvtk.writeUGrid(
ugrid=mesh,
- filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(n_zfill)+"."+sol_ext,
+ filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext,
verbose=verbose)
if (write_strains) or (write_strain_vs_radius):
- if (ref_mesh is not None):
- assert (mesh.GetCellData().HasArray("Strain_"+CYL_or_PPS))
- farray_strain = mesh.GetCellData().GetArray("Strain_"+CYL_or_PPS)
+ if (mesh_w_local_basis is not None):
+ assert (mesh.GetCellData().HasArray(strain_array_name+"_"+CYL_or_PPS))
+ farray_strain = mesh.GetCellData().GetArray(strain_array_name+"_"+CYL_or_PPS)
else:
- farray_strain = mesh.GetCellData().GetArray("Strain")
+ farray_strain = mesh.GetCellData().GetArray(strain_array_name)
if (write_strains):
if (n_sector_ids == 0):
@@ -122,8 +157,8 @@ def compute_strains(
strain_file.write("\n")
if (write_strain_vs_radius):
- assert (ref_mesh.GetCellData().HasArray("r"))
- farray_r = ref_mesh.GetCellData().GetArray("r")
+ assert (mesh_w_local_basis.GetCellData().HasArray("r"))
+ farray_r = mesh_w_local_basis.GetCellData().GetArray("r")
for k_cell in xrange(n_cells):
strain_vs_radius_file.write(" ".join([str(val) for val in [k_frame, farray_r.GetTuple1(k_cell)]+list(farray_strain.GetTuple(k_cell))]) + "\n")
strain_vs_radius_file.write("\n")
@@ -132,5 +167,12 @@ def compute_strains(
if (write_strains):
strain_file.close()
+ if (plot_strains):
+ ddic.plot_strains(
+ working_folder=sol_folder,
+ working_basenames=[sol_basename],
+ suffix=None,
+ verbose=verbose)
+
if (write_strain_vs_radius):
strain_vs_radius_file.close()
diff --git a/compute_unwarped_images.py b/compute_unwarped_images.py
new file mode 100644
index 0000000..2c39f86
--- /dev/null
+++ b/compute_unwarped_images.py
@@ -0,0 +1,95 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2012-2016 ###
+### ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import glob
+import numpy
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+def compute_unwarped_images(
+ images_folder,
+ images_basename,
+ sol_folder,
+ sol_basename,
+ sol_ext="vtu",
+ verbose=0):
+
+ ref_image_zfill = len(glob.glob(images_folder+"/"+images_basename+"_*.vti")[0].rsplit("_")[-1].split(".")[0])
+ ref_image_filename = images_folder+"/"+images_basename+"_"+str(0).zfill(ref_image_zfill)+".vti"
+ ref_image = myvtk.readImage(
+ filename=ref_image_filename)
+
+ image = vtk.vtkImageData()
+ image.SetOrigin(ref_image.GetOrigin())
+ image.SetSpacing(ref_image.GetSpacing())
+ image.SetExtent(ref_image.GetExtent())
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ image.AllocateScalars(vtk.VTK_FLOAT, 1)
+ else:
+ image.SetScalarTypeToFloat()
+ image.SetNumberOfScalarComponents(1)
+ image.AllocateScalars()
+ scalars = image.GetPointData().GetScalars()
+
+ sol_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
+ n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_"+"[0-9]"*sol_zfill+"."+sol_ext))
+ #n_frames = 1
+
+ X = numpy.empty(3)
+ U = numpy.empty(3)
+ x = numpy.empty(3)
+ I = numpy.empty(1)
+ m = numpy.empty(1)
+ for k_frame in xrange(n_frames):
+ mypy.my_print(verbose, "k_frame = "+str(k_frame))
+
+ def_image = myvtk.readImage(
+ filename=images_folder+"/"+images_basename+"_"+str(k_frame).zfill(ref_image_zfill)+".vti")
+
+ interpolator = myvtk.getImageInterpolator(
+ image=def_image)
+
+ mesh = myvtk.readUGrid(
+ filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext)
+
+ probe = vtk.vtkProbeFilter()
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ probe.SetInputData(image)
+ probe.SetSourceData(mesh)
+ else:
+ probe.SetInput(image)
+ probe.SetSource(mesh)
+ probe.Update()
+ probed_image = probe.GetOutput()
+ scalars_mask = probed_image.GetPointData().GetArray("vtkValidPointMask")
+ scalars_U = probed_image.GetPointData().GetArray("displacement")
+
+ for k_point in xrange(image.GetNumberOfPoints()):
+ scalars_mask.GetTuple(k_point, m)
+ if (m[0] == 0):
+ I[0] = 0.
+ else:
+ image.GetPoint(k_point, X)
+ scalars_U.GetTuple(k_point, U)
+ x = X + U
+ interpolator.Interpolate(x, I)
+ scalars.SetTuple(k_point, I)
+
+ myvtk.writeImage(
+ image=image,
+ filename=sol_folder+"/"+sol_basename+"-unwarped_"+str(k_frame).zfill(sol_zfill)+".vti")
diff --git a/compute_warped_images.py b/compute_warped_images.py
new file mode 100644
index 0000000..b454ab4
--- /dev/null
+++ b/compute_warped_images.py
@@ -0,0 +1,110 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2012-2016 ###
+### ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import glob
+import numpy
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+def compute_warped_images(
+ ref_image_folder,
+ ref_image_basename,
+ sol_folder,
+ sol_basename,
+ ref_frame=0,
+ sol_ext="vtu",
+ verbose=0):
+
+ ref_image_zfill = len(glob.glob(ref_image_folder+"/"+ref_image_basename+"_*.vti")[0].rsplit("_")[-1].split(".")[0])
+ ref_image_filename = ref_image_folder+"/"+ref_image_basename+"_"+str(ref_frame).zfill(ref_image_zfill)+".vti"
+ ref_image = myvtk.readImage(
+ filename=ref_image_filename)
+
+ interpolator = myvtk.getImageInterpolator(
+ image=ref_image)
+ #I = numpy.empty(1)
+ #interpolator.Interpolate([0.35, 0.25, 0.], I)
+
+ image = vtk.vtkImageData()
+ image.SetOrigin(ref_image.GetOrigin())
+ image.SetSpacing(ref_image.GetSpacing())
+ image.SetExtent(ref_image.GetExtent())
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ image.AllocateScalars(vtk.VTK_FLOAT, 1)
+ else:
+ image.SetScalarTypeToFloat()
+ image.SetNumberOfScalarComponents(1)
+ image.AllocateScalars()
+ scalars = image.GetPointData().GetScalars()
+
+ sol_zfill = len(glob.glob(sol_folder+"/"+sol_basename+"_*."+sol_ext)[0].rsplit("_")[-1].split(".")[0])
+ n_frames = len(glob.glob(sol_folder+"/"+sol_basename+"_"+"[0-9]"*sol_zfill+"."+sol_ext))
+ #n_frames = 1
+
+ X = numpy.empty(3)
+ U = numpy.empty(3)
+ x = numpy.empty(3)
+ I = numpy.empty(1)
+ m = numpy.empty(1)
+ for k_frame in xrange(n_frames):
+ mypy.my_print(verbose, "k_frame = "+str(k_frame))
+
+ mesh = myvtk.readUGrid(
+ filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext)
+ #print mesh
+
+ warp = vtk.vtkWarpVector()
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ warp.SetInputData(mesh)
+ else:
+ warp.SetInput(mesh)
+ warp.Update()
+ warped_mesh = warp.GetOutput()
+ #myvtk.writeUGrid(
+ #ugrid=warped_mesh,
+ #filename=sol_folder+"/"+sol_basename+"-warped_"+str(k_frame).zfill(sol_zfill)+"."+sol_ext)
+
+ probe = vtk.vtkProbeFilter()
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ probe.SetInputData(image)
+ probe.SetSourceData(warped_mesh)
+ else:
+ probe.SetInput(image)
+ probe.SetSource(warped_mesh)
+ probe.Update()
+ probed_image = probe.GetOutput()
+ scalars_mask = probed_image.GetPointData().GetArray("vtkValidPointMask")
+ scalars_U = probed_image.GetPointData().GetArray("displacement")
+ #myvtk.writeImage(
+ #image=probed_image,
+ #filename=sol_folder+"/"+sol_basename+"_"+str(k_frame).zfill(sol_zfill)+".vti")
+
+ for k_point in xrange(image.GetNumberOfPoints()):
+ scalars_mask.GetTuple(k_point, m)
+ if (m[0] == 0):
+ I[0] = 0.
+ else:
+ image.GetPoint(k_point, x)
+ scalars_U.GetTuple(k_point, U)
+ X = x - U
+ interpolator.Interpolate(X, I)
+ scalars.SetTuple(k_point, I)
+
+ myvtk.writeImage(
+ image=image,
+ filename=sol_folder+"/"+sol_basename+"-warped_"+str(k_frame).zfill(sol_zfill)+".vti")
diff --git a/compute_warped_mesh.py b/compute_warped_mesh.py
new file mode 100644
index 0000000..b5a803f
--- /dev/null
+++ b/compute_warped_mesh.py
@@ -0,0 +1,81 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2012-2016 ###
+### ###
+### University of California at San Francisco (UCSF), USA ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import numpy
+import os
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+def compute_warped_mesh(
+ mesh_folder,
+ mesh_basename,
+ images,
+ structure,
+ deformation,
+ evolution,
+ verbose=0):
+
+ mypy.my_print(verbose, "*** compute_warped_mesh ***")
+
+ mesh = myvtk.readUGrid(
+ filename=mesh_folder+"/"+mesh_basename+".vtk",
+ verbose=verbose-1)
+ n_points = mesh.GetNumberOfPoints()
+ n_cells = mesh.GetNumberOfCells()
+
+ if os.path.exists(mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk"):
+ ref_mesh = myvtk.readUGrid(
+ filename=mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk",
+ verbose=verbose-1)
+ else:
+ ref_mesh = None
+
+ farray_disp = myvtk.createFloatArray(
+ name="displacement",
+ n_components=3,
+ n_tuples=n_points,
+ verbose=verbose-1)
+ mesh.GetPointData().AddArray(farray_disp)
+
+ mapping = Mapping(images, structure, deformation, evolution)
+
+ X = numpy.empty(3)
+ x = numpy.empty(3)
+ U = numpy.empty(3)
+ if ("zfill" not in images.keys()):
+ images["zfill"] = len(str(images["n_frames"]))
+ for k_frame in xrange(images["n_frames"]):
+ t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
+ mapping.init_t(t)
+
+ for k_point in xrange(n_points):
+ mesh.GetPoint(k_point, X)
+ mapping.x(X, x)
+ U = x - X
+ farray_disp.SetTuple(k_point, U)
+
+ myvtk.addStrainsFromDisplacements(
+ mesh=mesh,
+ disp_array_name="displacement",
+ ref_mesh=ref_mesh,
+ verbose=verbose-1)
+
+ myvtk.writeUGrid(
+ ugrid=mesh,
+ filename=mesh_folder+"/"+mesh_basename+"_"+str(k_frame).zfill(images["zfill"])+".vtk",
+ verbose=verbose-1)
diff --git a/fedic.py b/fedic.py
index 6274373..3653935 100644
--- a/fedic.py
+++ b/fedic.py
@@ -16,10 +16,10 @@ import os
import shutil
import time
-import myVTKPythonLibrary as myVTK
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
-import myFEniCSPythonLibrary as myFEniCS
-from print_tools import *
+import dolfin_dic as ddic
########################################################################
@@ -34,7 +34,6 @@ def fedic(
working_basename,
images_folder,
images_basename,
- images_zfill=2,
images_n_frames=None,
images_ref_frame=0,
images_quadrature=None,
@@ -51,6 +50,7 @@ def fedic(
residual_type="Iref", # Iref, Iold, Iref-then-Iold
relax_type="constant", # constant, aitken, gss
relax_init=1.0,
+ initialize_DU_with_DUold=0,
tol_res=None,
tol_res_rel=None,
tol_dU=None,
@@ -61,17 +61,20 @@ def fedic(
tab = 0
- print_str(tab,"Checking number of frames…")
+ mypy.print_str(tab,"Checking number of frames…")
+ image_filenames = glob.glob(images_folder+"/"+images_basename+"_[0-9]*.vti")
+ images_zfill = len(image_filenames[0].rsplit("_",1)[-1].split(".",1)[0])
+ #mypy.print_var(tab+1,"images_zfill",images_zfill)
if (images_n_frames is None):
- images_n_frames = len(glob.glob(images_folder+"/"+images_basename+"_"+"[0-9]"*images_zfill+".vti"))
+ images_n_frames = len(image_filenames)
assert (images_n_frames > 1), "images_n_frames = "+str(images_n_frames)+" <= 1. Aborting."
- print_var(tab+1,"images_n_frames",images_n_frames)
+ mypy.print_var(tab+1,"images_n_frames",images_n_frames)
assert (abs(images_ref_frame) < images_n_frames), "abs(images_ref_frame) = "+str(images_ref_frame)+" >= images_n_frames. Aborting."
images_ref_frame = images_ref_frame%images_n_frames
- print_var(tab+1,"images_ref_frame",images_ref_frame)
+ mypy.print_var(tab+1,"images_ref_frame",images_ref_frame)
- print_str(tab,"Loading mesh…")
+ mypy.print_str(tab,"Loading mesh…")
assert (mesh is not None or ((mesh_folder is not None) and (mesh_basename is not None))), "Must provide a mesh (mesh = "+str(mesh)+") or a mesh file (mesh_folder = "+str(mesh_folder)+", mesh_basename = "+str(mesh_basename)+"). Aborting."
if (mesh is None):
mesh_filebasename = mesh_folder+"/"+mesh_basename
@@ -80,23 +83,26 @@ def fedic(
mesh = dolfin.Mesh(mesh_filename)
dX = dolfin.dx(mesh)
mesh_V0 = dolfin.assemble(dolfin.Constant(1)*dX)
- print_var(tab+1,"mesh_n_cells",len(mesh.cells()))
- print_sci(tab+1,"mesh_V0",mesh_V0)
+ mypy.print_var(tab+1,"mesh_n_cells",len(mesh.cells()))
+ mypy.print_sci(tab+1,"mesh_V0",mesh_V0)
- print_str(tab,"Computing quadrature degree for images…")
+ mypy.print_str(tab,"Computing quadrature degree for images…")
ref_image_filename = images_folder+"/"+images_basename+"_"+str(images_ref_frame).zfill(images_zfill)+".vti"
if (images_quadrature is None):
- images_quadrature = myFEniCS.compute_quadrature_degree_from_points_count(
+ images_quadrature = ddic.compute_quadrature_degree_from_points_count(
image_filename=ref_image_filename,
mesh_filebasename=mesh_filebasename,
verbose=1)
- print_var(tab+1,"images_quadrature",images_quadrature)
-
- print_str(tab,"Loading reference image…")
- images_dimension = myVTK.computeImageDimensionality(
- image_filename=ref_image_filename,
+ mypy.print_var(tab+1,"images_quadrature",images_quadrature)
+
+ mypy.print_str(tab,"Loading reference image…")
+ ref_image = myvtk.readImage(
+ filename=ref_image_filename,
+ verbose=verbose)
+ images_dimension = myvtk.getImageDimensionality(
+ image=ref_image,
verbose=0)
- print_var(tab+1,"images_dimension",images_dimension)
+ mypy.print_var(tab+1,"images_dimension",images_dimension)
assert (images_dimension in (2,3)), "images_dimension must be 2 or 3. Aborting."
fe = dolfin.FiniteElement(
family="Quadrature",
@@ -113,19 +119,19 @@ def fedic(
cell=mesh.ufl_cell(),
degree=images_quadrature,
quad_scheme="default")
- te._quad_scheme = "default" # shouldn't be needed
- for k in xrange(images_dimension**2): # shouldn't be needed
- te.sub_elements()[k]._quad_scheme = "default" # shouldn't be needed
+ te._quad_scheme = "default" # should not be needed
+ for k in xrange(images_dimension**2): # should not be needed
+ te.sub_elements()[k]._quad_scheme = "default" # should not be needed
if (images_expressions_type == "cpp"):
Iref = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="im",
im_is_def=0),
element=fe)
Iref.init_image(ref_image_filename)
DIref = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="grad",
im_is_def=0),
@@ -133,17 +139,17 @@ def fedic(
DIref.init_image(ref_image_filename)
elif (images_expressions_type == "py"):
if (images_dimension == 2):
- Iref = myFEniCS.ExprIm2(
+ Iref = ddic.ExprIm2(
filename=ref_image_filename,
element=fe)
- DIref = myFEniCS.ExprGradIm2(
+ DIref = ddic.ExprGradIm2(
filename=ref_image_filename,
element=ve)
elif (images_dimension == 3):
- Iref = myFEniCS.ExprIm3(
+ Iref = ddic.ExprIm3(
filename=ref_image_filename,
element=fe)
- DIref = myFEniCS.ExprGradIm3(
+ DIref = ddic.ExprGradIm3(
filename=ref_image_filename,
element=ve)
else:
@@ -151,10 +157,15 @@ def fedic(
Iref_int = dolfin.assemble(Iref * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V0
Iref_norm = (dolfin.assemble(Iref**2 * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V0)**(1./2)
assert (Iref_norm > 0.), "Iref_norm = "+str(Iref_norm)+" <= 0. Aborting."
- print_var(tab+1,"Iref_int",Iref_int)
- print_var(tab+1,"Iref_norm",Iref_norm)
+ mypy.print_var(tab+1,"Iref_int",Iref_int)
+ mypy.print_var(tab+1,"Iref_norm",Iref_norm)
+
+ file_error_basename = working_folder+"/"+working_basename+"-error"
+ file_error = open(file_error_basename+".dat", "w")
+ file_error.write("#k_frame err_im"+"\n")
+ file_error.write(" ".join([str(val) for val in [images_ref_frame, 0.]])+"\n")
- print_str(tab,"Defining functions…")
+ mypy.print_str(tab,"Defining functions…")
vfs = dolfin.VectorFunctionSpace(
mesh=mesh,
family="Lagrange",
@@ -169,28 +180,31 @@ def fedic(
name="previous displacement")
Uold.vector().zero()
Uold_norm = 0.
+ DUold = dolfin.Function(
+ vfs,
+ name="previous displacement increment")
dU = dolfin.Function(
vfs,
name="displacement correction")
dU_ = dolfin.TrialFunction(vfs)
dV_ = dolfin.TestFunction(vfs)
- print_str(tab,"Printing initial solution…")
+ mypy.print_str(tab,"Printing initial solution…")
if not os.path.exists(working_folder):
os.mkdir(working_folder)
- pvd_basename = working_folder+"/"+working_basename+"_"
- for vtu_filename in glob.glob(pvd_basename+"*.vtu"):
+ pvd_basename = working_folder+"/"+working_basename
+ for vtu_filename in glob.glob(pvd_basename+"_[0-9]*.vtu"):
os.remove(vtu_filename)
- file_pvd = dolfin.File(pvd_basename+".pvd")
+ file_pvd = dolfin.File(pvd_basename+"__.pvd")
file_pvd << (U, float(images_ref_frame))
- os.remove(pvd_basename+".pvd")
- shutil.move(pvd_basename+"".zfill(6)+".vtu", pvd_basename+str(images_ref_frame).zfill(6)+".vtu")
+ os.remove(pvd_basename+"__.pvd")
+ shutil.move(pvd_basename+"__"+"".zfill(6)+".vtu", pvd_basename+"_"+str(images_ref_frame).zfill(6)+".vtu")
if (print_iterations):
- for filename in glob.glob(working_folder+"/"+working_basename+"-frame=*.*"):
+ for filename in glob.glob(working_folder+"/"+working_basename+"-frame=[0-9]*.*"):
os.remove(filename)
- print_str(tab,"Defining mechanical energy…")
+ mypy.print_str(tab,"Defining mechanical energy…")
E = dolfin.Constant(1.0)
nu = dolfin.Constant(regul_poisson)
kappa = E/3/(1-2*nu) # = E/3 if nu = 0
@@ -229,100 +243,99 @@ def fedic(
DDpsi_m = dolfin.derivative(Dpsi_m, U, dU_)
mesh_V = dolfin.assemble(J*dX)
- print_sci(tab+1,"mesh_V",mesh_V)
+ mypy.print_sci(tab+1,"mesh_V",mesh_V)
- volume_basename = working_folder+"/"+working_basename+"-volume"
- file_volume = open(volume_basename+".dat","w")
+ file_volume_basename = working_folder+"/"+working_basename+"-volume"
+ file_volume = open(file_volume_basename+".dat","w")
file_volume.write("#k_frame mesh_V"+"\n")
file_volume.write(" ".join([str(val) for val in [images_ref_frame, mesh_V]])+"\n")
- print_str(tab,"Defining deformed image…")
+ mypy.print_str(tab,"Defining deformed image…")
scaling = numpy.array([1.,0.])
- #scaling = dolfin.Constant([1.,0.])
if (images_expressions_type == "cpp"):
Idef = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="im",
im_is_def=1),
element=fe)
- Idef.init_disp(U)
Idef.init_dynamic_scaling(scaling)
+ Idef.init_disp(U)
DIdef = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="grad",
im_is_def=1),
element=ve)
+ Idef.init_dynamic_scaling(scaling)
DIdef.init_disp(U)
- DIdef.init_dynamic_scaling(scaling)
if ("-wHess" in tangent_type):
assert (0), "ToDo"
Iold = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="im",
im_is_def=1),
element=fe)
+ Iold.init_dynamic_scaling(scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
Iold.init_disp(Uold)
- Iold.init_dynamic_scaling(scaling)
DIold = dolfin.Expression(
- cppcode=myFEniCS.get_ExprIm_cpp(
+ cppcode=ddic.get_ExprIm_cpp(
im_dim=images_dimension,
im_type="grad",
im_is_def=1),
element=ve)
+ DIold.init_dynamic_scaling(scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
DIold.init_disp(Uold)
- DIold.init_dynamic_scaling(scaling)
elif (images_expressions_type == "py"):
if (images_dimension == 2):
- Idef = myFEniCS.ExprDefIm2(
+ Idef = ddic.ExprDefIm2(
U=U,
scaling=scaling,
element=fe)
- DIdef = myFEniCS.ExprGradDefIm2(
+ DIdef = ddic.ExprGradDefIm2(
U=U,
scaling=scaling,
element=ve)
if ("-wHess" in tangent_type):
- DDIdef = myFEniCS.ExprHessDefIm2(
+ DDIdef = ddic.ExprHessDefIm2(
U=U,
scaling=scaling,
element=te)
- Iold = myFEniCS.ExprDefIm2(
+ Iold = ddic.ExprDefIm2(
U=Uold,
- scaling=scaling,
+ scaling=scaling, # 2016/07/25: ok, same scaling must apply to Idef & Iold…
element=fe)
- DIold = myFEniCS.ExprGradDefIm2(
+ DIold = ddic.ExprGradDefIm2(
U=Uold,
- scaling=scaling,
+ scaling=scaling, # 2016/07/25: ok, same scaling must apply to Idef & Iold…
element=ve)
elif (images_dimension == 3):
- Idef = myFEniCS.ExprDefIm3(
+ Idef = ddic.ExprDefIm3(
U=U,
scaling=scaling,
element=fe)
- DIdef = myFEniCS.ExprGradDefIm3(
+ DIdef = ddic.ExprGradDefIm3(
U=U,
scaling=scaling,
element=ve)
if ("-wHess" in tangent_type):
- DDIdef = myFEniCS.ExprHessDefIm3(
+ DDIdef = ddic.ExprHessDefIm3(
U=U,
- scaling=scaling,
+ scaling=scaling, # 2016/07/25: ok, same scaling must apply to Idef & Iold…
element=te)
- Iold = myFEniCS.ExprDefIm3(
+ Iold = ddic.ExprDefIm3(
U=Uold,
scaling=scaling,
element=fe)
- DIold = myFEniCS.ExprGradDefIm3(
+ DIold = ddic.ExprGradDefIm3(
U=Uold,
- scaling=scaling,
+ scaling=scaling, # 2016/07/25: ok, same scaling must apply to Idef & Iold…
element=ve)
else:
assert (0), "\"images_expressions_type\" (="+str(images_expressions_type)+") must be \"cpp\" or \"py\". Aborting."
- print_str(tab,"Defining correlation energy…")
+ mypy.print_str(tab,"Defining correlation energy…")
psi_c = (Idef - Iref)**2/2
Dpsi_c = (Idef - Iref) * dolfin.dot(DIdef, dV_)
DDpsi_c = dolfin.dot(DIdef, dU_) * dolfin.dot(DIdef, dV_)
@@ -338,23 +351,24 @@ def fedic(
B0 = dolfin.assemble(b0, form_compiler_parameters={'quadrature_degree':images_quadrature})
res_norm0 = B0.norm("l2")
assert (res_norm0 > 0.), "res_norm0 = "+str(res_norm0)+" <= 0. Aborting."
- print_var(tab+1,"res_norm0",res_norm0)
+ mypy.print_var(tab+1,"res_norm0",res_norm0)
#regul_level = dolfin.Constant(regul_level) * 5*(1+nu)*(1-2*nu)/(5-4*nu)
regul_level = dolfin.Constant(regul_level)
A = None
if (tangent_type == "Iref"):
- print_str(tab,"Assembly…")
+ mypy.print_str(tab,"Matrix assembly… (image term)")
A = dolfin.assemble((1.-regul_level) * DDpsi_c_ref * dX, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
- A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, add_values=True)
+ mypy.print_str(tab,"Matrix assembly… (regularization term)")
+ A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, form_compiler_parameters={'quadrature_degree':mesh_degree}, add_values=True)
B = None
- print_str(tab,"Looping over frames…")
+ mypy.print_str(tab,"Looping over frames…")
n_iter_tot = 0
global_success = True
for forward_or_backward in ["forward", "backward"]:
- print_var(tab,"forward_or_backward",forward_or_backward)
+ mypy.print_var(tab,"forward_or_backward",forward_or_backward)
if (forward_or_backward == "forward"):
k_frames_old = range(images_ref_frame , images_n_frames-1, +1)
@@ -362,18 +376,20 @@ def fedic(
elif (forward_or_backward == "backward"):
k_frames_old = range(images_ref_frame , 0, -1)
k_frames = range(images_ref_frame-1, -1, -1)
- print_var(tab,"k_frames",k_frames)
+ mypy.print_var(tab,"k_frames",k_frames)
if (forward_or_backward == "backward"):
U.vector().zero()
U_norm = 0.
Uold.vector().zero()
Uold_norm = 0.
+ DUold.vector().zero()
+ scaling[:] = [1.,0.]
tab += 1
success = True
for (k_frame,k_frame_old) in zip(k_frames,k_frames_old):
- print_var(tab-1,"k_frame",k_frame)
+ mypy.print_var(tab-1,"k_frame",k_frame)
if (print_iterations):
frame_basename = working_folder+"/"+working_basename+"-frame="+str(k_frame).zfill(images_zfill)
@@ -383,7 +399,7 @@ def fedic(
file_pvd_frame = dolfin.File(frame_basename+"_.pvd")
file_pvd_frame << (U, 0.)
- print_str(tab,"Loading image, image gradient and image hessian…")
+ mypy.print_str(tab,"Loading image, image gradient and image hessian…")
image_filename = images_folder+"/"+images_basename+"_"+str(k_frame).zfill(images_zfill)+".vti"
Idef.init_image(image_filename)
DIdef.init_image(image_filename)
@@ -395,12 +411,13 @@ def fedic(
# linear system: matrix
if (tangent_type == "Iold"):
- print_str(tab,"Assembly…")
+ mypy.print_str(tab,"Matrix assembly… (image term)")
A = dolfin.assemble((1.-regul_level) * DDpsi_c_old * dX, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
- A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, add_values=True)
- #print_var(tab,"A",A.array())
+ mypy.print_str(tab,"Matrix assembly… (regularization term)")
+ A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, form_compiler_parameters={'quadrature_degree':mesh_degree}, add_values=True)
+ #mypy.print_var(tab,"A",A.array())
#A_norm = A.norm("l2")
- #print_var(tab,"A_norm",A_norm)
+ #mypy.print_var(tab,"A_norm",A_norm)
if (print_iterations):
U.vector().zero()
@@ -412,7 +429,10 @@ def fedic(
err_im = im_diff/Iref_norm
file_dat_frame.write(" ".join([str(val) for val in [-1, None, None, None, None, None, None, im_diff, err_im, None]])+"\n")
- print_str(tab,"Running registration…")
+ if (initialize_DU_with_DUold):
+ U.vector().axpy(1., DUold.vector())
+
+ mypy.print_str(tab,"Running registration…")
tab += 1
k_iter = 0
if (residual_type.startswith("Iref")):
@@ -420,17 +440,18 @@ def fedic(
elif (residual_type.startswith("Iold")):
using_Iold_residual = True
while (True):
- print_var(tab-1,"k_iter",k_iter)
+ mypy.print_var(tab-1,"k_iter",k_iter)
n_iter_tot += 1
# linear system: matrix assembly
if (tangent_type.startswith("Idef")):
- print_str(tab,"Matrix assembly…")
+ mypy.print_str(tab,"Matrix assembly… (image term)")
A = dolfin.assemble((1.-regul_level) * DDpsi_c * dX, tensor=A, form_compiler_parameters={'quadrature_degree':images_quadrature})
- A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, add_values=True)
- #print_var(tab,"A",A.array())
+ mypy.print_str(tab,"Matrix assembly… (regularization term)")
+ A = dolfin.assemble(( regul_level) * DDpsi_m * dX, tensor=A, form_compiler_parameters={'quadrature_degree':mesh_degree}, add_values=True)
+ #mypy.print_var(tab,"A",A.array())
#A_norm = A.norm("l2")
- #print_sci(tab,"A_norm",A_norm)
+ #mypy.print_sci(tab,"A_norm",A_norm)
# linear system: residual assembly
if (k_iter > 0):
@@ -439,20 +460,23 @@ def fedic(
elif (k_iter > 1):
B_old[:] = B[:]
res_old_norm = res_norm
- print_str(tab,"Residual assembly…")
if (using_Iold_residual):
+ mypy.print_str(tab,"Residual assembly… (image term)")
B = dolfin.assemble(- (1.-regul_level) * Dpsi_c_old * dX, tensor=B, form_compiler_parameters={'quadrature_degree':images_quadrature})
- B = dolfin.assemble(- ( regul_level) * Dpsi_m * dX, tensor=B, add_values=True)
+ mypy.print_str(tab,"Residual assembly… (regularization term)")
+ B = dolfin.assemble(- ( regul_level) * Dpsi_m * dX, tensor=B, form_compiler_parameters={'quadrature_degree':mesh_degree}, add_values=True)
else:
+ mypy.print_str(tab,"Residual assembly… (image term)")
B = dolfin.assemble(- (1.-regul_level) * Dpsi_c * dX, tensor=B, form_compiler_parameters={'quadrature_degree':images_quadrature})
- B = dolfin.assemble(- ( regul_level) * Dpsi_m * dX, tensor=B, add_values=True)
- #print_var(tab,"B",B.array())
+ mypy.print_str(tab,"Residual assembly… (regularization term)")
+ B = dolfin.assemble(- ( regul_level) * Dpsi_m * dX, tensor=B, form_compiler_parameters={'quadrature_degree':mesh_degree}, add_values=True)
+ #mypy.print_var(tab,"B",B.array())
# residual error
res_norm = B.norm("l2")
- #print_sci(tab,"res_norm",res_norm)
+ #mypy.print_sci(tab,"res_norm",res_norm)
err_res = res_norm/res_norm0
- print_sci(tab,"err_res",err_res)
+ mypy.print_sci(tab,"err_res",err_res)
if (k_iter == 0):
err_res_rel = 1.
@@ -463,12 +487,12 @@ def fedic(
dB[:] = B[:] - B_old[:]
dres_norm = dB.norm("l2")
err_res_rel = dres_norm / res_old_norm
- print_sci(tab,"err_res_rel",err_res_rel)
+ mypy.print_sci(tab,"err_res_rel",err_res_rel)
# linear system: solve
- print_str(tab,"Solve…")
+ mypy.print_str(tab,"Solve…")
dolfin.solve(A, dU.vector(), B)
- #print_var(tab,"dU",dU.vector().array())
+ #mypy.print_var(tab,"dU",dU.vector().array())
# relaxation
if (relax_type == "constant"):
@@ -479,7 +503,7 @@ def fedic(
relax = relax_init
else:
relax *= (-1.) * B_old.inner(dB) / dres_norm**2
- print_sci(tab,"relax",relax)
+ mypy.print_sci(tab,"relax",relax)
elif (relax_type == "gss"):
phi = (1 + math.sqrt(5)) / 2
relax_a = (1-phi)/(2-phi)
@@ -492,43 +516,43 @@ def fedic(
tab += 1
relax_k = 0
while (True):
- print_var(tab-1,"relax_k",relax_k)
- print_sci(tab,"relax_a",relax_a)
- print_sci(tab,"relax_b",relax_b)
+ mypy.print_var(tab-1,"relax_k",relax_k)
+ mypy.print_sci(tab,"relax_a",relax_a)
+ mypy.print_sci(tab,"relax_b",relax_b)
if (need_update_c):
relax_c = relax_b - (relax_b - relax_a) / phi
relax_list.append(relax_c)
- print_sci(tab,"relax_c",relax_c)
+ mypy.print_sci(tab,"relax_c",relax_c)
U.vector().axpy(relax_c-relax_cur, dU.vector())
relax_cur = relax_c
relax_fc = dolfin.assemble((1.-regul_level) * psi_c * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})
- print_sci(tab,"relax_fc",relax_fc)
- relax_fc += dolfin.assemble(( regul_level) * psi_m * dX)
- print_sci(tab,"relax_fc",relax_fc)
- #print_sci(tab,"relax_fc",relax_fc)
+ #mypy.print_sci(tab,"relax_fc",relax_fc)
+ relax_fc += dolfin.assemble(( regul_level) * psi_m * dX, form_compiler_parameters={'quadrature_degree':mesh_degree})
+ #mypy.print_sci(tab,"relax_fc",relax_fc)
if (numpy.isnan(relax_fc)):
relax_fc = float('+inf')
- print_sci(tab,"relax_fc",relax_fc)
+ #mypy.print_sci(tab,"relax_fc",relax_fc)
+ mypy.print_sci(tab,"relax_fc",relax_fc)
relax_vals.append(relax_fc)
- #print_var(tab,"relax_list",relax_list)
- #print_var(tab,"relax_vals",relax_vals)
+ #mypy.print_var(tab,"relax_list",relax_list)
+ #mypy.print_var(tab,"relax_vals",relax_vals)
if (need_update_d):
relax_d = relax_a + (relax_b - relax_a) / phi
relax_list.append(relax_d)
- print_sci(tab,"relax_d",relax_d)
+ mypy.print_sci(tab,"relax_d",relax_d)
U.vector().axpy(relax_d-relax_cur, dU.vector())
relax_cur = relax_d
relax_fd = dolfin.assemble((1.-regul_level) * psi_c * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})
- print_sci(tab,"relax_fd",relax_fd)
- relax_fd += dolfin.assemble(( regul_level) * psi_m * dX)
- print_sci(tab,"relax_fd",relax_fd)
- #print_sci(tab,"relax_fd",relax_fd)
+ #mypy.print_sci(tab,"relax_fd",relax_fd)
+ relax_fd += dolfin.assemble(( regul_level) * psi_m * dX, form_compiler_parameters={'quadrature_degree':mesh_degree})
+ #mypy.print_sci(tab,"relax_fd",relax_fd)
if (numpy.isnan(relax_fd)):
relax_fd = float('+inf')
- print_sci(tab,"relax_fd",relax_fd)
+ #mypy.print_sci(tab,"relax_fd",relax_fd)
+ mypy.print_sci(tab,"relax_fd",relax_fd)
relax_vals.append(relax_fd)
- #print_var(tab,"relax_list",relax_list)
- #print_var(tab,"relax_vals",relax_vals)
+ #mypy.print_var(tab,"relax_list",relax_list)
+ #mypy.print_var(tab,"relax_vals",relax_vals)
if (relax_fc < relax_fd):
relax_b = relax_d
relax_d = relax_c
@@ -547,7 +571,7 @@ def fedic(
relax_k += 1
tab -= 1
U.vector().axpy(-relax_cur, dU.vector())
- #print_var(tab,"relax_vals",relax_vals)
+ #mypy.print_var(tab,"relax_vals",relax_vals)
if (print_iterations):
iter_basename = frame_basename+"-iter="+str(k_iter).zfill(3)
@@ -557,7 +581,7 @@ def fedic(
os.system("gnuplot -e \"set terminal pdf; set output '"+iter_basename+".pdf'; plot '"+iter_basename+".dat' u 1:2 w p title 'psi_int'\"")
relax = relax_list[numpy.argmin(relax_vals)]
- print_sci(tab,"relax",relax)
+ mypy.print_sci(tab,"relax",relax)
else:
assert (0), "relax_type must be \"constant\", \"aitken\" or \"gss\". Aborting."
@@ -566,7 +590,7 @@ def fedic(
U_norm = U.vector().norm("l2")
if (print_iterations):
- #print_var(tab,"U",U.vector().array())
+ #mypy.print_var(tab,"U",U.vector().array())
file_pvd_frame << (U, float(k_iter+1))
# displacement error
@@ -577,15 +601,15 @@ def fedic(
err_dU = dU_norm/U_norm
else:
err_dU = dU_norm/Uold_norm
- print_sci(tab,"err_dU",err_dU)
+ mypy.print_sci(tab,"err_dU",err_dU)
# image error
if (k_iter > 0):
im_diff_old = im_diff
im_diff = (dolfin.assemble(psi_c * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V0)**(1./2)
- #print_sci(tab,"im_diff",im_diff)
+ #mypy.print_sci(tab,"im_diff",im_diff)
err_im = im_diff/Iref_norm
- print_sci(tab,"err_im",err_im)
+ mypy.print_sci(tab,"err_im",err_im)
if (print_iterations):
file_dat_frame.write(" ".join([str(val) for val in [k_iter, res_norm, err_res, err_res_rel, relax, dU_norm, U_norm, err_dU, im_diff, err_im]])+"\n")
@@ -603,19 +627,19 @@ def fedic(
# exit
if (success):
- print_str(tab,"Nonlinear solver converged…")
+ mypy.print_str(tab,"Nonlinear solver converged…")
break
if (k_iter == n_iter_max-1):
if (residual_type=="Iref-then-Iold") and not (using_Iold_residual):
- print_str(tab,"Warning! Nonlinear solver failed to converge…using Iold instead of Iref. (k_frame = "+str(k_frame)+")")
+ mypy.print_str(tab,"Warning! Nonlinear solver failed to converge…using Iold instead of Iref. (k_frame = "+str(k_frame)+")")
using_Iold_residual = True
U.vector()[:] = Uold.vector()[:]
U_norm = Uold_norm
k_iter = 0
continue
else:
- print_str(tab,"Warning! Nonlinear solver failed to converge… (k_frame = "+str(k_frame)+")")
+ mypy.print_str(tab,"Warning! Nonlinear solver failed to converge… (k_frame = "+str(k_frame)+")")
global_success = False
break
@@ -633,18 +657,22 @@ def fedic(
break
# solution update
+ DUold.vector()[:] = U.vector()[:] - Uold.vector()[:]
Uold.vector()[:] = U.vector()[:]
Uold_norm = U_norm
mesh_V = dolfin.assemble(J*dX)
- print_sci(tab+1,"mesh_V",mesh_V)
+ mypy.print_sci(tab+1,"mesh_V",mesh_V)
file_volume.write(" ".join([str(val) for val in [k_frame, mesh_V]])+"\n")
- print_str(tab,"Printing solution…")
- file_pvd = dolfin.File(pvd_basename+".pvd")
+ mypy.print_str(tab,"Printing solution…")
+ file_pvd = dolfin.File(pvd_basename+"__.pvd")
file_pvd << (U, float(k_frame))
- os.remove(pvd_basename+".pvd")
- shutil.move(pvd_basename+"".zfill(6)+".vtu", pvd_basename+str(k_frame).zfill(6)+".vtu")
+ os.remove(
+ pvd_basename+"__.pvd")
+ shutil.move(
+ pvd_basename+"__"+"".zfill(6)+".vtu",
+ pvd_basename+"_"+str(k_frame).zfill(6)+".vtu")
if (images_dynamic_scaling):
p = numpy.empty((2,2))
@@ -656,20 +684,33 @@ def fedic(
q[0] = dolfin.assemble(Idef*Iref * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})
q[1] = dolfin.assemble(Iref * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})
scaling[:] = numpy.linalg.solve(p,q)
- print_var(tab,"scaling", scaling)
+ mypy.print_var(tab,"scaling", scaling)
+
+ if (images_expressions_type == "cpp"): # should not be needed
+ Idef.update_dynamic_scaling(scaling) # should not be needed
+ DIdef.update_dynamic_scaling(scaling) # should not be needed
+ if ("-wHess" in tangent_type): # should not be needed
+ DDIdef.update_dynamic_scaling(scaling) # should not be needed
+ Iold.update_dynamic_scaling(scaling) # should not be needed
+ DIold.update_dynamic_scaling(scaling) # should not be needed
im_diff = (dolfin.assemble(psi_c * dX, form_compiler_parameters={'quadrature_degree':images_quadrature})/mesh_V0)**(1./2)
err_im = im_diff/Iref_norm
- print_sci(tab,"err_im",err_im)
+ mypy.print_sci(tab,"err_im",err_im)
+
+ file_error.write(" ".join([str(val) for val in [k_frame, err_im]])+"\n")
tab -= 1
if not (success) and not (continue_after_fail):
break
- print_var(tab,"n_iter_tot",n_iter_tot)
+ mypy.print_var(tab,"n_iter_tot",n_iter_tot)
+
+ file_error.close()
+ os.system("gnuplot -e \"set terminal pdf; set output '"+file_error_basename+".pdf'; set grid; set yrange [0:1]; plot '"+file_error_basename+".dat' u 1:2 lw 3 notitle\"")
file_volume.close()
- #os.remove(pvd_basename+".pvd")
+ os.system("gnuplot -e \"set terminal pdf; set output '"+file_volume_basename+".pdf'; set grid; set yrange [0:*]; plot '"+file_volume_basename+".dat' u 1:2 lw 3 notitle\"")
return global_success
diff --git a/generate_images.py b/generate_images.py
new file mode 100644
index 0000000..2255094
--- /dev/null
+++ b/generate_images.py
@@ -0,0 +1,665 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2016 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import glob
+import math
+import numpy
+import os
+import random
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+#class ImagesInfo():
+ #def __init__(self, n_dim, L, n_voxels, n_integration, T, n_frames, data_type, images_folder, images_basename):
+ #assert (n_dim in (1,2,3))
+ #self.n_dim = n_dim
+
+ #if (type(L) == float):
+ #assert (L>0)
+ #self.L = numpy.array([L]*self.n_dim)
+ #elif (type(L) == int):
+ #assert (L>0)
+ #self.L = numpy.array([float(L)]*self.n_dim)
+ #else:
+ #assert (len(L) == self.n_dim)
+ #self.L = numpy.array(L)
+ #assert ((self.L>0).all())
+
+ #if (type(n_voxels) == int):
+ #assert (n_voxels>0)
+ #self.n_voxels = numpy.array([n_voxels]*self.n_dim)
+ #else:
+ #assert (len(n_voxels) == self.n_dim)
+ #self.n_voxels = numpy.array(n_voxels)
+ #assert ((self.n_voxels>0).all())
+
+ #if (type(n_integration) == int):
+ #assert (n_integration>0)
+ #self.n_integration = numpy.array([n_integration]*self.n_dim)
+ #else:
+ #assert (len(n_integration) == self.n_dim)
+ #self.n_integration = numpy.array(n_integration)
+ #assert ((self.n_integration>0).all())
+
+ #assert (T>0.)
+ #self.T = T
+
+ #assert (n_frames>0)
+ #self.n_frames = n_frames
+
+ #assert (data_type in ("int", "float", "unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float" "uint8", "uint16", "uint32", "uint64", "ufloat"))
+ #self.data_type = data_type
+
+ #self.images_folder = images_folder
+ #self.images_basename = images_basename
+
+#class StructureInfo():
+ #def __init__(self, images, type, **kwargs):
+ #assert (type in ("no", "heart"))
+ #self["type"] = type
+ #if (self["type"] == "heart"):
+ #self.Ri = kwargs["Ri"]
+ #self.Re = kwargs["Re"]
+ #if (images.n_dim == 3):
+ #self.Zmin = kwargs["Zmin"] if ("Zmin" in kwargs.keys()) else 0.
+ #self.Zmax = kwargs["Zmax"] if ("Zmax" in kwargs.keys()) else images.L[2]
+
+#class TextureInfo():
+ #def __init__(self, type, **kwargs):
+ #assert (type in ("no", "taggX", "taggY", "taggZ"))
+ #self["type"] = type
+
+#class NoiseInfo():
+ #def __init__(self, type, **kwargs):
+ #self["type"] = type
+
+#class DeformationInfo():
+ #def __init__(self, type, **kwargs):
+ #self["type"] = type
+
+#class EvolutionInfo():
+ #def __init__(self, type, **kwargs):
+ #self["type"] = type
+
+########################################################################
+
+class Image():
+ def __init__(self, images, structure, texture, noise):
+ self.L = images["L"]
+
+ # structure
+ if (structure["type"] == "no"):
+ self.I0_structure = self.I0_structure_no
+ elif (structure["type"] == "heart"):
+ if (images["n_dim"] == 2):
+ self.I0_structure = self.I0_structure_heart_2
+ self.R = float()
+ self.Ri = structure["Ri"]
+ self.Re = structure["Re"]
+ elif (images["n_dim"] == 3):
+ self.I0_structure = self.I0_structure_heart_3
+ self.R = float()
+ self.Ri = structure["Ri"]
+ self.Re = structure["Re"]
+ self.Zmin = structure.Zmin if ("Zmin" in structure.keys()) else 0.
+ self.Zmax = structure.Zmax if ("Zmax" in structure.keys()) else images["L"][2]
+ else:
+ assert (0), "n_dim must be \"2\" or \"3 for \"heart\" type structure. Aborting."
+ else:
+ assert (0), "structure type must be \"no\" or \"heart\". Aborting."
+
+ # texture
+ if (texture["type"] == "no"):
+ self.I0_texture = 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
+ else:
+ self.I0_texture = self.I0_texture_tagging_X
+ elif (images["n_dim"] == 2):
+ if ("-signed" in texture["type"]):
+ self.I0_texture = self.I0_texture_tagging_signed_XY
+ else:
+ self.I0_texture = self.I0_texture_tagging_XY
+ elif (images["n_dim"] == 3):
+ if ("-signed" in texture["type"]):
+ self.I0_texture = self.I0_texture_tagging_signed_XYZ
+ else:
+ self.I0_texture = self.I0_texture_tagging_XYZ
+ else:
+ assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
+ self.s = texture["s"]
+ elif (texture["type"].startswith("taggX")):
+ self.I0_texture = self.I0_texture_tagging_X
+ self.s = texture["s"]
+ elif (texture["type"].startswith("taggY")):
+ self.I0_texture = self.I0_texture_tagging_Y
+ self.s = texture["s"]
+ elif (texture["type"].startswith("taggZ")):
+ self.I0_texture = self.I0_texture_tagging_Z
+ self.s = texture["s"]
+ else:
+ assert (0), "texture type must be \"no\", \"tagging\", \"taggX\", \"taggY\" or \"taggZ\". Aborting."
+
+ # noise
+ if (noise["type"] == "no"):
+ self.I0_noise = self.I0_noise_no
+ elif (noise["type"] == "normal"):
+ self.I0_noise = self.I0_noise_normal
+ self.avg = noise.avg if ("avg" in noise.keys()) else 0.
+ self.std = noise.std
+ else:
+ assert (0), "noise type must be \"no\" or \"normal\". Aborting."
+
+ def I0(self, X, i, g=None):
+ self.I0_structure(X, i, g)
+ self.I0_texture(X, i, g)
+ self.I0_noise(i, g)
+
+ def I0_structure_no(self, X, i, g=None):
+ i[0] = 1.
+ #if (g is not None): g[:] = 1.
+
+ def I0_structure_heart_2(self, X, i, g=None):
+ self.R = ((X[0]-self.L[0]/2)**2 + (X[1]-self.L[1]/2)**2)**(1./2)
+ if (self.R >= self.Ri) and (self.R <= self.Re):
+ i[0] = 1.
+ #if (g is not None): g[:] = 1.
+ else:
+ i[0] = 0.
+ #if (g is not None): g[:] = 0.
+ #i[0] = 1.
+
+ def I0_structure_heart_3(self, X, i, g=None):
+ self.R = ((X[0]-self.L[0]/2)**2 + (X[1]-self.L[1]/2)**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.
+ #if (g is not None): g[:] = 1.
+ else:
+ i[0] = 0.
+ #if (g is not None): g[:] = 0.
+
+ def I0_texture_no(self, X, i, g=None):
+ i[0] *= 1.
+ #if (g is not None): g[:] *= 0.
+
+ def I0_texture_tagging_signed_X(self, X, i, g=None):
+ i[0] *= math.sin(math.pi*X[0]/self.s)
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_X(self, X, i, g=None):
+ i[0] *= abs(math.sin(math.pi*X[0]/self.s))
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_signed_Y(self, X, i, g=None):
+ i[0] *= math.sin(math.pi*X[1]/self.s)
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_Y(self, X, i, g=None):
+ i[0] *= abs(math.sin(math.pi*X[1]/self.s))
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_signed_Z(self, X, i, g=None):
+ i[0] *= math.sin(math.pi*X[2]/self.s)
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_Z(self, X, i, g=None):
+ i[0] *= abs(math.sin(math.pi*X[2]/self.s))
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_signed_XY(self, X, i, g=None):
+ i[0] *= (math.sin(math.pi*X[0]/self.s)
+ + math.sin(math.pi*X[1]/self.s))/2
+ assert (i[0] >= 0.)
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_XY(self, X, i, g=None):
+ i[0] *= (abs(math.sin(math.pi*X[0]/self.s))
+ + abs(math.sin(math.pi*X[1]/self.s)))/2
+ assert (i[0] >= 0.)
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_signed_XYZ(self, X, i, g=None):
+ i[0] *= (math.sin(math.pi*X[0]/self.s)
+ + math.sin(math.pi*X[1]/self.s)
+ + math.sin(math.pi*X[2]/self.s))/3
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_texture_tagging_XYZ(self, X, i, g=None):
+ 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
+ #if (g is not None):
+ #assert (0), "ToDo"
+
+ def I0_noise_no(self, i, g=None):
+ pass
+
+ def I0_noise_normal(self, i, g=None):
+ i[0] += random.normalvariate(self.avg, self.std)
+ #if (g is not None): g[k] += [2*random.normalvariate(self.avg, self.std) for k in xrange(len(g))]
+
+########################################################################
+
+class Mapping:
+ def __init__(self, images, structure, deformation, evolution):
+ 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"] == "trans"):
+ self.init_t = self.init_t_trans
+ self.X = self.X_trans
+ self.x = self.x_trans
+ self.D = numpy.empty(3)
+ elif (self.deformation["type"] == "rot"):
+ self.init_t = self.init_t_rot
+ self.X = self.X_rot
+ self.x = self.x_rot
+ 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
+ elif (self.deformation["type"] == "heart"):
+ assert (structure["type"] == "heart"), "structure type must be \"heart\" for \"heart\" type deformation. 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\", \"trans\", \"rot\", \"homogeneous\" or \"heart\". 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 ("+evolution["type"]+") type must be \"linear\" or \"sine\". 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_trans(self, t):
+ self.D[0] = self.deformation["Dx"]*self.phi(t) if ("Dx" in self.deformation.keys()) else 0.
+ self.D[1] = self.deformation["Dy"]*self.phi(t) if ("Dy" in self.deformation.keys()) else 0.
+ self.D[2] = self.deformation["Dz"]*self.phi(t) if ("Dz" in self.deformation.keys()) else 0.
+
+ def init_t_rot(self, t):
+ self.C[0] = self.deformation["Cx"] if ("Cx" in self.deformation.keys()) else 0.
+ self.C[1] = self.deformation["Cy"] if ("Cy" in self.deformation.keys()) else 0.
+ self.C[2] = self.deformation["Cz"] if ("Cz" in self.deformation.keys()) else 0.
+ Rx = self.deformation["Rx"]*math.pi/180*self.phi(t) if ("Rx" in self.deformation.keys()) else 0.
+ Ry = self.deformation["Ry"]*math.pi/180*self.phi(t) if ("Ry" in self.deformation.keys()) else 0.
+ Rz = self.deformation["Rz"]*math.pi/180*self.phi(t) if ("Rz" in self.deformation.keys()) 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):
+ Exx = self.deformation["Exx"]*self.phi(t) if ("Exx" in self.deformation.keys()) else 0.
+ Eyy = self.deformation["Eyy"]*self.phi(t) if ("Eyy" in self.deformation.keys()) else 0.
+ Ezz = self.deformation["Ezz"]*self.phi(t) if ("Ezz" in self.deformation.keys()) else 0.
+ Exy = self.deformation["Exy"]*self.phi(t) if ("Exy" in self.deformation.keys()) else 0.
+ Eyx = self.deformation["Eyx"]*self.phi(t) if ("Eyx" in self.deformation.keys()) else 0.
+ Exz = self.deformation["Exz"]*self.phi(t) if ("Exz" in self.deformation.keys()) else 0.
+ Ezx = self.deformation["Ezx"]*self.phi(t) if ("Ezx" in self.deformation.keys()) else 0.
+ Eyz = self.deformation["Eyz"]*self.phi(t) if ("Eyz" in self.deformation.keys()) else 0.
+ Ezy = self.deformation["Ezy"]*self.phi(t) if ("Ezy" in self.deformation.keys()) else 0.
+ self.F = numpy.array([[math.sqrt(1.+Exx), Exy , Exz ],
+ [ Eyx , math.sqrt(1.+Eyy), Eyz ],
+ [ Ezx , Ezy , math.sqrt(1.+Ezz)]])
+ 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.keys()) else 0.
+ self.dRe = self.deformation["dRi"]*self.phi(t) if ("dRi" in self.deformation.keys()) else 0.
+ self.dTi = self.deformation["dTi"]*self.phi(t) if ("dTi" in self.deformation.keys()) else 0.
+ self.dTe = self.deformation["dTe"]*self.phi(t) if ("dTe" in self.deformation.keys()) 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_trans(self, x, X, Finv=None):
+ X[:] = x - self.D
+ #if (Finv is not None): Finv[:,:] = numpy.identity(numpy.sqrt(numpy.size(Finv)))
+
+ def X_rot(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)
+ #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_trans(self, X, x, F=None):
+ x[:] = X + self.D
+ #if (F is not None): F[:,:] = numpy.identity(numpy.sqrt(numpy.size(F)))
+
+ def x_rot(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)
+ #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)
+
+########################################################################
+
+def generateImages(
+ images,
+ structure,
+ texture,
+ noise,
+ deformation,
+ evolution,
+ generate_image_gradient=False,
+ verbose=0):
+
+ mypy.my_print(verbose, "*** generateImages ***")
+
+ vtk_image = vtk.vtkImageData()
+
+ if (images["n_dim"] == 1):
+ vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0, 0, 0, 0])
+ elif (images["n_dim"] == 2):
+ vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, 0])
+ elif (images["n_dim"] == 3):
+ vtk_image.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, images["n_voxels"][2]-1])
+ else:
+ assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
+
+ spacing = numpy.array(images["L"])/numpy.array(images["n_voxels"])
+ if (images["n_dim"] == 1):
+ spacing = numpy.array([images["L"][0]/images["n_voxels"][0], 1., 1.])
+ elif (images["n_dim"] == 2):
+ spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], 1.])
+ elif (images["n_dim"] == 2):
+ spacing = numpy.array([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], images["L"][2]/images["n_voxels"][2]])
+ vtk_image.SetSpacing(spacing)
+
+ origin = numpy.array(vtk_image.GetSpacing())/2
+ if (images["n_dim"] == 1):
+ origin[1] = 0.
+ origin[2] = 0.
+ elif (images["n_dim"] == 2):
+ origin[2] = 0.
+ vtk_image.SetOrigin(origin)
+
+ n_points = vtk_image.GetNumberOfPoints()
+ vtk_image_scalars = myvtk.createFloatArray(
+ name="ImageScalars",
+ n_components=1,
+ n_tuples=n_points,
+ verbose=verbose-1)
+ vtk_image.GetPointData().SetScalars(vtk_image_scalars)
+ if (generate_image_gradient):
+ vtk_image_gradient = myvtk.createFloatArray(
+ name="ImageScalarsGradient",
+ n_components=images["n_dim"],
+ n_tuples=n_points,
+ verbose=verbose-1)
+ vtk_image.GetPointData().SetVectors(vtk_image_gradient)
+
+ if not os.path.exists(images["folder"]):
+ os.mkdir(images["folder"])
+
+ x0 = numpy.empty(3)
+ x = numpy.empty(3)
+ X = numpy.empty(3)
+ if (generate_image_gradient):
+ F = numpy.empty((3,3))
+ Finv = numpy.empty((3,3))
+ else:
+ F = None
+ Finv = None
+ dx = spacing[0:images["n_dim"]]/images["n_integration"][0:images["n_dim"]]
+ global_min = float("+Inf")
+ global_max = float("-Inf")
+ I = numpy.empty(1)
+ i = numpy.empty(1)
+ if (generate_image_gradient):
+ G = numpy.empty(images["n_dim"])
+ g = numpy.empty(images["n_dim"])
+ else:
+ G = None
+ g = None
+ image = Image(images, structure, texture, noise)
+ mapping = Mapping(images, structure, deformation, evolution)
+ if ("zfill" not in images.keys()):
+ images["zfill"] = len(str(images["n_frames"]))
+ for k_frame in xrange(images["n_frames"]):
+ t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
+ print "t = "+str(t)
+ mapping.init_t(t)
+ for k_point in xrange(n_points):
+ vtk_image.GetPoint(k_point, x0)
+ #print "x0 = "+str(x0)
+ x[:] = x0[:]
+ #print "x = "+str(x)
+ I[0] = 0.
+ #print "I = "+str(I)
+ #if (generate_image_gradient): G[:] = 0.
+ #print "G = "+str(G)
+ if (images["n_dim"] == 1):
+ for k_x in xrange(images["n_integration"][0]):
+ x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
+ mapping.X(x, X, Finv)
+ image.I0(X, i, g)
+ I += i
+ #if (generate_image_gradient): G += numpy.dot(g, Finv)
+ I /= images["n_integration"][0]
+ #if (generate_image_gradient): G /= images["n_integration"][0]
+ elif (images["n_dim"] == 2):
+ for k_y in xrange(images["n_integration"][1]):
+ x[1] = x0[1] - dx[1]/2 + (k_y+1./2)*dx[1]/images["n_integration"][1]
+ for k_x in xrange(images["n_integration"][0]):
+ x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
+ #print "x = "+str(x)
+ mapping.X(x, X, Finv)
+ #print "X = "+str(X)
+ #print "Finv = "+str(Finv)
+ image.I0(X, i, g)
+ #print "i = "+str(i)
+ #print "g = "+str(g)
+ I += i
+ #if (generate_image_gradient): G += numpy.dot(g, Finv)
+ I /= images["n_integration"][1]*images["n_integration"][0]
+ #if (generate_image_gradient):G /= images["n_integration"][1]*images["n_integration"][0]
+ elif (images["n_dim"] == 3):
+ for k_z in xrange(images["n_integration"][2]):
+ x[2] = x0[2] - dx[2]/2 + (k_z+1./2)*dx[2]/images["n_integration"][2]
+ for k_y in xrange(images["n_integration"][1]):
+ x[1] = x0[1] - dx[1]/2 + (k_y+1./2)*dx[1]/images["n_integration"][1]
+ for k_x in xrange(images["n_integration"][0]):
+ x[0] = x0[0] - dx[0]/2 + (k_x+1./2)*dx[0]/images["n_integration"][0]
+ mapping.X(x, X, Finv)
+ image.I0(X, i, g)
+ I += i
+ #if (generate_image_gradient): G += numpy.dot(g, Finv)
+ I /= images["n_integration"][2]*images["n_integration"][1]*images["n_integration"][0]
+ #if (generate_image_gradient): G /= images["n_integration"][2]*images["n_integration"][1]*images["n_integration"][0]
+ else:
+ assert (0), "n_dim must be \"1\", \"2\" or \"3\". Aborting."
+ vtk_image_scalars.SetTuple(k_point, I)
+ #if (generate_image_gradient): vtk_image_gradient.SetTuple(k_point, G)
+ if (I[0] < global_min): global_min = I[0]
+ if (I[0] > global_max): global_max = I[0]
+ myvtk.writeImage(
+ image=vtk_image,
+ filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+ verbose=verbose-1)
+
+ if (images["data_type"] in ("float")):
+ pass
+ elif (images["data_type"] in ("unsigned char", "unsigned short", "unsigned int", "unsigned long", "unsigned float", "uint8", "uint16", "uint32", "uint64", "ufloat")):
+ #print "global_min = "+str(global_min)
+ #print "global_max = "+str(global_max)
+ shifter = vtk.vtkImageShiftScale()
+ shifter.SetShift(-global_min)
+ if (images["data_type"] in ("unsigned char", "uint8")):
+ shifter.SetScale(float(2**8-1)/(global_max-global_min))
+ shifter.SetOutputScalarTypeToUnsignedChar()
+ elif (images["data_type"] in ("unsigned short", "uint16")):
+ shifter.SetScale(float(2**16-1)/(global_max-global_min))
+ shifter.SetOutputScalarTypeToUnsignedShort()
+ elif (images["data_type"] in ("unsigned int", "uint32")):
+ shifter.SetScale(float(2**32-1)/(global_max-global_min))
+ shifter.SetOutputScalarTypeToUnsignedInt()
+ elif (images["data_type"] in ("unsigned long", "uint64")):
+ shifter.SetScale(float(2**64-1)/(global_max-global_min))
+ shifter.SetOutputScalarTypeToUnsignedLong()
+ elif (images["data_type"] in ("unsigned float", "ufloat")):
+ shifter.SetScale(1./(global_max-global_min))
+ shifter.SetOutputScalarTypeToFloat()
+ for k_frame in xrange(images["n_frames"]):
+ vtk_image = myvtk.readImage(
+ filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+ verbose=verbose-1)
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ shifter.SetInputData(vtk_image)
+ else:
+ shifter.SetInput(vtk_image)
+ shifter.Update()
+ vtk_image = shifter.GetOutput()
+ myvtk.writeImage(
+ image=vtk_image,
+ filename=images["folder"]+"/"+images["basename"]+"_"+str(k_frame).zfill(images["zfill"])+".vti",
+ verbose=verbose-1)
+ else:
+ assert (0), "Wrong data type. Aborting."
diff --git a/generate_undersampled_images.py b/generate_undersampled_images.py
new file mode 100644
index 0000000..9ae62e5
--- /dev/null
+++ b/generate_undersampled_images.py
@@ -0,0 +1,164 @@
+#coding=utf8
+
+########################################################################
+### ###
+### Created by Martin Genet, 2016 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+########################################################################
+
+import numpy
+import vtk
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+########################################################################
+
+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"
+ fedic.generateImages(
+ 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"
+ fedic.generateImages(
+ 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"
+ fedic.generateImages(
+ 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 xrange(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 = vtk.vtkImageData()
+ imageXY.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, 0])
+ imageXY.SetSpacing([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], 1.])
+ imageXY.SetOrigin([images["L"][0]/images["n_voxels"][0]/2, images["L"][1]/images["n_voxels"][1]/2, 0.])
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ imageXY.AllocateScalars(vtk.VTK_FLOAT, 1)
+ else:
+ imageXY.SetScalarTypeToFloat()
+ imageXY.SetNumberOfScalarComponents(1)
+ imageXY.AllocateScalars()
+ scalars = imageXY.GetPointData().GetScalars()
+ x = numpy.empty(3)
+ for k_point in xrange(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 = vtk.vtkImageData()
+ imageXYZ.SetExtent([0, images["n_voxels"][0]-1, 0, images["n_voxels"][1]-1, 0, images["n_voxels"][2]-1])
+ imageXYZ.SetSpacing([images["L"][0]/images["n_voxels"][0], images["L"][1]/images["n_voxels"][1], images["L"][2]/images["n_voxels"][2]])
+ imageXYZ.SetOrigin([images["L"][0]/images["n_voxels"][0]/2, images["L"][1]/images["n_voxels"][1]/2, images["L"][2]/images["n_voxels"][2]/2])
+ if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
+ imageXYZ.AllocateScalars(vtk.VTK_FLOAT, 1)
+ else:
+ imageXYZ.SetScalarTypeToFloat()
+ imageXYZ.SetNumberOfScalarComponents(1)
+ imageXYZ.AllocateScalars()
+ scalars = imageXYZ.GetPointData().GetScalars()
+ x = numpy.empty(3)
+ for k_point in xrange(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")
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 1131bf4..424be96 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -46,7 +46,9 @@ class MyExpr : public Expression
{
vtkSmartPointer<vtkImageInterpolator> interpolator;
double static_scaling;'''+('''
- Array<double>* dynamic_scaling;
+ Array<double>* dynamic_scaling; // does not work
+ double dynamic_scaling_a; // should not be needed
+ double dynamic_scaling_b; // should not be needed
Function* U;
mutable Array<double> UX;
mutable Array<double> x;''')*(im_is_def)+('''
@@ -56,6 +58,8 @@ public:
MyExpr():
Expression('''+str(im_dim)*(im_type=="grad")+''')'''+(''',
+ dynamic_scaling_a(1.),
+ dynamic_scaling_b(0.),
UX('''+str(im_dim)+'''),
x(3)''')*(im_is_def)+(''',
X3D(3)''')*(not im_is_def)*(im_dim==2)+'''
@@ -65,10 +69,19 @@ public:
void init_dynamic_scaling(
const Array<double> &scaling)
{
- //dynamic_scaling = scaling;
- //std::cout << "dynamic_scaling = " << dynamic_scaling->str(1) << std::endl;
+ //dynamic_scaling = scaling; // does not work
+ // std::cout << "dynamic_scaling = " << dynamic_scaling->str(1) << std::endl; // does not work
+ dynamic_scaling_a = scaling[0]; // should not be needed
+ dynamic_scaling_b = scaling[1]; // should not be needed
}
+ void update_dynamic_scaling( // should not be needed
+ const Array<double> &scaling) // should not be needed
+ { // should not be needed
+ dynamic_scaling_a = scaling[0]; // should not be needed
+ dynamic_scaling_b = scaling[1]; // should not be needed
+ } // should not be needed
+
void init_disp(
Function* UU)
{
@@ -146,35 +159,47 @@ public:
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+('''
- expr[0] = expr[0]/static_scaling;''')*(im_type=="im")+('''
+ expr[0] /= static_scaling;''')*(im_type=="im")+('''
- expr[0] = expr[0]/static_scaling;
- expr[1] = expr[1]/static_scaling;''')*(im_type=="grad")*(im_dim==2)+('''
+ expr[0] /= static_scaling;
+ expr[1] /= static_scaling;''')*(im_type=="grad")*(im_dim==2)+('''
- expr[0] = expr[0]/static_scaling;
- expr[1] = expr[1]/static_scaling;
- expr[2] = expr[2]/static_scaling;''')*(im_type=="grad")*(im_dim==3)+('''
+ expr[0] /= static_scaling;
+ expr[1] /= static_scaling;
+ expr[2] /= static_scaling;''')*(im_type=="grad")*(im_dim==3)+('''
std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose)+(('''
- std::cout << "in (im)" << std::endl;
- std::cout << "expr = " << expr.str(1) << std::endl;
- std::cout << "dynamic_scaling = " << dynamic_scaling->str(1) << std::endl;
- expr[0] = expr[0] * (*dynamic_scaling)[0] + (*dynamic_scaling)[1];
- std::cout << "out (im)" << std::endl;''')*(im_type=="im")+('''
-
- expr[0] = expr[0] * (*dynamic_scaling)[0];
- expr[1] = expr[1] * (*dynamic_scaling)[0];''')*(im_type=="grad")*(im_dim==2)+('''
-
- std::cout << "in (grad)" << std::endl;
- std::cout << "expr = " << expr.str(1) << std::endl;
- std::cout << "dynamic_scaling = " << dynamic_scaling->str(1) << std::endl;
- expr[0] = expr[0] * (*dynamic_scaling)[0];
- expr[1] = expr[1] * (*dynamic_scaling)[0];
- expr[2] = expr[2] * (*dynamic_scaling)[0];
- std::cout << "out (grad)" << std::endl;''')*(im_type=="grad")*(im_dim==3)+('''
-
- std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose))*(im_is_def)*0+'''
+ // std::cout << "in (im)" << std::endl;
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // expr[0] *= (*dynamic_scaling)[0]; // does not work
+ // expr[0] += (*dynamic_scaling)[1]; // does not work
+ expr[0] *= dynamic_scaling_a; // should not be needed
+ expr[0] += dynamic_scaling_b; // should not be needed
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // std::cout << "out (im)" << std::endl;''')*(im_type=="im")+('''
+
+ // std::cout << "in (grad)" << std::endl;
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // expr[0] *= (*dynamic_scaling)[0]; // does not work
+ // expr[1] *= (*dynamic_scaling)[0]; // does not work
+ expr[0] *= dynamic_scaling_a; // should not be needed
+ expr[1] *= dynamic_scaling_a; // should not be needed
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // std::cout << "out (grad)" << std::endl;''')*(im_type=="grad")*(im_dim==2)+('''
+
+ // std::cout << "in (grad)" << std::endl;
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // expr[0] *= (*dynamic_scaling)[0]; // does not work
+ // expr[1] *= (*dynamic_scaling)[0]; // does not work
+ // expr[2] *= (*dynamic_scaling)[0]; // does not work
+ expr[0] *= dynamic_scaling_a; // should not be needed
+ expr[1] *= dynamic_scaling_a; // should not be needed
+ expr[2] *= dynamic_scaling_a; // should not be needed
+ // std::cout << "expr = " << expr.str(1) << std::endl;
+ // std::cout << "out (grad)" << std::endl;''')*(im_type=="grad")*(im_dim==3)+('''
+
+ std::cout << "expr = " << expr.str(1) << std::endl;''')*(verbose))*(im_is_def)+'''
}
};
diff --git a/image_expressions_py.py b/image_expressions_py.py
index f6cde4c..5724b8b 100644
--- a/image_expressions_py.py
+++ b/image_expressions_py.py
@@ -11,8 +11,10 @@
import dolfin
import numpy
-import myVTKPythonLibrary as myVTK
-from myVTKPythonLibrary.mat_vec_tools import *
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
########################################################################
@@ -32,12 +34,12 @@ class ExprIm2(dolfin.Expression):
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
@@ -58,12 +60,12 @@ class ExprIm3(dolfin.Expression):
self.init_image(filename=filename)
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
@@ -83,15 +85,15 @@ class ExprGradIm2(dolfin.Expression):
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageGradient(
+ self.image = myvtk.addImageGradient(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -114,15 +116,15 @@ class ExprGradIm3(dolfin.Expression):
self.init_image(filename=filename)
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageGradient(
+ self.image = myvtk.addImageGradient(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -144,15 +146,15 @@ class ExprHessIm2(dolfin.Expression):
self.X = numpy.array([float()]*2+[Z])
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageHessian(
+ self.image = myvtk.addImageHessian(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -172,15 +174,15 @@ class ExprHessIm3(dolfin.Expression):
self.init_image(filename=filename)
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageHessian(
+ self.image = myvtk.addImageHessian(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -203,12 +205,12 @@ class ExprDefIm2(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
@@ -238,12 +240,12 @@ class ExprDefIm3(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="linear",
out_value=0.,
@@ -273,15 +275,15 @@ class ExprGradDefIm2(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageGradient(
+ self.image = myvtk.addImageGradient(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -313,15 +315,15 @@ class ExprGradDefIm3(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageGradient(
+ self.image = myvtk.addImageGradient(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -353,15 +355,15 @@ class ExprHessDefIm2(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageHessian(
+ self.image = myvtk.addImageHessian(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
@@ -387,15 +389,15 @@ class ExprHessDefIm3(dolfin.Expression):
self.scaling = scaling
def init_image(self, filename):
- self.image = myVTK.readImage(
+ self.image = myvtk.readImage(
filename=filename,
verbose=0)
self.s = getScalingFactor(
scalar_type_as_string=self.image.GetScalarTypeAsString())
- self.image = myVTK.computeImageHessian(
+ self.image = myvtk.addImageHessian(
image=self.image,
verbose=0)
- self.interpolator = myVTK.createImageInterpolator(
+ self.interpolator = myvtk.getImageInterpolator(
image=self.image,
mode="nearest",
out_value=self.s,
diff --git a/plot_strains.py b/plot_strains.py
index f0d3521..eb1d8d5 100644
--- a/plot_strains.py
+++ b/plot_strains.py
@@ -15,13 +15,18 @@ import os
def plot_strains(
working_folder,
working_basenames,
+ working_scalings=None,
ref_folder=None,
ref_basename=None,
- components="all",
+ components="all", # all, circ-long, or rad-circ
suffix="",
verbose=1):
- assert (components in ("all", "circ-long", "rad-circ"))
+
+ if (working_scalings is None):
+ working_scalings = [1.] * len(working_basenames)
+ else:
+ assert (len(working_scalings) == len(working_basenames))
if (ref_folder is not None) and (ref_basename is not None):
lines = open(ref_folder+"/"+ref_basename+"-strains.dat").readlines()
@@ -32,11 +37,17 @@ def plot_strains(
#print "n_frames = " + str(n_frames)
#print "n_sectors = " + str(n_sectors)
- plotfile = open("plot_strains"+("-"+suffix)*(suffix!="")+".plt", "w")
+ assert (components in ("all", "circ-long", "rad-circ"))
+
+ 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")
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
@@ -92,13 +103,14 @@ plot 0 linecolor rgb "black" notitle,\\
'''))
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 ($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+'''-strains.dat" using ('''+str(working_scaling)+'''*($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 ('''+str(working_scaling)+'''*($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)*(''',\\
''')+(k_basename==len(working_basenames)-1)*('''
'''))
plotfile.close()
- os.system("gnuplot plot_strains"+("-"+suffix)*(suffix!="")+".plt")
+ os.system("gnuplot "+plotfile_basename+".plt")
diff --git a/print_tools.py b/print_tools.py
deleted file mode 100644
index 608624b..0000000
--- a/print_tools.py
+++ /dev/null
@@ -1,42 +0,0 @@
-#coding=utf8
-
-########################################################################
-### ###
-### Created by Martin Genet, 2016 ###
-### ###
-### École Polytechnique, Palaiseau, France ###
-### ###
-########################################################################
-
-########################################################################
-
-#class Printer:
- #def __init__(self):
- #self.cur_level = 0
-
- #def inc(self):
- #self.cur_level += 1
-
- #def dec(self):
- #self.cur_level -= 1
-
- #def print_str(self, string, var_level=0):
- #print " | "*(self.cur_level+var_level) + string
-
- #def print_var(self, name, val, var_level=0):
- #print " | "*(self.cur_level+var_level) + name + " = " + str(val)
-
- #def print_sci(self, name, val, var_level=0):
- #print " | "*(self.cur_level+var_level) + name.ljust(13) + " = " + format(val,".4e")
-
-########################################################################
-
-def print_str(tab, string):
- print " | "*tab + string
-
-def print_var(tab, name, val):
- print " | "*tab + name + " = " + str(val)
-
-def print_sci(tab, name, val):
- print " | "*tab + name.ljust(13) + " = " + format(val,".4e")
-
--
GitLab