Module pycpd.deformable_registration
Expand source code
from builtins import super
import numpy as np
import numbers
from .emregistration import EMRegistration
from .utility import gaussian_kernel, low_rank_eigen
class DeformableRegistration(EMRegistration):
"""
Deformable registration.
Attributes
----------
alpha: float (positive)
Represents the trade-off between the goodness of maximum likelihood fit and regularization.
beta: float(positive)
Width of the Gaussian kernel.
low_rank: bool
Whether to use low rank approximation.
num_eig: int
Number of eigenvectors to use in lowrank calculation.
"""
def __init__(self, alpha=None, beta=None, low_rank=False, num_eig=100, *args, **kwargs):
super().__init__(*args, **kwargs)
if alpha is not None and (not isinstance(alpha, numbers.Number) or alpha <= 0):
raise ValueError(
"Expected a positive value for regularization parameter alpha. Instead got: {}".format(alpha))
if beta is not None and (not isinstance(beta, numbers.Number) or beta <= 0):
raise ValueError(
"Expected a positive value for the width of the coherent Gaussian kerenl. Instead got: {}".format(beta))
self.alpha = 2 if alpha is None else alpha
self.beta = 2 if beta is None else beta
self.W = np.zeros((self.M, self.D))
self.G = gaussian_kernel(self.Y, self.beta)
self.low_rank = low_rank
self.num_eig = num_eig
if self.low_rank is True:
self.Q, self.S = low_rank_eigen(self.G, self.num_eig)
self.inv_S = np.diag(1./self.S)
self.S = np.diag(self.S)
self.E = 0.
def update_transform(self):
"""
Calculate a new estimate of the deformable transformation.
See Eq. 22 of https://arxiv.org/pdf/0905.2635.pdf.
"""
if self.low_rank is False:
A = np.dot(np.diag(self.P1), self.G) + \
self.alpha * self.sigma2 * np.eye(self.M)
B = self.PX - np.dot(np.diag(self.P1), self.Y)
self.W = np.linalg.solve(A, B)
elif self.low_rank is True:
# Matlab code equivalent can be found here:
# https://github.com/markeroon/matlab-computer-vision-routines/tree/master/third_party/CoherentPointDrift
dP = np.diag(self.P1)
dPQ = np.matmul(dP, self.Q)
F = self.PX - np.matmul(dP, self.Y)
self.W = 1 / (self.alpha * self.sigma2) * (F - np.matmul(dPQ, (
np.linalg.solve((self.alpha * self.sigma2 * self.inv_S + np.matmul(self.Q.T, dPQ)),
(np.matmul(self.Q.T, F))))))
QtW = np.matmul(self.Q.T, self.W)
self.E = self.E + self.alpha / 2 * np.trace(np.matmul(QtW.T, np.matmul(self.S, QtW)))
def transform_point_cloud(self, Y=None):
"""
Update a point cloud using the new estimate of the deformable transformation.
Attributes
----------
Y: numpy array, optional
Array of points to transform - use to predict on new set of points.
Best for predicting on new points not used to run initial registration.
If None, self.Y used.
Returns
-------
If Y is None, returns None.
Otherwise, returns the transformed Y.
"""
if Y is not None:
G = gaussian_kernel(X=Y, beta=self.beta, Y=self.Y)
return Y + np.dot(G, self.W)
else:
if self.low_rank is False:
self.TY = self.Y + np.dot(self.G, self.W)
elif self.low_rank is True:
self.TY = self.Y + np.matmul(self.Q, np.matmul(self.S, np.matmul(self.Q.T, self.W)))
return
def update_variance(self):
"""
Update the variance of the mixture model using the new estimate of the deformable transformation.
See the update rule for sigma2 in Eq. 23 of of https://arxiv.org/pdf/0905.2635.pdf.
"""
qprev = self.sigma2
# The original CPD paper does not explicitly calculate the objective functional.
# This functional will include terms from both the negative log-likelihood and
# the Gaussian kernel used for regularization.
self.q = np.inf
xPx = np.dot(np.transpose(self.Pt1), np.sum(
np.multiply(self.X, self.X), axis=1))
yPy = np.dot(np.transpose(self.P1), np.sum(
np.multiply(self.TY, self.TY), axis=1))
trPXY = np.sum(np.multiply(self.TY, self.PX))
self.sigma2 = (xPx - 2 * trPXY + yPy) / (self.Np * self.D)
if self.sigma2 <= 0:
self.sigma2 = self.tolerance / 10
# Here we use the difference between the current and previous
# estimate of the variance as a proxy to test for convergence.
self.diff = np.abs(self.sigma2 - qprev)
def get_registration_parameters(self):
"""
Return the current estimate of the deformable transformation parameters.
Returns
-------
self.G: numpy array
Gaussian kernel matrix.
self.W: numpy array
Deformable transformation matrix.
"""
return self.G, self.W
Classes
class DeformableRegistration (alpha=None, beta=None, low_rank=False, num_eig=100, *args, **kwargs)
-
Deformable registration.
Attributes
alpha
:float (positive)
- Represents the trade-off between the goodness of maximum likelihood fit and regularization.
beta
:float(positive)
- Width of the Gaussian kernel.
low_rank
:bool
- Whether to use low rank approximation.
num_eig
:int
- Number of eigenvectors to use in lowrank calculation.
Expand source code
class DeformableRegistration(EMRegistration): """ Deformable registration. Attributes ---------- alpha: float (positive) Represents the trade-off between the goodness of maximum likelihood fit and regularization. beta: float(positive) Width of the Gaussian kernel. low_rank: bool Whether to use low rank approximation. num_eig: int Number of eigenvectors to use in lowrank calculation. """ def __init__(self, alpha=None, beta=None, low_rank=False, num_eig=100, *args, **kwargs): super().__init__(*args, **kwargs) if alpha is not None and (not isinstance(alpha, numbers.Number) or alpha <= 0): raise ValueError( "Expected a positive value for regularization parameter alpha. Instead got: {}".format(alpha)) if beta is not None and (not isinstance(beta, numbers.Number) or beta <= 0): raise ValueError( "Expected a positive value for the width of the coherent Gaussian kerenl. Instead got: {}".format(beta)) self.alpha = 2 if alpha is None else alpha self.beta = 2 if beta is None else beta self.W = np.zeros((self.M, self.D)) self.G = gaussian_kernel(self.Y, self.beta) self.low_rank = low_rank self.num_eig = num_eig if self.low_rank is True: self.Q, self.S = low_rank_eigen(self.G, self.num_eig) self.inv_S = np.diag(1./self.S) self.S = np.diag(self.S) self.E = 0. def update_transform(self): """ Calculate a new estimate of the deformable transformation. See Eq. 22 of https://arxiv.org/pdf/0905.2635.pdf. """ if self.low_rank is False: A = np.dot(np.diag(self.P1), self.G) + \ self.alpha * self.sigma2 * np.eye(self.M) B = self.PX - np.dot(np.diag(self.P1), self.Y) self.W = np.linalg.solve(A, B) elif self.low_rank is True: # Matlab code equivalent can be found here: # https://github.com/markeroon/matlab-computer-vision-routines/tree/master/third_party/CoherentPointDrift dP = np.diag(self.P1) dPQ = np.matmul(dP, self.Q) F = self.PX - np.matmul(dP, self.Y) self.W = 1 / (self.alpha * self.sigma2) * (F - np.matmul(dPQ, ( np.linalg.solve((self.alpha * self.sigma2 * self.inv_S + np.matmul(self.Q.T, dPQ)), (np.matmul(self.Q.T, F)))))) QtW = np.matmul(self.Q.T, self.W) self.E = self.E + self.alpha / 2 * np.trace(np.matmul(QtW.T, np.matmul(self.S, QtW))) def transform_point_cloud(self, Y=None): """ Update a point cloud using the new estimate of the deformable transformation. Attributes ---------- Y: numpy array, optional Array of points to transform - use to predict on new set of points. Best for predicting on new points not used to run initial registration. If None, self.Y used. Returns ------- If Y is None, returns None. Otherwise, returns the transformed Y. """ if Y is not None: G = gaussian_kernel(X=Y, beta=self.beta, Y=self.Y) return Y + np.dot(G, self.W) else: if self.low_rank is False: self.TY = self.Y + np.dot(self.G, self.W) elif self.low_rank is True: self.TY = self.Y + np.matmul(self.Q, np.matmul(self.S, np.matmul(self.Q.T, self.W))) return def update_variance(self): """ Update the variance of the mixture model using the new estimate of the deformable transformation. See the update rule for sigma2 in Eq. 23 of of https://arxiv.org/pdf/0905.2635.pdf. """ qprev = self.sigma2 # The original CPD paper does not explicitly calculate the objective functional. # This functional will include terms from both the negative log-likelihood and # the Gaussian kernel used for regularization. self.q = np.inf xPx = np.dot(np.transpose(self.Pt1), np.sum( np.multiply(self.X, self.X), axis=1)) yPy = np.dot(np.transpose(self.P1), np.sum( np.multiply(self.TY, self.TY), axis=1)) trPXY = np.sum(np.multiply(self.TY, self.PX)) self.sigma2 = (xPx - 2 * trPXY + yPy) / (self.Np * self.D) if self.sigma2 <= 0: self.sigma2 = self.tolerance / 10 # Here we use the difference between the current and previous # estimate of the variance as a proxy to test for convergence. self.diff = np.abs(self.sigma2 - qprev) def get_registration_parameters(self): """ Return the current estimate of the deformable transformation parameters. Returns ------- self.G: numpy array Gaussian kernel matrix. self.W: numpy array Deformable transformation matrix. """ return self.G, self.W
Ancestors
Subclasses
Methods
def get_registration_parameters(self)
-
Return the current estimate of the deformable transformation parameters.
Returns
self.G: numpy array
- Gaussian kernel matrix.
self.W: numpy array
- Deformable transformation matrix.
Expand source code
def get_registration_parameters(self): """ Return the current estimate of the deformable transformation parameters. Returns ------- self.G: numpy array Gaussian kernel matrix. self.W: numpy array Deformable transformation matrix. """ return self.G, self.W
def transform_point_cloud(self, Y=None)
-
Update a point cloud using the new estimate of the deformable transformation.
Attributes
Y
:numpy array
, optional- Array of points to transform - use to predict on new set of points. Best for predicting on new points not used to run initial registration. If None, self.Y used.
Returns
If Y is None, returns None. Otherwise, returns the transformed Y.
Expand source code
def transform_point_cloud(self, Y=None): """ Update a point cloud using the new estimate of the deformable transformation. Attributes ---------- Y: numpy array, optional Array of points to transform - use to predict on new set of points. Best for predicting on new points not used to run initial registration. If None, self.Y used. Returns ------- If Y is None, returns None. Otherwise, returns the transformed Y. """ if Y is not None: G = gaussian_kernel(X=Y, beta=self.beta, Y=self.Y) return Y + np.dot(G, self.W) else: if self.low_rank is False: self.TY = self.Y + np.dot(self.G, self.W) elif self.low_rank is True: self.TY = self.Y + np.matmul(self.Q, np.matmul(self.S, np.matmul(self.Q.T, self.W))) return
def update_transform(self)
-
Calculate a new estimate of the deformable transformation. See Eq. 22 of https://arxiv.org/pdf/0905.2635.pdf.
Expand source code
def update_transform(self): """ Calculate a new estimate of the deformable transformation. See Eq. 22 of https://arxiv.org/pdf/0905.2635.pdf. """ if self.low_rank is False: A = np.dot(np.diag(self.P1), self.G) + \ self.alpha * self.sigma2 * np.eye(self.M) B = self.PX - np.dot(np.diag(self.P1), self.Y) self.W = np.linalg.solve(A, B) elif self.low_rank is True: # Matlab code equivalent can be found here: # https://github.com/markeroon/matlab-computer-vision-routines/tree/master/third_party/CoherentPointDrift dP = np.diag(self.P1) dPQ = np.matmul(dP, self.Q) F = self.PX - np.matmul(dP, self.Y) self.W = 1 / (self.alpha * self.sigma2) * (F - np.matmul(dPQ, ( np.linalg.solve((self.alpha * self.sigma2 * self.inv_S + np.matmul(self.Q.T, dPQ)), (np.matmul(self.Q.T, F)))))) QtW = np.matmul(self.Q.T, self.W) self.E = self.E + self.alpha / 2 * np.trace(np.matmul(QtW.T, np.matmul(self.S, QtW)))
def update_variance(self)
-
Update the variance of the mixture model using the new estimate of the deformable transformation. See the update rule for sigma2 in Eq. 23 of of https://arxiv.org/pdf/0905.2635.pdf.
Expand source code
def update_variance(self): """ Update the variance of the mixture model using the new estimate of the deformable transformation. See the update rule for sigma2 in Eq. 23 of of https://arxiv.org/pdf/0905.2635.pdf. """ qprev = self.sigma2 # The original CPD paper does not explicitly calculate the objective functional. # This functional will include terms from both the negative log-likelihood and # the Gaussian kernel used for regularization. self.q = np.inf xPx = np.dot(np.transpose(self.Pt1), np.sum( np.multiply(self.X, self.X), axis=1)) yPy = np.dot(np.transpose(self.P1), np.sum( np.multiply(self.TY, self.TY), axis=1)) trPXY = np.sum(np.multiply(self.TY, self.PX)) self.sigma2 = (xPx - 2 * trPXY + yPy) / (self.Np * self.D) if self.sigma2 <= 0: self.sigma2 = self.tolerance / 10 # Here we use the difference between the current and previous # estimate of the variance as a proxy to test for convergence. self.diff = np.abs(self.sigma2 - qprev)
Inherited members