reaarrange files
This commit is contained in:
@@ -0,0 +1,210 @@
|
||||
import numpy
|
||||
import pygeodesic.geodesic as geodesic
|
||||
|
||||
__all__ = ['find_coefficients', 'transform']
|
||||
|
||||
|
||||
def cdist(K: numpy.ndarray, B: numpy.ndarray) -> numpy.ndarray:
|
||||
"""Calculate Euclidean distance between K[i, :] and B[j, :].
|
||||
|
||||
Arguments
|
||||
---------
|
||||
K : numpy.array
|
||||
B : numpy.array
|
||||
"""
|
||||
K = numpy.atleast_2d(K)
|
||||
B = numpy.atleast_2d(B)
|
||||
assert K.ndim == 2
|
||||
assert B.ndim == 2
|
||||
|
||||
K = numpy.expand_dims(K, 1)
|
||||
B = numpy.expand_dims(B, 0)
|
||||
D = K - B
|
||||
return numpy.linalg.norm(D, axis=2)
|
||||
|
||||
|
||||
# def cdist(K: numpy.ndarray, B: numpy.ndarray, faces: numpy.ndarray) -> numpy.ndarray:
|
||||
# geoalg = geodesic.PyGeodesicAlgorithmExact(points, faces)
|
||||
|
||||
|
||||
|
||||
def pairwise_radial_basis(K: numpy.ndarray, B: numpy.ndarray) -> numpy.ndarray:
|
||||
"""Compute the TPS radial basis function phi(r) between every row-pair of K
|
||||
and B where r is the Euclidean distance.
|
||||
|
||||
Arguments
|
||||
---------
|
||||
K : numpy.array
|
||||
n by d vector containing n d-dimensional points.
|
||||
B : numpy.array
|
||||
m by d vector containing m d-dimensional points.
|
||||
|
||||
Return
|
||||
------
|
||||
P : numpy.array
|
||||
n by m matrix where.
|
||||
P(i, j) = phi( norm( K(i,:) - B(j,:) ) ),
|
||||
where phi(r) = r^2*log(r), if r >= 1
|
||||
r*log(r^r), if r < 1
|
||||
"""
|
||||
# r_mat(i, j) is the Euclidean distance between K(i, :) and B(j, :).
|
||||
r_mat = cdist(K, B)
|
||||
|
||||
pwise_cond_ind1 = r_mat >= 1
|
||||
pwise_cond_ind2 = r_mat < 1
|
||||
r_mat_p1 = r_mat[pwise_cond_ind1]
|
||||
r_mat_p2 = r_mat[pwise_cond_ind2]
|
||||
|
||||
# P correcponds to the matrix K from [1].
|
||||
P = numpy.empty(r_mat.shape)
|
||||
P[pwise_cond_ind1] = (r_mat_p1**2) * numpy.log(r_mat_p1)
|
||||
P[pwise_cond_ind2] = r_mat_p2 * numpy.log(numpy.power(r_mat_p2, r_mat_p2))
|
||||
|
||||
return P
|
||||
|
||||
|
||||
def find_coefficients(control_points: numpy.ndarray,
|
||||
target_points: numpy.ndarray,
|
||||
lambda_: float = 0.,
|
||||
solver: str = 'exact') -> numpy.ndarray:
|
||||
"""Given a set of control points and their corresponding points, compute the
|
||||
coefficients of the TPS interpolant deforming surface.
|
||||
|
||||
Arguments
|
||||
---------
|
||||
control_points : numpy.array
|
||||
p by d vector of control points
|
||||
target_points : numpy.array
|
||||
p by d vector of corresponding target points on the deformed
|
||||
surface
|
||||
lambda_ : float
|
||||
regularization parameter
|
||||
solver : str
|
||||
the solver to get the coefficients. default is 'exact' for the exact
|
||||
solution. Or use 'lstsq' for the least square solution.
|
||||
|
||||
Return
|
||||
------
|
||||
coef : numpy.ndarray
|
||||
the coefficients
|
||||
|
||||
.. seealso::
|
||||
|
||||
http://cseweb.ucsd.edu/~sjb/pami_tps.pdf
|
||||
"""
|
||||
# ensure data type and shape
|
||||
control_points = numpy.atleast_2d(control_points)
|
||||
target_points = numpy.atleast_2d(target_points)
|
||||
if control_points.shape != target_points.shape:
|
||||
raise ValueError(
|
||||
'Shape of and control points {cp} and target points {tp} are not the same.'.
|
||||
format(cp=control_points.shape, tp=target_points.shape))
|
||||
|
||||
p, d = control_points.shape
|
||||
|
||||
# The matrix
|
||||
K = pairwise_radial_basis(control_points, control_points)
|
||||
P = numpy.hstack([numpy.ones((p, 1)), control_points])
|
||||
|
||||
# Relax the exact interpolation requirement by means of regularization.
|
||||
K = K + lambda_ * numpy.identity(p)
|
||||
|
||||
# Target points
|
||||
M = numpy.vstack([
|
||||
numpy.hstack([K, P]),
|
||||
numpy.hstack([P.T, numpy.zeros((d + 1, d + 1))])
|
||||
])
|
||||
Y = numpy.vstack([target_points, numpy.zeros((d + 1, d))])
|
||||
|
||||
# solve for M*X = Y.
|
||||
# At least d+1 control points should not be in a subspace; e.g. for d=2, at
|
||||
# least 3 points are not on a straight line. Otherwise M will be singular.
|
||||
solver = solver.lower()
|
||||
if solver == 'exact':
|
||||
X = numpy.linalg.solve(M, Y)
|
||||
elif solver == 'lstsq':
|
||||
X, _, _, _ = numpy.linalg.lstsq(M, Y, None)
|
||||
else:
|
||||
raise ValueError('Unknown solver: ' + solver)
|
||||
|
||||
return X
|
||||
|
||||
|
||||
def transform(source_points: numpy.ndarray, control_points: numpy.ndarray,
|
||||
coefficient: numpy.ndarray) -> numpy.ndarray:
|
||||
"""Transform the source points form the original surface to the destination
|
||||
(deformed) surface.
|
||||
|
||||
Arguments
|
||||
---------
|
||||
source_points : numpy.array
|
||||
n by d array of source points to be transformed
|
||||
control_points : numpy.array
|
||||
the control points used in the function `find_coefficients`
|
||||
coefficient : numpy.array
|
||||
the computed coefficients
|
||||
|
||||
Return
|
||||
------
|
||||
deformed_points : numpy.array
|
||||
n by d array of the transformed point on the target surface
|
||||
"""
|
||||
source_points = numpy.atleast_2d(source_points)
|
||||
control_points = numpy.atleast_2d(control_points)
|
||||
if source_points.shape[-1] != control_points.shape[-1]:
|
||||
raise ValueError(
|
||||
'Dimension of source points ({sd}D) and control points ({cd}D) are not the same.'.
|
||||
format(sd=source_points.shape[-1], cd=control_points.shape[-1]))
|
||||
|
||||
n = source_points.shape[0]
|
||||
|
||||
A = pairwise_radial_basis(source_points, control_points)
|
||||
K = numpy.hstack([A, numpy.ones((n, 1)), source_points])
|
||||
|
||||
deformed_points = numpy.dot(K, coefficient)
|
||||
return deformed_points
|
||||
|
||||
|
||||
|
||||
|
||||
class TPS:
|
||||
"""The thin plate spline deformation warpping.
|
||||
"""
|
||||
|
||||
def __init__(self,
|
||||
control_points: numpy.ndarray,
|
||||
target_points: numpy.ndarray,
|
||||
lambda_: float = 0.,
|
||||
solver: str = 'exact'):
|
||||
"""Create a instance that preserve the TPS coefficients.
|
||||
|
||||
Arguments
|
||||
---------
|
||||
control_points : numpy.array
|
||||
p by d vector of control points
|
||||
target_points : numpy.array
|
||||
p by d vector of corresponding target points on the deformed
|
||||
surface
|
||||
lambda_ : float
|
||||
regularization parameter
|
||||
solver : str
|
||||
the solver to get the coefficients. default is 'exact' for the
|
||||
exact solution. Or use 'lstsq' for the least square solution.
|
||||
"""
|
||||
self.control_points = control_points
|
||||
self.coefficient = find_coefficients(
|
||||
control_points, target_points, lambda_, solver)
|
||||
|
||||
def __call__(self, source_points):
|
||||
"""Transform the source points form the original surface to the
|
||||
destination (deformed) surface.
|
||||
|
||||
Arguments
|
||||
---------
|
||||
source_points : numpy.array
|
||||
n by d array of source points to be transformed
|
||||
"""
|
||||
return transform(source_points, self.control_points,
|
||||
self.coefficient)
|
||||
|
||||
transform = __call__
|
||||
Reference in New Issue
Block a user