diff --git a/.gitignore b/.gitignore index e8b6601..6d82b88 100644 --- a/.gitignore +++ b/.gitignore @@ -25,4 +25,7 @@ trying_stuff docs/build # Tox -.tox \ No newline at end of file +.tox + +# pytest +.coverage \ No newline at end of file diff --git a/README.rst b/README.rst index 02896d5..d94db72 100644 --- a/README.rst +++ b/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. diff --git a/README_for_docs.rst b/README_for_docs.rst new file mode 100644 index 0000000..4a57d1e --- /dev/null +++ b/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) `, + :meth:`remove_scale(lmk_sets) ` +* Rigid Rotation, Ordinary and Generalized Procrustes alignment: \ + :meth:`rotate(src_sets,tar_sets) `, + :meth:`opa(src_set,tar_set) `, + :meth:`gpa(all_sets) ` +* Thin-plate spline warping: \ + :meth:`tps_warp(X, Y, pts) ` +* Reading from and writing to \*.dta files: \ + :meth:`read_dta(fn) `, + :meth:`write_dta(fn,lmk_sets,names) ` + +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. diff --git a/docs/source/_autosummary/morphops.tps.rst b/docs/source/_autosummary/morphops.tps.rst new file mode 100644 index 0000000..13c8bd9 --- /dev/null +++ b/docs/source/_autosummary/morphops.tps.rst @@ -0,0 +1,10 @@ +morphops.tps +============ + +.. contents:: + :local: + +.. automodule:: morphops.tps + +.. Members +.. ======= \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py index 16568a0..dc75692 100644 --- a/docs/source/conf.py +++ b/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 --------------------------------------------------- diff --git a/docs/source/index.rst b/docs/source/index.rst index 06cf69f..124a6e3 100644 --- a/docs/source/index.rst +++ b/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 =============== diff --git a/morphops/__init__.py b/morphops/__init__.py index b26a12c..2c95e5a 100644 --- a/morphops/__init__.py +++ b/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 diff --git a/morphops/_version.py b/morphops/_version.py index bad32e9..0876fc4 100644 --- a/morphops/_version.py +++ b/morphops/_version.py @@ -1 +1 @@ -__version__ = '0.1.4' \ No newline at end of file +__version__ = '0.1.7' \ No newline at end of file diff --git a/morphops/lmk_util.py b/morphops/lmk_util.py index da1a168..298ce74 100644 --- a/morphops/lmk_util.py +++ b/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 \ No newline at end of file + 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) \ No newline at end of file diff --git a/morphops/tps.py b/morphops/tps.py new file mode 100644 index 0000000..7595546 --- /dev/null +++ b/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) + + diff --git a/setup.py b/setup.py index f51da5f..2910890 100644 --- a/setup.py +++ b/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 ) \ No newline at end of file diff --git a/tests/test_lmk_util.py b/tests/test_lmk_util.py index c41e446..fd12672 100644 --- a/tests/test_lmk_util.py +++ b/tests/test_lmk_util.py @@ -42,6 +42,12 @@ "when X is the unit square + it's pi/4 rotated form, as " "4*[(1-cos(pi/4))^2 + cos(pi/4)^2]")] +distance_matrix_data = [ + (np.zeros((4,2)), np.zeros((5,2)), np.zeros((4,5)), + "when X is (p1,k) zeros and Y is (p2,k) zeros, as (pi,p2) zeros"), + ([[0,0],[0,1]],[[1,1],[1,0]],[[np.sqrt(2),1],[1,np.sqrt(2)]], + "when X and Y are oppo sides of a unit square, as 1s and sqrt(2)s")] + @pytest.mark.parametrize("X, scn", num_lmk_sets_fail_data) def test_num_lmk_sets_fail(X, scn): print("num_lmk_sets should fail -", scn) @@ -88,4 +94,11 @@ def test_ssqd_fail(X, err_msg, scn): def test_ssqd(X, ans, scn): print("ssqd should give the sum of squared differences " "between all pairs of matrices -", scn) - assert np.isclose(lmk_util.ssqd(X), ans) \ No newline at end of file + assert np.isclose(lmk_util.ssqd(X), ans) + +@pytest.mark.parametrize("X, Y, ans, scn", distance_matrix_data) +def test_distance_matrix(X, Y, ans, scn): + print("distance_matrix should give the distance between each lmk in X " + "and Y when -", scn) + assert np.allclose(lmk_util.distance_matrix(X,Y), ans) + diff --git a/tests/test_tps.py b/tests/test_tps.py new file mode 100644 index 0000000..5bc1072 --- /dev/null +++ b/tests/test_tps.py @@ -0,0 +1,74 @@ +import pytest +import numpy as np +import morphops.tps as tps +from .helpers import make_haus, make_ngon, get_2d_rot + +pentagon = make_ngon(5) +pentagon_45 = make_ngon(5, np.pi/4.0) +fourgon = make_ngon(4) +fourgon_45 = make_ngon(4, np.pi/4.0) + +K_matrix_2_data = [(np.zeros((5,2)), None, np.zeros((5,5)), + "X is (p,2) zeros. The result should be (p,p) zeros."), + (np.ones((5,2)), None, np.zeros((5,5)), + "X is (p,2) ones. The result should be (p,p) zeros."), + (fourgon, None, + [[0, 2*np.log(2), 8*np.log(2), 2*np.log(2)], + [2*np.log(2), 0, 2*np.log(2), 8*np.log(2)], + [8*np.log(2), 2*np.log(2), 0, 2*np.log(2)], + [2*np.log(2), 8*np.log(2), 2*np.log(2), 0]], + "X is square of side sqrt(2). Result should be same as in 89 " + "paper, sec E.")] + +K_matrix_3_data = [(np.column_stack((np.ones(4),fourgon)), None, + [[0, np.sqrt(2), 2, np.sqrt(2)], + [np.sqrt(2), 0, np.sqrt(2), 2], + [2, np.sqrt(2), 0, np.sqrt(2)], + [np.sqrt(2), 2, np.sqrt(2), 0]], + "X is square of side sqrt(2). Result should be the distance " + "matrix.")] + +tps_coefs_affine_data = [(fourgon, fourgon + 1, + np.zeros((4,2)), [[1,1],[1,0],[0,1]], "X is a square of side " + "sqrt(2), Y = X + (1,1). Result should be that W is all " + "zeros and A is a row of ones stacked on identity." + ), + (fourgon, fourgon_45, + np.zeros((4,2)),np.row_stack(([0,0],get_2d_rot(-np.pi/4))), + "X is a square of side sqrt(2), Y = XR. Result should be that " + "W is all zeros and A is a row of zeros stacked on R.T." + )] + +tps_warp_affine_data = [(fourgon, fourgon + 1, + fourgon, fourgon + 1, + "X is a square of side sqrt(2), Y = X + (1,1), pts = X. " + "Result should Y."), + (fourgon, fourgon_45, + fourgon, np.dot(fourgon, get_2d_rot(-np.pi/4)), + "X is a square of side sqrt(2), Y = XR, pts = square. Result " + "should be pts*R.T." + )] + +@pytest.mark.parametrize("X, Y, ans, scn", K_matrix_2_data) +def test_K_matrix_2(X, Y, ans, scn): + print("K_matrix should evaluate the rbf at the distance between each " + "lmk in X and Y when -", scn) + assert np.allclose(tps.K_matrix(X,Y), ans) + +@pytest.mark.parametrize("X, Y, ans, scn", K_matrix_3_data) +def test_K_matrix_3(X, Y, ans, scn): + print("K_matrix should evaluate the rbf at the distance between each " + "lmk in X and Y when -", scn) + assert np.allclose(tps.K_matrix(X,Y), ans) + +@pytest.mark.parametrize("X, Y, W_ans, A_ans, scn", tps_coefs_affine_data) +def test_tps_coefs(X, Y, W_ans, A_ans, scn): + print("tps_coefs should evaluate spline weights W, A when -", scn) + W, A = tps.tps_coefs(X, Y) + assert np.allclose(W, W_ans) + assert np.allclose(A, A_ans) + +@pytest.mark.parametrize("X, Y, pts, warped_pts_ans, scn", tps_warp_affine_data) +def test_tps_warp(X, Y, pts, warped_pts_ans, scn): + print("tps_warp should given warped pts when -", scn) + assert np.allclose(tps.tps_warp(X,Y,pts), warped_pts_ans) diff --git a/tox.ini b/tox.ini index 4f3e2a5..b232658 100644 --- a/tox.ini +++ b/tox.ini @@ -1,12 +1,21 @@ [tox] -envlist = py35,py36,python3.7 +envlist = py36,python3.7 [travis] python = - 3.5: py35 3.6: py36 3.7: python3.7 +[testenv:py36] +deps = + pytest + pytest-cov +commands = + pytest -q -s tests + pytest --cov=morphops tests/ + [testenv] -deps = pytest -commands = pytest -q -s tests \ No newline at end of file +deps = + pytest +commands = + pytest -q -s tests \ No newline at end of file