The 3d_dtcwt_directionality.py script in the examples directory shows how one may demonstrate the directional sensitivity of the 3D DT-CWT complex subband coefficients. It computes empirically the maximally sensitive directions for each subband and plots them in an interactive figure using matplotlib. A screenshot is reproduced below:
There are some points to note about this diagram. Each subband is labeled sich that ‘1’ refers to the first subband, ‘5’ the fifth and so forth. On this diagram the subbands are all four apart reflecting the fact that, for example, subbands 2, 3 and 4 are positioned in the other four quadrants of the upper hemisphere reflecting the position of subband 1. There are seven visible subband directions in the +ve quadrant of the hemisphere and hence there are 28 directions in total over all four quadrants.
The source for the script is shown below:
#!/bin/python """ An example of the directional selectivity of 3D DT-CWT coefficients. This example creates a 3D array holding an image of a sphere and performs the 3D DT-CWT transform on it. The locations of maxima (and their images about the mid-point of the image) are determined for each complex coefficient at level 2. These maxima points are then shown on a single plot to demonstrate the directions in which the 3D DT-CWT transform is selective. """ # Import the libraries we need from matplotlib.pyplot import * import numpy as np from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection from dtcwt import dtwavexfm3, dtwaveifm3, biort, qshift # Specify details about sphere and grid size GRID_SIZE = 128 SPHERE_RAD = int(0.45 * GRID_SIZE) + 0.5 # Compute an image of the sphere grid = np.arange(-(GRID_SIZE>>1), GRID_SIZE>>1) X, Y, Z = np.meshgrid(grid, grid, grid) r = np.sqrt(X*X + Y*Y + Z*Z) sphere = (0.5 + np.clip(SPHERE_RAD-r, -0.5, 0.5)).astype(np.float32) # Specify number of levels and wavelet family to use nlevels = 2 b = biort('near_sym_a') q = qshift('qshift_a') # Form the DT-CWT of the sphere. We use discard_level_1 since we're # uninterested in the inverse transform and this saves us some memory. Yl, Yh = dtwavexfm3(sphere, nlevels, b, q, discard_level_1=True) # Plot maxima figure(figsize=(8,8)) ax = gcf().add_subplot(1,1,1, projection='3d') ax.set_aspect('equal') ax.view_init(35, 75) # Plot unit sphere +ve octant thetas = np.linspace(0, np.pi/2, 10) phis = np.linspace(0, np.pi/2, 10) def sphere_to_xyz(r, theta, phi): st, ct = np.sin(theta), np.cos(theta) sp, cp = np.sin(phi), np.cos(phi) return r * np.asarray((st*cp, st*sp, ct)) tris =  rad = 0.99 # so that points plotted latter are not z-clipped for t1, t2 in zip(thetas[:-1], thetas[1:]): for p1, p2 in zip(phis[:-1], phis[1:]): tris.append([ sphere_to_xyz(rad, t1, p1), sphere_to_xyz(rad, t1, p2), sphere_to_xyz(rad, t2, p2), sphere_to_xyz(rad, t2, p1), ]) sphere = Poly3DCollection(tris, facecolor='w', edgecolor=(0.6,0.6,0.6)) ax.add_collection3d(sphere) locs =  scale = 1.1 for idx in xrange(Yh[-1].shape): Z = Yh[-1][:,:,:,idx] C = np.abs(Z) max_loc = np.asarray(np.unravel_index(np.argmax(C), C.shape)) - np.asarray(C.shape)*0.5 max_loc /= np.sqrt(np.sum(max_loc * max_loc)) # Only record directions in the +ve octant (or those from the -ve quadrant # which can be flipped). if np.all(np.sign(max_loc) == 1): locs.append(max_loc) ax.text(max_loc * scale, max_loc * scale, max_loc * scale, str(idx+1)) elif np.all(np.sign(max_loc) == -1): locs.append(-max_loc) ax.text(-max_loc * scale, -max_loc * scale, -max_loc * scale, str(idx+1)) # Plot all directions as a scatter plot locs = np.asarray(locs) ax.scatter(locs[:,0], locs[:,1], locs[:,2], c=np.arange(locs.shape)) w = 1.1 ax.auto_scale_xyz([0, w], [0, w], [0, w]) legend() title('3D DT-CWT subband directions for +ve hemisphere quadrant') tight_layout() show() # vim:sw=4:sts=4:et