Add fit estimate to tps (#123)

* clunky, but working copy of javafor computing dMtxDat

* working with tests, still clunky

* added stand-alone test, not dependent on pre-computed json from java client

Author: djkapner

renderapi/transform/leaf/thin_plate_spline.py

@@ -180,6 +180,101 @@ def inverse_tform(
             newpts.append(npt)
         return np.array(newpts)
 
+    @staticmethod
+    def fit(A, B, computeAffine=True):
+        """function to fit this transform given the corresponding sets of points A & B
+
+        Parameters
+        ----------
+        A : numpy.array
+            a Nx2 matrix of source points
+        B : numpy.array
+            a Nx2 matrix of destination points
+
+        Returns
+        -------
+        dMatrix : numpy.array
+            ndims x nLm
+        aMatrix : numpy.array
+            ndims x ndims, affine matrix
+        bVector : numpy.array
+            ndims x 1, translation vector
+        """
+
+        if not all([A.shape[0] == B.shape[0], A.shape[1] == B.shape[1] == 2]):
+            raise EstimationError(
+                'shape mismatch! A shape: {}, B shape {}'.format(
+                    A.shape, B.shape))
+
+        # build displacements
+        ndims = B.shape[1]
+        nLm = B.shape[0]
+        y = (B - A).flatten()
+
+        # compute K
+        # tempting to matricize this, but, nLm x nLm can get big
+        # settle for vectorize
+        kMatrix = np.zeros((ndims * nLm, ndims * nLm))
+        for i in range(nLm):
+            r = np.linalg.norm(A[i, :] - A, axis=1)
+            nrm = np.zeros_like(r)
+            ind = np.argwhere(r > 1e-8)
+            nrm[ind] = r[ind] * r[ind] * np.log(r[ind])
+            kMatrix[i * ndims, 0::2] = nrm
+            kMatrix[(i * ndims + 1)::2, 1::2] = nrm
+
+        # compute L
+        lMatrix = kMatrix
+        if computeAffine:
+            pMatrix = np.tile(np.eye(ndims), (nLm, ndims + 1))
+            for d in range(ndims):
+                pMatrix[0::2, d*ndims] = A[:, d]
+                pMatrix[1::2, d*ndims + 1] = A[:, d]
+            lMatrix = np.zeros(
+                    (ndims * (nLm + ndims + 1), ndims * (nLm + ndims + 1)))
+            lMatrix[
+                    0: pMatrix.shape[0],
+                    kMatrix.shape[1]: kMatrix.shape[1] + pMatrix.shape[1]] = \
+                pMatrix
+            pMatrix = np.transpose(pMatrix)
+            lMatrix[
+                    kMatrix.shape[0]: kMatrix.shape[0] + pMatrix.shape[0],
+                    0: pMatrix.shape[1]] = pMatrix
+            lMatrix[0: ndims * nLm, 0: ndims * nLm] = kMatrix
+            y = np.append(y, np.zeros(ndims * (ndims + 1)))
+
+        wMatrix = np.linalg.solve(lMatrix, y)
+
+        dMatrix = np.reshape(wMatrix[0: ndims * nLm], (ndims, nLm), order='F')
+        aMatrix = None
+        bVector = None
+        if computeAffine:
+            aMatrix = np.reshape(
+                    wMatrix[ndims * nLm: ndims * nLm + ndims * ndims],
+                    (ndims, ndims),
+                    order='F')
+            bVector = wMatrix[ndims * nLm + ndims * ndims:]
+
+        return dMatrix, aMatrix, bVector
+
+    def estimate(self, A, B, computeAffine=True):
+        """method for setting this transformation with the best fit
+        given the corresponding points A,B
+        Parameters
+        ----------
+        A : numpy.array
+            a Nx2 matrix of source points
+        B : numpy.array
+            a Nx2 matrix of destination points
+        computeAffine: boolean
+            whether to include an affine computation
+        """
+
+        self.dMtxDat, self.aMtx, self.bVec = self.fit(
+                A, B, computeAffine=computeAffine)
+        (self.nLm, self.ndims) = B.shape
+        self.srcPts = np.transpose(A)
+
     @property
     def dataString(self):
         header = 'ThinPlateSplineR2LogR {} {}'.format(self.ndims, self.nLm)

test/rendersettings.py

@@ -9,7 +9,7 @@
     'owner': 'renderowner',
     'project': 'renderproject',
     'client_scripts': '/path/to/client_scripts'
-    }
+        }
 
 DEFAULT_RENDER_CLIENT = dict(DEFAULT_RENDER, **{
     'client_script': '/path/to/client_scripts/run_ws_client.sh',
@@ -21,7 +21,7 @@
     'RENDER_OWNER': DEFAULT_RENDER['owner'],
     'RENDER_PROJECT': DEFAULT_RENDER['project'],
     'RENDER_CLIENT_SCRIPTS': DEFAULT_RENDER['client_scripts']
-    }
+        }
 
 DEFAULT_RENDER_CLIENT_ENVIRONMENT_VARIABLES = dict(
     DEFAULT_RENDER_ENVIRONMENT_VARIABLES, **{
@@ -33,8 +33,12 @@
 
 TEST_TILESPECS_FILE = os.path.join(TEST_FILES_DIR, 'tilespecs.json')
 
-TEST_THINPLATESPLINE_FILE = os.path.join(TEST_FILES_DIR,
-                                         'thin_plate_spline.json')
+TEST_THINPLATESPLINE_FILE = os.path.join(
+        TEST_FILES_DIR,
+        'thin_plate_spline.json')
+TEST_THINPLATESPLINEAFFINE_FILE = os.path.join(
+        TEST_FILES_DIR,
+        'thin_plate_spline_affine.json')
 
 INTERPOLATED_TRANSFORM_TILESPEC = os.path.join(
     TEST_FILES_DIR, 'tilespec_interpolated.json')

test/test_files/thin_plate_spline_affine.json

@@ -0,0 +1,5 @@
+{
+  "type" : "leaf",
+  "className" : "mpicbg.trakem2.transform.ThinPlateSplineTransform",
+  "dataString" : "ThinPlateSplineR2LogR 2 100 @P4R64SLPrGq8kudjBEERFDykx/MHaYpjP2iTdKla/W/AWp5BXNVtaMBXQhhumRdE eJx9VH1UjHkULps6ks+S1IqOWkQpW9k4Xu3aPUeFxIxEUdmtlTgVq0+bmSaqpVCRYWoxanvLUFrK8qYimzmkqCWl8bErsaeRb0527/x+ve85c+3Z+ec59/bc5zy/e583AwP9n1ciy5YaOUwV6kyWZbf+HSXUuYgvh7/XDvER6p+hLkxaJdRKxC8Wi8UhsmNCzSK+yuA/f7yvj3DAn1DnIv6AP6Ee8CfUSsQf8CfULOKr/n9PH/nCmIv4vD++5v3xtRLxeX98zSK+Sp+P7/bRnrAvjLw/vub98bUS8Xl/fM0ivkqfL/jDd8J7wT74GiPvj6+ViM/742sW8VX6fMEfzg2+E94L9sHXGJWIz/vjaxbxVfp8wR/OMc4NvhPeC/ahRDweeX98zSK+Sp8v+MPfFc4xzg2+E94L9sHXGFnEV+nzBX/4O8ffFc4xzg2+E94L9sEiHo8qfb7gD/3fwd85/q5wjnFu8J3wXrAPoUbInDRbmS9yDOZKG/peG/YXcNK6ZwFmtm7c7rqQ9L7wOdyu+mh3qXk3YKjMdZsX9IMdzarKOGm9lXGo7C3MPWp0C6oGnRUvG236mZxV0fcmKGO5/Nm/jOc8zJmcs/3q8zfdmO2SpwVlRWrmpwW3mGcWpoCaMvtQF+i/+GrFowfA6321cpotzMktPH3DQWed9rvTrcz+COtI81RLTpFaVd8U/RsTnXvpfYN/IRd7LPvJ1e9TuIhp2Qb3ZH2A6+6rk5yhL9k7UbsdeDvT14+8DXMlnhXLI0HH4sozwwgm762ka8LcGk7+IO2a7/FgJlnatqSjOJTzLhiUvTU2kfGLE7VHdY4D7J4bJcrhvBVZ1W3nKoGnGO1zRg5zW47WvksGneQJ85vvMHkFL5r67vlx8vphnTuHHGYSwpTDPRPVXKBoQ3fTLAfGbprGasySU4yds/Pwyph50HcbX/bSGniRd2c4tXDyug8nJy8TgU6vw9rRIVSvPZjqbe4nel+0hxO93xM8iJ5vxWCil6UOJXolylSi5zKpiOr5Xqd6Trb0vYP30/cGBJL3dpYGkPfuan1D3rtBWkbeu/7QbvLejsqZ5L2+K03Ie+tS/iDv/drImN6jfAS9x2Jnco+L5iJyD7V9MrnHg/mL6D2sbOk9ao3pPdKayT3Kn04n99Aeiqd5SUqhedl4muSFW9NJ8lL6IZPmZfQkmpfKQJqXGzdIXgJbltG8SCJJXsJnLaN5dmqlee7R0Dwfy6B5XhBE8/xlB82zVwPNs+YHkueQC440z5vH0jybx+j0TFQ/toD+xsv3vRLBt1VVS10x7GPbJCNtOez5uU37egVF732kP9jahPIaUumciyHVkY4B/V7toahR4FvZZf3wLOyj8ta56puwZ+mI8MO6u5vaqFoZijtCaL8mk/LudNO5G3KqE+8A73utmHG3CPbxbsPjmHjYc3XFUI9CuF/j9fx5NZCL/QVrehIJrm4bR/oHZuwhPNP0q3RunTXVOVWm28dn5Seu6fY8e+2sBN39ikdG+EIusl2VDrm6vIm0cwwJ9sRYkn7FRcobnuNP5yyXU52E97BniWLq3G/gfur3l+3EkIv6vVODNjN+W2SyKxWvGLspigsLJxYSdP02lvRbaw8Q3vSwE3TO3onqaP4CXPJppuM26D/tXfSwCXju/ilB8YxffFWOYXczY+caHZ61NIGg7LAd6ZsF5BGepNKNzC3es4noZHyw0uVi5yKfJ7q8FR3J6NLluK3LyVj3fYwtHTZZ992V7ziYSlCS5kn658WTCe+uhQedc59OdfwzIEejLjDiX2Gvz0XO8XGw5/wWhwYx7L1kb9ZjJ7hD3GKTkmaK6W9If7f/Lco7uJDOxZlSnejVcOeepgN/7oO7K5Z2yM5ADo7bX3LfBLlIDCsqGAo5GRRrVD2e4tI+0i8OHUp4FzX/0LlNNlSn6xPIYeD2jXVvIJdRQ45OuQQ5tUye+XkN5HbrneC8PMix9nZSXhjFgym0/0pDeK79OWTuyHMXqtNo8S+ZhThG"
+}

test/test_transform.py

@@ -7,6 +7,7 @@
 import pytest
 
 EPSILON = 0.0000000001
+EPSILON2 = 0.000000001
 
 
 def cross_py23_reload(module):
@@ -804,6 +805,104 @@ def test_thinplatespline():
     dst_pts = t.tform(src_pts)
     assert(dst_pts.shape == src_pts.shape)
 
+    # check load from json
+    jt = t.to_dict()
+    nt = renderapi.transform.ThinPlateSplineTransform(
+            json=jt)
+    assert nt == t
+
+
+def test_thinplatespline_apply():
+    # tests some copied behavior from trakem2
+    j = json.load(open(rendersettings.TEST_THINPLATESPLINE_FILE, 'r'))
+    t = renderapi.transform.ThinPlateSplineTransform(
+            dataString=j['dataString'])
+    del t.dMtxDat
+    pt = np.array([1.234, 45.678])
+    npt = t.apply(pt)
+    assert np.all(pt == npt)
+
+
+def estimate_test(jpath, computeAffine=True):
+    # test that the estimate method can produce same results
+    with open(jpath, 'r') as f:
+        j = json.load(f)
+
+    t = renderapi.transform.ThinPlateSplineTransform(
+            dataString=j['dataString'])
+    # exact points
+    src1 = np.transpose(t.srcPts)
+    # some in-between points
+    x = np.linspace(
+            src1[:, 0].min(),
+            src1[:, 0].max(),
+            int(np.sqrt(t.nLm)) * 3)
+    y = np.linspace(
+            src1[:, 1].min(),
+            src1[:, 1].max(),
+            int(np.sqrt(t.nLm)) * 3)
+    xt, yt = np.meshgrid(x, y)
+    src2 = np.transpose(np.vstack((xt.flatten(), yt.flatten())))
+    dst1_a = t.tform(src1)
+    dst2_a = t.tform(src2)
+
+    # estimate
+    t.estimate(src1, dst1_a, computeAffine=computeAffine)
+    dst1_b = t.tform(src1)
+    dst2_b = t.tform(src2)
+    delta1 = np.linalg.norm(dst1_b - dst1_a, axis=1)
+    delta2 = np.linalg.norm(dst2_b - dst2_a, axis=1)
+
+    assert delta1.max() < EPSILON2
+    assert delta2.max() < EPSILON2
+
+    with pytest.raises(renderapi.errors.EstimationError):
+        t.estimate(src1, dst1_a[1:, :])
+
+
+def test_thinplatespline_estimate():
+    estimate_test(
+            rendersettings.TEST_THINPLATESPLINE_FILE,
+            computeAffine=False)
+    estimate_test(
+            rendersettings.TEST_THINPLATESPLINEAFFINE_FILE,
+            computeAffine=True)
+    thinplate_estimate_nojson(computeAffine=True)
+    thinplate_estimate_nojson(computeAffine=False)
+
+
+def thinplate_estimate_nojson(computeAffine=True):
+    # an estimate test that does not depend on pre-computed json
+    x = np.linspace(0, 1000, 40)
+    xt, yt = np.meshgrid(x, x)
+    src = np.transpose(np.vstack((xt.flatten(), yt.flatten())))
+
+    def poly2d(src):
+        dst = np.zeros_like(src)
+        dx = (src[:, 0] - src[:, 0].mean()) / (src[:, 0].ptp())
+        dy = (src[:, 1] - src[:, 1].mean()) / (src[:, 1].ptp())
+        dst[:, 0] = src[:, 0] + 0.5 * dx * dy + 7.0 * dx * dy * dy
+        dst[:, 1] = src[:, 1] + 0.7 * dy * dy - 4.0 * dx * dx * dy
+        a = np.array([[1.1, -0.04], [-0.02, 0.95]])
+        b = np.array([3.0, 4.0])
+        dst = np.transpose(a.dot(np.transpose(dst))) + b
+        return dst
+
+    dst = poly2d(src)
+    t = renderapi.transform.ThinPlateSplineTransform()
+    t.estimate(src, dst, computeAffine=computeAffine)
+
+    x = np.linspace(0, 1000, 120)
+    xt, yt = np.meshgrid(x, x)
+    test_src = np.transpose(np.vstack((xt.flatten(), yt.flatten())))
+    p_dst = poly2d(test_src)
+    t_dst = t.tform(test_src)
+
+    delta = np.linalg.norm(p_dst - t_dst, axis=1)
+    # can it match the polynomial within a pixel?
+    # for low-N, it won't
+    assert delta.max() < 1.0
+
 
 def test_encode64():
     # case for Stephan's '@' character