Deformation fields with thin-plates
import unittest | |
import numpy | |
from numpy.testing import * | |
set_package_path() | |
from transformation import * | |
restore_path() | |
class TestDeformationCreation(NumpyTestCase): | |
def test_1DCreation(self): | |
size = (10,) | |
points = numpy.array(((2,), (8, ))) | |
displacements = [0, (-4,)] | |
field = denseDeformationFieldFromSparse(size, points, displacements) | |
for point,displacement in zip(points, displacements): | |
assert_almost_equal(field[tuple(point.tolist())], displacement) | |
def test_2DCreation(self): | |
size = (11, 11) | |
points = numpy.array(((5., 6.), (5., 4.), (6., 5.), (4., 5.))) | |
displacements = [(0, 1), (0, 1) , (0, -1), (0, -1)] | |
field = denseDeformationFieldFromSparse(size, points, displacements) | |
for point,displacement in zip(points, displacements): | |
assert_almost_equal(field[tuple(point.tolist())], displacement) | |
def test_2DCreation_bis(self): | |
size = (10, 10) | |
points = numpy.array(((5., 5.), (2., 8.), (8., 2.))) | |
displacements = [0, (5, -4), (0, 2)] | |
field = denseDeformationFieldFromSparse(size, points, displacements) | |
for point,displacement in zip(points, displacements): | |
assert_almost_equal(field[tuple(point.tolist())], displacement) | |
def test_3DCreation(self): | |
size = (10, 10, 10) | |
points = numpy.array(((5., 5., 5.), (2., 8., 2.), (8., 2., 8.), (2., 2., 8.))) | |
displacements = [0, (5, -4, 3), (0, 2, 0), -2] | |
field = denseDeformationFieldFromSparse(size, points, displacements) | |
for point,displacement in zip(points, displacements): | |
assert_almost_equal(field[tuple(point.tolist())], displacement) | |
if __name__ == "__main__": | |
unittest.main() | |
__all__ = ['denseDeformationFieldFromSparse'] | |
import numpy | |
small_ = 1e-100 | |
Ufunction = [None, (lambda x: numpy.abs(x)**3), (lambda x: x**2 * numpy.where(x<small_, 0, numpy.log(x**2))), (lambda x: numpy.abs(x))] | |
def dist2hd(x,y): | |
""" | |
Generate a 'coordinate' of the solution at a time | |
""" | |
d = numpy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype) | |
for i in xrange(x.shape[1]): | |
diff2 = x[:,i,None] - y[:,i] | |
diff2 **= 2 | |
d += diff2 | |
numpy.sqrt(d,d) | |
return d | |
def denseDeformationFieldFromSparse(size, points, displacements): | |
""" | |
Creates a smooth deformation field based on : | |
- its size | |
- the points where the deformation is known | |
- the displacement at those points | |
size must be of the same size than the number of coordinates of points. | |
This function uses Thin Plates to interpolate the field (cf Bookstein 1989 doi:10.1109/34.24792) | |
""" | |
#Computation of the coefficients that will be stored in function | |
points=numpy.asfarray(points) | |
distances = dist2hd(points, points) | |
totalsize = sum(points.shape) | |
L = numpy.zeros((totalsize+1, totalsize+1)) | |
L[:len(points),:len(points)] = Ufunction[len(size)](distances) | |
L[len(points),:len(points)] = 1 | |
L[:len(points),len(points)] = 1 | |
L[len(points)+1:,:len(points)] = points.T | |
L[:len(points),len(points)+1:] = points | |
Y = numpy.zeros((points.shape[0]+len(size)+1, points.shape[1])) | |
for i in range(0, len(displacements)): | |
Y[i] = displacements[i] | |
function = numpy.linalg.solve(L, Y) | |
#Creation of the dense deformation field | |
field = numpy.empty(size+(len(size),)) | |
for indice in range(len(size)): | |
coeffs = numpy.arange(0, size[indice], dtype=int) | |
newaxises = [1] * (len(size)-indice-1) | |
coeffs.shape = [coeffs.shape[0]]+newaxises | |
field[...,indice] = coeffs | |
view = field.reshape(-1, len(size)) | |
distances = dist2hd(view, points) | |
distances = Ufunction[len(size)](distances) | |
view[:] = function[len(points)] + numpy.dot(distances, function[:len(points)]) + numpy.dot(view, function[len(points)+1:]) | |
return field | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment