Advice on plotting a large number of filled polygons

[mforets]

Do you have any suggestion or example about using the library to plot filled polygons? I'm interested in large numbers (~1e6) -- in which case my current workflow with Plots+Jupyter notebook is to save to a separate file though that's a bit cumbersome and the plot quality is not really great; otherwise, like trying to visualize inside the notebook, it just crashes.

Previously for another project i've succesfully used this script, but i'd like to know if it's possible to use the julia object (a Vector{VPolygon{..}}) and send it to Gnuplot.jl directly.

cc @dfcaporale

Thanks in advance.

[lazarusA]

You can send polygons directly to gnuplot. In Gallery, https://lazarusa.github.io/gnuplot-examples/gallery/ there is one example, the 5th, where you just send the points to form the polygons.

[gcalderone]

There are several possible approaches here:

• draw the polygon as gnuplot objects (see examples by @lazarusA above);
• convert all polygon coordinates to a suitable Vector{String} and pass that straight to @gp. This is the simplest method, but it is going to be slow if the number of polygons is too high;
• write a binary file containing all coordinates, and read it in gnuplot.

The following code shows examples for the last two approaches:

using Gnuplot

# A structure to store polygon coordinates
struct Polygon
    x::Vector{Float64}
    y::Vector{Float64}
end

# Generate random polygons
function generate_polygons(n)
    data = Vector{Polygon}()
    for i in 1:Int(n)
        n = mod(i, 5) + 3 # number of vertex
        ρ =  rand(n)
        θ = (rand(n) .+  (0:(n-1))) .* 2π / n
        x = ρ .* cos.(θ) .+ (rand() - 0.5) * 10
        y = ρ .* sin.(θ) .+ (rand() - 0.5) * 10
        push!(data, Polygon(x, y))
    end
    return data
end

#=
Approach 1:

Create a Vector{String}() with all polygon coordinates.
Insert an empty string after each polygon.
=#
function approach1(data)
    tmp = Vector{String}()
    for p in data
        push!(tmp, join(string.(p.x) .* " ".* string.(p.y) .* " " .*
                      string(length(p.x)), "\n") * "\n")
    end
    @gp cblabel="N vertex" "set style fill noborder" tmp "w filledcurves notit lc pal"
end


#=
Approach 2:

Write a binary file with all coordinates.
=#
function approach2(data)
    (path, io) = mktemp()

    n = sum(length.(getfield.(data, :x))) + length(data)
    gpsource = " '$path'binary record=$n format='%float%float%int' u 1:2:3 "
    for polygon in data
        for i in 1:length(polygon.x)
            write(io, convert(Float32, polygon.x[i]))
            write(io, convert(Float32, polygon.y[i]))
            write(io, convert(Int32  , length(polygon.x)))
        end
        write(io, convert(Float32, NaN))
        write(io, convert(Float32, NaN))
        write(io, convert(Int32, 0))
    end
    close(io)
    @gp cblabel="N vertex" "set style fill noborder" "plot $gpsource w filledcurves notit lc pal"
end

And here's the benchmarks:

# Start with just 10^4 polygons
julia> data = generate_polygons(1e4);

julia> @time approach1(data)
  3.629743 seconds (518.90 k allocations: 67.459 MiB, 0.31% gc time)

julia> @time approach2(data)
  0.133400 seconds (181.19 k allocations: 2.881 MiB)

# Try with 10^6 polygons
julia> data = generate_polygons(1e6);

julia> @time approach2(data)
  2.883426 seconds (18.00 M allocations: 282.346 MiB)

The window refresh (e.g. with the qt terminal) may requires a few more seconds.

[gcalderone]

and here's the results for 10^6 polygons:

[mforets]

This is really great, thanks a lot! I look forward to adding Gnuplot as an optional plotting backend to our library LazySets (JuliaReach/LazySets.jl#2108).