Maximum of multiple images or arrays in Julia

Question

I want to find a maximum of several images: load them into an array and find a maximum along first dimension.

Python code for example:

import cv2
import sys
import numpy as np

imgs_paths = sys.argv[1:]
imgs = list(map(cv2.imread, imgs_paths))
imgs_arr = np.array(imgs, dtype=np.float32)
imgs_max = np.max(imgs_arr, 0)

What I did is the following:

using Colors, Images

function im_to_array(im)
    img_array = permutedims(channelview(im), (2,3,1)) 
    img_array = Float32.(img_array)
    return img_array
end


imgs = map(Images.load, imgs_paths)
imgs_arr = map(im_to_array, imgs)
a = imgs_arr
b = reshape(cat(a..., dims=1), tuple(length(a), size(a[1])...))
imgs_max = maximum(b, dims=1)

But it's ugly.

I found more simple way to get a maximum (code is below) but it's performance is awful. May be it does not what I'm expecting.

function im_to_array(im)
    img_array = permutedims(channelview(im), (2,3,1)) 
    img_array = Float32.(img_array)
    return img_array
end

imgs = map(Images.load, imgs_paths)
imgs_arr = map(im_to_array, imgs)
imgs_max = max.(imgs_arr...)

Run time of the first method on 120 FHD images is ~5 seconds on my laptop. And I can't figure out run time of the second method because I was waiting for ~30 minutes and it didn't stop. I'm testing it on Julia 1.4.1

Is there a better way to find a maximum of multiple images?

UPD: here is simple case of what I want:

a = [zeros(Int8, 8, 8, 3), zeros(Int8, 8, 8, 3), zeros(Int8, 8, 8, 3)] # 3 black images with shape 8x8
max.(a) #doesn't work
max.(a...) #works with this simple input but when I test it on 120 FHD images it's extremely slow 

UPD2: I tested both methods on a smaller number of images.

function max1(imgs_arr)
    a = imgs_arr
    b = reshape(cat(a..., dims=1), tuple(length(a), size(a[1])...))
    imgs_max = maximum(b, dims=1)
    return imgs_max
end

function max2(imgs_arr)
    return max.(imgs_arr...)
end

imgs_arr = my_imgs_arrays[1:5]

@time max1(imgs_arr)
@time max2(imgs_arr)

  0.247060 seconds (5.29 k allocations: 142.657 MiB)
  0.154158 seconds (44.85 k allocations: 26.388 MiB)

imgs_arr = my_imgs_arrays[1:15]

@time max1(imgs_arr)

@time max2(imgs_arr)

  0.600093 seconds (72.38 k allocations: 382.923 MiB)
  0.769446 seconds (1.24 M allocations: 71.374 MiB)

imgs_arr = my_imgs_arrays[1:25]

@time max1(imgs_arr)

@time max2(imgs_arr)

  1.057548 seconds (23.08 k allocations: 618.309 MiB)
  5.270050 seconds (151.52 M allocations: 2.329 GiB, 4.77% gc time)

So, more images I use - more slowly it works.

Answer 1

It seems like you're looking to do a pairwise max reduction across a number of images. First, here's a function to generate random "images":

rand_images(k, dims...) = [rand(UInt8, dims...) for _ = 1:k]

I'll generate a vector of three random 10x12 images:

julia> images = rand_images(3, 10, 12)
3-element Array{Array{UInt8,2},1}:
 [0x51 0xdc … 0xf7 0x1e; 0xe1 0x10 … 0xd8 0x98; … ; 0x54 0x45 … 0x7a 0xaf; 0x7b 0xfc … 0x0a 0x81]
 [0xc8 0xa5 … 0xa8 0x81; 0x92 0x89 … 0x9f 0xbe; … ; 0x6a 0x03 … 0xb1 0xfd; 0x34 0xa9 … 0xa3 0x50]
 [0x26 0x9b … 0x2a 0x7c; 0x5c 0x7d … 0x8d 0x2b; … ; 0x32 0x1b … 0x57 0xdf; 0x96 0xa1 … 0x2a 0xc9]

One straightforward way to do this is to do a pairwise max reduction:

julia> using BenchmarkTools

julia> @btime reduce(images) do a, b
           max.(a, b)
       end
  400.485 ns (2 allocations: 416 bytes)
10×12 Array{UInt8,2}:
 0xc8  0xdc  0x82  0xa7  0xa6  0xce  0xcd  0xb2  0x6e  0xba  0xf7  0x81
 0xe1  0x89  0x9f  0xeb  0x89  0xdf  0xd2  0xd2  0xab  0xea  0xd8  0xbe
 0xeb  0xdd  0x9e  0xe2  0xf5  0x4b  0xd2  0xe8  0xe4  0xf8  0xb9  0xf8
 0x63  0xa3  0xd7  0xea  0xf0  0x93  0xed  0xf7  0xfb  0xfb  0x9f  0xbb
 0xf2  0x51  0xf0  0xd4  0xfc  0xcf  0xf4  0xdd  0xeb  0xc3  0xe9  0xf9
 0xf8  0x72  0xfa  0x92  0x72  0xaa  0xa2  0xed  0xa1  0xdf  0xf1  0xd0
 0xef  0xe6  0x64  0xb3  0xd0  0x6a  0xce  0x9e  0x96  0xba  0xed  0xf9
 0xdb  0xc5  0x52  0xb3  0xf7  0xd1  0xdd  0xba  0xac  0xbc  0xd3  0xa1
 0x6a  0x45  0x88  0xda  0xf5  0xc6  0xcf  0x64  0xbc  0xf9  0xb1  0xfd
 0x96  0xfc  0xb1  0xc0  0xc4  0xcf  0x89  0xb4  0xe8  0xad  0xa3  0xc9

That's pretty fast: 400ns. I would time it on images of size comparable to what you're doing, but you didn't mention images sizes that I can see (the code isn't data dependent, so the data in the images shouldn't matter).

The reduction computes a maximal slice, reducing that with an image at a time, which may not be the fastest way to do this. It seems like it may be faster to compute each maximal "pixel" one at a time across all the images, which is a bit more complicated but can also be done:

function max_images(images::Vector{<:Array})
    M = copy(images[1])
    for i = 1:length(M)
        for j = 2:length(images)
            M[i] = max(M[i], images[j][i])
        end
    end
    return M
end

This works but it takes 421 nanoseconds which is slower than the array reduce version! Oops. One of the reasons is that there's no guarantee that the images are all the same size so there's bounds checking in the inner loop indexing into each image. We can skip that at our own risk by putting an inbounds annotation on @inbounds M[i] = max(M[i], images[j][i]). That brings the time down to 282 ns. There's a bit more speed that can be gained by telling the compiler that it can safely reorder both of the loops to take advantage of instruction-level parallelism by putting the @simd macro on each for loop. That brings the time down to 240 ns. The final version of the code is:

function max_images(images::Vector{<:Array})
    M = copy(images[1])
    @simd for i = 1:length(M)
        @simd for j = 2:length(images)
            @inbounds M[i] = max(M[i], images[j][i])
        end
    end
    return M
end