diff --git a/Energy.py b/Energy.py
index a645ef42bb0639043bb6ac83b2253ad86360a0c6..064cca612fd9e3ac23e70fb2139937635c6e1ada 100644
--- a/Energy.py
+++ b/Energy.py
@@ -29,6 +29,14 @@ class Energy():
+ def call_before_assembly(self,
+ *args,
+ **kwargs):
+
+ pass
+
+
+
def call_before_solve(self,
*args,
**kwargs):
diff --git a/Energy_GeneratedImage.py b/Energy_GeneratedImage.py
index 8757b91291f3569ba936bf521b92c4f7c6f2f5a1..d22e75058906ad08963f1eaac8cc7e426f42f2c4 100644
--- a/Energy_GeneratedImage.py
+++ b/Energy_GeneratedImage.py
@@ -29,8 +29,7 @@ class WarpedImageEnergy(Energy):
quadrature_degree,
name="im",
w=1.,
- ref_frame=0,
- dynamic_scaling=False):
+ ref_frame=0):
self.problem = problem
self.printer = self.problem.printer
@@ -39,7 +38,6 @@ class WarpedImageEnergy(Energy):
self.name = name
self.w = w
self.ref_frame = ref_frame
- self.dynamic_scaling = dynamic_scaling
self.printer.print_str("Defining warped image correlation energy…")
self.printer.inc()
@@ -78,6 +76,7 @@ class WarpedImageEnergy(Energy):
self.printer.print_str("Defining measure…")
+ # dV
self.form_compiler_parameters = {
"quadrature_degree":self.quadrature_degree,
"quadrature_scheme":"default"}
@@ -86,50 +85,33 @@ class WarpedImageEnergy(Energy):
domain=self.problem.mesh,
metadata=self.form_compiler_parameters)
- self.printer.print_str("Loading reference image…")
- self.printer.inc()
-
# ref_frame
assert (abs(self.ref_frame) < self.image_series.n_frames),\
"abs(ref_frame) = "+str(abs(self.ref_frame))+" >= "+str(self.image_series.n_frames)+" = image_series.n_frames. Aborting."
self.ref_frame = self.ref_frame%self.image_series.n_frames
+ self.ref_image_filename = self.image_series.get_image_filename(self.ref_frame)
self.printer.print_var("ref_frame",self.ref_frame)
- # Iref
- self.Iref = dolfin.Expression(
- cppcode=ddic.get_ExprIm_cpp(
+ self.printer.dec()
+ self.printer.print_str("Defining deformed image…")
+
+ # Igen
+ self.Igen = dolfin.Expression(
+ cppcode=ddic.get_ExprGenIm_cpp(
im_dim=self.image_series.dimension,
- im_type="im",
- im_is_def=0),
+ im_type="im"),
element=self.fe)
- self.ref_image_filename = self.image_series.get_image_filename(self.ref_frame)
- self.Iref.init_image(self.ref_image_filename)
-
- self.Iref_int = dolfin.assemble(self.Iref * self.dV)/self.problem.mesh_V0
- self.printer.print_var("Iref_int",self.Iref_int)
+ self.Igen.init_image(self.ref_image_filename)
+ self.Igen.init_disp(self.problem.U)
- self.Iref_norm = (dolfin.assemble(self.Iref**2 * self.dV)/self.problem.mesh_V0)**(1./2)
- assert (self.Iref_norm > 0.),\
- "Iref_norm = "+str(self.Iref_norm)+" <= 0. Aborting."
- self.printer.print_var("Iref_norm",self.Iref_norm)
-
- # DIref
- self.DIref = dolfin.Expression(
- cppcode=ddic.get_ExprIm_cpp(
+ # DIgen
+ self.DIgen = dolfin.Expression(
+ cppcode=ddic.get_ExprGenIm_cpp(
im_dim=self.image_series.dimension,
- im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
- im_is_def=0),
+ im_type="grad"),
element=self.ve)
- self.ref_image_grad_filename = self.image_series.get_image_grad_filename(self.ref_frame)
- self.DIref.init_image(self.ref_image_grad_filename)
-
- self.printer.dec()
- self.printer.print_str("Defining deformed image…")
-
- self.scaling = numpy.array([1.,0.])
- if (self.dynamic_scaling):
- self.p = numpy.empty((2,2))
- self.q = numpy.empty(2)
+ self.Igen.init_image(self.ref_image_filename)
+ self.Igen.init_disp(self.problem.U)
# Idef
self.Idef = dolfin.Expression(
@@ -140,7 +122,6 @@ class WarpedImageEnergy(Energy):
element=self.fe)
self.Idef.init_image(self.ref_image_filename)
self.Idef.init_disp(self.problem.U)
- self.Idef.init_dynamic_scaling(self.scaling)
# DIdef
self.DIdef = dolfin.Expression(
@@ -151,32 +132,9 @@ class WarpedImageEnergy(Energy):
element=self.ve)
self.DIdef.init_image(self.ref_image_filename)
self.DIdef.init_disp(self.problem.U)
- self.DIdef.init_dynamic_scaling(self.scaling)
self.printer.print_str("Defining previous image…")
- # Iold
- self.Iold = dolfin.Expression(
- cppcode=ddic.get_ExprIm_cpp(
- im_dim=self.image_series.dimension,
- im_type="im",
- im_is_def=1),
- element=self.fe)
- self.Iold.init_image(self.ref_image_filename)
- self.Iold.init_disp(self.problem.Uold)
- self.Iold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
-
- # DIold
- self.DIold = dolfin.Expression(
- cppcode=ddic.get_ExprIm_cpp(
- im_dim=self.image_series.dimension,
- im_type="grad" if (self.image_series.grad_basename is None) else "grad_no_deriv",
- im_is_def=1),
- element=self.ve)
- self.DIold.init_image(self.ref_image_filename)
- self.DIold.init_disp(self.problem.Uold)
- self.DIold.init_dynamic_scaling(self.scaling) # 2016/07/25: ok, same scaling must apply to Idef & Iold…
-
self.printer.print_str("Defining correlation energy…")
# Phi_ref
@@ -187,15 +145,9 @@ class WarpedImageEnergy(Energy):
self.Phi_Iref.init_image(self.ref_image_filename)
# Psi_c
- self.Psi_c = self.Phi_Iref * (self.Idef - self.Iref)**2/2
- self.DPsi_c = self.Phi_Iref * (self.Idef - self.Iref) * dolfin.dot(self.DIdef, self.problem.dU_test)
-
- self.DDPsi_c = self.Phi_Iref * dolfin.dot(self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIdef, self.problem.dU_test)
- self.DDPsi_c_old = self.Phi_Iref * dolfin.dot(self.DIold, self.problem.dU_trial) * dolfin.dot(self.DIold, self.problem.dU_test)
- self.DDPsi_c_ref = self.Phi_Iref * dolfin.dot(self.DIref, self.problem.dU_trial) * dolfin.dot(self.DIref, self.problem.dU_test)
-
- self.Psi_c_old = self.Phi_Iref * (self.Idef - self.Iold)**2/2
- self.DPsi_c_old = self.Phi_Iref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
+ self.Psi_c = self.Phi_Iref * (self.Igen - self.Idef)**2/2
+ self.DPsi_c = self.Phi_Iref * (self.Igen - self.Idef) * dolfin.dot(self.DIgen - self.DIdef, self.problem.dU_test)
+ self.DDPsi_c = self.Phi_Iref * dolfin.dot(self.DIgen - self.DIdef, self.problem.dU_trial) * dolfin.dot(self.DIgen - self.DIdef, self.problem.dU_test)
# forms
self.ener_form = self.Psi_c * self.dV
@@ -208,7 +160,14 @@ class WarpedImageEnergy(Energy):
def reinit(self):
- self.scaling[:] = [1.,0.]
+ pass
+
+
+
+ def call_before_assembly(self):
+
+ self.Igen.generate_image()
+ self.DIgen.generate_image()
@@ -226,16 +185,6 @@ class WarpedImageEnergy(Energy):
self.def_grad_image_filename = self.image_series.get_image_grad_filename(k_frame)
self.DIdef.init_image(self.def_grad_image_filename)
- self.printer.print_str("Loading previous image for correlation energy…")
-
- # Iold
- self.old_image_filename = self.image_series.get_image_filename(k_frame_old)
- self.Iold.init_image(self.old_image_filename)
-
- # DIold
- self.old_grad_image_filename = self.image_series.get_image_grad_filename(k_frame_old)
- self.DIold.init_image(self.old_grad_image_filename)
-
def call_after_solve(self):
@@ -246,7 +195,7 @@ class WarpedImageEnergy(Energy):
def get_qoi_names(self):
- return [self.name+"_ener", self.name+"_err"]
+ return [self.name+"_ener"]
@@ -254,7 +203,5 @@ class WarpedImageEnergy(Energy):
self.ener = (dolfin.assemble(self.ener_form)/self.problem.mesh_V0)**(1./2)
self.printer.print_sci(self.name+"_ener",self.ener)
- self.err = self.ener/self.Iref_norm
- self.printer.print_sci(self.name+"_err",self.err)
- return [self.ener, self.err]
+ return [self.ener]
diff --git a/NonlinearSolver.py b/NonlinearSolver.py
index 7bba64cf5312416c168934068df05ea97f256378..c983d8dbc25274030d3f1af97b4953c3a2f2ea15 100644
--- a/NonlinearSolver.py
+++ b/NonlinearSolver.py
@@ -179,6 +179,8 @@ class NonlinearSolver():
self.res_old_vec = self.res_vec.copy()
self.res_old_norm = self.res_norm
+ self.problem.call_before_assembly()
+
# linear system: residual assembly
self.printer.print_str("Residual assembly…",newline=False)
timer = time.time()
@@ -281,6 +283,7 @@ class NonlinearSolver():
self.printer.print_var("relax_k",relax_k,-1)
# self.printer.print_sci("relax_a",relax_a)
# self.printer.print_sci("relax_b",relax_b)
+ self.problem.call_before_assembly()
if (need_update_c):
relax_c = relax_b - (relax_b - relax_a) / phi
relax_list.append(relax_c)
diff --git a/Problem_ImageRegistration.py b/Problem_ImageRegistration.py
index fdf6396a8dc1e490f6abd6d46089f2d902bfd691..d0898f9edd6cf8346c5add70472f96fb5b00ebb5 100644
--- a/Problem_ImageRegistration.py
+++ b/Problem_ImageRegistration.py
@@ -179,6 +179,17 @@ class ImageRegistrationProblem(Problem):
+ def call_before_assembly(self,
+ *kargs,
+ **kwargs):
+
+ for energy in self.energies:
+ energy.call_before_assembly(
+ *kargs,
+ **kwargs)
+
+
+
def assemble_ener(self):
ener_form = 0.
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 4b98af09470cc064994c0f6b46ebe0fc485e80c5..361fe93094eb34320c0ead07f6bf9cec2b04de09 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -14,6 +14,7 @@ def get_ExprIm_cpp(
im_dim,
im_type="im",
im_is_def=0,
+ disp_type="fenics", # "vtk"
verbose=0):
assert (im_dim in (2,3))