highpass sampling question

[scf32]

Is there any particular reason for the highpass sampling functions to return M x N x 6 x 6 arrays instead of M x N x 6?

I am only interested because my code is calling upon this quite a lot and having to then slice [: , :, sb, sb] rather than just [:, :, sb]. If there is unnecessary computation going on I could do without it as my SLP function spends 60% of its time on resampling!

I wonder if the upsample_highpass() function could be given an offset argument? That way it would sample from a regular grid which didn't have to be fixed between actual sample locations, and I'd only have to call it 12 times with a separate offset for each, and wouldn't even need to define a grid of sampling locations first.

Of course, this is quite an application specific request but I don't want to go changing it myself if there is a good reason for having it this way.

[rjw57]

Is there any particular reason for the highpass sampling functions to return M x N x 6 x 6 arrays instead of M x N x 6?

I'm not sure I understand this. As far as I'm aware it doesn't return MxNx6x6 arrays, only ever MxNx6 ones:

import numpy as np
import dtcwt
import dtcwt.sampling

# Load "Lena" and transform
lena = np.load('lena.npz')['lena']
lena_t = dtcwt.Transform2d().forward(lena)

# What is the coarsest highpass shape?
print('Coarsest highpass shape: {0}'.format(lena_t.highpasses[-1].shape))

# Rescale it to the size of Lena
resampled_hp = dtcwt.sampling.rescale_highpass(lena_t.highpasses[-1],
    lena.shape, method='bilinear')
print('Re-sampled to: {0}'.format(resampled_hp.shape))

# Sample from two arbitrary locations
xs, ys = np.array(((4, 7),)), np.array(((8, 9),))
sampled_hp = dtcwt.sampling.sample_highpass(lena_t.highpasses[-1],
    xs, ys, method='bilinear')
print('Sampled highpass has shape: {0}'.format(sampled_hp.shape))

results in:

Coarsest highpass shape: (64, 64, 6)
Re-sampled to: (512, 512, 6)
Sampled highpass has shape: (1, 2, 6)

Am I missing something? (Beyond the annoying wrinkle that xs and ys need to be two-dimensional.)

[rjw57]

... or is this a request for a {rescale,sample}_subband function which acts on individual MxN subbands? Something like:

sbidx = 4 # valid values are integers on the interval [0, 6)
rescaled_sb = dtcwt.sampling.rescale_subband(lena_t.highpasses[-1][:,:,sbidx],
    sbidx, lena.shape)

Or, even better, taking a sequence of subband indices corresponding to slices of the input. That way {rescale,sample}_highpass become convenience wrappers.

Is that the sort of thing you're getting at?

[rjw57]

... oh... or is it that it is non-obvious that the shape argument to {rescale,sample}_highpass takes two elements, not three? That isn't obvious from the documentation. If so, do please suggest a clarification. Doing dtcwt.sampling.rescale_highpass(..., (M,N)) is correct. Doing dtcwt.sampling.rescale_highpass(..., (M,N,6)) is not.

[rjw57]

I wonder if the upsample_highpass() function could be given an offset argument? That way it would sample from a regular grid which didn't have to be fixed between actual sample locations, and I'd only have to call it 12 times with a separate offset for each, and wouldn't even need to define a grid of sampling locations first.

I'm not sure I follow this. Do you mean something like (not-tested):

hp = lena_t.highpasses[-1]

# Set up a sample grid with 8 pixel pitch
xs, ys = np.meshgrid(np.arange(0, lena.shape[1], 8),
    np.arange(0, lena.shape[0], 8))

# Sample from the grid offset by (dx, dy)
for dx in range(8):
    for dy in range(8):
        sampled_hp = dtcwt.sampling.sample_highpass(hp, xs + dx, ys + dy)
        # ... do something with sampled_hp

Although in that case it's probably a bit faster to resample hp to lena.shape to begin with and then sample via slice notation (again, untested):

hp = lena_t.highpasses[-1]

# Upscale highpass
upscale_hp = dtcwt.sampling.rescale_highpass(hp, lena.shape)

# Sample from the 8 pixel pitch grid offset by (dx, dy)
for dx in range(8):
    for dy in range(8):
        sampled_hp = upscale_hp[dx::8, dy::8]
        # ... do something with sampled_hp

[rjw57]

@scf32 If either of #91 or #92 are related to this issue, please weigh in on them too if you have an opinion.