Deformation fields with thin-plates

test_thin_plates.py

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()
  

thin_plates.py

__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