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: https://github.com/JuliaAstro/AstroImages.jl/issues/10#issuecomment-480535886. 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: https://github.com/JuliaPlots/Makie.jl/issues/335.
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))
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…
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…)
@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.
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
imgfirst.
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.)
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).
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).
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.
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)
