Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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