shape broadcast error with odd shaped image arrays

[pr4deepr]

Hi Thanks for making this library. I'm writing single z planes using append_plane function and I've been running into issues with image shapes that are odd.

If my array is odd shape and I try subsampling, I get a shape broadcast error.

For example, here is a minimal example:

import numpy as np 
import npy2bdv
save_path = "test.h5"

#initial image
img = np.zeros((147,1001,1001))
nz,ny,nx = img.shape

bdv_writer = npy2bdv.BdvWriter(save_path, nchannels=1,subsamp=((1,2,2),(1, 8, 8),(1, 16, 16)),
                            blockdim=((8, 32, 32), (8, 32, 32), (8, 32, 32)),
                            compression='gzip', overwrite=True)

bdv_writer.append_view(stack=None, virtual_stack_dim=(nz, ny, nx), channel=0)
bdv_writer.append_plane(plane=img[0], channel=0,z=0)

I get the error:

\npy2bdv\npy2bdv.py:409, in BdvWriter.append_plane(self, plane, z, time, illumination, channel, tile, angle)
    407 print(plane.shape)
    408 print(self.subsamp[ilevel])
--> 409 dataset[z, :, :] = self._subsample_plane(plane, self.subsamp[ilevel]).astype('int16')

File h5py\_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py\_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File ~\AppData\Roaming\Python\Python39\site-packages\h5py\_hl\dataset.py:997, in Dataset.__setitem__(self, args, val)
    994     mshape = val.shape
    996 # Perform the write, with broadcasting
--> 997 mspace = h5s.create_simple(selection.expand_shape(mshape))
    998 for fspace in selection.broadcast(mshape):
    999     self.id.write(mspace, fspace, val, mtype, dxpl=self._dxpl)
...
    267     # All dimensions from target_shape should either have been popped
    268     # to match the selection shape, or be 1.
    269     raise TypeError("Can't broadcast %s -> %s" % (source_shape, self.array_shape))  # array shape

TypeError: Can't broadcast (501, 501) -> (500, 500)

I believe its do with how npy2bdv and the skimage.transform.downscale_local_mean function calculates image shapes... Latter returns 501,501, whereas npy2bdv calculates it as (500,500)

npy2bdv calculates it here: https://github.com/nvladimus/npy2bdv/blob/4d1ffd356b0ec99f45e26fe8327275602fd12095/npy2bdv/npy2bdv.py#L512

downscale_local_mean function uses this block function: https://github.com/scikit-image/scikit-image/blob/441fe68b95a86d4ae2a351311a0c39a4232b6521/skimage/measure/block.py#L78

I've essentially modified your code to round up the shape which seems to solve this:

grp.create_dataset('cells', chunks=self.chunks[ilevel],
                                   shape = np.ceil(virtual_stack_dim / self.subsamp[ilevel]),
                                   compression=self.compression, dtype='int16')

Not sure if this is the best way to do this.

Cheers Pradeep

[nvladimus]

Hi, Pradeep, Thanks for reporting this issue. Can you please create a PR? I will test it and merge into the main if all tests pass. Many thanks, Nikita

[nvladimus]

Fixed and merged PR #15, thanks!