Assuming you can easily compute the infimum and supremum of the density function $f$ on a given square with axis-aligned sides, you can use a quadtree to perform an approximate sampling. Let's assume that $f$ is bounded.
Work in 2D Cartesian coordinates and write $S_{l, b, \delta}$ for the axis-aligned square containing all points $(x,y)$ with $l \leq x < l+\delta$ and $b \leq y < b+\delta$, where $\delta > 0$. Assume that you start with a square domain $S_{l_0, b_0, \delta_0}$, within which you want to sample points. Fix an $\epsilon > 0$ and proceed to build a quadtree as follows, starting from the quadtree with a single node $S_{l_0, b_0, \delta_0}$ (a leaf).
For every node $S_{l, b, \delta}$ of the quadtree, compute the infimum $m_{l, b, \delta} := \inf\{f(x,y) | (x,y) \in S_{l, b, \delta} \}$ and supremum $M_{l, b, \delta} := \sup\{f(x,y) | (x,y) \in S_{l, b, \delta} \}$ of the density function $f$ on the node. These two quantities give bounds for the integral of $f$ over the node:
$$
m_{l, b, \delta}\delta^2
\leq
\int_{S_{l, b, \delta}} f(x,y) \,d(x,y)
\leq
M_{l, b, \delta} \delta^2
$$
If $(M_{l, b, \delta} - m_{l, b, \delta})\delta^2 > \epsilon$ for any leaf $S_{l, b, \delta}$ of the quadtree, proceed to subdivide the leaf into four new leaves, each with side $\delta/2$:
$$
S_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}
\hspace{5mm} \text{for} \hspace{5mm}
\lambda, \beta \in \{0, 1\}
$$
Because $f$ is bounded, this process is guaranteed to terminate and the maximum depth for the quadtree is $O\left(\log(\epsilon^{-1})\right)$, or more specifically it is bounded above by the smallest $d \geq 0$ such that:
$$
(M_{l_0, b_0, \delta_0} - m_{l_0, b_0, \delta_0}) \delta_0^2 2^{-2d} \leq \epsilon
$$
The number of nodes in the quadtree can be $O\left(\epsilon^{-1}\right)$ in specific ill-behaved cases, but it can be guaranteed to be $O\left(\log(\epsilon^{-1})\right)$ if some additional regularity requirements are imposed on $f$ (I seem to recall that $f$ being Lipschitz continuous with a fixed upper bound on the Lipschitz constant is enough).
When subdivision is terminated, define an approximate integral $I_{l, b, \delta}$ for each node as follows:
- if $S_{l, b, \delta}$ is a leaf, set the approximate integral to $I_{l, b, \delta} := M_{l, b, \delta} \delta^2$;
- if $S_{l, b, \delta}$ is an internal node, set $I_{l, b, \delta}$ to be the sum of the approximate integrals over its children:
$$
I_{l, b, \delta}
:=
\sum_{\lambda=0}^1 \sum_{\beta=0}^1
I_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}
$$
Because of the termination condition for quadtree subdivision, we now have the following approximation for every node $S_{l, b, \delta}$ in the quadtree:
$$
\left| I_{l, b, \delta} - \int_{S_{l, b, \delta}} f(x,y) \,d(x,y)\right| \leq \epsilon \delta^2
$$
Once the quadtree is built, $\epsilon$-approximate sampling of a point according to the density function $f$ can be performed efficiently by descending the quadtree. Starting from the root, at every internal node $S_{l, b, \delta}$, randomly choose one of the four children $S_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}$ with probability given by:
$$
\mathbb{P}(\lambda, \beta) := \frac{I_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}}{I_{l, b, \delta}}
$$
Once a leaf is reached, sample a uniformly distributed point within it.
This produces independent and identically distributed samples for an $\epsilon$-approximation to $f$, each sample taking $O(\log(\epsilon^{-1}))$ time and using $O(\log(\epsilon^{-1}))$ samples from a uniform probability distribution on $[0,1]$ (one sample for each step of the descent from root to leaf plus two samples for the leaf).
