You can try fitting a multi-dimensional polynomial regression. It seems that for your data, a two dimensional regression model should be fine:
$$ (a)x^2 + (b)x + (c)y^2 + (d)y + (e)xy + (f) $$
In python for example, you can fit a proper model using:
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
X = np.array([ [x0,y0], [x1,y1], ....
I would expect most of the works use Generative Adversarial Networks (GANs) for this because they are powerful generative models capable of learning the complex underlying probability distribution.
In this amazing work the authors used a Conditional GAN, in which they can generate an image conditioned on semantic segmentation map. In your case, you might ...