Large Array Data

Neuroinformatics Open Software Summer School 2026

Introduction

Setting up

Open these slides on your laptop

neuroinformatics.dev/course-large-array-data

They contain links to exercises, polls and further references.

Also, activate your (conda) environment and install the helper functions for the course

conda activate large-arrays

uv pip install git+https://github.com/neuroinformatics-unit/course-large-array-data

Download the exercises

The instructor team

Alessandro Felder

Igor Tatarnikov

Kimberly Meechan

Harry Carey

Jim Bednar

Mentimeter

Mentimeter

Overall objectives

• Understand computational challenges related to large array data and key strategies to address them
• Gain hands-on Python experience and increase confidence in handling large array data
• Appreciate the diversity of the open-source Python ecosystem for large array data and how you are part of it

Schedule

• Day 1: fundamentals, motivation, reading and writing large arrays
• Day 2: processing large arrays in parallel, standards for large imaging data
• Day 3: advanced topics and case studies: compression, benchmarking, visualisation, real-life applications
• Day 4-5: collaboration

What is an array?

• Data arranged in a grid
• Data all have the same type

Why does data come in array form?

• Naturally lends itself to what the data represents (e.g. microscopy)
• Arrays can be indexed efficiently in memory

Data has a “size”

Computers store data as bits , arranged in bytes of 8 bits.

• Kilobyte (KB) = 1000 bytes
• Megabyte (MB) = 1000 kilobytes
• Gigabyte (GB) = 1000 megabytes
• Terabyte (TB) = 1000 gigabytes

Data has a “size” (sidenote)

Computers store data as bits , arranged in bytes of 8 bits. There are also binary system equivalents

• Kibibyte (KiB) = 1024 bytes
• Mebibyte (MiB) = 1024 kibibytes
• …

We treat KB and KiB as equivalent for our purposes.

Why do we care about arrays?

Lots of scientific data comes in (large) array form.

Electrophysiological probes

Electrophysiological probe recordings

Brain Wide Map

Brain Wide Map

Brain Wide Map

49.7 TB across 2048 files and 139 subjects

Week-long recordings

Week-long recordings

In a typical experiment (for example the ones described in the section Recording and quantification of social behaviours), we record from over 30 unique hardware devices and over 60 unique data streams, and when including electrophysiology data from two Neuropixels probes, the total acquisition throughput can reach up to 200 GB per hour.

Two-photon microscopy

Whole brain imaging

A whole marine worm

Volume electron microscopy

A whole marine worm

(~2TB)

Geography

Astronomy

Mentimeter

Mentimeter

Mentimeter

Starting small

Exploring array size and shape with numpy

numpy arrays

numpy is the default Python library to process array data.

Common data types

NumPy supports a very wide range of data types, but there are a few that are most common in scientific applications:

NumPy datatype	Full name	Range of values
bool	Boolean (8-bit)	[False, True]
uint8	Unsigned integer 8-bit	0…255
uint16	Unsigned integer 16-bit	0…65535
float32	Float 32-bit	\(-3.4 \times 10^{38}...+3.4 \times 10^{38}\)
float64	Float 64-bit	\(-1.7 \times 10^{308}...+1.7 \times 10^{308}\)

Exploring a numpy array

import numpy as np

array_shape = (5,6)
dtype = np.uint16
array = np.random.randint(
  low=0,
  high=100,
  size=array_shape,
  dtype=dtype)

print(array)
print("Array shape: "+str(array.shape))

Exploring a numpy array

[[32 43 65 95 37 88]
 [ 4 34 82 30 26 22]
 [69 78 42 87 56 81]
 [ 5 41 48 63 92 53]
 [72 56 62 23 33 51]]
Array shape: (5, 6)

Exploring a numpy array

Extract its bytes.

array_bytes = array.tobytes()

print(f"{len(array_bytes)} bytes total")
print()
print("First five elements:")
for i in range(0, 10, 2):
    word = int.from_bytes(array_bytes[i:i+2], byteorder='little')
    print(f'{word:016b}', end=' ')

print("...")

Exploring the array size

Extract its bytes

60 bytes total

First five elements:
0000000000100000 0000000000101011 0000000001000001 0000000001011111 0000000000100101 ...

numpy is great!

It works well for many applications and is very stable.

Using numpy keeps your code simple, more accessible, and easy-to-read.

So what’s the problem?

A mental model

Data being processed on a computer.

A mental model

Data being processed on a computer.

A mental model

Data being processed on a computer.

A mental model

Data being processed on a computer.

Memory allocation

When we make a numpy array, the memory it needs is allocated eagerly (i.e. immediately).

This is usually fine, but…

Large Array Data

…vastly exceeds an average laptop’s memory (8-32 GB) in size.

Therefore, allocating memory eagerly is a bad idea!

A mental model

Data being processed on a computer.

Exercise 1 (Eager memory allocation)

• notebooks/1_memory_cap.ipynb
• Activate your conda environment

conda activate large-arrays

• Run jupyter

jupyter lab

• Open http://localhost:8888/lab?
• Navigate to the downloaded notebook and work through it.

Exercise 1 (demo)

Coffee break?

A mental model

Eager memory allocation is a problem for large arrays.

We need more control

“Lazy” operations:

Hold off allocating memory/doing computations for as long as we can…

…and only operate on the minimum amount of in-memory data needed.

We get more control by being lazy :)

An example

Our goal is to plot a slice of a 3D brain image.

Eager strategy:

• Load the whole image into memory
• Access the slice
• Plot the slice

An example

Our goal is to plot a slice of a 3D brain image on disk.

Lazy strategy:

• Store the index of the slice in-memory
• Figure out where on disk to get the slice from
• Only when plotting, access the slice

numpy can do this

A memory map lets us treat a file on disk as if it were in memory.

The operating system loads parts of it from disk into memory only as they’re actually accessed.

numpy can do this

import numpy as np

shape = (100, 512, 512)
dtype = np.uint16

arr = np.memmap("volume.raw", dtype=dtype, mode="r", shape=shape)

slice_index = 50
print(arr[50, :, :]) # only necessary "pages" loaded into memory

Limitations of numpy memory maps

• Data has to be on local disk (not cloud)
• No compression
• No parallelisation

How could we improve on this?

How could we improve on this?

We would benefit from increased flexibility like

• Allowing large arrays to be cut into smaller, customisable “chunks” that can live anywhere
• Allow control over when operations are executed

How could we improve on this?

dask has solved this problem for us.

The dask array

An array of numpy arrays

that also behaves like a lazy numpy array.

What does this mean?

What does this mean?

The dask array

An evolution of the numpy array

• Same interface (for most functionality)
• Allows “chunking” of the array into smaller chunks
• Gives control over when data is actually loaded into memory

dask keeps track of lazy operations in a “task graph”

The dask array

import dask.array as da
rng = da.random.default_rng()
x = rng.random((10000, 10000), chunks=(1000, 1000))
print(x)

The dask array

dask.array<random, shape=(10000, 10000), dtype=float64, chunksize=(1000, 1000), chunktype=numpy.ndarray>

Lazy evaluation means print(x) doesn’t load array data into memory.

The dask array

import dask.array as da
rng = da.random.default_rng()
x = rng.random((10000, 10000), chunks=(1000, 1000))
print(x[0:100, 0:100])

Lazy evaluation means print(x) and slicing still doesn’t load array data into memory.

The dask array

dask.array<getitem, shape=(100, 100), dtype=float64, chunksize=(100, 100), chunktype=numpy.ndarray>

The dask array

import dask.array as da
rng = da.random.default_rng()
x = rng.random((10000, 10000), chunks=(1000, 1000))
print(x[0:100, 0:100].compute())

.compute method gives us control over when to actually load into memory.

The dask array

[[0.03253803 0.2614327  0.64681631 ... 0.49464121 0.46261251 0.17876863]
 [0.01697999 0.6643703  0.58008443 ... 0.01420531 0.90376143 0.31415469]
 [0.95924724 0.16851846 0.09130982 ... 0.35808623 0.82692309 0.88373208]
 ...
 [0.68139916 0.6499788  0.85691498 ... 0.96517229 0.01258917 0.35944043]
 [0.27329427 0.40688852 0.66911353 ... 0.17843952 0.24625911 0.67287817]
 [0.359546   0.85340175 0.92695889 ... 0.3839975  0.44808929 0.62934189]]

.compute method gives us control over when to actually load into memory.

dask delayed

import numpy as np
from dask import delayed

def for_loop_sum(array):
  result = array[0]
  for a in array[1:]:
    result += a
  return result

array = np.random.randint(
  low=0,
  high=100,
  size=(100,100),
)
print(delayed(for_loop_sum(x[0:100, 0:100]))) # lazy
print(for_loop_sum(x[0:100, 0:100])) # eager

delayed function allows us to execute our own functions lazily.

dask delayed

Delayed('ndarray-6639df84-fc12-4b55-b394-c2452a27e7d0')
[10713 10070  9443 10340  8764 10698  9925  9931  9548  9109  9852  9778
  9391  9987  9454  9959 10309 11093 10239 10865  8978  8706  9623  9754
  9861  9302  9413  8719  9348 11095 10199 11380 10790  9714 10914  8905
  9777 10032 10765  9924  9883 10885 10268  9401  9282  9391  9267 10121
  9658  9377  9686 10151  9805  9766 10398 10098  9550  9668  9545  9832
  9912  9328  9486  9309 10835  9481  9830  9937 10057  9875  9994 10111
 10204  9341  9334  9370  9709 10305  9290  9037 10736 10070  9729  8456
  9310  9280  9567  9878 10244  8613  9541 10346 10905 10115 10486  9773
 10633  9783 10230 10062]

delayed function allows us to execute our own functions lazily.

Exercise 2 (Joint)

Loading of large array slices with

• numpy (eager, via tifffile.imread)
• numpy.memmap (lazy, via tifffile.memmap)
• dask

showcasing the increased flexibility of dask.

Exercise 2 (Joint)

If you want to run on your own data

• notebooks/2_numpy_to_dask.ipynb
• Activate your conda environment

conda activate large-arrays

• Run jupyter

jupyter lab

• Open http://localhost:8888/lab?
• Navigate to the downloaded notebook and work through it, make it point to your own data

Reading and writing dask arrays

So far, we have only created dask arrays in-memory (lazily). But how do we store them on disk?

• Having chunks helps us here too.
• We can store each chunk in a separate file.
• We can store metadata in the same (super-)file.

How do we store them on disk?

zarr has solved this problem for us.

Zarr files

• open-source file format specification
• n-dimensional arrays
• split into chunks
• Each chunk is a separate file

Note: Zarr and zarr-python

Zarr is not the same as zarr-python.

• Zarr is a file format specification
• zarr-python is a Python library to work with Zarr files
  • Confusingly, used like

import zarr

...

zarr_array = zarr.open(path)

Exercise 3

Exercise 3

Showcases

• dask reading a file from a public cloud archive
• dask processing the data and storing part of it to zarr

Exercise 3

• notebooks/3_dask_remote_to_zarr_io.ipynb
• Activate your conda environment

conda activate large-arrays

• Run jupyter

jupyter lab

• Open http://localhost:8888/lab?
• Navigate to the downloaded notebook and work through it, make it point to your own data

Summary

• numpy eager memory allocation is a bad idea for large arrays
• numpy.memmap helps but has limitations

Summary

• dask is more flexible
  • Lazy by default
  • Customise chunks
  • Control function execution with delayed
  • Works with remote data
• zarr
  • Flexible file format for chunked arrays
  • Each chunk = one file on disk

OME-Zarr introduction

Recap on Zarr

• open-source file format specification
• n-dimensional arrays
• split into chunks
• Each chunk is a separate file
• Enables parallel processing of massive datasets

Working with Zarr

• One of the main libraries for interacting with Zarr in python is zarr-python.

• Quick demo of using zarr-python to make arrays and groups: notebooks/zarr_python_example.ipynb

What’s in a Zarr?

The zarr specification dictates how a Zarr file should be structured.

As we saw in the demo, the main components of a Zarr file are:

• groups (which you can think of as containers) and arrays (our actual numeric data)
• groups can contain any combination of groups and arrays, allowing a tree of datasets to be made
• each level of the tree (a group or array) contains a zarr.json file with basic metadata

Tree / hierarchy

Zarr is general (pros)

One of Zarr’s great strengths is that it can be applied to any kind of array data across many fields.

Their website’s dataset section includes:

• climate model outputs
• sea surface temperature
• bio-imaging data
• …

Zarr is general (cons)

• Equally, one of Zarr’s great weaknesses is that it is general!

• When we use it for a specific kind of data (say biological images), there is no consistent place for related metadata to be stored.

• For example:

  • pixel size?
  • units?
  • meaning of dimensions? E.g. time / z / y / x / channels

This makes it harder to share datasets with others, and to process datasets that weren’t created by yourself.

OME-Zarr

Specific fields can get around this limitation by extending the Zarr specification.

• OME-Zarr = Zarr + bio-imaging related metadata
• It also introduces conventions on how to layout the Zarr files to handle common use cases:
  • multi-resolution datasets (we’ll come back to this!)
  • high-content screening
  • storing label images alongside raw images

Working with OME-Zarr

Let’s take a look at some OME-Zarr metadata using ome-zarr-models.

• OME-Zarr models is a python package specialised to validating and loading OME-Zarr metadata.

• Quick demo of accessing metadata: notebooks/ome-zarr-models_example.ipynb

Multi-scales

Part of the OME-Zarr specification adds multiscales metadata.

• It allows storing the same dataset at multiple resolutions
• May also hear this called a ‘pyramidal’ file format

Why multi-scale?

Viewing a whole slice from a multi-terabyte dataset can still be hundreds of GB!

Specialised viewers

Specialised viewers dynamically load different resolution levels.

• Web examples:

  • neuroglancer - human organ atlas kidney
  • webknossos - segmented connectomics dataaset

• Local examples:

  • Napari
  • BigDataViewer (Fiji)
  • …

Versioning

Spec versions can get complicated.

• The main zarr spec is currently v3
• The OME-Zarr spec is currently v0.5
• OME-Zarr changes more often than Zarr

News from OME-Zarr land

• follow the Announcements for the ome-ngff tag on image.sc
• Join a community call or office hours - community pages
• ngff docs

Ecosystem

The ecosystem of tools / packages around OME-Zarr is rapidly evolving.

Don’t be afraid to try different options, and see what works best for your data.

• The resources section has useful lists of tools, training materials etc for OME-Zarr
• Tools list
• Training materials

Try it out!

• notebooks/creating_ome_zarr.ipynb: creating OME-Zarr step by step
• notebooks/converting_ome_zarr.ipynb: libraries for converting to OME-Zarr
• notebooks/visualising_ome_zarr.ipynb: viewing OME-Zarr in napari

Parallelising array operations with dask

Parallelising array operations with dask

How should we coordinate array operations so Python can execute them at the same time?

• “In parallel”, “concurrently”

Concurrency

Two types of concurrency built into Python

• “Multi-threading”
• “Multi-processing”

Multi-threading

• Single Python process
• Threads share resources (e.g. memory)

Multi-threading example

import threading

counter = 0
lock = threading.Lock()

def increment(n):
    global counter
    for _ in range(n):
        with lock: # make sure counter doesn't get updated simultaneously
            counter += 1

threads = [threading.Thread(target=increment, args=(100_000,)) for _ in range(4)]

with ThreadPoolExecutor(max_workers=4) as pool:
    futures = [pool.submit(increment, 100_000) for _ in range(4)]
    for f in futures:
        f.result()

print(counter)

Multi-threading example

400000

Multi-threading with Queue

import threading
import queue

q = queue.Queue()

def worker():
    while True:
        item = q.get()
        print(f'Working on {item}')
        print(f'Finished {item}')
        q.task_done()

# Turn-on the worker thread.
threading.Thread(target=worker, daemon=True).start()

# Send thirty task requests to the worker.
for item in range(30):
    q.put(item)

# Block until all tasks are done.
q.join()
print('All work completed')

https://docs.python.org/3/library/queue.html

Multi-threading

• Single Python process
• Threads share resources (e.g. memory)
• Need for coordination of sharing memory may actually hurt performance

Multi-processing

• Main Python process spawns other processes
• No shared resources
  • Each process has its own memory space
  • Memory needs to be copied to the subprocess
  • Overhead of launching extra processes

Multi-processing example

import multiprocessing as mp

def worker(xs):
    print(len(xs))
    return sum(xs)

chunks = [[1,2],[3,4],[5,6],[7,8]]
ctx = mp.get_context("fork")
with ctx.Pool(4) as pool:
   partial_sums = pool.map(worker, chunks)

print(partial_sums)

total = sum(partial_sums)
print(total)

Multi-processing example

2222



[3, 7, 11, 15]
36

Python concurrency

Python has concurrency built-in.

• Threads share memory
• Processes don’t

Use Queue for safe communication.

But writing concurrent code is complex and error-prone!

Dask provides some abstractions so you don’t have to worry about the intricacies.

Processing chunks

We can split operations on arrays into two types:

1. The output element depends only on its input element (e.g threshold)
2. The output element also depends on some neighbouring elements (e.g. Gaussian filter)

Processing chunks

We can split operations on arrays into two types:

1. The output element depends only on its input element (e.g threshold)
2. The output element also depends on some neighbouring elements (e.g. Gaussian filter)

Neighbouring elements may live in a different chunk!

Operating on dask/numpy arrays

Some extra dask/numpy syntax that will come in handy

• where
• Boolean comparison of arrays
• nanmean

Operating on dask chunks

You can tell dask to apply a function lazily on each chunk

• map_block (each chunk has all the info)
• map_overlap (you need info from neighbouring chunks)

Exercise 4

Exercise 4 will guide you through a realistic analysis of multiphoton imaging data that uses chunkwise dask operations

• Without overlap
• with overlap

Also showing how to choose between multi-threading and multi-processing.

Exercise 4

• notebooks/4_parallel_processing_with_dask.ipynb
• Activate your conda environment

conda activate large-arrays

• Run jupyter

jupyter lab

• Open http://localhost:8888/lab?
• Navigate to the downloaded notebook and work through it.

More notes on dask

• Scales to multiple compute nodes
• Various extra libraries for GPU support
• Very general

More notes on performance

Favour robustness and simplicity over speed

• Performance may differ depending on data, hardware…
• Optimising performance is very tricky

Other Python array libraries

• Xarray: label arrays with e.g units, dimensions, …
• pytorch: aimed at machine-learning, native GPU support
• JAX: aimed at machine-learning
• cupy: array library for Python on nvidia GPU
• …

Array API

Goal: abstract the array backend

def unit_vector_array_api(x: AnyArray) -> AnyArray:
    """Works with any supported array backend."""
    # NumPy, CuPy, PyTorch, JAX, ...
    xp = x.__array_namespace__()
    return x / xp.sqrt(xp.sum(x * x))

Array API

Array API

Goal: work with different array backends

def unit_vector_array_api(x: AnyArray) -> AnyArray:
    """Works with any supported array backend."""
    # NumPy, CuPy, PyTorch, JAX, ...
    xp = x.__array_namespace__()
    return x / xp.sqrt(xp.sum(x * x))

Array API

• Work-in-progress
• Some difficulties in unifying array libraries
• Choose for flexibility over performance

Summary

• Python has threads and processes
• Use map_blocks to apply chunkwise operations
• Use map_overlap if you need information from neighbouring chunks,
  • you will need to think about the depth parameter
• By default, dask Arrays operate concurrently using multi-threading
• Favour robustness and simplicity over performance
• There are other open-source Python array libraries

Compression

Data has a “size”

Computers store data as bits.

60 bytes total

First five elements:
0000000000100000 0000000000101011 0000000001000001 0000000001011111 0000000000100101 ...

Data compression

We can be smart about how to store this data with fewer bits.

60 bytes total

First five elements:
0000000000100000 0000000000101011 0000000001000001 0000000001011111 0000000000100101 ...

Data compression

We could “encode” something like “assume the first 8 bits are always 0” and halve the file size.

60 bytes total

First five elements:
0000000000100000 0000000000101011 0000000001000001 0000000001011111 0000000000100101 ...

Lossless and lossy compression

Some compression algorithms preserve all the information.

Others are lossy, e.g. “.jpg” files

What is a codec?

Codec stands for encoder/decoder.

If defines how data is written to/read from disk.

Different codecs will work better/worse for different data.

Compression trade-offs

Generally, we expect a smaller size on disk to result in slower read/write.

But compression algorithms can be clever about this too.

Case study

Two audio codecs work well on electrophysiology data.

Better compression ratio, slower decompression speed

Using Codecs in zarr

Zarr provides a number of Codecs that you can specify when writing files.

Exercise 5

Compares write speeds and disk size for different compression levels.

Exercise 5

• notebooks/5_remote_file_compression.ipynb
• Activate your conda environment

conda activate large-arrays

• Run jupyter

jupyter lab

• Open http://localhost:8888/lab?
• Navigate to the downloaded notebook and work through it, make it point to your own data

Zarr benchmarks

Recap on Zarr

• open-source file format specification
• n-dimensional arrays
• split into chunks
• Each chunk is a separate file
• Enables parallel processing of massive datasets
• OME-Zarr extends this with bioimaging-related metadata + multiscales

Zarr settings

There are lots of settings that can be adjusted when writing Zarr and OME-Zarr images:

• Chunk size
• Compressor
• Compression level
• Other compressor options like ‘shuffle’
• …

How to choose?

These settings affect:

• Write time

• Read time

• Size on disk

• Number of files

• Memory usage for different operations

• Choosing settings is a trade-off between these different factors. You have to decide what you want to prioritise for your data.

HEFTIE Project

Benchmarks were made as part of the HEFTIE project:

• David Stansby, Ruaridh Gollifer and Kimberly Meechan (UCL ARC)
• Also made the online OME-Zarr book that we used some exercises from yesterday.

Benchmarks

• Benchmarks repository

• Benchmarks report

• 3 images: heart (335 MB), dense segmentation, sparse segmentation

• Images used are on zenodo

• All sized 806 x 629 x 629

• All 16-bit unsigned integer

Images

Heart: HiP-CT scan of a heart from the Human Organ Atlas

Images

Dense: segmented neurons from electron microscopy

Images

Sparse: A few select segmented neurons from electron microscopy

Benchmarks

• pytest-benchmark
• Vary different settings and record:
  • Time to write the whole array (lower is better)
  • Time to read the whole array (lower is better)
  • Compression ratio: data size in memory / data size on disk (higher is better)
• Each run 5 times and mean values used

Write time

Higher compression level = higher compression ratio = longer write time

Read time

Higher compression level doesn’t usually result in longer read times. For many compressors this is a feature of their design, with a large one-off cost of compressing the data, but no slow down when reading it

Chunk size

Larger chunks = lower compression ratio (but only slightly)

Chunk size

Larger chunks = faster write and read times (mainly due to fewer total files)

Image type

Heart image Dense segmentation Sparse segmentation

Segmentations compress much more!

Software library

Which software library you use affects read / write times. Tensorstore had the best performance.

Benchmarking Zarr conclusions

Some general recommendations:

• Use tensorstore library for fastest read/write

• blosc-zstd is a good default compressor

• Use a high compression level to get the smallest file size (usually means longer write times, but not much effect on read time)

• Smaller chunks = smaller overall file size + less memory usage

• Larger chunks = faster read/write times + fewer files

• It’s a balance! People often use 64x64x64 or 128x128x128

• This is different for different data! Worth testing some small samples of your own data with different settings

Sharding

One of the settings we didn’t benchmark was ‘sharding’:

• Shards contain multiple chunks
• Each shard is a file on disk (rather than each individual chunk)
• A chunk is the smallest unit you can read; a shard is the smallest unit you can write
• Main purpose is to reduces number of files
• Some file systems may have performance issues storing and accessing a large number of small files.

Benchmarking sharding

Let’s experiment with sharding options:

• notebooks/sharding.ipynb

Large Array Data

Case studies from the Neuroinformatics Unit

Photon-mosaic

Two-photon microscopy

Live multi-photon imaging data

• Very long arrays in time
• Can be 2D or 3D in space
• Can have several channels

Context

• Some existing open-source software packages (e.g. suite2p, CaImAn)
• Not particularly interoperable
• Experiments becoming more complex:
  • Match videos from several days
  • More 3D data
  • From deeper in the brain

We need flexibility to customise!

Modularity

Vision: Mix and match like

import photon_mosaic as pm

pm.correct_motion(method="suite2p")
pm.rois_caiman = extract_rois(method="caiman")

# for comparison
pm.rois_suite2p = extract_rois(method="suite2p")

compare(rois_caiman, rois_suite2p)

A project in its very early stages

Modularity

Demo time!

Mental models

• Shared representation of arrays and processing steps
• Inspired by spikeinterface

Array handling strategy

• Only chunk along time
• Rely on ROIExtractors to read (lazily)
• We support native writing to Zarr and binary.

photon-mosaic-pipeline for batch-processing on HPC

• Built on top of snakemake.
• Configure with a single file.

Community

• Collaboration with Allen Institute for Neural Dynamics and International Brain Laboratory
• Community-driven: open meetings every two weeks!

Alessandro Felder

Laura Porta

Alessio Buccino

Arielle Leon

Johannes Friedrich

Sean McCulloch

Georg Raiser

Photon-mosaic

• photon-mosaic.org
• github.com/photon-mosaic
• Zulip channel

Photon-mosaic

Collaboration day ideas? Open an issue.

brainglobe-atlasapi

BrainGlobe Initiative

Established 2020 with three aims:

1. Develop general-purpose tools to help others build interoperable software for computational neuroanatomy.
2. Develop specialist software for specific analysis and visualisation needs.
3. Reduce barriers of entry, and facilitate the building of an ecosystem of computational neuroanatomy tools.

BrainGlobe atlases

Problem - the data analysis tools are fragmented

• Model species
• Imaging modality
• Anatomical focus
• Developmental stage

BrainGlobe atlases

brainglobe-atlasapi

BrainGlobe atlases

BrainGlobe atlases

Atlases are made up of 4 main components:

1. Common coordinate space - text based metadata
2. Template - uint16 image
3. Annotations - uint32 image and meshes for each region
4. Terminology - text based metadata

brainglobe-atlasapi V2

• Images stored as .tiff files
• Meshes stored as .obj files
• Distributed as a single .tar.gz file
• No reuse of data between atlases

brainglobe-atlasapi V3

• Images stored as OME-Zarr files
• Meshes stored in a precomputed formats
• Each component stored and downloaded lazily
• Data is reused between atlases

brainglobe-atlasapi V3

brainglobe-atlasapi annotations

• Multi-resolution 3D image for per voxel annotation
• Mask for each region (798, 1320, 800, 1140)
  • ~900 GB uncompressed for 10um resolution!!
• Compressed with Blosc zstd codec
  • Only ~180 MB for 10um, 25um, 50um, and 100um resolutions

Cellfinder

Cellfinder

• part of the BrainGlobe ecosystem
• real-world example of many concepts demonstrated here

Cellfinder

Finding centre coordinates of cells (bright blobs) in large arrays.

Serial two-photon tomography

cellfinder input data

cellfinder input data

Cellfinder

Finding centre coordinates of cells (bright blobs) in large arrays.

Arrays are ~100GB / channel, stored as 2D tiffs (40 MB each)

How to find cells

• A series of 2D and 3D filtering and thresholding operations

• Store centres of bright blobs (“cell candidates”)

• Use a classification algorithm

  • Works on small cube around candidate centre

Strategy

• Sweep along first array axis
  • Keep a “slab” in memory
• Filter and threshold
  • Store centres of blobs separately (small)
• Sweep again
  • Extract cubes and add classification metadata

First sweep (filter+threshold)

First sweep (filter+threshold)

• 3D operations in main thread
• Communicate both ways with sub-threads/processes
  • Using two queues
• Custom exception handling

Second sweep (classify)

Second sweep (classify)

• cache whole slab
  • for reuse across cell centres

Code examples

• Taken without edit directly from cellfinder codebase
• See if you understand them!

Lazy loading

lazy_imread = delayed(tifffile.imread)  # lazy reader

def read_with_dask(path: str | Path) -> da.Array:

    ...

    shape, dtype = get_tiff_meta(filenames[0])
    lazy_arrays = [
        lazy_imread(fn, is_ome=False)
        for fn in get_sorted_file_paths(filenames)
    ]
    dask_arrays = [
        da.from_delayed(delayed_reader, shape=shape, dtype=dtype)
        for delayed_reader in lazy_arrays
    ]
    stack = da.stack(dask_arrays, axis=0)
    return stack

2D filtering in a sub-process

@inference_wrapper
def _plane_filter(
    process: ProcessWithException,
    tile_processor: TileProcessor,
    n_threads: int,
    buffers: List[Tuple[torch.Tensor, torch.Tensor]],
):
    ...

    while True:
        msg = process.get_msg_from_mainthread()
        if msg == EOFSignal:
            return
        # with torch multiprocessing, tensors are shared in memory - so
        # just update in place
        token, i = msg
        tensor, masks = buffers[token]

        plane, mask = tile_processor.get_tile_mask(tensor[i : i + 1, :, :])
        tensor[i : i + 1, :, :] = plane
        masks[i : i + 1, :, :] = mask

        # tell the main thread we processed all the planes for this tensor
        process.send_msg_to_mainthread(None)

Distributing data across subprocesses

def start_dataset_thread(self, num_workers: int) -> None:

    # include queue for host thread
    ctx = mp.get_context("spawn")
    # we use maxsize=0 to prevent potential locking issues. But, we never
    # actually request more than one data batch at a time over a queue
    queues = [ctx.Queue(maxsize=0) for _ in range(num_workers + 1)]
    self._worker_queues = queues

    self._dataset_thread = ThreadWithExceptionMPSafe(
        target=_read_data_send_cuboids,
        args=(self.src_image_data, queues),
        pass_self=True,
    )
    self._dataset_thread.start()

Cellfinder summary

• dask for lazy loading
• pytorch for processing
  • suited to ML
  • suited to GPU
• Coordinate via Python concurrency

Cellfinder efficiently processes large arrays.

Cellfinder summary

Cellfinder efficiently processes large arrays.

Cellfinder

Collaboration day ideas? Open an issue.

Overall summary

Overall objectives

• Understand computational challenges related to large array data and key strategies to address them
• Gain hands-on Python experience and increase confidence in handling large array data
• Appreciate the diversity of the open-source Python ecosystem for large array data and how you are part of it

Understand computational challenges and how to address them

Large arrays do not fit into memory

We can deal with this through

• Chunking
• Lazy evaluation
• Parallelising across chunks
• Compressing chunks on disk

Handling large array data with Python

• zarr-python and dask are generalist libraries for processing large arrays (in parallel, with compression, anywhere).
• Use zarr files, ideally following a community standard (e.g. OME-Zarr).

The open-source ecosystem

The open-source ecosystem

• There are many array libraries and related open-source projects out there.
  • They also depend on layers of generalist open-source libraries.
• Each have their niche.
• Including in (neuro)science.
• By using them, you are part of it!
• Hopefully, interact more during collaboration days!

Feedback

Please provide your honest feedback :)

Thank you!

References

Aeon. 2026. Neuropixels Recording Guide. https://aeon.swc.ucl.ac.uk/getting_started/npx_recording/#id1.

Buccino, Alessio P, Olivier Winter, David Bryant, David Feng, Karel Svoboda, and Joshua H Siegle. 2023. “Compression Strategies for Large-Scale Electrophysiology Data.” Journal of Neural Engineering 20 (5): 056009. https://doi.org/10.1088/1741-2552/acf5a4.

Claudi, Federico, Luigi Petrucco, Adam L. Tyson, Tiago Branco, Troy W. Margrie, and Ruben Portugues. 2020. “BrainGlobe Atlas API: A Common Interface for Neuroanatomical Atlases.” Journal of Open Source Software 5 (54): 2668. https://doi.org/10.21105/joss.02668.

HP Tech Takes. n.d. How Much RAM Do I Need? A Capacity Guide. HP. Accessed July 14, 2026. https://www.hp.com/us-en/shop/tech-takes/how-much-ram-do-i-need-every-use.

International Brain Laboratory. 2026. IBL - Brain Wide Map. https://doi.org/10.48324/dandi.000409/0.260309.1324.

Juavinett, Ashley L., George Bekheet, and Anne K. Churchland. 2019. “Chronically Implanted Neuropixels Probes Enable High-Yield Recordings in Freely Moving Mice.” eLife 8: e47188. https://doi.org/10.7554/eLife.47188.

Jun, James J., Nicholas A. Steinmetz, Joshua H. Siegle, et al. 2017. “Fully Integrated Silicon Probes for High-Density Recording of Neural Activity.” Nature 551 (7679): 232–36. https://doi.org/10.1038/nature24636.

Osten, Pavel, and Troy W. Margrie. 2013. “Mapping Brain Circuitry with a Light Microscope.” Nature Methods 10 (6): 515–23. https://doi.org/10.1038/nmeth.2477.

Sirmpilatze, Nikoloz, Alessandro Felder, Dinora Abdulazhanova, et al. 2025. Mapping the Magnetoreceptive Brain: A 3D Digital Atlas of the Migratory Bird Eurasian Blackcap ( Sylvia Atricapilla ). https://doi.org/10.1101/2025.03.04.641293.

Stansby, David, Kimberly Meechan, and Ruaridh Gollifer. 2026. An Introduction to OME-Zarr for Big Bioimaging Data. https://ome-zarr-book.readthedocs.io/.

Vergara, Hernando M., Constantin Pape, Kimberly I. Meechan, et al. 2021. “Whole-Body Integration of Gene Expression and Single-Cell Morphology.” Cell 184 (18): 4819–4837.e22. https://doi.org/10.1016/j.cell.2021.07.017.