How are regression lines calculated for 2D histograms? I've been looking into how to fit a line to a 2D histogram and have come across different pieces of information which I can't quite piece together. To show what I mean, I want to achieve something like the image below which is taken from a paper:

In the paper they refer to the line simply as the "linear regression result" but it is evidently not a regular linear regression as we are not dealing with just points in a plane as usual. They do not give further details on this that I could catch. I've seen plots like this in other places too even with confidence intervals for the regression line.
Looking around elsewhere I found this question in stack overflow. They advise OP to use the first principal component for the line which makes sense but this should be possible with regression as that is what they use in the paper. In the comments there is also some talk about "weighted least squares" which is a variant of ordinary least sqaures that is useful when the residuals show heteroscedasticity. However, I do not understand why this regression would be relevant to the problem or if thats what they are doing in the paper.
 A: The paper doesn’t seem to say that the regression line is estimated from the histogram. The usual approach, that they likely took, is to calculate regression on the same raw data that was used for generating histogram and show  them both on same plot.
A: For future reference, if anyone struggles (like I did) to see that you can map the regression line from the raw data directly into your histogram plot, here is a small example in R. Props to Tim and the commenters for the insight:
# simulate the raw data
set.seed(123)
nsim <- 1000
noise <- 0.5
p <- rnorm(nsim)
x <- p + rnorm(nsim, sd = noise)
y <- p + rnorm(nsim, sd = noise)

# define the breaks
binSize <- 1
xbrks <- seq(floor(min(x)), ceiling(max(x)), binSize)
ybrks <- seq(floor(min(y)), ceiling(max(y)), binSize)

# calculate 2D histogram and plot
H <- table(findInterval(x, xbrks), findInterval(y, ybrks))
image(xbrks, ybrks, H, col = rev(topo.colors(max(H))))
abline(lm(x ~ y, cbind.data.frame(x = x, y = y)))


The regression is based on the original unbinned data:
points(x, y, col = 2)


