Tps init

.gitignore

@@ -25,4 +25,7 @@ trying_stuff
 docs/build
 
 # Tox
-.tox
+.tox
+
+# pytest
+.coverage

README.rst

@@ -7,6 +7,13 @@ Welcome to Morphops!
 Morphops implements common operations and algorithms for geometric
 morphometrics, in python 3.
 
+Some high-level operations in the current version are
+
+* Centering, rescaling data
+* Rigid Rotation, Ordinary and Generalized Procrustes alignment
+* Thin-plate spline warping
+* Reading from and writing to \*.dta files
+
 Dependencies
 ------------
 
@@ -17,8 +24,8 @@ Installation
 
 :code:`pip install morphops`
 
-Usage
------
+Usage Examples
+--------------
 
 .. code-block:: python
 
@@ -31,6 +38,10 @@ Usage
    # Perform Generalized Procrustes alignment to align A, B, C.
    # :func:`gpa` is in the procrustes module.
    res = mops.gpa([A,B,C])
-   
+
    # res['X0_ald'] contains the aligned A, B, C.
    # res['X0_mu'] contains the mean of the aligned A, B, C.
+
+   # Create a Thin-plate Spline warp from A to B and warp C.
+   warped_C = mops.tps_warp(A, B, C)
+   # warped_C contains the image of the pts in C under the TPS warp.

README_for_docs.rst

@@ -0,0 +1,62 @@
+.. image:: https://travis-ci.com/vaipatel/morphops.svg?branch=master
+    :target: https://travis-ci.com/vaipatel/morphops
+
+|
+
+.. contents::
+    :local:
+
+|
+
+Welcome to Morphops!
+====================
+
+Morphops implements common operations and algorithms for geometric
+morphometrics, in python 3.
+
+Some high-level operations in the current version are
+
+* Centering, rescaling data: \
+  :meth:`remove_position(lmk_sets) <morphops.procrustes.remove_position>`,
+  :meth:`remove_scale(lmk_sets) <morphops.procrustes.remove_scale>`
+* Rigid Rotation, Ordinary and Generalized Procrustes alignment: \
+  :meth:`rotate(src_sets,tar_sets) <morphops.procrustes.rotate>`,
+  :meth:`opa(src_set,tar_set) <morphops.procrustes.opa>`,
+  :meth:`gpa(all_sets) <morphops.procrustes.gpa>`
+* Thin-plate spline warping: \
+  :meth:`tps_warp(X, Y, pts) <morphops.tps.tps_warp>`
+* Reading from and writing to \*.dta files: \
+  :meth:`read_dta(fn) <morphops.io.read_dta>`,
+  :meth:`write_dta(fn,lmk_sets,names) <morphops.io.write_dta>`
+
+Dependencies
+------------
+
+* numpy
+
+Installation
+------------
+
+:code:`pip install morphops`
+
+Usage Examples
+--------------
+
+.. code-block:: python
+
+   import morphops as mops
+   # Create 3 landmark sets, each having 5 landmarks in 2 dimensions.
+   A = [[0,0],[2,0],[2,2],[1,3],[0,2]]
+   B = [[0.1,-0.1],[2,0],[2.3,1.8],[1,3],[0.4,2]]
+   C = [[-0.1,-0.1],[2.1,0],[2,1.8],[0.9,3.1],[-0.4,2.1]]
+
+   # Perform Generalized Procrustes alignment to align A, B, C.
+   # :func:`gpa` is in the procrustes module.
+   res = mops.gpa([A,B,C])
+
+   # res['X0_ald'] contains the aligned A, B, C.
+   # res['X0_mu'] contains the mean of the aligned A, B, C.
+
+   # Create a Thin-plate Spline warp from A to B and warp C.
+   warped_C = mops.tps_warp(A, B, C)
+   # warped_C contains the image of the pts in C under the TPS warp.

docs/source/_autosummary/morphops.tps.rst

@@ -0,0 +1,10 @@
+morphops.tps
+============
+
+.. contents::
+    :local:
+
+.. automodule:: morphops.tps
+
+.. Members
+.. =======

docs/source/conf.py

@@ -28,7 +28,7 @@
 # The short X.Y version
 version = '0.1'
 # The full version, including alpha/beta/rc tags
-release = '0.1.4'
+release = '0.1.7'
 
 # -- General configuration ---------------------------------------------------
 

docs/source/index.rst

@@ -3,7 +3,7 @@
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
 
-.. include:: ../../README.rst
+.. include:: ../../README_for_docs.rst
 
 morphops module
 ===============

morphops/__init__.py

@@ -13,16 +13,18 @@
 .. autosummary::
     :toctree: _autosummary
 
-    lmk_util
-    io
     procrustes
+    tps
+    io
+    lmk_util
 """
 
 # Not sure if prepending morphops to module is necessary, but makes it easier to import into jupyter for testing.
 from .lmk_util import \
-    transpose, num_coords, num_lmks, num_lmk_sets, ssqd
+    transpose, num_coords, num_lmks, num_lmk_sets, ssqd, distance_matrix
 from .io import \
     MopsFileReadError, MopsFileWriteError, read_dta, write_dta
+from .tps import *
 from .procrustes import \
     get_position, get_scale, remove_position, remove_scale, \
     rotate, opa, gpa

morphops/_version.py

@@ -1 +1 @@
-__version__ = '0.1.4'
+__version__ = '0.1.7'

morphops/lmk_util.py

@@ -75,4 +75,21 @@ def ssqd(X):
     ssq = 0
     for i in np.arange(n_lmk_sets - 1):
         ssq += np.sum(np.square(X[i:] - X[i]))
-    return ssq*1.0/n_lmk_sets
+    return ssq*1.0/n_lmk_sets
+
+def distance_matrix(X,Y):
+    """For (p1,k)-shaped X and (p2,k)-shaped Y, returns the (p1,p2) matrix 
+    where the element at [i,j] is the distance between X[i,:] and Y[j,:].
+
+    This is basically a matrix version of the following code.
+
+    .. code-block:: python
+
+        for i in range(len(X)):
+            for j in range(len(Y)):
+                d[i,j] = np.sqrt(np.sum(np.square(np.array(X[i]) - np.array(Y[j]))))
+    """
+    XX = np.tile(np.sum(np.square(X),axis=1),(len(Y),1)).T
+    YY = np.tile(np.sum(np.square(Y),axis=1),(len(X),1))
+    XY = np.dot(X, transpose(Y))
+    return np.sqrt(XX + YY - 2*XY)

morphops/tps.py

@@ -0,0 +1,190 @@
+"""Provides thin-plate splines related operations and algorithms.
+
+Given two sets of points, the thin-plate spline can interpolate from one to the 
+other in a manner that minimizes the "integral bending norm"[bookstein89]_.
+
+Importantly, it has a remarkable connection to Kendall's shape space in the 
+following way: The non-zero eigenvectors of the bending energy matrix form an 
+orthonormal basis in the tangent space of shape coordinates [bookstein96]_.
+
+References
+----------
+.. [bookstein89] Bookstein, F.L., 1989. Principal warps: Thin-plate splines 
+    and the decomposition of deformations. IEEE Transactions on pattern 
+    analysis and machine intelligence, 11(6), pp.567-585.
+.. [bookstein96] Bookstein, F.L., 1996. Biometrics, biomathematics and the 
+    morphometric synthesis. Bulletin of mathematical biology, 58(2), p.313.
+"""
+
+import numpy as np
+import math
+import morphops.lmk_util as lmk_util
+import warnings
+
+def K_matrix(X, Y=None):
+    """Calculates the upper-right (p,p) submatrix of the (p+k+1,p+k+1)-shaped 
+    L matrix.
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+
+        A (p,k) array of p points in k=2 or k=3 dimensions.
+
+    Y : (m,2) or (m,3) shaped array-like, optional
+        A (m,k) array of p points in k=2 or k=3 dimensions. `Y` must have the 
+        same k as `X`.
+
+        If `Y` is `None`, it is just set to `X`.
+
+    Returns
+    -------
+    K : np.ndarray
+        A (p,p) array where the element at [i,j] is :math:`U(\|X_i - Y_j\|)`. The definition of U depends on k. 
+
+        In particular, if k = 2, then  :math:`U(r) = r^2 \log(r^2)`, else 
+        :math:`U(r) = r`.
+
+        Note that using :math:`\\alpha U(r)` instead of :math:`U(r)` for some constant :math:`\\alpha \in \mathbb{R}` will not matter when warping, as it will become inverted when calculating :math:`L^{-1}`, and get multiplied by itself when calculating the non-uniform part of the warp.
+    """
+    num_coords = lmk_util.num_coords(X)
+    if (num_coords != 2) and (num_coords != 3):
+        raise ValueError("The input matrix must have landmarks with "
+                         "coordinates in either 2 or 3 dimensions.")
+    if Y is None:
+        Y = X
+    r = lmk_util.distance_matrix(X, Y)
+    if (num_coords == 2):
+        r_sqd = np.square(r)
+        # Make a copy of r_sqd where 0->1. This copy will be passed to log.
+        # This way log(1) will be 0 and we wont get NaN and warnings. 
+        r_sqd_cl = np.copy(r_sqd)
+        r_sqd_cl[np.isclose(r_sqd_cl,0)] = 1
+        return np.multiply(r_sqd, np.log(r_sqd_cl))
+    # else num_coords is 3
+    return r
+
+def P_matrix(X):
+    """Makes the minor diagonal submatrix P of the (p+k+1,p+k+1)-shaped L 
+    matrix.
+
+    Basically just stacks a column of 1s before the coordinate columns in `X`.
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p points in k=2 or k=3 dimensions.
+
+    Returns
+    -------
+    P : np.ndarray
+        A (p,k+1) array, which is 1 in the first column, and exactly `X` in the 
+        remaining columns.
+    """
+    ones = np.ones(lmk_util.num_lmks(X))
+    return np.column_stack((ones, X))
+
+def L_matrix(X):
+    """Makes the (p+k+1,p+k+1)-shaped L matrix that gets inverted when 
+    calculating the thin-plate spline "from" `X`.
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p landmarks in k=2 or k=3 dimensions for one specimen.
+
+    Returns
+    -------
+    L : np.ndarray
+        A (p+k+1,p+k+1) array of the form [[K | P][P.T | 0]].
+    """
+    n_coords = lmk_util.num_coords(X)
+    n_lmks = lmk_util.num_lmks(X)
+    K = K_matrix(X)
+    P = P_matrix(X)
+    L = np.zeros((n_lmks + n_coords + 1, n_lmks + n_coords + 1))
+    L[0:n_lmks,0:n_lmks] = K
+    L[0:n_lmks,n_lmks:] = P
+    L[n_lmks:,0:n_lmks] = np.transpose(P)
+    return L
+
+def bending_energy_matrix(X):
+    """Returns the upper right (pxp) submatrix of L^(-1).
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p landmarks in k=2 or k=3 dimensions for one specimen.
+
+    Returns
+    -------
+    L_inv : np.ndarray
+        The upper right (p,p) submatrix of the inverse of the `L_matrix` of `X`.
+    """
+    n_lmks = lmk_util.num_lmks(X)
+    L = L_matrix(X)
+    L_inv = np.linalg.inv(L)
+    return L_inv[0:n_lmks,0:n_lmks]
+
+def tps_coefs(X, Y):
+    """Finds the thin-plate spline coefficients for the thin-plate spline 
+    function that interpolates from X to Y.
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p points in k=2 or k=3 dimensions.
+
+    Y : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p points in k=2 or k=3 dimensions. `Y` must have the 
+        same shape as `X`.
+
+    Returns
+    -------
+    W : np.ndarray
+        A (p,k) array of weights for the non-affine part of the spline.
+
+    A : np.ndarray
+        A (k+1,k) array of weights for the affine part of the spline.
+    """
+    n_coords = lmk_util.num_coords(X)
+    n_lmks = lmk_util.num_lmks(X)
+    Y_0 = np.row_stack((Y, np.zeros((n_coords+1,n_coords))))
+    L = L_matrix(X)
+    Q = np.linalg.solve(L, Y_0)
+    if np.any(np.isnan(Q)):
+        raise ValueError("The result of L_inv*Y contained NaN values.")
+    # return W and A.
+    return Q[0:n_lmks], Q[n_lmks:]
+
+def tps_warp(X, Y, pts):
+    """Maps points `pts` to their image under the thin-plate spline function generated by :func:`tps_coefs` of `X` and `Y`.
+
+    A point :math:`(a,b)` gets mapped to
+    .. math:: f(a,b)
+
+    Parameters
+    ----------
+    X : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p points in k=2 or k=3 dimensions.
+
+    Y : (p,2) or (p,3) shaped array-like
+        A (p,k) array of p points in k=2 or k=3 dimensions. `Y` must have the 
+        same shape as `X`.
+
+    pts : (m,2) or (m,3) shaped array-like, optional
+        A (m,k) array of m points in k=2 or k=3 dimensions. `pts` must have the 
+        same coordinate dimensions k as `X`.
+
+    Returns
+    -------
+    warped_pts : (m,2) or (m,3) shaped array-like, optional
+        A (m,k) array of points corresponding to the image of `pts` under the thin-plate spline produced by `X`, `Y`.
+    """
+    W, A = tps_coefs(X, Y)
+    U = K_matrix(pts, X)
+    P = P_matrix(pts)
+    # The warped pts are the affine part + the non-uniform part
+    return np.matmul(P,A) + np.matmul(U,W)
+
+

setup.py

@@ -24,6 +24,8 @@
         'License :: OSI Approved :: MIT License',
         'Topic :: Scientific/Engineering'
     ],
-    install_requires= ['numpy'],
+    setup_requires= ['numpy >= 1.10.0'],
+    install_requires= ['numpy >= 1.10.0'],
+    python_requires='>=3.5',
     include_package_data=True
 )