Best practices for handling multi-channel movies in Images.jl discussion

[tlnagy]

I'm not sure that this is the best place to ask this, but I generally take movies that are broken up along channel info and time stored as separate TIFF files in a directory (using MicroManager by the way). Currently I'm just using a bunch of looping load calls for everything in readdir. Is there a better way of handling this? I would prefer to have the time dimension and channel info in the image metadata if possible.

[timholy]

TIFF is in many ways a terrible format. I like NRRD, which is (decently?) supported by NRRD.jl. It supports mmapping a large file and should load metadata properly, but obviously please do report any issues you discover.

[timholy]

(NRRD lets you concatenate all your image files into a single file.)

[tlnagy]

Is that what your lab uses? I agree TIFF kinda sucks for this use-case. Also, is it possible to tell if a file that is loaded is being mmapped or loaded into memory? There's not too much documentation on it.

[timholy]

We have a custom acquisition program that outputs in our own format, ImagineFormat. But the binary data are in the same format as NRRD (which is to say a raw dump); it's only the header that's different.

I don't know of a way to tell whether an array is mmapped. IIUC, it's possible that most of your arrays are, as malloc uses mmap internally for large requests.

[tlnagy]

NRRD.jl has some issues, I have some fixes here: JuliaIO/NRRD.jl#3. It appears that Fiji's NRRD writer has some weird behavior because it encodes a XYT image as

Gray Images.Image with:
  data: 2048×2044×96 Array{UInt16,3}
  properties:
    spaceunits:  pixel pixel pixel
    colorspace: Gray
    spatialorder:  x y z
    spacedirections:  [1.0,0.0,0.0] [0.0,1.0,0.0] [0.0,0.0,1.0]
    comments:  Created by Nrrd_Writer at Sun Sep 11 13:43:43 PDT 2016
    pixelspacing:  1.0 1.0 1.0

This has a couple problems. First, the spatialorder has z not t and the spaceunits are all wrong (this actually gives SIUnits some issues since it pixel isn't a SIUnit. Is it possible to create a XYT image in Images.jl and save it? with the correct spaceunits?

[timholy]

Indeed let's fix whatever needs doing on the Julia side. Fiji's issues obviously would be ideal to fix too, but worst-case scenario one can edit the header, esp. if Fiji has the option to save as separate .nhdr and .raw files.

Overall I'd say NRRD.jl needs to do a lot more to make sure it can read/write NRRD files that are compatible with other software, ideally even when other software writes buggy headers. The implementation is just taken directly from the description of the format, but sometimes I find the description slightly confusing. And there's nothing like a test to discover you've gotten something wrong.

For creating an image in Julia, this seems to work:

julia> using Images

julia> img = grayim(collect(reshape(0x01:0x30, 4, 6, 2)))
Gray Images.Image with:
  data: 4×6×2 Array{FixedPointNumbers.UFixed{UInt8,8},3}
  properties:
    colorspace: Gray
    spatialorder:  x y z

julia> img["timedim"] = 3
3

julia> img["spatialorder"] = ["x", "y"]
2-element Array{String,1}:
 "x"
 "y"

julia> img
Gray Images.Image with:
  data: 4×6×2 Array{FixedPointNumbers.UFixed{UInt8,8},3}
  properties:
    timedim: 3
    colorspace: Gray
    spatialorder:  x y

julia> save("/tmp/test.nrrd", img)

which dumps out with header

NRRD0001
type: uint8
dimension: 3
sizes: 4 6 2
kinds: space space time
encoding: raw
endian: little
spacings: 1.0 1.0 nan

which seems plausible to me.

[timholy]

What does the actual header written by Fiji look like?

[tlnagy]

This what Fiji created from a virtual stack of TIFF images I had:

NRRD0004
# Created by Nrrd_Writer at Sun Sep 11 13:43:43 PDT 2016
type: uint16
encoding: raw
endian: big
dimension: 3
sizes: 2048 2044 96
space dimension: 3
space directions: (1.0,0,0) (0,1.0,0) (0,0,1.0)
space units: "pixel" "pixel" "pixel"

which looks like a bug in Fiji moreso than in NRRD.jl

[timholy]

Thanks.

Assuming you're on a little endian machine, that endian: big poses an interesting challenge. Formerly, it would have obviated the use mmap on the input array, which would be rather unfortunate if you recorded hours of data and had a 1TB file to deal with. In the past you would have been essentially forced to copy the file to little-endian format to work with it in Julia.

But now we have MappedArrays, and we can do on-the-fly transformation with essentially no "extra" overhead:

julia> using MappedArrays

julia> a = [0x0102, 0x0304]
2-element Array{UInt16,1}:
 0x0102
 0x0304

julia> aswap = mappedarray((ntoh,hton), a)
2-element MappedArrays.MappedArray{UInt16,1,Array{UInt16,1},Base.#ntoh,Base.#hton}:
 0x0201
 0x0403

julia> aswap[1]
0x0201

Performance:

julia> a = rand(UInt16, 10^4);

julia> aswap = mappedarray((ntoh,hton), a);

julia> using BenchmarkTools

julia> @benchmark sum($a)
BenchmarkTools.Trial: 
  samples:          10000
  evals/sample:     60
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  0.00 bytes
  allocs estimate:  0
  minimum time:     866.00 ns (0.00% GC)
  median time:      874.00 ns (0.00% GC)
  mean time:        884.48 ns (0.00% GC)
  maximum time:     1.45 μs (0.00% GC)

julia> @benchmark sum($aswap)
BenchmarkTools.Trial: 
  samples:          10000
  evals/sample:     10
  time tolerance:   5.00%
  memory tolerance: 1.00%
  memory estimate:  0.00 bytes
  allocs estimate:  0
  minimum time:     1.17 μs (0.00% GC)
  median time:      1.22 μs (0.00% GC)
  mean time:        1.23 μs (0.00% GC)
  maximum time:     15.93 μs (0.00% GC)

Julia is so awesome, sometimes I just have to cheer.

Note that MappedArrays requires Julia 0.5, so to do this automatically we'll have to bump the Julia version of this package another notch. Probably best to do this in conjunction with release of the rewritten Images?

[timholy]

Oh, and I'm slow to realize this, but perhaps you're just using Fiji to convert a bunch of TIFFs into one NRRD? You should be able to do that in pure Julia just fine. The sequence would be (1) write the NRRD header to a .nhdr file, (2) open a .raw file, (3) loop over the TIFFs, reading them in one-by-one and appending them to the .raw file, (4) close the .raw file.

[timholy]

This section of code:

https://github.com/JuliaIO/NRRD.jl/blob/5849ec0c998a7e2ad2b79acb9e1767c808406c53/src/NRRD.jl#L373-L383

might help when it comes to writing the header.

[tlnagy]

Oh, and I'm slow to realize this, but perhaps you're just using Fiji to convert a bunch of TIFFs into one NRRD? You should be able to do that in pure Julia just fine. The sequence would be (1) write the NRRD header to a .nhdr file, (2) open a .raw file, (3) loop over the TIFFs, reading them in one-by-one and appending them to the .raw file, (4) close the .raw file.

Yeah, this the case. Thanks for the intro to MappedArrays though. They look super cool. My machine is little-endian so that's a bug in the Fiji exporter. Might try and file a bug report with them. But more importantly, I would like to get proper NRRD support in micromanager itself since it would be able to populate all those fields of the NRRD file with the original data.

This section of code:

https://github.com/JuliaIO/NRRD.jl/blob/5849ec0c998a7e2ad2b79acb9e1767c808406c53/src/NRRD.jl#L373-L383

might help when it comes to writing the header.

I might do this in the mean time to test my optical flow algorithm.

[tlnagy]

I opened up an issue with micro-manager here: micro-manager/micro-manager#446 We'll see what they think.

[tlnagy]

I haven't seen any examples of NRRD storing channel information, do you happen to know what that looks like? I want to make sure it can handle my data needs. Thanks!

[timholy]

I played with this a bit yesterday. At least some files have an appalling lack of concrete information about color or hints how to interpret the dimensions. 😦 Assuming that a dimension of size 3 means color is oh-so-fragile, because I can't represent a 3d grayscale image with 3 pixels along the fastest dimension.

However, starting with format version 3 one can use the kinds field to indicate color & colorspace. See http://teem.sourceforge.net/nrrd/format.html.

[tlnagy]

At least some files have an appalling lack of concrete information about color or hints how to interpret the dimensions.

That file is scary 😱

How are channels represented in Images/ImageCore? I'm interested because I generally have images with distinct channels, i.e. I imaged cells with a 405nm laser and brightfield.

[timholy]

You might be interested in Overlay in Images. Something like it will certainly still be around in the new Images, even if I do decide to implement it slightly differently.

[tlnagy]

Overlay is super cool, it's what I was looking for. I don't really understand how that type relates to other images though. I would think that a channel is simply another dimension that could be handled similarly to time. The only difference between a multichannel and channel-less image would be during display where we can collapse the multichannel image into a single image with multiple colors. This is different than timeseries images where a such a consistent mapping does not exist (i.e. why I brought up https://github.com/timholy/Images.jl/issues/557)

[timholy]

The only useful thing an Overlay currently does is to blend images for display.

Are you saying you'd like to see an interface that provides access like img[c, y, x, t] where c is the "channel" number? Or perhaps img[y, x, t].c, where in that case you could give names to your channels?

[tlnagy]

I was initially thinking the former, but I'm intrigued by the latter. How would the image be stored in memory with the latter?

With regards to the first option, I was thinking a 4D matrix* where you would have img[x,y,t,c] and whenever you sliced say img[:, :, 1, :] it would show this as a RGB image where the c dimension would be collapsed for viewing. You would still be able access the channels separately via img[:, :, 1, 1], for example.

• if we also had a z dimension then we couldn't do this until one z-slice is selected

[timholy]

Closed by #577