How to fit intensity peaks from a image?

Question

I have an image that I can convert to a text file / table of intensities. In this image, several regions show higher intensity, i.e., a peak. A peak may be stored in a 10x10 matrix with rows and columns corresponding to the x1, x2 coordinates in the image. Basically my data looks like the two-dimensional normal distribution in the first plot of this pdf. I want to know the center of the peak and the full width half maximum for each of those matrices.

What should be the approach assuming that I can approximate this matrices to a multidimensional Gaussian? I have z, x1 and x2, but I would like to get the regression parameters.

Answer 1

What you need to do is fit a nonlinear regression model where the intensity is the response as a function of the x1, x2 coordinates. The functional form will be the PDF of a multivariate normal. In R, you can use ?nls to do this. You can write the function yourself, but for convenience you may prefer to use ?dmvnorm from the mvtnorm package. Here is a simplified example:

library(mvtnorm)  # for dmvnorm
set.seed(7124)    # this will make the code exactly reproducible
d = expand.grid(x1=1:10, x2=1:10)
i = dmvnorm(as.matrix(d), mean=c(2.5, 3), sigma=rbind(c( 1.0, -0.6),
                                                      c(-0.6,  1.0) ) )
e = rnorm(100, mean=0, sd=0.1)
y = 5 + 40*i + e;  d$y = y

round(matrix(i, nrow=10, ncol=10, byrow=F), 2)  # output omitted
round(matrix(y, nrow=10, ncol=10, byrow=F), 2)  # output omitted

mvn = function(m1, m2, v1, v2, cv12){  # you need this function for nls()
  y = dmvnorm(as.matrix(d[,1:2]), mean=c(m1, m2), sigma=rbind(c(v1,  cv12),
                                                              c(cv12,  v2) ))
  return(y)
}
mod = nls(y~b0+b1*mvn(m1, m2, v1, v2, cv12), d, algorithm="port",
          start=list(b0=mean(y), b1=diff(range(y)), 
                     m1=d[which.max(d$y),1], m2=d[which.max(d$y),2], 
                     v1=10, v2=10, cv12=0))
summary(mod)
# 
# Parameters:
#      Estimate Std. Error t value Pr(>|t|)    
# b0    4.99398    0.01112  449.02   <2e-16 ***
# b1   39.27467    0.50379   77.96   <2e-16 ***
# m1    2.49315    0.01129  220.87   <2e-16 ***
# m2    3.00054    0.01143  262.58   <2e-16 ***
# v1    0.96444    0.02354   40.97   <2e-16 ***
# v2    0.99235    0.02332   42.55   <2e-16 ***
# cv12 -0.58553    0.01913  -30.60   <2e-16 ***
# ---
# Residual standard error: 0.1012 on 93 degrees of freedom