This question is interesting because it concerns the general problem of regression in the sense of characterizing (or estimating based on data) the conditional distribution of one variable, $X_1,$ based on values of another variable $X_2.$ Their lack of statistical independence implies there is something to be gained from this.
We may take our cues from the theory of Generalized Linear Models (GLMs), which suppose the conditional distribution of $X_1$ is a member of a finitely parameterized family of distributions and the parameters depend first on a linear function of $X_2$ and, usually, by a nonlinear transformation of that result (aka the "link function"). In the present question no such parametric assumption is made. We should therefore work with the conditional distribution of $X_1$ directly.
Suppose, then, that $X_2=x_2$ is a given value and that we can make sense of the conditional distribution of $X_1$ given $X_2=x_2.$ Let this distribution function be $F_{x_2}.$ (For the moment, let's drop the $x_2$ subscript for brevity.) By definition, for any number $x_1,$
$$\Pr(X_1 \le x_1\mid X_2=x_2) = F(x_1).$$
The axioms of probability imply $F$ has a right inverse $F^{-1}$ characterized by
$$F\left(F^{-1}(p)\right) = p$$
for all $0\lt p \lt 1.$
Let $Z$ be a random variable with a uniform distribution on the interval $(0,1).$ This means that for all $0\le p\le 1,$ $\Pr(Z \le p) = p$. Note that for any number $x,$
$$\Pr\left(F^{-1}(Z) \le x\right) = \Pr\left(F(F^{-1}(Z)) \le F(x)\right) = \Pr\left(Z \le F(x)\right) = F(x).$$
This shows that the distribution function of the random variable $F^{-1}(Z)$ is $F.$
Suppose it is possible to find such a uniform $Z$ that is independent of $X_2.$ (One can always do this by enlarging the original probability space if necessary.) In this case the foregoing result holds for every possible value $x_2.$ Consequently, if we define the function $f:\mathbb{R}^2\to \mathbb{R}$ by
$$f(x_2, z) = F_{x_2}^{-1}(z),$$
then $(f(X_2,Z), X_2)$ is a bivariate random variable for which
Its second component, $X_2,$ is identical to the second component of the original variable $(X_1,X_2);$ and
The conditional distribution of its first component equals the conditional distribution of $X_1$ given $X_2.$
Consequently,
The random variables $(X_1,X_2)$ and $(f(X_2,Z), X_2)$ are equal in distribution.
A fortiori, $f(X_2,Z)$ and $X_1$ are equal in distribution.
Examples
As the first example, consider the intriguing illustration presented in the reply by Dave at https://stats.stackexchange.com/a/443150/919/. It displays a $42$-element dataset consisting of a sequence of ordered pairs
$$(x_1,x_2) = (i/10, \pm i/10)$$
as $i$ ranges from $-20$ to $20$ inclusive. (The value $(0,0)$ is included twice, even though that is not apparent from the scatterplot.) We may conceive of this in terms of its empirical distribution, which assigns the probability $1/42$ to each element of the dataset, thereby defining a bivariate random variable $(X_1,X_2).$
In this case we will need to define $F_{x_2}$ for each value $x_2$ that appears as a second coordinate in the dataset (it doesn't matter how we define $F_{x_2}$ elsewhere, because those values have no probability). With this in mind, compute
$$F_{i}(x) = \left\{\matrix{0 & x \lt -|i| \\ 1/2 & -|i| \le x \lt |i| \\ 1 & \text{otherwise}}\right.$$
for any number $i.$ One right inverse of $F_{i}$ maps the interval $[0,1/2]$ to $-|i|$ and the interval $(1/2,1]$ to $|i|.$ You can see this implemented in the code below.
The claim is that
the distribution of $(X_1, X_2)$ is the same as that of $(F_{X_2}^{-1}(Z), X_2).$
Rather than do a direct calculation, let's have the computer simulate the latter using R
, to be compared to the scatterplot of $(X_1, X_2)$ in Dave's reply:
set.seed(17)
n <- 1e5 # Simulation size
f <- function(i, p) ifelse (p <= 1/2, -1, 1) * abs(i) # Inverse conditional distribution
X.2 <- sample((-20):20 / 10, n, replace=TRUE) # Sample of X[2]
Z <- runif(n) # Independent sample of Z
X.1 <- f(X.2, Z) # Construction of X[1]
This scatterplot confirms it works correctly:

Another example suggested in comments is a uniform distribution on a circle:

Generating these data required only two changes to the code: $(1+X_2)/2$ has a Beta$(1/2,1/2)$ distribution and $F^{-1}$ now reflects the equation of the circle:
X.2 <- rbeta(n, 1/2, 1/2) * 2 - 1
F.inv <- function(p, i) ifelse (p <= 1/2, -1, 1) * sqrt(1 - i^2)
Finally, the more difficult situations occur when the scatterplot cannot be conceived of as exhibiting either of the $X_i$ as a function of the other, as in this uniform distribution on a square:

Here, the marginal distribution of $X_2$ is a mixture of a uniform distribution on the interval $[-1,1]$ and equal point masses at $\pm 1,$ which can be simulated via
X.2 <- pmax(-1, pmin(1, runif(n, -2, 2)))
The function $F^{-1}_{i}$ depends on whether its argument is in the interval $(-1,1)$ or equal to $\pm1;$ it can be implemented as
F.inv <- function(p, i) ifelse (abs(i)==1, 2*p - 1, ifelse (p <= 1/2, -1, 1))