WarpedView or any other suggestion to warp large arrays (Tiff stacks) efficiently?

[gltchr]

I was suggested to re-post from JuliaDiscourse here, hoping to get a more specific help.

Hello,

In order to avoid the XY-problem I will try to provide a context and a desired result first, and only later some particular questions.

Context: I don't have much experience in Julia and (partially in order to fix that) I am making an application for analysis of the fast camera footage now. The raw files are collections of Tiffs or Tiff-stacks with file size around several GB. Camera captures 3 views of the same object filtered through different spectral bandpass filters. These images are combined side-by-side (via optics) on a single camera sensor. The analysis is based on relation of spectral intensities captured from the object. This implies that 3 ROIs have to be selected from raw footage (displaying the same image captured at different wavelengths), precisely aligned and pixel-wise divided one by another.

Algorithm:

1. Load raw footage in a camera-native format (12bit) -> N4f12 to minimize memory footprint.
   Output: Array{N4f12,3}
2. Crop out and align 3 ROIs -> WarpedView.
   Expected output: 3 of view(::Array{N4f12,3})
3. Apply same transformations to a whole-sensor calibration image (this time with allocation) -> warp. Expected output: 3 Array{N4f12,2}
4. Divide (pixel-wise) (for calibration) cropped views of stacks (2) by the corresponding calibration matrices (3)
   Output: 3 of Array{Float32,3}. This will be an intermediate result that should be saved, hence the allocation.
5. Perform pixel-wise division of stacks yielding 2 ratio (e.g. for stack a,b,c -> a/b, b/c)
   Output: 2 of Array{Float32,3}

I was successful implementing a MWE which did allocations at every stage of the process and worked on a rather small test-stack (30 frames). Now the intention is to substitute views where appropriate in order to make it to work with the full size files, (preferably without freezing my laptop 🙃). Ordering was intended to minimize unnecessary allocations and make all the transformations and most importantly, cropping before doing any multiplications or divisions.

Initial naive implementation of (2) was:

for i = 1:size(normed_stack,3)
      L1[:,:,i] = warp(normed_stack[:,:,i], L1trfm, crop_size);
      L2[:,:,i] = warp(normed_stack[:,:,i], L2trfm, crop_size);
      L3[:,:,i] = warp(normed_stack[:,:,i], L3trfm, crop_size);
end

However as this the first time working with views in Julia I am not quite sure how I can return a reference to it after iterating? In order to iterate over it and/or append to it is should be declared in some way, but I don't quite get how does one declare an empty view? 🫠
This is what I am trying to do:

using ImageTransformtations, Rotations, CoordinateTransformations

rawstack = rand(UInt16, (100,600,10)) # Just 10 frames for now 
trfm1=  recenter(RotMatrix(pi/128), center(rawstack[:,:,1]))  ∘ Translation(-10, 20) 

function warpncrop(rawstack::AbstractArray, trfm::Transformation, cropsize::NTuple{2,Int})
      for (i,slc) in enumerate(eachslice(rawstack, dims=3))
              wv[:,:,i] = WarpedView(slc, trfm, cropsize) # This assignment is wrong, but it is trying to convey an idea
      end
      return wv
end

I am not insisting on that particular implementation and will welcome any comment on the general approach as well.

Thanks in advance and all the best,
Igor

[mkitti]

Could you discuss the nature of the 12-bit camera output? Are the pixels initially packed?

https://github.com/JaneliaSciComp/UInt12Arrays.jl

You may want to take a look at the JuliaArrays Github organization. For example, MappedArrays allows you to apply an arbitrary function on get. Also see LazyArrays.

[mkitti]

I'm missing some details of the last example. What is wv? How was it constructed?

[gltchr]

Hello, thank you for taking your time to reply!

I will definitely take a look on MappedArrays and LazyArrays, thanks for a pointer!

Speaking of the raw data it is a stack of tiff images, recorded with a Photron fast camera. The sensor bit-depth is 12bits, however Tiff files are saved as UInt16. I am unsure if the benefits of exporting to packed UInt12 Array will be noticeable in this case, as (compared to further promotions to Float32) the difference between UInt16 (N4f12) and UInt12 are not that great.

With regards to wv that is exactly the culprit - It is undefined and not declared as I don't know how to do it.
I would like it to be a 3d Array of Views. E.g. appending a WarpedView of a new frame on each iteration. Hope that makes sense.