Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion muDIC/mesh/meshUtilities.py
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ def ok(event):
reset_button = Button(but_ax2, 'Reset')

rectangle = RectangleSelector(overview, line_select_callback,
drawtype='box', useblit=True,
useblit=True,
button=[1, 3], # don't use middle button
minspanx=5, minspany=5,
spancoords='pixels')
Expand Down
151 changes: 140 additions & 11 deletions muDIC/vlab/deformation_fields.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import numpy as np


def harmonic_x(xs, _, amp=1.1, omega=0.05 * np.pi, frame=1):
def harmonic_x(xs, _, amp=1.1, omega=0.05 * np.pi, Cx=0, Cy=0, frame=1):
"""
Displacement field being sinusoidal along the x-axis

Expand All @@ -14,6 +14,10 @@ def harmonic_x(xs, _, amp=1.1, omega=0.05 * np.pi, frame=1):
The y-coordinates
omega : float
The angular frequency in radians/pixel
Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation


Returns
Expand All @@ -26,11 +30,44 @@ def harmonic_x(xs, _, amp=1.1, omega=0.05 * np.pi, frame=1):
# Center around x and y
xs = xs.astype(float)

xs_mapped = amp * np.sin(omega * xs) * float(frame)
xs_mapped = amp * np.sin(omega * (xs-Cx)) * float(frame)
return xs_mapped, np.zeros_like(xs_mapped)


def linear_x(xs, _, slope=0.001, frame=1):
def harmonic_y(_, ys, amp=1.1, omega=0.05 * np.pi, Cx=0, Cy=0, frame=1):
"""
Displacement field being sinusoidal along the y-axis


Parameters
----------
xs : array
The x-coodinates
ys : array
The y-coordinates
omega : float
The angular frequency in radians/pixel
Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation


Returns
-------
u_x,u_y : array
The displacement field values in each direction

"""

# Center around x and y
ys = ys.astype(float)

ys_mapped = amp * np.sin(omega * (ys-Cy)) * float(frame)
return np.zeros_like(ys_mapped), ys_mapped


def linear_x(xs, _, slope=0.001, Cx=np.nan, Cy=np.nan, frame=1):
"""
Displacement field being linearly increasig along x with zero in the center

Expand All @@ -43,7 +80,11 @@ def linear_x(xs, _, slope=0.001, frame=1):
The y-coordinates
slope : float
Displacement increment per pixel along X

Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation


Returns
-------
Expand All @@ -54,13 +95,54 @@ def linear_x(xs, _, slope=0.001, frame=1):

# Center around x and y
xs = xs.astype(float)
center = (xs.max() - xs.min()) / 2.
if np.isnan(Cx):
center = (xs.max() - xs.min()) / 2.
else:
center = Cx

xs_mapped = float(frame) * slope * xs - float(frame) * slope * center
xs_mapped = float(frame) * slope * (xs - center)
return xs_mapped, np.zeros_like(xs_mapped)


def harmonic_bilat(xs, ys, amp=1.1, omega=0.05 * np.pi, frame=1):
def linear_y(_, ys, slope=0.001, Cx=np.nan, Cy=np.nan, frame=1):
"""
Displacement field being linearly increasig along y with zero in the center


Parameters
----------
xs : array
The x-coodinates
ys : array
The y-coordinates
slope : float
Displacement increment per pixel along X
Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation

Returns
-------
u_x,u_y : array
The displacement field values in each direction

"""

# Center around x and y
ys = ys.astype(float)
if np.isnan(Cy):
center = (ys.max() - ys.min()) / 2.
else:
center = Cy

ys_mapped = float(frame) * slope * (ys - center)
return np.zeros_like(ys_mapped), ys_mapped


def harmonic_bilat(xs, ys, ampx=1.1, ampy=1.1,
omegax=0.05 * np.pi, omegay=0.05 * np.pi,
Cx=0, Cy=0, frame=1):
"""
Displacement field being sinusoidal along the both axes

Expand All @@ -71,10 +153,51 @@ def harmonic_bilat(xs, ys, amp=1.1, omega=0.05 * np.pi, frame=1):
The x-coodinates
ys : array
The y-coordinates
omega : float
The angular frequency
ampx : fload
Amplitude in x
ampy : fload
Amplitude in y
omegax : float
The angular frequency in x
omegay : float
The angular frequency in y
img_shape : tuple
The shape of the field
Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation

Returns
-------
u_x,u_y : array
The displacement field values in each direction

"""

xs_mapped = ampx * np.sin(omegax * (xs-Cx)) * np.sin(omegay * (ys-Cy)) * \
float(frame)
ys_mapped = ampy * np.sin(omegax * (xs-Cx)) * np.sin(omegay * (ys-Cy)) * \
float(frame)
return xs_mapped, ys_mapped


def rigid_rotation(xs, ys, theta=np.pi/50, Cx=0, Cy=0, frame=1):
"""
Displacement field from rigid rotation about the out of plane axis.

Parameters
----------
xs : array
The x-coodinates
ys : array
The y-coordinates
theta : float
(default np.pi/50) Rotation angle
Cx : float
X Location of the center for the deformation
Cy : float
Y Location of the center for the deformation


Returns
Expand All @@ -84,6 +207,12 @@ def harmonic_bilat(xs, ys, amp=1.1, omega=0.05 * np.pi, frame=1):

"""

xs_mapped = amp * np.sin(omega * xs) * np.sin(omega * ys) * float(frame)
ys_mapped = amp * np.sin(omega * xs) * np.sin(omega * ys) * float(frame)
xs_mapped = np.cos(theta * float(frame))*(xs-Cx) - \
np.sin(theta * float(frame))*(ys-Cy) - (xs-Cx)
ys_mapped = np.sin(theta * float(frame))*(xs-Cx) + \
np.cos(theta * float(frame))*(ys-Cy) - (ys-Cy)

# Warning: inversion starts zooming out for large rotations.
# xs_mapped = -(ys-Cy)*(theta*float(frame))
# ys_mapped = (xs-Cx)*(theta*float(frame))
return xs_mapped, ys_mapped
75 changes: 47 additions & 28 deletions muDIC/vlab/image_deformer.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,17 @@ def num_diff(xs, ys, func, component=(1, 1)):
raise ValueError("Incompatible component request")


def inverse(xs, ys, function, tol=1e-6,frame=1):
def inverse(xs, ys, function, tol=1e-6, frame=1):
"""Calculates the inverse of a displacement function

Returns the arguments corresponding to the function values given.
This is done by a Newton solver which minimizes:

g(X) = X + f(X,Y) - x = 0
where x is the coordinate of the displacement material point, having the coordinate X in its undeformed state.
The coordinate Y is kept constant thoughout the optimization, and the routine is therefore not guarantied to
where x is the coordinate of the displacement material point, having the
coordinate X in its undeformed state.
The coordinate Y is kept constant thoughout the optimization, and the
routine is therefore not guarantied to
converge.

Parameters
Expand All @@ -56,11 +58,11 @@ def inverse(xs, ys, function, tol=1e-6,frame=1):
logger = logging.getLogger()

def f_x(x, y, a_x):
u_x, _ = function(x, y,frame=frame)
u_x, _ = function(x, y, frame=frame)
return x + u_x - a_x

def f_y(y, x, a_y):
_, u_y = function(x, y,frame=frame)
_, u_y = function(x, y, frame=frame)
return y + u_y - a_y

x_deformed = xs.flatten()
Expand All @@ -69,23 +71,27 @@ def f_y(y, x, a_y):
y_deformed = ys.flatten()
y_undeformed_guess = y_deformed.copy()

x = optimize.newton(f_x, x_undeformed_guess, args=(y_undeformed_guess, x_deformed))
y = optimize.newton(f_y, y_undeformed_guess, args=(x_undeformed_guess, y_deformed))
x = optimize.newton(f_x, x_undeformed_guess,
args=(y_undeformed_guess, x_deformed))
y = optimize.newton(f_y, y_undeformed_guess,
args=(x_undeformed_guess, y_deformed))
# A COUPLED APPROACH CAN'T BE THAT HARD TO IMPLEMENT!!

error_x = np.max(np.abs(f_x(np.array(x), y_undeformed_guess, x_deformed)))
error_y = np.max(np.abs(f_y(np.array(y), x_undeformed_guess, y_deformed)))

if max(error_x, error_y) > tol:
raise ValueError("The displacement function was not inverted successfully")
raise ValueError("The displacement function was not inverted "
"successfully")

logger.info("Larges error in y is: %f" % error_y)
logger.info("Larges error in x is: %f" % error_x)
logger.info("Largest error in y is: %f" % error_y)
logger.info("Largest error in x is: %f" % error_x)

return np.reshape(x, xs.shape), np.reshape(y, ys.shape)


class ImageDeformer(object):
def __init__(self, coordinate_mapper, multiplicative,order=4):
def __init__(self, coordinate_mapper, multiplicative, order=4):
"""Image deformer

Deforms an image according to the settings given at instantiation.
Expand All @@ -94,27 +100,32 @@ def __init__(self, coordinate_mapper, multiplicative,order=4):
Parameters
----------
coordinate_mapper : func
A function taking two coordinates and returns two mapped coordinates
A function taking two coordinates and returns two mapped
coordinates
multiplicative : bool
Whether the coordinate mapping is multiplicative or not. Use True for deformation gradients
and use False for displacement fields.
Whether the coordinate mapping is multiplicative or not. Use True
for deformation gradients and use False for displacement fields.
order : int
The interpolation order used

Example
-------
Make a ImageDeformed which deforms the image according to a deformation gradient.
Make a ImageDeformed which deforms the image according to a
deformation gradient.
NOTE: You should rather use the factory imageDeformer_from_defGrad.

>>> import numpy as np
>>> from muDIC import vlab
>>> F = np.array([[1.1,.0], [0., 1.0]], dtype=float)
>>> coordinate_mapper = vlab.image_deformer.CoordinateMapper(partial(map_coords_by_defgrad, F=F))
>>> vlab.image_deformer.ImageDeformer(coordinate_mapper,multiplicative=True)
>>> coordinate_mapper = vlab.image_deformer.CoordinateMapper(
partial(map_coords_by_defgrad, F=F))
>>> vlab.image_deformer.ImageDeformer(coordinate_mapper,
multiplicative=True)

NOTE
----
The image deformed returns a list of deformed images when called, the first being undeformed
The image deformed returns a list of deformed images when called,
the first being undeformed

"""
self.multiplicative = multiplicative
Expand All @@ -125,7 +136,6 @@ def __call__(self, img, steps=2):
n, m = np.shape(img)
def_imgs = []


xn, yn = np.meshgrid(np.arange(m), np.arange(n))
xn_mapped, yn_mapped = xn, yn

Expand All @@ -134,11 +144,14 @@ def __call__(self, img, steps=2):
Ic = img
else:
if self.multiplicative:
xn_mapped, yn_mapped = self.coodinate_mapper(xn_mapped, yn_mapped)
xn_mapped, yn_mapped = self.coodinate_mapper(xn_mapped,
yn_mapped)
else:
xn_mapped, yn_mapped = self.coodinate_mapper(xn, yn,frame=i)
xn_mapped, yn_mapped = self.coodinate_mapper(xn, yn,
frame=i)

Ic = nd.map_coordinates(img, np.array([yn_mapped, xn_mapped]), order=self.order, cval=0)
Ic = nd.map_coordinates(img, np.array([yn_mapped, xn_mapped]),
order=self.order, cval=0)

def_imgs.append(Ic)

Expand All @@ -147,7 +160,8 @@ def __call__(self, img, steps=2):

def map_coords_by_defgrad(xs, ys, F):
"""Maps coordinates based on a deformation gradient
NOTE. The coordinates are first centered, implying that a rotation is about the middle of the image.
NOTE. The coordinates are first centered, implying that a rotation is about
the middle of the image.

Parameters
----------
Expand Down Expand Up @@ -192,14 +206,16 @@ def imageDeformer_from_defGrad(F):

Examples
--------
First, we define the deformation gradient and then make an ImageDeformer object:
First, we define the deformation gradient and then make an
ImageDeformer object:
>>> import numpy as np
>>> from muDIC import vlab
>>> F = np.array([[1.1,.0], [0., 1.0]], dtype=float)
>>> image_deformer = vlab.imageDeformer_from_defGrad(F)

"""
return ImageDeformer(partial(map_coords_by_defgrad, F=F), multiplicative=True)
return ImageDeformer(partial(map_coords_by_defgrad, F=F),
multiplicative=True)


def imageDeformer_from_uFunc(uFunc, **kwargs):
Expand All @@ -213,13 +229,16 @@ def imageDeformer_from_uFunc(uFunc, **kwargs):
Returns
-------
imageDeformer
An ImageDeformer object which deformes images according to the displacement function
An ImageDeformer object which deformes images according to the
displacement function

Examples
--------
First, we define the deformation gradient and then make an ImageDeformer object:
First, we define the deformation gradient and then make an
ImageDeformer object:
>>> from muDIC import vlab
>>> image_deformer = vlab.imageDeformer_from_uFunc(vlab.deformation_fields.bilat_harmonic,omega=4*np.pi)
>>> image_deformer = vlab.imageDeformer_from_uFunc(
vlab.deformation_fields.bilat_harmonic,omega=4*np.pi)

"""
coordinate_generator = partial(uFunc, **kwargs)
Expand Down