MLE for Poisson conditioned on multivariate Gaussian?

I am writing some Python code to fit 2D Gaussians to fluorescent emitters on a dark background to determine the subpixel-resolution (x, y) position of the fluorescent emitter. The crude, pixel-resolution (x, y) locations of the pixels are stored in a list xy. The height of the Gaussian represents the predicted pixel intensity at that location. Each 2D Gaussian has 5 parameters, and my end goal is to find the optimal value of those 5 parameters for each peak using MLE.

The number of photons that arrive at each pixel can be assumed to be Poisson-distributed (discrete random variable). The parameter of this Poisson random variable is equal to the height of the 2D Gaussian at that pixel, i.e. the Poisson random variable is conditioned on the multivariate Gaussian. I've been thinking for a few days now, but I can't figure it out - what would an MLE algorithm for a Poisson conditioned on a multivariate Gaussian look like? My main point of confusion is that the Poisson random variable has some optimal parameter, but that just gives me the height of the 2D Gaussian at that pixel, and I want the Gaussian parameters that would give me that optimal height.

def twoDGauss(xy, x_guess, y_guess, size_guess, count_guess, background_guess):
    x, y = xy[:, 0], xy[:, 1]
    psf = count_guess/(2 * np.pi * size_guess ** 2)
    psf *= np.exp(-((x - x_guess) ** 2 + (y - y_guess) ** 2)/(2 * size_guess ** 2))
    psf += background_guess
    return psf

That's my code defining the 2D Gaussian so far. I'm following the ideas in Ovesny et al, 2014, specifically in the supplementary section, if anyone is interested.