From 709d437ec6119f0669c07276e1a1ae35eb73822f Mon Sep 17 00:00:00 2001 From: Kaidong Chai Date: Thu, 7 Nov 2024 15:03:49 -0500 Subject: [PATCH] Cleanup --- linc_convert/__init__.py | 1 + linc_convert/cli.py | 1 + linc_convert/modalities/__init__.py | 1 + linc_convert/modalities/df/__init__.py | 1 + linc_convert/modalities/df/cli.py | 1 + linc_convert/modalities/df/multi_slice.py | 88 +++--- linc_convert/modalities/df/single_slice.py | 56 ++-- linc_convert/modalities/lsm/__init__.py | 1 + linc_convert/modalities/lsm/cli.py | 1 + linc_convert/modalities/lsm/mosaic.py | 294 +++++++++++---------- linc_convert/utils/j2k.py | 10 +- linc_convert/utils/math.py | 1 + linc_convert/utils/orientation.py | 6 +- linc_convert/utils/zarr.py | 7 +- scripts/utils.py | 31 +-- tests/helper.py | 1 + 16 files changed, 270 insertions(+), 231 deletions(-) diff --git a/linc_convert/__init__.py b/linc_convert/__init__.py index 3f242ef..1d5a710 100644 --- a/linc_convert/__init__.py +++ b/linc_convert/__init__.py @@ -1,3 +1,4 @@ """Data conversion tools for the LINC project.""" + __all__ = ["modalities", "utils"] from . import modalities, utils diff --git a/linc_convert/cli.py b/linc_convert/cli.py index ea9a2d0..aaeca28 100644 --- a/linc_convert/cli.py +++ b/linc_convert/cli.py @@ -1,4 +1,5 @@ """Root command line entry point.""" + from cyclopts import App help = "Collection of conversion scripts for LINC datasets" diff --git a/linc_convert/modalities/__init__.py b/linc_convert/modalities/__init__.py index 04518af..7ae40c8 100644 --- a/linc_convert/modalities/__init__.py +++ b/linc_convert/modalities/__init__.py @@ -1,3 +1,4 @@ """Converters for all imaging modalities.""" + __all__ = ["df", "lsm"] from . import df, lsm diff --git a/linc_convert/modalities/df/__init__.py b/linc_convert/modalities/df/__init__.py index 3649993..a6d9d6f 100644 --- a/linc_convert/modalities/df/__init__.py +++ b/linc_convert/modalities/df/__init__.py @@ -1,4 +1,5 @@ """Dark Field microscopy converters.""" + __all__ = ["cli", "multi_slice", "single_slice"] from . import cli, multi_slice, single_slice diff --git a/linc_convert/modalities/df/cli.py b/linc_convert/modalities/df/cli.py index 3c51f07..9ca3fdc 100644 --- a/linc_convert/modalities/df/cli.py +++ b/linc_convert/modalities/df/cli.py @@ -1,4 +1,5 @@ """Entry-points for Dark Field microscopy converter.""" + from cyclopts import App from linc_convert.cli import main diff --git a/linc_convert/modalities/df/multi_slice.py b/linc_convert/modalities/df/multi_slice.py index 9c3bf74..d66136f 100644 --- a/linc_convert/modalities/df/multi_slice.py +++ b/linc_convert/modalities/df/multi_slice.py @@ -4,6 +4,7 @@ We do not recompute the image pyramid but instead reuse the JPEG2000 levels (obtained by wavelet transform). """ + # stdlib import ast import json @@ -130,7 +131,7 @@ def convert( new_width = jp2.shape[1] new_size = (new_height, new_width) if has_channel: - new_size += (3.,) + new_size += (3.0,) print(len(inp), new_size, nblevel, has_channel) # Prepare chunking options @@ -158,15 +159,25 @@ def convert( vxw, vxh = get_pixelsize(j2k) subdat = WrappedJ2K(j2k, level=level) subdat_size = subdat.shape - print("Convert level", level, "with shape", shape, - "for slice", idx, "with size", subdat_size) + print( + "Convert level", + level, + "with shape", + shape, + "for slice", + idx, + "with size", + subdat_size, + ) # offset while attaching x = floordiv(shape[-2] - subdat_size[-2], 2) y = floordiv(shape[-1] - subdat_size[-1], 2) if max_load is None or (shape[-2] < max_load and shape[-1] < max_load): - array[..., idx, x:x + subdat_size[1], y:y + subdat_size[2]] = subdat[...] + array[..., idx, x : x + subdat_size[1], y : y + subdat_size[2]] = ( + subdat[...] + ) else: ni = ceildiv(shape[-2], max_load) @@ -175,8 +186,14 @@ def convert( for i in range(ni): for j in range(nj): print(f"\r{i+1}/{ni}, {j+1}/{nj}", end=" ") - start_x, end_x = i*max_load, min((i+1)*max_load, shape[-2]) - start_y, end_y = j*max_load, min((j+1)*max_load, shape[-1]) + start_x, end_x = ( + i * max_load, + min((i + 1) * max_load, shape[-2]), + ) + start_y, end_y = ( + j * max_load, + min((j + 1) * max_load, shape[-1]), + ) if end_x <= x or end_y <= y: continue @@ -186,28 +203,30 @@ def convert( array[ ..., idx, - x + start_x:x + min(end_x, subdat_size[-2]), - y + start_y:y + min(end_y, subdat_size[-1]), + x + start_x : x + min(end_x, subdat_size[-2]), + y + start_y : y + min(end_y, subdat_size[-1]), ] = subdat[ ..., - start_x:min((i+1)*max_load, subdat_size[-2]), - start_y:min((j+1)*max_load, subdat_size[-1]), + start_x : min((i + 1) * max_load, subdat_size[-2]), + start_y : min((j + 1) * max_load, subdat_size[-1]), ] print("") # Write OME-Zarr multiscale metadata print("Write metadata") - multiscales = [{ - "version": "0.4", - "axes": [ - {"name": "z", "type": "space", "unit": "micrometer"}, - {"name": "y", "type": "distance", "unit": "micrometer"}, - {"name": "x", "type": "space", "unit": "micrometer"} - ], - "datasets": [], - "type": "jpeg2000", - "name": "", - }] + multiscales = [ + { + "version": "0.4", + "axes": [ + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "distance", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"}, + ], + "datasets": [], + "type": "jpeg2000", + "name": "", + } + ] if has_channel: multiscales[0]["axes"].insert(0, {"name": "c", "type": "channel"}) @@ -225,26 +244,25 @@ def convert( level["coordinateTransformations"] = [ { "type": "scale", - "scale": [1.0] * has_channel + [ + "scale": [1.0] * has_channel + + [ 1.0, - (shape0[0]/shape[0])*vxh, - (shape0[1]/shape[1])*vxw, - ] + (shape0[0] / shape[0]) * vxh, + (shape0[1] / shape[1]) * vxw, + ], }, { "type": "translation", - "translation": [0.0] * has_channel + [ + "translation": [0.0] * has_channel + + [ 0.0, - (shape0[0]/shape[0] - 1)*vxh*0.5, - (shape0[1]/shape[1] - 1)*vxw*0.5, - ] - } + (shape0[0] / shape[0] - 1) * vxh * 0.5, + (shape0[1] / shape[1] - 1) * vxw * 0.5, + ], + }, ] multiscales[0]["coordinateTransformations"] = [ - { - "scale": [1.0] * (3 + has_channel), - "type": "scale" - } + {"scale": [1.0] * (3 + has_channel), "type": "scale"} ] omz.attrs["multiscales"] = multiscales @@ -278,7 +296,7 @@ def convert( # Write sidecar .json file json_name = os.path.splitext(out)[0] - json_name += '.json' + json_name += ".json" dic = {} dic["PixelSize"] = json.dumps([vxw, vxh]) dic["PixelSizeUnits"] = "um" diff --git a/linc_convert/modalities/df/single_slice.py b/linc_convert/modalities/df/single_slice.py index 566c5c1..7bc7782 100644 --- a/linc_convert/modalities/df/single_slice.py +++ b/linc_convert/modalities/df/single_slice.py @@ -4,6 +4,7 @@ It does not recompute the image pyramid but instead reuse the JPEG2000 levels (obtained by wavelet transform). """ + # stdlib import ast import os @@ -134,27 +135,29 @@ def convert( print(f"\r{i+1}/{ni}, {j+1}/{nj}", end="") array[ ..., - i*max_load:min((i+1)*max_load, shape[-2]), - j*max_load:min((j+1)*max_load, shape[-1]), + i * max_load : min((i + 1) * max_load, shape[-2]), + j * max_load : min((j + 1) * max_load, shape[-1]), ] = subdat[ ..., - i*max_load:min((i+1)*max_load, shape[-2]), - j*max_load:min((j+1)*max_load, shape[-1]), + i * max_load : min((i + 1) * max_load, shape[-2]), + j * max_load : min((j + 1) * max_load, shape[-1]), ] print("") # Write OME-Zarr multiscale metadata print("Write metadata") - multiscales = [{ - "version": "0.4", - "axes": [ - {"name": "y", "type": "space", "unit": "micrometer"}, - {"name": "x", "type": "space", "unit": "micrometer"} - ], - "datasets": [], - "type": "jpeg2000", - "name": "", - }] + multiscales = [ + { + "version": "0.4", + "axes": [ + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"}, + ], + "datasets": [], + "type": "jpeg2000", + "name": "", + } + ] if has_channel: multiscales[0]["axes"].insert(0, {"name": "c", "type": "channel"}) @@ -172,24 +175,23 @@ def convert( level["coordinateTransformations"] = [ { "type": "scale", - "scale": [1.0] * has_channel + [ - (shape0[0]/shape[0])*vxh, - (shape0[1]/shape[1])*vxw, - ] + "scale": [1.0] * has_channel + + [ + (shape0[0] / shape[0]) * vxh, + (shape0[1] / shape[1]) * vxw, + ], }, { "type": "translation", - "translation": [0.0] * has_channel + [ - (shape0[0]/shape[0] - 1)*vxh*0.5, - (shape0[1]/shape[1] - 1)*vxw*0.5, - ] - } + "translation": [0.0] * has_channel + + [ + (shape0[0] / shape[0] - 1) * vxh * 0.5, + (shape0[1] / shape[1] - 1) * vxw * 0.5, + ], + }, ] multiscales[0]["coordinateTransformations"] = [ - { - "scale": [1.0] * (2 + has_channel), - "type": "scale" - } + {"scale": [1.0] * (2 + has_channel), "type": "scale"} ] omz.attrs["multiscales"] = multiscales diff --git a/linc_convert/modalities/lsm/__init__.py b/linc_convert/modalities/lsm/__init__.py index 3c60643..1ce6617 100644 --- a/linc_convert/modalities/lsm/__init__.py +++ b/linc_convert/modalities/lsm/__init__.py @@ -1,3 +1,4 @@ """Light Sheet Microscopy converters.""" + __all__ = ["cli", "mosaic"] from . import cli, mosaic diff --git a/linc_convert/modalities/lsm/cli.py b/linc_convert/modalities/lsm/cli.py index 5c7e8d1..2026322 100644 --- a/linc_convert/modalities/lsm/cli.py +++ b/linc_convert/modalities/lsm/cli.py @@ -1,4 +1,5 @@ """Entry-points for Dark Field microscopy converter.""" + from cyclopts import App from linc_convert.cli import main diff --git a/linc_convert/modalities/lsm/mosaic.py b/linc_convert/modalities/lsm/mosaic.py index fcd463b..9cbfe2c 100644 --- a/linc_convert/modalities/lsm/mosaic.py +++ b/linc_convert/modalities/lsm/mosaic.py @@ -4,6 +4,7 @@ Example input files can be found at https://lincbrain.org/dandiset/000004/0.240319.1924/files?location=derivatives%2F """ + # stdlib import ast import os @@ -29,18 +30,18 @@ @mosaic.default def convert( - inp: str, - out: str = None, - *, - chunk: int = 128, - compressor: str = 'blosc', - compressor_opt: str = "{}", - max_load: int = 512, - nii: bool = False, - orientation: str = 'coronal', - center: bool = True, - thickness: float | None = None, - voxel_size: list[float] = (1, 1, 1), + inp: str, + out: str = None, + *, + chunk: int = 128, + compressor: str = "blosc", + compressor_opt: str = "{}", + max_load: int = 512, + nii: bool = False, + orientation: str = "coronal", + center: bool = True, + thickness: float | None = None, + voxel_size: list[float] = (1, 1, 1), ) -> None: """ Convert a collection of tiff files generated by the LSM pipeline into ZARR. @@ -97,13 +98,10 @@ def convert( max_load += 1 CHUNK_PATTERN = re.compile( - r'^(?P\w*)' - r'_z(?P[0-9]+)' - r'_y(?P[0-9]+)' - r'(?P\w*)$' + r"^(?P\w*)" r"_z(?P[0-9]+)" r"_y(?P[0-9]+)" r"(?P\w*)$" ) - all_chunks_dirnames = list(sorted(glob(os.path.join(inp, '*_z*_y*')))) + all_chunks_dirnames = list(sorted(glob(os.path.join(inp, "*_z*_y*")))) all_chunks_info = dict( dirname=[], prefix=[], @@ -124,70 +122,70 @@ def convert( # parse all directory names for dirname in all_chunks_dirnames: parsed = CHUNK_PATTERN.fullmatch(os.path.basename(dirname)) - all_chunks_info['dirname'].append(dirname) - all_chunks_info['prefix'].append(parsed.group('prefix')) - all_chunks_info['suffix'].append(parsed.group('suffix')) - all_chunks_info['z'].append(int(parsed.group('z'))) - all_chunks_info['y'].append(int(parsed.group('y'))) + all_chunks_info["dirname"].append(dirname) + all_chunks_info["prefix"].append(parsed.group("prefix")) + all_chunks_info["suffix"].append(parsed.group("suffix")) + all_chunks_info["z"].append(int(parsed.group("z"))) + all_chunks_info["y"].append(int(parsed.group("y"))) # default output name if not out: - out = all_chunks_info['prefix'][0] + all_chunks_info['suffix'][0] - out += '.nii.zarr' if nii else '.ome.zarr' - nii = nii or out.endswith('.nii.zarr') + out = all_chunks_info["prefix"][0] + all_chunks_info["suffix"][0] + out += ".nii.zarr" if nii else ".ome.zarr" + nii = nii or out.endswith(".nii.zarr") # parse all individual file names - nchunkz = max(all_chunks_info['z']) - nchunky = max(all_chunks_info['y']) + nchunkz = max(all_chunks_info["z"]) + nchunky = max(all_chunks_info["y"]) allshapes = [[(0, 0, 0) for _ in range(nchunky)] for _ in range(nchunkz)] nchannels = 0 dtype = None for zchunk in range(nchunkz): for ychunk in range(nchunky): - for i in range(len(all_chunks_info['dirname'])): - if all_chunks_info['z'][i] == zchunk + 1 \ - and all_chunks_info['y'][i] == ychunk + 1: + for i in range(len(all_chunks_info["dirname"])): + if ( + all_chunks_info["z"][i] == zchunk + 1 + and all_chunks_info["y"][i] == ychunk + 1 + ): break - dirname = all_chunks_info['dirname'][i] - planes_filenames \ - = list(sorted(glob(os.path.join(dirname, '*.tiff')))) + dirname = all_chunks_info["dirname"][i] + planes_filenames = list(sorted(glob(os.path.join(dirname, "*.tiff")))) PLANE_PATTERN = re.compile( - os.path.basename(dirname) + - r'_plane(?P[0-9]+)' - r'_c(?P[0-9]+)' - r'.tiff$' + os.path.basename(dirname) + r"_plane(?P[0-9]+)" + r"_c(?P[0-9]+)" + r".tiff$" ) for fname in planes_filenames: parsed = PLANE_PATTERN.fullmatch(os.path.basename(fname)) - all_chunks_info['planes'][i]['fname'] += [fname] - all_chunks_info['planes'][i]['z'] += [int(parsed.group('z'))] - all_chunks_info['planes'][i]['c'] += [int(parsed.group('c'))] + all_chunks_info["planes"][i]["fname"] += [fname] + all_chunks_info["planes"][i]["z"] += [int(parsed.group("z"))] + all_chunks_info["planes"][i]["c"] += [int(parsed.group("c"))] f = TiffFile(fname) dtype = f.pages[0].dtype yx_shape = f.pages[0].shape - all_chunks_info['planes'][i]['yx_shape'].append(yx_shape) + all_chunks_info["planes"][i]["yx_shape"].append(yx_shape) - nplanes = max(all_chunks_info['planes'][i]['z']) - nchannels = max(nchannels, max(all_chunks_info['planes'][i]['c'])) + nplanes = max(all_chunks_info["planes"][i]["z"]) + nchannels = max(nchannels, max(all_chunks_info["planes"][i]["c"])) - yx_shape = set(all_chunks_info['planes'][i]['yx_shape']) + yx_shape = set(all_chunks_info["planes"][i]["yx_shape"]) if not len(yx_shape) == 1: - raise ValueError('Incompatible chunk shapes') + raise ValueError("Incompatible chunk shapes") yx_shape = list(yx_shape)[0] allshapes[zchunk][ychunk] = (nplanes, *yx_shape) # check that all chink shapes are compatible for zchunk in range(nchunkz): if len(set(shape[1] for shape in allshapes[zchunk])) != 1: - raise ValueError('Incompatible Y shapes') + raise ValueError("Incompatible Y shapes") for ychunk in range(nchunky): if len(set(shape[ychunk][0] for shape in allshapes)) != 1: - raise ValueError('Incompatible Z shapes') + raise ValueError("Incompatible Z shapes") if len(set(shape[2] for subshapes in allshapes for shape in subshapes)) != 1: - raise ValueError('Incompatible X shapes') + raise ValueError("Incompatible X shapes") # compute full shape fullshape = [0, 0, 0] @@ -201,30 +199,36 @@ def convert( # Prepare chunking options opt = { - 'chunks': [nchannels] + [chunk] * 3, - 'dimension_separator': r'/', - 'order': 'F', - 'dtype': np.dtype(dtype).str, - 'fill_value': None, - 'compressor': make_compressor(compressor, **compressor_opt), + "chunks": [nchannels] + [chunk] * 3, + "dimension_separator": r"/", + "order": "F", + "dtype": np.dtype(dtype).str, + "fill_value": None, + "compressor": make_compressor(compressor, **compressor_opt), } # write first level - omz.create_dataset('0', shape=[nchannels, *fullshape], **opt) - array = omz['0'] - print('Write level 0 with shape', [nchannels, *fullshape]) - for i, dirname in enumerate(all_chunks_info['dirname']): - chunkz = all_chunks_info['z'][i] - 1 - chunky = all_chunks_info['y'][i] - 1 - planes = all_chunks_info['planes'][i] - for j, fname in enumerate(planes['fname']): - subz = planes['z'][j] - 1 - subc = planes['c'][j] - 1 - yx_shape = planes['yx_shape'][j] + omz.create_dataset("0", shape=[nchannels, *fullshape], **opt) + array = omz["0"] + print("Write level 0 with shape", [nchannels, *fullshape]) + for i, dirname in enumerate(all_chunks_info["dirname"]): + chunkz = all_chunks_info["z"][i] - 1 + chunky = all_chunks_info["y"][i] - 1 + planes = all_chunks_info["planes"][i] + for j, fname in enumerate(planes["fname"]): + subz = planes["z"][j] - 1 + subc = planes["c"][j] - 1 + yx_shape = planes["yx_shape"][j] zstart = sum(shape[0][0] for shape in allshapes[:chunkz]) - ystart = sum(shape[1] for subshapes in allshapes for shape in subshapes[:chunky]) - print(f'Write plane ({subc}, {zstart + subz}, {ystart}:{ystart + yx_shape[0]})', end='\r') + ystart = sum( + shape[1] for subshapes in allshapes for shape in subshapes[:chunky] + ) + print( + f"Write plane " + f"({subc}, {zstart + subz}, {ystart}:{ystart + yx_shape[0]})", + end="\r", + ) slicer = ( subc, zstart + subz, @@ -234,7 +238,7 @@ def convert( f = TiffFile(fname) array[slicer] = f.asarray() - print('') + print("") # build pyramid using median windows level = 0 @@ -243,10 +247,10 @@ def convert( prev_shape = prev_array.shape[-3:] level += 1 - new_shape = list(map(lambda x: max(1, x//2), prev_shape)) + new_shape = list(map(lambda x: max(1, x // 2), prev_shape)) if all(x < chunk for x in new_shape): break - print('Compute level', level, 'with shape', new_shape) + print("Compute level", level, "with shape", new_shape) omz.create_dataset(str(level), shape=[nchannels, *new_shape], **opt) new_array = omz[str(level)] @@ -258,74 +262,81 @@ def convert( for cz in range(ncz): for cy in range(ncy): for cx in range(ncx): - - print(f'chunk ({cz}, {cy}, {cx}) / ({ncz}, {ncy}, {ncx})', - end='\r') + print(f"chunk ({cz}, {cy}, {cx}) / ({ncz}, {ncy}, {ncx})", end="\r") dat = prev_array[ - ..., - cz*max_load:(cz+1)*max_load, - cy*max_load:(cy+1)*max_load, - cx*max_load:(cx+1)*max_load, - ] + ..., + cz * max_load: (cz + 1) * max_load, + cy * max_load: (cy + 1) * max_load, + cx * max_load: (cx + 1) * max_load, + ] crop = [0 if x == 1 else x % 2 for x in dat.shape[-3:]] slicer = [slice(-1) if x else slice(None) for x in crop] dat = dat[(Ellipsis, *slicer)] pz, py, px = dat.shape[-3:] - dat = dat.reshape([ - nchannels, - max(pz//2, 1), min(pz, 2), - max(py//2, 1), min(py, 2), - max(px//2, 1), min(px, 2), - ]) + dat = dat.reshape( + [ + nchannels, + max(pz // 2, 1), + min(pz, 2), + max(py // 2, 1), + min(py, 2), + max(px // 2, 1), + min(px, 2), + ] + ) dat = dat.transpose([0, 1, 3, 5, 2, 4, 6]) - dat = dat.reshape([ - nchannels, - max(pz//2, 1), - max(py//2, 1), - max(px//2, 1), - -1, - ]) + dat = dat.reshape( + [ + nchannels, + max(pz // 2, 1), + max(py // 2, 1), + max(px // 2, 1), + -1, + ] + ) dat = np.median(dat, -1) new_array[ - ..., - cz*max_load//2:(cz+1)*max_load//2, - cy*max_load//2:(cy+1)*max_load//2, - cx*max_load//2:(cx+1)*max_load//2, + ..., + cz * max_load // 2: (cz + 1) * max_load // 2, + cy * max_load // 2: (cy + 1) * max_load // 2, + cx * max_load // 2: (cx + 1) * max_load // 2, ] = dat - print('') + print("") nblevel = level # Write OME-Zarr multiscale metadata - print('Write metadata') - multiscales = [{ - 'version': '0.4', - 'axes': [ - {"name": "z", "type": "space", "unit": "micrometer"}, - {"name": "y", "type": "space", "unit": "micrometer"}, - {"name": "x", "type": "space", "unit": "micrometer"} - ], - 'datasets': [], - 'type': 'median window 2x2x2', - 'name': '', - }] - multiscales[0]['axes'].insert(0, {"name": "c", "type": "channel"}) + print("Write metadata") + multiscales = [ + { + "version": "0.4", + "axes": [ + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"}, + ], + "datasets": [], + "type": "median window 2x2x2", + "name": "", + } + ] + multiscales[0]["axes"].insert(0, {"name": "c", "type": "channel"}) voxel_size = list(map(float, reversed(voxel_size))) factor = [1] * 3 for n in range(nblevel): shape = omz[str(n)].shape[-3:] - multiscales[0]['datasets'].append({}) - level = multiscales[0]['datasets'][-1] + multiscales[0]["datasets"].append({}) + level = multiscales[0]["datasets"][-1] level["path"] = str(n) # We made sure that the downsampling level is exactly 2 # However, once a dimension has size 1, we stop downsampling. if n > 0: - shape_prev = omz[str(n-1)].shape[-3:] + shape_prev = omz[str(n - 1)].shape[-3:] if shape_prev[0] != shape[0]: factor[0] *= 2 if shape_prev[1] != shape[1]: @@ -336,57 +347,56 @@ def convert( level["coordinateTransformations"] = [ { "type": "scale", - "scale": [1.0] + [ - factor[0] * voxel_size[0], - factor[1] * voxel_size[1], - factor[2] * voxel_size[2], - ] + "scale": [1.0] + + [ + factor[0] * voxel_size[0], + factor[1] * voxel_size[1], + factor[2] * voxel_size[2], + ], }, { "type": "translation", - "translation": [0.0] + [ - (factor[0] - 1) * voxel_size[0] * 0.5, - (factor[1] - 1) * voxel_size[1] * 0.5, - (factor[2] - 1) * voxel_size[2] * 0.5, - ] - } + "translation": [0.0] + + [ + (factor[0] - 1) * voxel_size[0] * 0.5, + (factor[1] - 1) * voxel_size[1] * 0.5, + (factor[2] - 1) * voxel_size[2] * 0.5, + ], + }, ] multiscales[0]["coordinateTransformations"] = [ - { - "scale": [1.0] * 4, - "type": "scale" - } + {"scale": [1.0] * 4, "type": "scale"} ] omz.attrs["multiscales"] = multiscales if not nii: - print('done.') + print("done.") return # Write NIfTI-Zarr header # NOTE: we use nifti2 because dimensions typically do not fit in a short # TODO: we do not write the json zattrs, but it should be added in # once the nifti-zarr package is released - shape = list(reversed(omz['0'].shape)) - shape = shape[:3] + [1] + shape[3:] # insert time dimension + shape = list(reversed(omz["0"].shape)) + shape = shape[:3] + [1] + shape[3:] # insert time dimension affine = orientation_to_affine(orientation, *voxel_size) if center: affine = center_affine(affine, shape[:3]) header = nib.Nifti2Header() header.set_data_shape(shape) - header.set_data_dtype(omz['0'].dtype) + header.set_data_dtype(omz["0"].dtype) header.set_qform(affine) header.set_sform(affine) - header.set_xyzt_units(nib.nifti1.unit_codes.code['micron']) - header.structarr['magic'] = b'nz2\0' - header = np.frombuffer(header.structarr.tobytes(), dtype='u1') + header.set_xyzt_units(nib.nifti1.unit_codes.code["micron"]) + header.structarr["magic"] = b"nz2\0" + header = np.frombuffer(header.structarr.tobytes(), dtype="u1") opt = { - 'chunks': [len(header)], - 'dimension_separator': r'/', - 'order': 'F', - 'dtype': '|u1', - 'fill_value': None, - 'compressor': None, + "chunks": [len(header)], + "dimension_separator": r"/", + "order": "F", + "dtype": "|u1", + "fill_value": None, + "compressor": None, } - omz.create_dataset('nifti', data=header, shape=shape, **opt) - print('done.') + omz.create_dataset("nifti", data=header, shape=shape, **opt) + print("done.") diff --git a/linc_convert/utils/j2k.py b/linc_convert/utils/j2k.py index f4c0309..c61805b 100644 --- a/linc_convert/utils/j2k.py +++ b/linc_convert/utils/j2k.py @@ -1,5 +1,5 @@ - """Utilities for JPEG2000 files.""" + # stdlib import uuid from dataclasses import dataclass @@ -87,7 +87,7 @@ def __getitem__(self, index: tuple[slice] | slice) -> np.ndarray: if Ellipsis not in index: index += (Ellipsis,) if any(idx is None for idx in index): - raise TypeError('newaxis not supported') + raise TypeError("newaxis not supported") # substitute ellipses new_index = [] @@ -99,13 +99,13 @@ def __getitem__(self, index: tuple[slice] | slice) -> np.ndarray: if not has_seen_ellipsis: new_index += [slice(None)] * nb_ellipsis elif not last_was_ellipsis: - raise ValueError('Multiple ellipses should be contiguous') + raise ValueError("Multiple ellipses should be contiguous") has_seen_ellipsis = True last_was_ellipsis = True elif not isinstance(idx, slice): - raise TypeError('Only slices are supported') + raise TypeError("Only slices are supported") elif idx.step not in (None, 1): - raise ValueError('Striding not supported') + raise ValueError("Striding not supported") else: last_was_ellipsis = False new_index += [idx] diff --git a/linc_convert/utils/math.py b/linc_convert/utils/math.py index 53e2a96..b8f009d 100644 --- a/linc_convert/utils/math.py +++ b/linc_convert/utils/math.py @@ -1,4 +1,5 @@ """Math utilities.""" + import math from numbers import Number diff --git a/linc_convert/utils/orientation.py b/linc_convert/utils/orientation.py index 336c537..0d26bde 100644 --- a/linc_convert/utils/orientation.py +++ b/linc_convert/utils/orientation.py @@ -1,4 +1,5 @@ """Orientation of an array of voxels with respect to world space.""" + import numpy as np @@ -32,10 +33,7 @@ def orientation_ensure_3d(orientation: str) -> str: def orientation_to_affine( - orientation: str, - vxw: float = 1, - vxh: float = 1, - vxd: float = 1 + orientation: str, vxw: float = 1, vxh: float = 1, vxd: float = 1 ) -> np.ndarray: """ Build an affine matrix from an orientation string and voxel size. diff --git a/linc_convert/utils/zarr.py b/linc_convert/utils/zarr.py index 227fb18..a2c031c 100644 --- a/linc_convert/utils/zarr.py +++ b/linc_convert/utils/zarr.py @@ -1,4 +1,5 @@ """Zarr utilities.""" + import numcodecs import numcodecs.abc @@ -9,10 +10,10 @@ def make_compressor(name: str, **prm: dict) -> numcodecs.abc.Codec: if not isinstance(name, str): return name name = name.lower() - if name == 'blosc': + if name == "blosc": Compressor = numcodecs.Blosc - elif name == 'zlib': + elif name == "zlib": Compressor = numcodecs.Zlib else: - raise ValueError('Unknown compressor', name) + raise ValueError("Unknown compressor", name) return Compressor(**prm) diff --git a/scripts/utils.py b/scripts/utils.py index d0f54f2..6bf0f01 100644 --- a/scripts/utils.py +++ b/scripts/utils.py @@ -1,21 +1,22 @@ import math + import numcodecs import numpy as np def orientation_ensure_3d(orientation): orientation = { - 'coronal': 'LI', - 'axial': 'LP', - 'sagittal': 'PI', + "coronal": "LI", + "axial": "LP", + "sagittal": "PI", }.get(orientation.lower(), orientation).upper() if len(orientation) == 2: - if 'L' not in orientation and 'R' not in orientation: - orientation += 'R' - if 'P' not in orientation and 'A' not in orientation: - orientation += 'A' - if 'I' not in orientation and 'S' not in orientation: - orientation += 'S' + if "L" not in orientation and "R" not in orientation: + orientation += "R" + if "P" not in orientation and "A" not in orientation: + orientation += "A" + if "I" not in orientation and "S" not in orientation: + orientation += "S" return orientation @@ -25,9 +26,9 @@ def orientation_to_affine(orientation, vxw=1, vxh=1, vxd=1): vx = np.asarray([vxw, vxh, vxd]) for i in range(3): letter = orientation[i] - sign = -1 if letter in 'LPI' else 1 - letter = {'L': 'R', 'P': 'A', 'I': 'S'}.get(letter, letter) - index = list('RAS').index(letter) + sign = -1 if letter in "LPI" else 1 + letter = {"L": "R", "P": "A", "I": "S"}.get(letter, letter) + index = list("RAS").index(letter) affine[index, i] = sign * vx[i] return affine @@ -52,10 +53,10 @@ def make_compressor(name, **prm): if not isinstance(name, str): return name name = name.lower() - if name == 'blosc': + if name == "blosc": Compressor = numcodecs.Blosc - elif name == 'zlib': + elif name == "zlib": Compressor = numcodecs.Zlib else: - raise ValueError('Unknown compressor', name) + raise ValueError("Unknown compressor", name) return Compressor(**prm) diff --git a/tests/helper.py b/tests/helper.py index 242f3d8..eb865d7 100644 --- a/tests/helper.py +++ b/tests/helper.py @@ -39,6 +39,7 @@ def _cmp_zarr_archives(path1: str, path2: str) -> bool: print(f"Mismatch found in dataset: {key}") return False if zarr1[key].attrs != zarr2[key].attrs: + print("attrs mismatch") return False # If all checks pass