Plotting image data in Julia is much slower than MATLAB

I often visualize 2D raster grids of elevation data, where the grids can be 20,000x20,000 pixels, and each pixel contains an elevation value. In MATLAB, imagesc(x,y,Z) is lighting fast for this sort of thing, but the same task in Julia slows my computer to a crawl.

Following the suggestions here, I’ve tried heatmap and imshow with Plots.jl and ImageView, but every Julia solution I’ve found is extremely slow. Is there a fast solution for plotting large data grids in Julia?

Last time I looked at the problem of plotting large images (in 2019), all packages I tried failed pretty badly: cannot plot large images · Issue #10 · JuliaAstro/AstroImages.jl · GitHub. I don’t know if situation improved since then. I’m not very hopeful as an issue I opened about this for Makie was closed last week without resolution: Error plotting a *large* matrix · Issue #335 · JuliaPlots/Makie.jl · GitHub.

1 Like

Presumably Matlab is downsampling the image to plot it? This is pretty quick (0.007 seconds) with Images:

using Images, FixedPointNumbers
img = rand(N0f8, 2*10^4,2*10^4);
@time Gray.(imresize(img, 500,500))
1 Like

That’s a reasonable guess, @stevengj, but I don’t think Matlab downsamples. Here’s MATLAB:

>> tic;imagesc(Z);toc
Elapsed time is 0.056218 seconds.
>> tic;imresize(Z,.03);toc
Elapsed time is 0.196041 seconds.

That means that it takes 4x as long to resize the grid as it takes to plot it. Of course, plotting every Nth row and column would be much faster, but I’ve never seen any signs of aliasing, which would be easy to spot.

I will note that my data is Float64 precision, and plotting it at full resolution takes more like 13 seconds in Julia. Also, it appears that plotting with the Images package in Julia doesn’t allow zooming?

Unless you have a 20000x20000 pixel display, something is downsampling…

9 Likes

Good point.

My experience with Matlab is that imagesc() downsamples beyond some unspecified size, and going through contortions to try to get a true 1 pix-per-element plot is actually pretty difficult. (I found I had to use pcolor() to bypass the downsampling logic and run Matlab in headless mode to create larger-than-screen figures — specifically on Linux, since that probably matters given how unpredictable/unrepeatable Matlab’s plotting is…)

3 Likes

@time Gray.(img) is the time taken to map Gray onto the img, no? If you do g = reinterpret(Gray{Float64}, img) instead, that time goes away (and you aren’t making a copy of the data). It still takes time to display the image, however.

If you have Float64 data (normalized in [0,1]), you could also do e.g.

imresize!(Array{Gray{N0f8}}(undef, 500, 500), img)

to resize directly into a lower-precision array (since you’ll only have 8-bit precision when plotting anyway).

@stevengj I appreciate your help on this. Can you point me to the documentation for the Gray function? The syntax not very intuitive for a newbie like me, and it’s unclear how I can use other colormaps or adjust the color axis limits.

Gray is a type defined in Colors.jl. Gray(x), where x is a floationg-point number in [0,1], just returns an object of this type representing the corresponding gray color (0 = black, 1 = white). The syntax Gray.(img) is an instance of Julia’s “dot syntax” for vectorization: it just applies the Gray function elementwise to each entry of the matrix img to produce an array of Gray values. If you have any array of colors, then it will be displayed as an image (thanks to the ImageShow.jl package, which is loaded by the Images.jl package).

More generally, you can use any method you want to generate an array of colors from an image. If you want a variety of different colormaps etcetera, one option is to use the ColorSchemes.jl package. For example, here is an imagesc function implemented with ColorSchemes.jl, which automatically downsamples (with imresize) any data bigger than a specified maxsize:

using ColorSchemes, Images

function imagesc(data::AbstractMatrix{<:Real};
                 colorscheme::ColorScheme=ColorSchemes.viridis,
                 maxsize::Integer=500, rangescale=:extrema)
    s = maximum(size(data))
    if s > maxsize
        return imagesc(imresize(data, ratio=maxsize/s);   # imresize from Images.jl
                       colorscheme, maxsize, rangescale)
    end
    return get(colorscheme, data, rangescale) # get(...) from ColorSchemes.jl
end

This is pretty quick for me. e.g.

x = range(0,1,length=2*10^4)
y = x'
img = atan.(x, y); # a 20000x20000 Float64 array
@time imagesc(img)

reports a time of < 0.02 seconds.

Even more simply and generally, you can just call imresize to downsample the data, and then pass it to Plots.jl or Makie.jl or any other plotting package you want. Once the data is downsampled, plotting should be fast. Probably my imagesc function is quicker because it just shows a raw PNG image of the colormapped data, but on the other hand by using Plots.jl (and similar) you can add annotations and other niceties, and have many more options.

5 Likes

Do you think it’s important to low pass filter the image first? It’s mentioned in the docstring, and I’ve seen some aliasing artifacts before when resizing spectrograms to plot.

@stevengj That’s awesome! Thanks for explaining the Gray function and pointing me to the docs.

Following your MWE, I’m able to replicate your function, getting

@time imagesc(img)
0.065483 seconds (14.19 k allocations: 29.380 MiB, 30.42% compilation time)

That’s essentially the same as the performance I’m getting from Matlab’s imagesc:

[X,Y] = meshgrid(linspace(0,10,20e3),linspace(0,10,20e3)); 
Z = sin(X)+cos(Y); 

tic
imagesc(Z)
toc
Elapsed time is 0.054676 seconds.

However, from what you’re saying, it sounds like in Julia I need to choose between this implementation if I want to achieve Matlab’s initial plotting speed, or Plots.jl (etc) if I want Matlab’s ability to plot in x,y coordinates, zoom, pan, annotate, dynamically change colormaps, etc?

@ericphanson I don’t see any mention of antialiasing in the docs for Julia’s imresize, but I know that it’s the default behavior when downscaling with Matlab’s imresize. Would be good to have some clarity in the Julia docs about whether it’s antialiasing or not.

It says

If you are shrinking the image, you risk aliasing unless you low-pass filter img first.

No, you can first downsample the data with imresize and then pass to Plots.jl (etc), which should hopefully be only a little slower and still give you all of the features.

Depends on what the data looks like. If it’s relatively smooth then there might not be much high-frequency data to filter out.

If you have to do low-pass filtering, that imposes a significant slowdown AFAICT. (But I doubt Matlab’s imagesc function is doing that.)

1 Like

It would be great to have more Julia packages for downsampling huge datasets for visualization, since there are a variety of strategies for different kinds of data.

For inspiration, see e.g. the Python datashader package for visualizing point clouds, this thesis on downsampling huge timeseries, this paper on downsampling volumetric data giving a “new low-cost downsampling filter” (and also nicely reviewing a lot of related algorithms), and the Julia GigaSOM.jl package for visualizing huge cytometry datasets. There is also this discussion thread on software for visualizing huge medical images (hundreds of GiB).

6 Likes

In case anyone is interested — after I came across that paper a few months ago, I wrote this implementation to play with it a bit. Bottom panel shows what the reference paper can do to a noisy timestream like the top panel (versus simple binned averaging in the middle panel).

12 Likes

GitHub - org-arl/InteractiveViz.jl: Interactive visualization tools for Julia is built on top of Makie and performs adaptive resampling when zooming. Perhaps it can be of use to readers of this thread? Maybe additional downsampling algorithms can be added here as well.

The goals of the package are (among others)

  • To render perceptually accurate summaries at large scale, allowing drill down to individual data points.
  • To allow generation of data points on demand through a graphics pipeline, requiring computation only at a level of detail appropriate for display at the viewing resolution. Additional data points can be generated on demand when zooming or panning.
6 Likes

Gray is part of Images

julia> using Images

julia> @time save("test.png", rand(Gray, 20, 20))
  2.869541 seconds (1.63 M allocations: 94.804 MiB, 5.60% gc time, 85.79% compilation time)

julia> @time save("test.png", rand(Gray, 20, 20))
0.000943 seconds (29 allocations: 6.203 KiB)

julia> @time save("test.png", rand(Gray, 2_000, 2_000))
  0.328472 seconds (32 allocations: 38.164 MiB, 18.80% gc time)

julia> @time save("test.png", rand(Gray, 10_000, 10_000))
 10.394720 seconds (33 allocations: 953.753 MiB, 8.13% gc time)

julia> @time save("test.png", rand(Gray, 20_000, 20_000))
 46.454625 seconds (33 allocations: 3.725 GiB, 0.32% gc time)

RGB is also available

julia> @time save("test.png", rand(RGB, 20, 20))
  1.728649 seconds (1.23 M allocations: 69.212 MiB, 89.54% compilation time)

julia> @time save("test.png", rand(RGB, 20, 20))
0.086617 seconds (29 allocations: 14.141 KiB)

julia> @time save("test.png", rand(RGB, 2_000, 2_000))
  0.762473 seconds (32 allocations: 114.458 MiB)

julia> @time save("test.png", rand(RGB, 10_000, 10_000))
 28.764319 seconds (33 allocations: 2.794 GiB, 2.18% gc time)

I couldn’t do 20_000 x 20_000 in RGB as the little 8GB PC I’m typing this on OOMs and kills Julia

julia> @time save("test.png", rand(RGB, 200, 200))
  0.011166 seconds (32 allocations: 1.148 MiB)

test