Handling memory-mapped images

[tlnagy]

Hi @mkitti, thanks for the great package! I was wondering if you had an tips for getting memory-mapped TiffImages.jl images to work with Napari.jl. It would be great if I could pass a memory-mapped TIFF object to Napari.jl and it could request data per-slice instead of all at once.

[mkitti]

What is the current difficulty? What would you like it to look like from the Python?

I would take a look at the approach I took in https://github.com/mkitti/NumPyArrays.jl. Generally, if you have a memory mapped TIFF image, you should be able to obtain a pointer to the data that can be then wrapped into a NumPy-like array.

[mkitti]

Taking a closer look, what you probably want to do is override Napari.view_image, Napari.add_image, and/or anNapari.Image.primitive_array for your type. See
https://github.com/mkitti/Napari.jl/blob/main/src/image.jl

[tlnagy]

Given a XYZ image, the tricky part is that TiffImages loads each XY slice independently so I would need to be able to overload a function that lets me load each XY slice as it is needed. Is that possible?

E.g. As the user drags the Z slider, call a function that loads the corresponding XY slice. Looks like Napari supports this in Python with dask/zarr

Reinterpreting and permuting the dims should be fine.

[mkitti]

That might require more support on the Python side if you want to make the I/O completely lazy. I have not looked at what the update mechanism is that they are using for dask or zarr, but we would basically need a Python class that forwards indexing requests to Julia's getindex. A "jlarray". We would need to look into what API Napari is using to work with dask or zarr.

[mkitti]

I'm still confused what the current challenge is though. You have a memory mapped array. This is already lazy.

So you are memory mapping each XY slice. For each XY slice you essentially have a pointer, and that pointer can be passed to NumPy which will represent it as an array. What is currently preventing that from happening?

[tlnagy]

Sorry, I'm returning back to this now that tlnagy/TiffImages.jl#79 has been merged and TiffImages.jl has true mmap capabilities. This should make a lot of the stuff I was talking about above irrelevant. A memory-mapped TIFF is now a vector of chunks which are the different TIFF slices on disk.

The type of these chunks are (for an example 1024x1024 Gray{N0f16} image)

julia> typeof(diskimg.chunks[1])
Base.ReinterpretArray{Gray{N0f16}, 2, UInt8, Base.ReshapedArray{UInt8, 3, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{UInt64}}, true}, Tuple{}}, true}

which is basically from inner to outer -> the raw bytes on disk -> reshaped to a (2, 1024, 1024) array -> reinterpreted as Gray{N0f16}. This should definitely work with Napari.jl.

Currently, trying to call @view_image on the chunks gives me:

julia> @view_image diskimg.chunks[1]
ERROR: MethodError: no method matching NumPyArray(::Base.ReinterpretArray{Gray{N0f16}, 2, UInt8, Base.ReshapedArray{UInt8, 3, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{UInt64}}, true}, Tuple{}}, true})
Closest candidates are:
  NumPyArray(::StridedArray{T}, ::Bool) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:88
  NumPyArray(::PyCall.PyArray{T, N}) where {T, N} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:75
  NumPyArray(::AbstractArray{T}) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:87
  ...
Stacktrace:
 [1] view_image(::Base.ReinterpretArray{Gray{N0f16}, 2, UInt8, Base.ReshapedArray{UInt8, 3, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{UInt64}}, true}, Tuple{}}, true}; kwargs::Base.Pairs{Symbol, String, Tuple{Symbol}, NamedTuple{(:name,), Tuple{String}}})
   @ Napari.Layers.Image ~/.julia/packages/Napari/vaVWX/src/image.jl:220
 [2] top-level scope
   @ REPL[29]:1

Which looks like it isn't properly recursing down into the wrapped images (it tries to immediately cast the NumPyArray).

julia> NumPyArray(diskimg.chunks[4].parent)
ERROR: conversion to pointer not defined for SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{UInt64}}, true}
<snip>

julia> NumPyArray(diskimg.chunks[4].parent.parent)
ERROR: conversion to pointer not defined for SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{UInt64}}, true}
<snip>

julia> NumPyArray(diskimg.chunks[4].parent.parent.parent)
11334290172-element NumPyArray{UInt8, 1}:
 0x49
 0x49
 0x2b
 0x00
 0x08
 0x00
 0x00
 0x00
 0x10
 0x00
 0x20
 0x00
    ⋮
 0x2e
 0x6a
 0x6c
 0x20
 0x76
 0x30
 0x2e
 0x36
 0x2e
 0x31
 0x00

The last conversion appears to create a NumPy pointer to the entire TIFF image on disk. I'm having a hard time following why exactly the NumPyArray conversions are failing given that SubArrays, ReinterpretArrays, ReshapedArrays are all supported. Any help would be appreciated.

[tlnagy]

A quick MWE using tlnagy/exampletiffs

julia> using Napari, TiffImages

julia> img = TiffImages.load("exampletiffs/4D-series.ome.tif", mmap = true);

julia> @view_image img.chunks[1]
ERROR: MethodError: no method matching NumPyArrays.NumPyArray(::Base.ReinterpretArray{ColorTypes.Gray{FixedPointNumbers.Q0f7}, 2, UInt8, Base.ReshapedArray{UInt8, 2, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, true})
Closest candidates are:
  NumPyArrays.NumPyArray(::StridedArray{T}, ::Bool) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:88
  NumPyArrays.NumPyArray(::PyCall.PyArray{T, N}) where {T, N} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:75
  NumPyArrays.NumPyArray(::AbstractArray{T}) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:87
  ...
Stacktrace:
 [1] view_image(::Base.ReinterpretArray{ColorTypes.Gray{FixedPointNumbers.Q0f7}, 2, UInt8, Base.ReshapedArray{UInt8, 2, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, true}; kwargs::Base.Pairs{Symbol, String, Tuple{Symbol}, NamedTuple{(:name,), Tuple{String}}})
   @ Napari.Layers.Image ~/.julia/packages/Napari/vaVWX/src/image.jl:220
 [2] top-level scope
   @ REPL[16]:1

[mkitti]

The original problem is that you are using UInt64 rather Int64 to index to the array. Try using Int64 ranges instead to subindex instead.

Issue:
JuliaLang/julia#42120

Fix by me:
JuliaLang/julia#43262

[tlnagy]

Okay force casting the array indexing to use Int64 helps! Trying with the chunk still fails, but now we're a level up and the parent of the chunk now loads:

julia> @view_image diskimg.chunks[4]
ERROR: MethodError: no method matching NumPyArray(::Base.ReinterpretArray{Gray{N0f16}, 2, UInt8, Base.ReshapedArray{UInt8, 3, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, true})
Closest candidates are:
  NumPyArray(::StridedArray{T}, ::Bool) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:88
  NumPyArray(::PyCall.PyArray{T, N}) where {T, N} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:75
  NumPyArray(::AbstractArray{T}) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64} at ~/.julia/packages/NumPyArrays/vqSo5/src/NumPyArrays.jl:87
  ...
Stacktrace:
 [1] view_image(::Base.ReinterpretArray{Gray{N0f16}, 2, UInt8, Base.ReshapedArray{UInt8, 3, SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, true}; kwargs::Base.Pairs{Symbol, String, Tuple{Symbol}, NamedTuple{(:name,), Tuple{String}}})
   @ Napari.Layers.Image ~/.julia/packages/Napari/vaVWX/src/image.jl:220
 [2] top-level scope
   @ REPL[9]:1

julia> @view_image diskimg.chunks[4].parent
PyObject Viewer(axes=Axes(visible=False, labels=True, colored=True, dashed=False, arrows=True), camera=Camera(center=(0.0, 511.5, 511.5), zoom=0.9936035156249999, angles=(0.0, 0.0, 90.0), perspective=0.0, interactive=True), cursor=Cursor(position=(1.0, 1.0, 0.0), scaled=True, size=1, style=<CursorStyle.STANDARD: 'standard'>), dims=Dims(ndim=3, ndisplay=2, last_used=0, range=((0.0, 2.0, 1.0), (0.0, 1024.0, 1.0), (0.0, 1024.0, 1.0)), current_step=(1, 512, 512), order=(0, 1, 2), axis_labels=('0', '1', '2')), grid=GridCanvas(stride=1, shape=(-1, -1), enabled=False), layers=[<Image layer '(diskimg.chunks[4]).parent' at 0x7f6d1284ac40>], scale_bar=ScaleBar(visible=False, colored=False, ticks=True, position=<Position.BOTTOM_RIGHT: 'bottom_right'>, font_size=10.0, unit=None), text_overlay=TextOverlay(visible=False, color=array([0.5, 0.5, 0.5, 1. ]), font_size=10.0, position=<TextOverlayPosition.TOP_LEFT: 'top_left'>, text=''), overlays=Overlays(interaction_box=InteractionBox(points=None, show=False, show_handle=False, show_vertices=False, selection_box_drag=None, selection_box_final=None, transform_start=<napari.utils.transforms.transforms.Affine object at 0x7f6cfeb25d90>, transform_drag=<napari.utils.transforms.transforms.Affine object at 0x7f6cfeb25df0>, transform_final=<napari.utils.transforms.transforms.Affine object at 0x7f6cfeb25e50>, transform=<napari.utils.transforms.transforms.Affine object at 0x7f6cfeb25eb0>, allow_new_selection=True, selected_vertex=None)), help='', status='Ready', tooltip=Tooltip(visible=False, text=''), theme='dark', title='napari', mouse_move_callbacks=[<function InteractionBoxMouseBindings.initialize_mouse_events.<locals>.mouse_move at 0x7f6d1285a310>], mouse_drag_callbacks=[<function InteractionBoxMouseBindings.initialize_mouse_events.<locals>.mouse_drag at 0x7f6d12844790>], mouse_double_click_callbacks=[], mouse_wheel_callbacks=[<function dims_scroll at 0x7f6d271d7940>], _persisted_mouse_event={}, _mouse_drag_gen={}, _mouse_wheel_gen={}, keymap={'Shift': <function InteractionBoxMouseBindings.initialize_key_events.<locals>.hold_to_lock_aspect_ratio at 0x7f6d128443a0>, 'Control-Shift-R': <function InteractionBoxMouseBindings._reset_active_layer_affine at 0x7f6d12844ee0>, 'Control-Shift-A': <function InteractionBoxMouseBindings._transform_active_layer at 0x7f6d12844ca0>})

[tlnagy]

Forcibly calling Napari.Layers.Image.primitive_type seems to fix the problem:

julia> @view_image Napari.Layers.Image.primitive_array(diskimg.chunks[4])
PyObject Viewer(axes=Axes(visible=False, labels=True, colored=True, dashed=False, arrows=True), camera=Camera(center=(0.0, 511.5, 511.5), zoom=0.9936035156249999, angles=(0.0, 0.0, 90.0), perspective=0.0, interactive=True), cursor=Cursor(position=(1.0, 1.0, 0.0), scaled=True, size=1, style=<CursorStyle.STANDARD: 'standard'>), dims=Dims(ndim=3, ndisplay=2, last_used=0, range=((0.0, 1.0, 1.0), (0.0, 1024.0, 1.0), (0.0, 1024.0, 1.0)), current_step=(0, 512, 512), order=(0, 1, 2), axis_labels=('0', '1', '2')), grid=GridCanvas(stride=1, shape=(-1, -1), enabled=False), layers=[<Image layer 'Napari.Layers.Image.primitive_array(diskimg.chunks[4])' at 0x7f6d0d9d8a60>], scale_bar=ScaleBar(visible=False, colored=False, ticks=True, position=<Position.BOTTOM_RIGHT: 'bottom_right'>, font_size=10.0, unit=None), text_overlay=TextOverlay(visible=False, color=array([0.5, 0.5, 0.5, 1. ]), font_size=10.0, position=<TextOverlayPosition.TOP_LEFT: 'top_left'>, text=''), overlays=Overlays(interaction_box=InteractionBox(points=None, show=False, show_handle=False, show_vertices=False, selection_box_drag=None, selection_box_final=None, transform_start=<napari.utils.transforms.transforms.Affine object at 0x7f6cfea46b20>, transform_drag=<napari.utils.transforms.transforms.Affine object at 0x7f6cfea46b80>, transform_final=<napari.utils.transforms.transforms.Affine object at 0x7f6cfea46be0>, transform=<napari.utils.transforms.transforms.Affine object at 0x7f6cfea46c40>, allow_new_selection=True, selected_vertex=None)), help='', status='Ready', tooltip=Tooltip(visible=False, text=''), theme='dark', title='napari', mouse_move_callbacks=[<function InteractionBoxMouseBindings.initialize_mouse_events.<locals>.mouse_move at 0x7f6d0d99c700>], mouse_drag_callbacks=[<function InteractionBoxMouseBindings.initialize_mouse_events.<locals>.mouse_drag at 0x7f6d0d99c040>], mouse_double_click_callbacks=[], mouse_wheel_callbacks=[<function dims_scroll at 0x7f6d271d7940>], _persisted_mouse_event={}, _mouse_drag_gen={}, _mouse_wheel_gen={}, keymap={'Shift': <function InteractionBoxMouseBindings.initialize_key_events.<locals>.hold_to_lock_aspect_ratio at 0x7f6d0d9899d0>, 'Control-Shift-R': <function InteractionBoxMouseBindings._reset_active_layer_affine at 0x7f6d0d9898b0>, 'Control-Shift-A': <function InteractionBoxMouseBindings._transform_active_layer at 0x7f6d0d989820>})

[mkitti]

I guess that supposed I should use primitive_array before trying to use NumPyArray here:

Napari.jl/src/image.jl

Lines 216 to 223 in b2c1469

	view_image(img::PermutedDimsArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
	add_image(viewer, img::PermutedDimsArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
	view_image(img::SubArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
	add_image(viewer, img::SubArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
	view_image(img::Base.ReinterpretArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
	add_image(viewer, img::Base.ReinterpretArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
	view_image(img::Base.ReshapedArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
	add_image(viewer, img::Base.ReshapedArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)

I'm curious what the output of @which Napari.Layers.Image.primitive_array(diskimg.chunks[4]) is

[tlnagy]

julia> @which Napari.Layers.Image.primitive_array(diskimg.chunks[4])
primitive_array(A::AbstractArray{C}) where {F, T<:(FixedPoint{F}), C<:Colorant{T, 1}} in Napari.Layers.Image at /home/tlnagy/.julia/packages/Napari/vaVWX/src/image.jl:48

I guess now the big problem is how to overload getindex or similar so that we can load the right chunk from the disk. I'm not sure what the best approach this would be.

[mkitti]

That's just

Napari.jl/src/image.jl

Line 51 in b2c1469

reinterpret(F,channelview(A))

function primitive_array(A::AbstractArray{C}) where C <: Colorant{T,1} where T <: FixedPoint{F} where F
    # 1. Grab channel view which gets us Array{T,3}. This is a Base.ReinterpretArray
    # 2. Convert to Array{F,3}, still a Base.ReinterpretArray
    reinterpret(F,channelview(A))
end

[tlnagy]

Yeah, it looks like it's just converting it to a type that's compatible with Python since the FixedPointNumber is a Julia specific thing:

julia> NumPyArrays.PYARR_TYPES
Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, Ptr{PyCall.PyObject_struct}, UInt16, UInt32, UInt64, UInt8, PyCall.PyObject, ComplexF32, ComplexF64}

What are your thoughts on the getindex problem?