[R-sig-Geo] Find a circle center with spatial points

[Edzer Pebesma]

Here's a toy R script that fits circle parameters to noisy data, using
optim:

x0 = 10
y0 = 5
r = 4

n = 50
sd = .3

rad = runif(n) * 2 * pi
x = x0 + r * cos(rad) + rnorm(n, sd = sd)
y = y0 + r * sin(rad) + rnorm(n, sd = sd)
plot(x, y, asp = 1)

SSdist2circ = function(par, x, y)  {
	x0 = par[1]
	y0 = par[2]
	r = par[3]
	sum(abs(sqrt((x-x0)^2+(y-y0)^2) - r))
}


plotCirc = function(x0, y0, r, ...) {
	sq = seq(0, 2 * pi, length.out = 100)
	lines(x0 + r * cos(sq), y0 + r * sin(sq), ...)
}

# model:
plotCirc(x0, y0, r, col = 'blue')
# fitted:
p = optim(c(x0 = 0, y0 = 0, r = 1), SSdist2circ, x = x, y = y)
plotCirc(p$par[1], p$par[2], p$par[3], col = 'red')


On 18/03/16 23:28, Chris Reudenbach wrote:
> Oops :-[ , Barry thanks for clarification , I got it now.
> 
> I don't understand the math either but as a kind of "reparation"
> for not reading well enough I found the authors java code.
> https://www.spaceroots.org/documents/circle/CircleFitter.java
> After download you should for standalone compilation delete line 37
> "package org.spaceroots;"
> compile it with:
> 
>  javac CircleFitter.java -Xlint:unchecked
> 
> then you can run it with:
> 
> java CircleFitter input.file
> 
> The input data should be formated like:
> 
> # input
> # x y
> 0 5
> 1 4.5
> 2.5 4
> 3 3.5
> 4 2
> 5 0
> 
> the above input yields
> 
> initial circle: -001.69467803 -000.69446643 006.23578103
> converged after 7 iterations
> final circle: -001.34339845 -001.34426151 006.44308386
> with the format x,y,radius
> 
> which at least make sense with respect to the data.
> 
> You easily can  run it from R by system() and grab the output in a file
> using a OS depending pipe
> or you may run it using the rJava package which could be more complex.
> 
> cheers Chris
> 
> if you want you can call it from R with javaR
> 
> Am 18.03.2016 um 20:11 schrieb Barry Rowlingson:
>> On Fri, Mar 18, 2016 at 6:49 PM, Chris Reudenbach
>> <reudenbach at uni-marburg.de> wrote:
>>> Because it seems to be an arc and not a circle issue that you can
>>> solve the
>>> problem by
>>> picking arbitrary two points of your assumed "arc" then construct
>>> (calculate)  the perpendicular bisector of
>>> the line between them and do so for another arbitrary two points of the
>>> assumed "arc".
>>>
>>> The intersection of the perpendicular lines is the assumed center of the
>>> arc.
>>>
>>> If you iterate over all points this should be a pretty good
>>> estimation of
>>> the real center.
>>   This is the  "sample 3 points and find the fitted circle" idea, you
>> are likely to get massive "outliers" and if you take the mean
>> coordinate it could fail horribly. See the paper I linked to for an
>> example. They use the median to get an initial "robust" estimate of
>> x,y,R, and then use some specialised optimisation to improve the
>> estimate - you can't just throw it into "optim"!
>>
>>   Not sure I understand the maths in it yet though....
>>
>> Barry
>>
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160319/ac7a330a/attachment.bin>