Get representative values ​​from 2d density

Question

How do you get N pairs of values ​​that represent the joint probability (2d density) of much larger pairs of values?

I am sampling MCMC by function parameters, and I want to visualize the back density of this function by building, say, 20 translucent lines that visualize the density of the function. I know that I can just make a large enough sample to approximate the density ( like this ), but it can clutter the graph. Rather, I would think that something like percentiles would work. For instance. a pair for every 2% change in percentile should accurately represent the density using 50 pairs.

In particular, I take a sample from a Bayesian model with a two-parameter geometric series. The density does not necessarily increase monotonously (there may be several peaks). It's a little complicated, but to get started, multi-dimensional normal work will be more minimal:

library(MASS)
pairs = mvrnorm(10000, c(1,3), rbind(c(0.2, 0.1), c(0.1, 0.2)))

# Easy to just sample the rows to get an approximate representation:
pairs[sample(nrow(pairs), size=50, replace=FALSE)

# But how to get more certainty that the samples are really representative?
# This would probably start with the density
dens = kde2d(pairs[,1], pairs[,2], n=100)

Here dens$zis a 100 * 100 matrix containing the density for each of the combinations dens$xand dens$y, i.e. (coded) pairs in pairs. Here's a image(dens)visualization: