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

Barry Rowlingson b.rowlingson at lancaster.ac.uk
Sat Mar 19 15:45:07 CET 2016


And here's a package that implements that, along with finding an
initial point via the 3-point sampling circle method and a function
for making test data:

https://github.com/barryrowlingson/circlefitter/tree/master

no documentation, of course, except the README... And no license...
Edzer wrote a chunk of it... I don't think he'll mind :)

the least-squares method works pretty well for fairly point data that
defines the circle pretty well, but if you do points over a small arc
then I've had cases where the centre is on the wrong side of the
circle. I'm not sure if the optimiser has found the global minimum is
over there or whether it has just missed a minimum at the true x,y,r
values. But data with a large radius variance on a short arc section
is hard to tell which way it is arcing...

At some point I still might implement the methods in that paper...

Barry



On Sat, Mar 19, 2016 at 12:54 PM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
> 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
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



More information about the R-sig-Geo mailing list