Tps adaptive mesh refinement (#132)

* adding adaptive mesh estimate for TPS size reduction

* fixing className caps in json

* added some documentation header

* trivial change to retrigger build

* trivial change to retrigger build

* separating recursion from setup, raising max_iter exception

* Update thin_plate_spline.py

trivial change to retrigger build

* another trivial change for rebuild

* added doc string, default pass self to recursion, custom exception

* subclassing custom exception as EstimationError

* removing custom class __str__ method

* removing test of exception __str__

* reverting to simpler estimation error exception

Author: djkapner

renderapi/transform/leaf/thin_plate_spline.py

@@ -3,9 +3,15 @@
 from renderapi.utils import encodeBase64, decodeBase64
 from .transform import Transform
 import scipy.spatial
+import logging
+import sys
 __all__ = ['ThinPlateSplineTransform']
 
 
+logger = logging.getLogger(__name__)
+logger.addHandler(logging.StreamHandler(sys.stdout))
+
+
 class ThinPlateSplineTransform(Transform):
     """
     render-python class that can hold a dataString for
@@ -281,3 +287,127 @@ def dataString(self):
         b64_2 = encodeBase64(blk2)
 
         return '{} {} {}'.format(header, b64_1, b64_2)
+
+    @staticmethod
+    def mesh_refine(
+            new_src,
+            old_src,
+            old_dst,
+            old_tf=None,
+            computeAffine=True,
+            tol=1.0,
+            max_iter=50,
+            nworst=10,
+            niter=0):
+        """recursive kernel for adaptive_mesh_estimate()
+        Parameters
+        ----------
+        new_src : numpy.array
+            Nx2 array of new control source points. Adapts during recursion.
+            Seeded by adaptive_mesh_estimate.
+        old_src : numpy.array
+            Nx2 array of orignal control source points.
+        old_dst : numpy.array
+            Nx2 array of orignal control destination points.
+        old_tf : ThinPlateSplineTransform
+            transform constructed from old_src and old_dst, passed through
+            recursion iterations. Created if None.
+        computeAffine : boolean
+            whether returned transform will have aMtx
+        tol : float
+            in units of pixels, how close should the points match
+        max_iter: int
+            some limit on how many recursive attempts
+        nworst : int
+            per iteration, the nworst matching srcPts will be added
+        niter : int
+            passed through the recursion for stopping criteria
+
+        Returns
+        -------
+        ThinPlateSplineTransform
+        """
+
+        if old_tf is None:
+            old_tf = ThinPlateSplineTransform()
+            old_tf.estimate(old_src, old_dst, computeAffine=computeAffine)
+
+        new_tf = ThinPlateSplineTransform()
+        new_tf.estimate(
+                new_src,
+                old_tf.tform(new_src),
+                computeAffine=computeAffine)
+        new_dst = new_tf.tform(old_src)
+
+        delta = np.linalg.norm(new_dst - old_dst, axis=1)
+        ind = np.argwhere(delta > tol).flatten()
+
+        if ind.size == 0:
+            return new_tf
+
+        if niter == max_iter:
+            raise EstimationError(
+                    "Max number of iterations ({}) reached in"
+                    " ThinPlateSplineTransform.mesh_refine()".format(
+                        max_iter))
+
+        sortind = np.argsort(delta[ind])
+        new_src = np.vstack((new_src, old_src[ind[sortind[0: nworst]]]))
+
+        return ThinPlateSplineTransform.mesh_refine(
+            new_src,
+            old_src,
+            old_dst,
+            old_tf=old_tf,
+            computeAffine=computeAffine,
+            tol=tol,
+            max_iter=max_iter,
+            nworst=nworst,
+            niter=(niter + 1))
+
+    def adaptive_mesh_estimate(
+            self,
+            starting_grid=7,
+            computeAffine=True,
+            tol=1.0,
+            max_iter=50,
+            nworst=10):
+        """method for creating a transform with fewer control points
+        that matches the original transfom within some tolerance.
+        Parameters
+        ----------
+        starting_grid : int
+            estimate will start with an n x n grid
+        computeAffine : boolean
+            whether returned transform will have aMtx
+        tol : float
+            in units of pixels, how close should the points match
+        max_iter: int
+            some limit on how many recursive attempts
+        nworst : int
+            per iteration, the nworst matching srcPts will be added
+
+        Returns
+        -------
+        ThinPlateSplineTransform
+        """
+
+        mn = self.srcPts.min(axis=1)
+        mx = self.srcPts.max(axis=1)
+        xt, yt = np.meshgrid(
+                np.linspace(mn[0], mx[0], starting_grid),
+                np.linspace(mn[1], mx[1], starting_grid))
+        new_src = np.vstack((xt.flatten(), yt.flatten())).transpose()
+        old_src = self.srcPts.transpose()
+        old_dst = self.tform(old_src)
+
+        return ThinPlateSplineTransform.mesh_refine(
+                new_src,
+                old_src,
+                old_dst,
+                old_tf=self,
+                computeAffine=computeAffine,
+                tol=tol,
+                max_iter=max_iter,
+                nworst=nworst,
+                niter=0)

test/test_transform.py

@@ -170,13 +170,13 @@ def test_invert_Affine():
 def test_polynomial_scale():
     p = np.transpose(np.array([[0, 0]]))
     t = renderapi.transform.Polynomial2DTransform(params=p)
-    assert(t.order==0)
-    assert(t.scale==(1.0, 1.0))
-    
+    assert(t.order == 0)
+    assert(t.scale == (1.0, 1.0))
+
     p = np.transpose(np.array([[0, 0], [1.1, 0], [0, 0.9]]))
     t = renderapi.transform.Polynomial2DTransform(params=p)
-    assert(t.order==1)
-    assert(t.scale==(1.1, 0.9))
+    assert(t.order == 1)
+    assert(t.scale == (1.1, 0.9))
 
 
 def test_Polynomial_estimation(use_numpy=False):
@@ -835,12 +835,12 @@ def test_thinplatespline():
     with pytest.raises(renderapi.errors.RenderError):
         s = t2.dataString.split(' ')
         s[1] = str(int(s[1])+1)
-        t3 = renderapi.transform.ThinPlateSplineTransform(
+        renderapi.transform.ThinPlateSplineTransform(
                 dataString=" ".join(s))
     with pytest.raises(renderapi.errors.RenderError):
         s = t2.dataString.split(' ')
         s[2] = str(int(s[2])-4)
-        t3 = renderapi.transform.ThinPlateSplineTransform(
+        renderapi.transform.ThinPlateSplineTransform(
                 dataString=" ".join(s))
 
     x = np.linspace(0, 3840, 10)
@@ -963,3 +963,51 @@ def test_encode64():
     s = renderapi.utils.encodeBase64(x)
     y = renderapi.utils.decodeBase64(s)
     assert(np.all(x == y))
+
+
+def test_adaptive_estimate():
+    with open(rendersettings.TEST_THINPLATESPLINE_FILE, 'r') as f:
+        j = json.load(f)
+
+    tf = renderapi.transform.ThinPlateSplineTransform(
+            dataString=j['dataString'])
+
+    tol = 1.0
+    ntf = tf.adaptive_mesh_estimate(tol=1.0)
+
+    src = ntf.srcPts.transpose()
+    dsta = tf.tform(src)
+    dstb = ntf.tform(src)
+    assert(np.linalg.norm(dsta - dstb, axis=1).max() <= tol)
+
+    src = tf.srcPts.transpose()
+    dsta = tf.tform(src)
+    dstb = ntf.tform(src)
+    nover = np.argwhere(np.linalg.norm(dsta - dstb, axis=1) >= tol).size
+    assert(nover == 0)
+
+    with pytest.raises(renderapi.errors.EstimationError):
+        tf.adaptive_mesh_estimate(max_iter=1)
+
+    # invoke the recursion directly, without passing self
+    mn = tf.srcPts.min(axis=1)
+    mx = tf.srcPts.max(axis=1)
+    xt, yt = np.meshgrid(
+            np.linspace(mn[0], mx[0], 5),
+            np.linspace(mn[1], mx[1], 5))
+    new_src = np.vstack((xt.flatten(), yt.flatten())).transpose()
+    old_src = tf.srcPts.transpose()
+    old_dst = tf.tform(old_src)
+    ntf = tf.mesh_refine(
+        new_src,
+        old_src,
+        old_dst,
+        old_tf=None,
+        computeAffine=True,
+        tol=1.0,
+        max_iter=50,
+        nworst=10,
+        niter=0)
+    dsta = tf.tform(src)
+    dstb = ntf.tform(src)
+    assert(np.linalg.norm(dsta - dstb, axis=1).max() <= tol)