diff --git a/Energy_WarpedImage.py b/Energy_WarpedImage.py
index 94113039c4a2d02de32a012b73895701e88df701..6ad2162dcdaa50ce6b97e6bab886b670bfdf7be4 100644
--- a/Energy_WarpedImage.py
+++ b/Energy_WarpedImage.py
@@ -30,6 +30,7 @@ class WarpedImageEnergy(Energy):
name="im",
w=1.,
ref_frame=0,
+ im_is_cone=0,
dynamic_scaling=False):
self.problem = problem
@@ -182,20 +183,32 @@ class WarpedImageEnergy(Energy):
# Phi_ref
self.Phi_Iref = dolfin.Expression(
cppcode=ddic.get_ExprCharFuncIm_cpp(
- im_dim=self.image_series.dimension),
+ im_dim=self.image_series.dimension,
+ im_is_def=0,
+ im_is_cone=im_is_cone),
element=self.fe)
self.Phi_Iref.init_image(self.ref_image_filename)
+ # Phi_def
+ self.Phi_Idef = dolfin.Expression(
+ cppcode=ddic.get_ExprCharFuncIm_cpp(
+ im_dim=self.image_series.dimension,
+ im_is_def=1,
+ im_is_cone=im_is_cone),
+ element=self.fe)
+ self.Phi_Idef.init_image(self.ref_image_filename)
+ self.Phi_Idef.init_disp(self.problem.U)
+
# 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.Psi_c = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iref)**2/2
+ self.DPsi_c = self.Phi_Idef * 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.DDPsi_c = self.Phi_Idef * 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_Idef * 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_Idef * 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_old = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iold)**2/2
+ self.DPsi_c_old = self.Phi_Idef * self.Phi_Iref * (self.Idef - self.Iold) * dolfin.dot(self.DIdef, self.problem.dU_test)
# forms
self.ener_form = self.Psi_c * self.dV
diff --git a/fedic2.py b/fedic2.py
index cd2e049990f8040a8bf3230f0cb8b033f5da8b6d..a3b5f2a00b164cf9045b00e8278441aec7e29086 100644
--- a/fedic2.py
+++ b/fedic2.py
@@ -28,6 +28,7 @@ def fedic2(
images_quadrature_from="points_count", # points_count, integral
images_expressions_type="cpp", # cpp, py
images_dynamic_scaling=1,
+ images_is_cone=0,
mesh=None,
mesh_folder=None,
mesh_basename=None,
@@ -113,7 +114,8 @@ def fedic2(
quadrature_degree=images_quadrature,
w=1.-regul_level,
ref_frame=images_ref_frame,
- dynamic_scaling=images_dynamic_scaling)
+ dynamic_scaling=images_dynamic_scaling,
+ im_is_cone=images_is_cone)
problem.add_image_energy(warped_image_energy)
regularization_energy = ddic.RegularizationEnergy(
diff --git a/image_expressions_cpp.py b/image_expressions_cpp.py
index 541e41d8ae2f766a5bd96100abd4e9976ce70200..57592540161684f5e9fe7fb7206b6606a4574d09 100644
--- a/image_expressions_cpp.py
+++ b/image_expressions_cpp.py
@@ -245,11 +245,16 @@ public:
def get_ExprCharFuncIm_cpp(
im_dim,
+ im_is_def=0,
+ im_is_cone=0,
verbose=0):
assert (im_dim in (2,3))
+ if (im_is_cone):
+ assert (im_dim == 3)
ExprCharFuncIm_cpp = '''\
+#include <math.h>
#include <string.h>
#include <vtkSmartPointer.h>
@@ -262,13 +267,35 @@ namespace dolfin
class MyExpr : public Expression
{
double xmin, xmax, ymin, ymax, zmin, zmax;
+ mutable Array<double> UX;
+ mutable Array<double> x;'''+('''
+ Function* U;''')*(im_is_def)+('''
+ Array<double> O;
+ Array<double> n1;
+ Array<double> n2;
+ Array<double> n3;
+ Array<double> n4;
+ mutable double d1, d2, d3, d4;''')*(im_is_cone)+'''
public:
MyExpr():
- Expression()
+ Expression(),
+ UX('''+str(im_dim)+'''),
+ x('''+str(im_dim)+''')'''+(''',
+ O(3),
+ n1(3),
+ n2(3),
+ n3(3),
+ n4(3)''')*(im_is_cone)+'''
{
- }
+ }'''+('''
+
+ void init_disp(
+ Function* UU)
+ {
+ U = UU;
+ }''')*(im_is_def)+'''
void init_image(
const char* filename)
@@ -291,16 +318,71 @@ public:
std::cout << "ymin = " << ymin << std::endl;
std::cout << "ymax = " << ymax << std::endl;
std::cout << "zmin = " << zmin << std::endl;
- std::cout << "zmax = " << zmax << std::endl;''')*(verbose)+'''
+ std::cout << "zmax = " << zmax << std::endl;''')*(verbose)+('''
+
+ O[0] = (xmin+xmax)/2;
+ O[1] = (ymin+ymax)/2;
+ O[2] = zmax;
+
+ n1[0] = +cos(35. * M_PI/180.);
+ n1[1] = 0.;
+ n1[2] = -sin(35. * M_PI/180.);
+
+ n2[0] = -cos(35. * M_PI/180.);
+ n2[1] = 0.;
+ n2[2] = -sin(35. * M_PI/180.);
+
+ n3[0] = 0.;
+ n3[1] = +cos(40. * M_PI/180.);
+ n3[2] = -sin(40. * M_PI/180.);
+
+ n4[0] = 0.;
+ n4[1] = -cos(40. * M_PI/180.);
+ n4[2] = -sin(40. * M_PI/180.);''')*(im_is_cone)+'''
}
void eval(Array<double>& expr, const Array<double>& X) const
{'''+('''
std::cout << "X = " << X.str(1) << std::endl;''')*(verbose)+('''
- if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax))''')*(im_dim==2)+('''
+ U->eval(UX, X);''')*(im_is_def)+('''
+ std::cout << "UX = " << UX.str(1) << std::endl;''')*(verbose)+('''
- if ((X[0] >= xmin) && (X[0] <= xmax) && (X[1] >= ymin) && (X[1] <= ymax) && (X[2] >= zmin) && (X[2] <= zmax))''')*(im_dim==3)+'''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];''')*(im_dim==2)+('''
+ x[0] = X[0] + UX[0];
+ x[1] = X[1] + UX[1];
+ x[2] = X[2] + UX[2];''')*(im_dim==3)+('''
+ std::cout << "x = " << x.str(1) << std::endl;''')*(verbose)+(('''
+
+ if ((x[0] >= xmin)
+ && (x[0] <= xmax)
+ && (x[1] >= ymin)
+ && (x[1] <= ymax))''')*(im_dim==2)+('''
+ if ((x[0] >= xmin)
+ && (x[0] <= xmax)
+ && (x[1] >= ymin)
+ && (x[1] <= ymax)
+ && (x[2] >= zmin)
+ && (x[2] <= zmax))''')*(im_dim==3))*(not im_is_cone)+('''
+
+ d1 = (x[0]-O[0])*n1[0]
+ + (x[1]-O[1])*n1[1]
+ + (x[2]-O[2])*n1[2];
+ d2 = (x[0]-O[0])*n2[0]
+ + (x[1]-O[1])*n2[1]
+ + (x[2]-O[2])*n2[2];
+ d3 = (x[0]-O[0])*n3[0]
+ + (x[1]-O[1])*n3[1]
+ + (x[2]-O[2])*n3[2];
+ d4 = (x[0]-O[0])*n4[0]
+ + (x[1]-O[1])*n4[1]
+ + (x[2]-O[2])*n4[2];
+
+ if ((d1 >= 0.)
+ && (d2 >= 0.)
+ && (d3 >= 0.)
+ && (d4 >= 0.))''')*(im_is_cone)+'''
{
expr[0] = 1.;
}
@@ -313,5 +395,5 @@ public:
};
}'''
- #print ExprCharFuncIm_cpp
+ # print ExprCharFuncIm_cpp
return ExprCharFuncIm_cpp