Merge branch 'generate_image_using_imaging_model' into 'master'

Generate image using imaging model

See merge request mgenet/dolfin_dic!9

compute_unwarped_images.py

@@ -2,10 +2,10 @@
 
 ################################################################################
 ###                                                                          ###
-### Created by Martin Genet, 2012-2016                               ###
+### Created by Martin Genet, 2012-2016                                       ###
 ###                                                                          ###
-### University of California at San Francisco (UCSF), USA            ###
-### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### University of California at San Francisco (UCSF), USA                    ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland         ###
 ### École Polytechnique, Palaiseau, France                                   ###
 ###                                                                          ###
 ################################################################################

compute_warped_images.py

@@ -2,10 +2,10 @@
 
 ################################################################################
 ###                                                                          ###
-### Created by Martin Genet, 2012-2016                               ###
+### Created by Martin Genet, 2012-2016                                       ###
 ###                                                                          ###
-### University of California at San Francisco (UCSF), USA            ###
-### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### University of California at San Francisco (UCSF), USA                    ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland         ###
 ### École Polytechnique, Palaiseau, France                                   ###
 ###                                                                          ###
 ################################################################################
@@ -27,7 +27,10 @@ def compute_warped_images(
         working_folder,
         working_basename,
         ref_frame=0,
+        ref_image_model=None,
         working_ext="vtu",
+        working_displacement_field_name="displacement",
+        print_warped_mesh=0,
         verbose=0):
 
     ref_image_zfill = len(glob.glob(ref_image_folder+"/"+ref_image_basename+"_*.vti")[0].rsplit("_")[-1].split(".")[0])
@@ -35,10 +38,11 @@ def compute_warped_images(
     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)
+    if (ref_image_model is None):
+        ref_image_interpolator = myvtk.getImageInterpolator(
+            image=ref_image)
+        #I = numpy.empty(1)
+        #ref_image_interpolator.Interpolate([0.35, 0.25, 0.], I)
 
     image = vtk.vtkImageData()
     image.SetOrigin(ref_image.GetOrigin())
@@ -54,7 +58,7 @@ def compute_warped_images(
 
     working_zfill = len(glob.glob(working_folder+"/"+working_basename+"_*."+working_ext)[0].rsplit("_")[-1].split(".")[0])
     n_frames = len(glob.glob(working_folder+"/"+working_basename+"_"+"[0-9]"*working_zfill+"."+working_ext))
-    #n_frames = 1
+    # n_frames = 1
 
     X = numpy.empty(3)
     U = numpy.empty(3)
@@ -66,7 +70,9 @@ def compute_warped_images(
 
         mesh = myvtk.readUGrid(
             filename=working_folder+"/"+working_basename+"_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
-        #print mesh
+        # print mesh
+        assert (mesh.GetPointData().HasArray(working_displacement_field_name)), "no array '" + working_displacement_field_name + "' in mesh"
+        mesh.GetPointData().SetActiveVectors(working_displacement_field_name)
 
         warp = vtk.vtkWarpVector()
         if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
@@ -75,9 +81,10 @@ def compute_warped_images(
             warp.SetInput(mesh)
         warp.Update()
         warped_mesh = warp.GetOutput()
-        #myvtk.writeUGrid(
-            #ugrid=warped_mesh,
-            #filename=working_folder+"/"+working_basename+"-warped_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
+        if print_warped_mesh:
+            myvtk.writeUGrid(
+                ugrid=warped_mesh,
+                filename=working_folder+"/"+working_basename+"-warped_"+str(k_frame).zfill(working_zfill)+"."+working_ext)
 
         probe = vtk.vtkProbeFilter()
         if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
@@ -89,7 +96,7 @@ def compute_warped_images(
         probe.Update()
         probed_image = probe.GetOutput()
         scalars_mask = probed_image.GetPointData().GetArray("vtkValidPointMask")
-        scalars_U = probed_image.GetPointData().GetArray("displacement")
+        scalars_U = probed_image.GetPointData().GetArray(working_displacement_field_name)
         #myvtk.writeImage(
             #image=probed_image,
             #filename=working_folder+"/"+working_basename+"_"+str(k_frame).zfill(working_zfill)+".vti")
@@ -102,7 +109,10 @@ def compute_warped_images(
                 image.GetPoint(k_point, x)
                 scalars_U.GetTuple(k_point, U)
                 X = x - U
-                interpolator.Interpolate(X, I)
+                if (ref_image_model is None):
+                    ref_image_interpolator.Interpolate(X, I)
+                else:
+                    I[0] = ref_image_model(X)
             scalars.SetTuple(k_point, I)
 
         myvtk.writeImage(