We need to be able to calculate the likelihood of a candidate circle (i.e. a
center location and a radius) for a given, known point.
We assume that the observed point is distributed with an isotropic Gaussian
noise around its true location, with a candidate standard deviation. The
desired likelihood of the candidate circle and the candidate standard deviation
can be calculated by summing over all the positions on the circle the
likelihood of observing the given point under a Gaussian noise with the
candidate standard deviation.
The following code calculates this (log-)likelihood by performing an approximate
integration along the circle:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 | #' Calculate the log-likelihood of one circle (center and radius) for one
#' point with isotropic Gaussian noise around it
#'
#' See https://en.wikipedia.org/wiki/Gaussian_function#Two-dimensional_Gaussian_function
#'
#' @param x x coordinate of the observed point
#' @param y y coordinate of the observed point
#' @param sd Spread of the Gaussian noise (sd = sigma_x = sigma_y)
#' @param x0 x coordinate of circle center
#' @param x0 y coordinate of circle center
#' @param radius Radius of the circle
#' @param nSteps Number of steps used for approximate integration along the
#' circle
#'
loglikOnePoint = function(x, y, sd, x0, y0, radius, nSteps = 24) {
gridAngles = seq(0, 2 * pi, length.out = nSteps + 1)[1:nSteps]
gridX = x0 + cos(gridAngles) * radius
gridY = y0 + sin(gridAngles) * radius
lik = 0
A = 1 / (2 * pi * sd^2) # Cf. volume under two-dimensional Gaussian
# function on Wikipedia page
for (i in seq_along(gridX)) {
dsquare = (x-gridX[i])^2 + (y-gridY[i])^2
lik = lik + A * exp(- dsquare / (2*sd^2)) * 2 * pi * radius / nSteps
}
return(log(lik))
}
|