[python-dtcwt] 266/497: registration.qtildematrices: operate on transforms

[Ghislain Vaillant]

This is an automated email from the git hooks/post-receive script.

ghisvail-guest pushed a commit to branch debian/sid
in repository python-dtcwt.

commit 5ca079ff91ca4c4558d6aa009b6b0a7cdf75f64e
Author: Rich Wareham <rjw57 at cam.ac.uk>
Date:   Thu Jan 30 11:21:44 2014 +0000

    registration.qtildematrices: operate on transforms
    
    Instead of taking subband tuples directly, refactor qtildematrices to
    take transforms and warp them via an warptransform function. This is an
    effort to tidy up the public API to make it easier to fit OpenCL in
    later. Make warphighpass internal-only.
---
 dtcwt/registration.py | 74 +++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 57 insertions(+), 17 deletions(-)

diff --git a/dtcwt/registration.py b/dtcwt/registration.py
index 8e2443f..aac7dce 100644
--- a/dtcwt/registration.py
+++ b/dtcwt/registration.py
@@ -15,6 +15,7 @@ import itertools
 from six.moves import xrange
 
 import dtcwt
+from dtcwt.backend import backend_numpy
 import dtcwt.sampling
 import dtcwt.utils
 import numpy as np
@@ -34,7 +35,7 @@ __all__ = [
     'solvetransform',
     'velocityfield',
     'warp',
-    'warphighpass',
+    'warptransform',
 ]
 
 # Horizontal and vertical expected phase shifts for each subband
@@ -107,11 +108,27 @@ def confidence(sb1, sb2):
 Q_TRIU_INDICES = list(zip(*np.triu_indices(6)))
 Q_TRIU_FLAT_INDICES = np.ravel_multi_index(np.triu_indices(6), (6,6))
 
-def qtildematrices(Yh1, Yh2, levels):
+def qtildematrices(t_ref, t_target, levels):
     r"""
     Compute :math:`\tilde{Q}` matrices for given levels.
 
+    :param t_ref: the transformed reference image
+    :param t_target: the transformed target image
+    :param levels: a sequence of indices specifying which levels to examine
+    :returns: a tuple of :math:`\tilde{Q}` matrices for each index in *levels*
+
+    Both *t_ref* and *t_target* should be
+    :py:class:`dtcwt.backend.base.TransformDomainSignal`-compatible objects.
+    Indices in *levels* are 0-based.
+
+    The returned matrices are NxMx27 where NxM is the shape of the
+    corresponding level's highpass subbands.
+
     """
+    # Extract subbands
+    Yh1 = t_ref.subbands
+    Yh2 = t_target.subbands
+
     # A list of arrays of \tilde{Q} matrices for each level
     Qt_mats = []
 
@@ -303,6 +320,33 @@ def affinewarp(Yh, a, method=None):
     vxs, vys = affinevelocityfield(a, xs, ys)
     return normsample(Yh, xs+vxs, ys+vys, method=method)
 
+def warptransform(t, avecs, levels, method='bilinear'):
+    """
+    Return a warped version of a transformed image acting only on specified levels.
+
+    :param t: a transformed image
+    :param avecs: an array of affine distortion parameters
+    :param levels: a sequence of 0-based indices specifying which levels to act on
+
+    *t* should be a
+    :pyclass:`dtcwt.backend.base.TransformDomainSignal`-compatible instance.
+
+    .. note::
+
+        This function will clone the transform *t* but it is a shallow clone
+        where possible. Only the levels specified in *levels* will be
+        deep-copied and warped.
+
+    """
+    warped_subbands = list(t.subbands)
+
+    # Warp specified levels
+    for l in levels:
+        warped_subbands[l] = warphighpass(warped_subbands[l], avecs, method=method)
+
+    # Clone the transform
+    return backend_numpy.TransformDomainSignal(t.lowpass, tuple(warped_subbands), t.scales)
+
 def estimatereg(reference, target):
     """
     Estimate registration from reference image to target.
@@ -321,19 +365,16 @@ def estimatereg(reference, target):
     this function into a velocity field.
 
     """
-    # Extract highpass
-    reference_h = reference.subbands
-    target_h = target.subbands
-
-    # Make a copy of reference_h for warping
-    warped_h = list(reference_h)
+    # Extract number of levels and shape of level 3 subband
+    nlevels = len(reference.subbands)
+    avecs_shape = reference.subbands[2].shape[:2] + (6,)
 
     # Initialise matrix of 'a' vectors
-    avecs = np.zeros((warped_h[2].shape[0], warped_h[2].shape[1], 6))
+    avecs = np.zeros(avecs_shape)
 
     # Compute initial global transform
-    levels = list(x for x in xrange(len(warped_h)-1, len(warped_h)-3, -1) if x>=0)
-    Qt_mats = list(np.sum(np.sum(x, axis=0), axis=0) for x in qtildematrices(warped_h, target_h, levels))
+    levels = list(x for x in xrange(nlevels-1, nlevels-3, -1) if x>=0)
+    Qt_mats = list(np.sum(np.sum(x, axis=0), axis=0) for x in qtildematrices(reference, target, levels))
     Qt = np.sum(Qt_mats, axis=0)
 
     a = solvetransform(Qt)
@@ -341,18 +382,17 @@ def estimatereg(reference, target):
         avecs[:,:,idx] = a[idx]
 
     # Refine estimate
-    for it in xrange(2 * (len(warped_h)-3) - 1):
-        s, e = len(warped_h) - 2, len(warped_h) - 2 - (it+1)//2
-        levels = list(x for x in xrange(s, e-1, -1) if x>=2 and x<len(warped_h))
+    for it in xrange(2 * (nlevels-3) - 1):
+        s, e = nlevels - 2, nlevels - 2 - (it+1)//2
+        levels = list(x for x in xrange(s, e-1, -1) if x>=2 and x<nlevels)
         if len(levels) == 0:
             continue
 
         # Warp the levels we'll be looking at with the current best-guess transform
-        for l in levels:
-            warped_h[l] = warphighpass(reference_h[l], avecs, method='bilinear')
+        warped = warptransform(reference, avecs, levels, method='bilinear')
 
         qts = np.sum(list(dtcwt.sampling.rescale(_boxfilter(x, 3), avecs.shape[:2], method='bilinear')
-                          for x in qtildematrices(warped_h, target_h, levels)),
+                          for x in qtildematrices(warped, target, levels)),
                      axis=0)
 
         avecs += solvetransform(qts)

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/python-dtcwt.git