[python-dtcwt] 337/497: fixed 3d DTCWT numpy backend to match matlab implementation

Ghislain Vaillant ghisvail-guest at moszumanska.debian.org
Tue Jul 21 18:06:23 UTC 2015


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 99ffbc1472889b9d01bd88b681bc1c1cccb68efb
Author: timseries <timothy.daniel.roberts at gmail.com>
Date:   Thu Feb 6 21:45:09 2014 +0000

    fixed 3d DTCWT numpy backend to match matlab implementation
---
 dtcwt/backend/backend_numpy/transform3d.py |  79 ++++++++++-------------------
 dtcwt/transform3d.py                       |  43 ++++++----------
 tests/qbgn.mat                             | Bin 1193947 -> 826597 bytes
 tests/testagainstmatlab.py                 |   4 +-
 tests/verification.npz                     | Bin 4364786 -> 4510508 bytes
 5 files changed, 46 insertions(+), 80 deletions(-)

diff --git a/dtcwt/backend/backend_numpy/transform3d.py b/dtcwt/backend/backend_numpy/transform3d.py
index 2e5da6c..adf0454 100644
--- a/dtcwt/backend/backend_numpy/transform3d.py
+++ b/dtcwt/backend/backend_numpy/transform3d.py
@@ -97,7 +97,6 @@ class Transform3d(Transform3dBase):
             h0o, g0o, h1o, g1o, h2o, g2o = self.biort
         else:
             raise ValueError('Biort wavelet must have 6 or 4 components.')
-        
         # If qshift has 12 elements instead of 8, then it's a modified
         # rotationally symmetric wavelet
         # FIXME: there's probably a nicer way to do this
@@ -125,16 +124,14 @@ class Transform3d(Transform3dBase):
             # Transform
             if level == 0 and discard_level_1:
                 Yl = _level1_xfm_no_highpass(Yl, h0o, h1o, self.ext_mode)
-                if include_scale:
-                    Yscale[0] = Yl
             elif level == 0 and not discard_level_1:
                 Yl, Yh[level] = _level1_xfm(Yl, h0o, h1o, self.ext_mode)
-                if include_scale:
-                    Yscale[0] = Yl
             else:
                 Yl, Yh[level] = _level2_xfm(Yl, h0a, h0b, h1a, h1b, self.ext_mode)
-                if include_scale:
-                    Yscale[level] = Yl
+            if include_scale:
+                Yscale[level] = Yl.copy()
+        
+                #Yh[nlevels+1]=1 #to throw an error for debugging in nose
         if include_scale:
             return TransformDomainSignal(Yl, tuple(Yh), tuple(Yscale))
         else: 
@@ -272,7 +269,6 @@ def _level1_xfm(X, h0o, h1o, ext_mode):
     for f in xrange(work.shape[1] >> 1):
         # extract slice
         y = work[s0a, f, x2a].T
-
         # Do odd top-level filters on 3rd dim. The order here is important
         # since the second filtering will modify the elements of y as well
         # since y is merely a view onto work.
@@ -288,6 +284,8 @@ def _level1_xfm(X, h0o, h1o, ext_mode):
         # Do odd top-level filters on columns.
         work[s0a, :, f] = _BACKEND.colfilter(y2, h0o)
         work[s0b, :, f] = _BACKEND.colfilter(y2, h1o)
+        #if f==2:
+        #work[:,:,work.shape[2]+1]=1 #to throw an error so we can inspect y in the unit test
 
     # Return appropriate slices of output
     return (
@@ -299,7 +297,7 @@ def _level1_xfm(X, h0o, h1o, ext_mode):
             cube2c(work[x0a, x1a, x2b]),    # LLH
             cube2c(work[x0a, x1b, x2b]),    # HLH
             cube2c(work[x0b, x1a, x2b]),    # LHH
-            cube2c(work[x0b, x1b, x2b]),    # HLH
+            cube2c(work[x0b, x1b, x2b]),    # HHH
             ), axis=3)
         )
 
@@ -567,36 +565,14 @@ def cube2c(y):
     # IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 1, JANUARY 2012
     # eqs. (6) to (9)
 
-    A = y[1::2, 1::2, 1::2]
-    B = y[1::2, 1::2, 0::2]
-    C = y[1::2, 0::2, 1::2]
-    D = y[1::2, 0::2, 0::2]
-    E = y[0::2, 1::2, 1::2]
-    F = y[0::2, 1::2, 0::2]
-    G = y[0::2, 0::2, 1::2]
-    H = y[0::2, 0::2, 0::2]
-
-    # TODO: check if the above should be the below and, if so, fix c2cube
-    #
-    #A = y[0::2, 0::2, 0::2]
-    #B = y[0::2, 0::2, 1::2]
-    #C = y[0::2, 1::2, 0::2]
-    #D = y[0::2, 1::2, 1::2]
-    #E = y[1::2, 0::2, 0::2]
-    #F = y[1::2, 0::2, 1::2]
-    #G = y[1::2, 1::2, 0::2]
-    #H = y[1::2, 1::2, 1::2]
-
-    # TODO: check if the above should be the below and, if so, fix c2cube
-    #
-    #A = y[0::2, 0::2, 0::2]
-    #B = y[0::2, 1::2, 0::2]
-    #C = y[1::2, 0::2, 0::2]
-    #D = y[1::2, 1::2, 0::2]
-    #E = y[0::2, 0::2, 1::2]
-    #F = y[0::2, 1::2, 1::2]
-    #G = y[1::2, 0::2, 1::2]
-    #H = y[1::2, 1::2, 1::2]
+    A = y[0::2, 0::2, 0::2]
+    B = y[0::2, 1::2, 0::2]
+    C = y[1::2, 0::2, 0::2]
+    D = y[1::2, 1::2, 0::2]
+    E = y[0::2, 0::2, 1::2]
+    F = y[0::2, 1::2, 1::2]
+    G = y[1::2, 0::2, 1::2]
+    H = y[1::2, 1::2, 1::2]
 
     # Combine to form subbands
     p = ( A-G-D-F) * j2[0] + ( B-H+C+E) * j2[1]
@@ -604,6 +580,7 @@ def cube2c(y):
     r = ( A+G+D-F) * j2[0] + ( B+H-C+E) * j2[1]
     s = ( A+G-D+F) * j2[0] + (-B-H-C+E) * j2[1]
 
+    #j2[2]=1 #to throw an error
     # Form the 2 subbands in z.
     z = np.concatenate((
         p[:,:,:,np.newaxis],
@@ -635,22 +612,22 @@ def c2cube(z):
     r = z[:,:,:,2]
     s = z[:,:,:,3]
 
-    pr, pi = p.real, p.imag
-    qr, qi = q.real, q.imag
-    rr, ri = r.real, r.imag
-    sr, si = s.real, s.imag
+    pr, pi = p.real, p.imag #A,E
+    qr, qi = q.real, q.imag #B,F
+    rr, ri = r.real, r.imag #C,G
+    sr, si = s.real, s.imag #D,H
 
     y = np.zeros(np.asanyarray(z.shape[:3])*2, dtype=z.real.dtype)
 
-    y[1::2, 1::2, 1::2] = ( pr+qr+rr+sr)
-    y[0::2, 0::2, 1::2] = (-pr-qr+rr+sr)
-    y[1::2, 0::2, 0::2] = (-pr+qr+rr-sr)
-    y[0::2, 1::2, 0::2] = (-pr+qr-rr+sr)
+    y[0::2, 0::2, 0::2] = ( pr+qr+rr+sr)
+    y[1::2, 0::2, 1::2] = (-pr-qr+rr+sr)
+    y[1::2, 1::2, 0::2] = (-pr+qr+rr-sr)
+    y[0::2, 1::2, 1::2] = (-pr+qr-rr+sr)
 
-    y[1::2, 1::2, 0::2] = ( pi-qi+ri-si)
-    y[0::2, 0::2, 0::2] = (-pi+qi+ri-si)
-    y[1::2, 0::2, 1::2] = ( pi+qi-ri-si)
-    y[0::2, 1::2, 1::2] = ( pi+qi+ri+si)
+    y[0::2, 1::2, 0::2] = ( pi-qi+ri-si)
+    y[1::2, 1::2, 1::2] = (-pi+qi+ri-si)
+    y[1::2, 0::2, 0::2] = ( pi+qi-ri-si)
+    y[0::2, 0::2, 1::2] = ( pi+qi+ri+si)
 
     return y * scale
 
diff --git a/dtcwt/transform3d.py b/dtcwt/transform3d.py
index 63f6651..4af88ba 100644
--- a/dtcwt/transform3d.py
+++ b/dtcwt/transform3d.py
@@ -508,25 +508,14 @@ def cube2c(y):
     # IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 21, NO. 1, JANUARY 2012
     # eqs. (6) to (9)
 
-    A = y[1::2, 1::2, 1::2]
-    B = y[1::2, 1::2, 0::2]
-    C = y[1::2, 0::2, 1::2]
-    D = y[1::2, 0::2, 0::2]
-    E = y[0::2, 1::2, 1::2]
-    F = y[0::2, 1::2, 0::2]
-    G = y[0::2, 0::2, 1::2]
-    H = y[0::2, 0::2, 0::2]
-
-    # TODO: check if the above should be the below and, if so, fix c2cube
-    #
-    # A = y[0::2, 0::2, 0::2]
-    # B = y[0::2, 0::2, 1::2]
-    # C = y[0::2, 1::2, 0::2]
-    # D = y[0::2, 1::2, 1::2]
-    # E = y[1::2, 0::2, 0::2]
-    # F = y[1::2, 0::2, 1::2]
-    # G = y[1::2, 1::2, 0::2]
-    # H = y[1::2, 1::2, 1::2]
+    A = y[0::2, 0::2, 0::2]
+    B = y[0::2, 1::2, 0::2]
+    C = y[1::2, 0::2, 0::2]
+    D = y[1::2, 1::2, 0::2]
+    E = y[0::2, 0::2, 1::2]
+    F = y[0::2, 1::2, 1::2]
+    G = y[1::2, 0::2, 1::2]
+    H = y[1::2, 1::2, 1::2]
 
     # Combine to form subbands
     p = ( A-G-D-F) * j2[0] + ( B-H+C+E) * j2[1]
@@ -572,15 +561,15 @@ def c2cube(z):
 
     y = np.zeros(np.asanyarray(z.shape[:3])*2, dtype=z.real.dtype)
 
-    y[1::2, 1::2, 1::2] = ( pr+qr+rr+sr)
-    y[0::2, 0::2, 1::2] = (-pr-qr+rr+sr)
-    y[1::2, 0::2, 0::2] = (-pr+qr+rr-sr)
-    y[0::2, 1::2, 0::2] = (-pr+qr-rr+sr)
+    y[0::2, 0::2, 0::2] = ( pr+qr+rr+sr)
+    y[1::2, 0::2, 1::2] = (-pr-qr+rr+sr)
+    y[1::2, 1::2, 0::2] = (-pr+qr+rr-sr)
+    y[0::2, 1::2, 1::2] = (-pr+qr-rr+sr)
 
-    y[1::2, 1::2, 0::2] = ( pi-qi+ri-si)
-    y[0::2, 0::2, 0::2] = (-pi+qi+ri-si)
-    y[1::2, 0::2, 1::2] = ( pi+qi-ri-si)
-    y[0::2, 1::2, 1::2] = ( pi+qi+ri+si)
+    y[0::2, 1::2, 0::2] = ( pi-qi+ri-si)
+    y[1::2, 1::2, 1::2] = (-pi+qi+ri-si)
+    y[1::2, 0::2, 0::2] = ( pi+qi-ri-si)
+    y[0::2, 0::2, 1::2] = ( pi+qi+ri+si)
 
     return y * scale
 
diff --git a/tests/qbgn.mat b/tests/qbgn.mat
index 1d9445d..10edef1 100644
Binary files a/tests/qbgn.mat and b/tests/qbgn.mat differ
diff --git a/tests/testagainstmatlab.py b/tests/testagainstmatlab.py
index 9173f64..287b9a8 100644
--- a/tests/testagainstmatlab.py
+++ b/tests/testagainstmatlab.py
@@ -116,10 +116,10 @@ def test_rescale_highpass():
     assert_percentile_almost_equal_to_summary(Xrescale, verif['lena_upsample'], 60, tolerance=TOLERANCE)
 
 def test_transform3d_numpy():
-    transform = Transform3d(biort='near_sym_a',qshift='qshift_b')
+    transform = Transform3d(biort='near_sym_b',qshift='qshift_b')
     td_signal = transform.forward(qbgn, nlevels=3, include_scale=True, discard_level_1=False)
     Yl, Yh, Yscale = td_signal.lowpass, td_signal.subbands, td_signal.scales
-    #assert_almost_equal_to_summary_cube(Yl, verif['qbgn_Yl'], tolerance=TOLERANCE)
+    assert_almost_equal_to_summary_cube(Yl, verif['qbgn_Yl'], tolerance=TOLERANCE)
     for idx, a in enumerate(Yh):
         assert_almost_equal_to_summary_cube(a, verif['qbgn_Yh_{0}'.format(idx)], tolerance=TOLERANCE)
 
diff --git a/tests/verification.npz b/tests/verification.npz
index fa586ce..847c35a 100644
Binary files a/tests/verification.npz and b/tests/verification.npz differ

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



More information about the debian-science-commits mailing list