Regarding best fit circle

[Vaibhav Kumar]

Hi,

Can someone guide me regarding how to create best fit circle using given set of 2D points ?
i want to ask the approach for the same ?



Regards
Vaibhav

[Arjen Markus]

That is not a question about Tcl ;), but it has intrigued me in the past - at least for some definition of the best fit.

Do you mean: the smallest circle that will enclose all points?

Regards,

Arjen

[Vaibhav Kumar]

Yes i meant that.
I want to know the approach and then i will implement it in tcl.

[Arjen Markus]

On Tuesday, February 26, 2019 at 10:27:24 AM UTC+1, Vaibhav Kumar wrote:
> Yes i meant that.
> I want to know the approach and then i will implement it in tcl.

A quick search brought up this article: https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf

It should not be too hard to implement one of the algorithms in there in Tcl.

Regards,

Arjen

[Vaibhav Kumar]

So i will start implementing the full least squares method.
Any guidance or tips ?

[Arjen Markus]

On Tuesday, February 26, 2019 at 10:41:43 AM UTC+1, Vaibhav Kumar wrote:
> So i will start implementing the full least squares method.
> Any guidance or tips ?

Well, have a look at the Tcllib library, more specifically the math module in there. The linearalgebra package may have all the stuff you need. The Wiki has more information - https://wiki.tcl-lang.org/page/Tcllib

Regards,

Arjen

[Vaibhav Kumar]

Thanks i will try this

[Vaibhav Kumar]

I had coded the least squares method but it is giving an error. I tried to figured that out as well. My some values are becoming zeros.

The code is as below

proc CircleFit { X Y } {
set n [llength $X]
catch {
unset XX ;
unset YY ;
unset XY ;
unset XXY ;
unset YYY ;
unset XXX ;
unset XYY
}
set sumX 0 ;
set sumY 0 ;
set sumXX 0 ;
set sumYY 0 ;
set sumXY 0 ;
set sumXXY 0 ;
set sumYYY 0 ;
set sumXXX 0 ;
set sumXYY 0 ;
set sumXXplusYY 0 ;
set sumXXYplusYYY 0 ;
set sumXXXplusXYY 0;
for {set i 0} {$i < $n} {incr i} {
set XX [lappend XX [expr {[lindex $X $i] * [lindex $X $i]}]]
set YY [lappend YY [expr {[lindex $Y $i] * [lindex $Y $i]}]]
set XY [lappend XY [expr {[lindex $X $i] * [lindex $Y $i]}]]
set XXY [lappend XXY [expr {[lindex $X $i] * [lindex $X $i] * [lindex $Y $i]}]]
set YYY [lappend YYY [expr {[lindex $Y $i] * [lindex $Y $i] * [lindex $Y $i]}]]
set XXX [lappend XXX [expr {[lindex $X $i] * [lindex $X $i] * [lindex $X $i]}]]
set XYY [lappend XYY [expr {[lindex $X $i] * [lindex $Y $i] * [lindex $Y $i]}]]
set sumX [expr {$sumX + [lindex $X $i]}]
set sumY [expr {$sumY + [lindex $Y $i]}]
set sumXX [expr {$sumXX + [lindex $XX $i]}]
set sumYY [expr {$sumYY + [lindex $YY $i]}]
set sumXY [expr {$sumXY + [lindex $XY $i]}]
set sumXXY [expr {$sumXXY + [lindex $XXY $i]}]
set sumYYY [expr {$sumYYY + [lindex $YYY $i]}]
set sumXXX [expr {$sumXXX + [lindex $XXX $i]}]
set sumXYY [expr {$sumXYY + [lindex $XYY $i]}]
}
set nsumXXplusYY -[expr {$sumXX + $sumYY}]
set nsumXXYplusYYY -[expr {$sumXXY + $sumYYY}]
set nsumXXXplusXYY -[expr {$sumXXX + $sumXYY}]
set A(0,0) $sumX ;
set A(0,1) $sumY ;
set A(0,2) $n;
set A(1,0) $sumXY ;
set A(1,1) $sumYY ;
set A(1,2) $sumY;
set A(2,0) $sumXX ;
set A(2,1) $sumXY ;
set A(2,2) $sumX;
set B(0) $nsumXXplusYY ;
set B(1) $nsumXXYplusYYY ;
set B(2) $nsumXXXplusXYY ;

for {set i 0} {$i <= 1} {incr i 1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} {
puts $A($k,$i)
puts $A($i,$i)
puts $i\ $k
puts "-------------------"
set temp [expr {double($A($k,$i)) / double($A($i,$i))}]
for {set l [expr {$i + 1}]} {$l <= 2} {incr l 1} {
set A($k,$l) [expr {$A($k,$l) - ($A($i,$l) * $temp)}]
}
set B($k) [expr {$B($k) - (($B($i)) * $temp)}]
}
}
set B(2) [expr {$B(2) / $A(2,2)}]
for {set i 1} {$i >= 0} {incr i -1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} {
set B($i) [expr {$B($i) - $A($i,$k) * $B($k)}]
}
set B($i) [expr {$B($i) / $A($i,$i)}]
}
## Calculate the center and radius of the circle
set XC [expr {-(0.5 * $B(0))}]
set YC [expr {-(0.5 * $B(1))}]
set R [expr {sqrt((pow($B(0),2) + pow($B(1),2)) / 4 - $B(2))}]
return [list $XC $YC $R]
}


puts [CircleFit {10 20 30 40} {10 20 30 40}]




The log was as follows

3000
100
0 1
-------------------
3000
100
0 2
-------------------
0.0
0.0
1 2
-------------------
domain error: argument not in valid range
while executing
"expr {double($A($k,$i)) / double($A($i,$i))}"
(procedure "CircleFit" line 64)
invoked from within
"CircleFit {10 20 30 40} {10 20 30 40}"
invoked from within
"puts [CircleFit {10 20 30 40} {10 20 30 40}]"
(file "testTCL.tcl" line 86)

[Arjen Markus]

I have no time to investigate this problem at the moment, but just a bit of advice on the encoding: why append the results of various expressions to variables like XX and YY? You do not need them as such. Just using the result in the subsequent expressions will make it easier to read.

Also, you might want to use a [foreach] loop:

foreach x $X y $Y {
...
}

this is in general (in my experience) easier to read (and get right) than [for] with an index.

Regards,

Arjen

[Rich]

Vaibhav Kumar <vk53...@gmail.com> wrote:
> proc CircleFit { X Y } {
> set n [llength $X]
> catch {
> unset XX ;
> unset YY ;
> unset XY ;
> unset XXY ;
> unset YYY ;
> unset XXX ;
> unset XYY
> }

Since you are using a proc, none of these variables exist until you set
them, so 'unsetting' them here is redundant.

Additionally, there is no need for a catch, read about the -nocomplain
option to unset in the documentation. Also, note that unset can unset
plural names at once, so all of the above in the catch, plus the catch,
can become one single line.

> set sumX 0 ;
> set sumY 0 ;
> set sumXX 0 ;
> set sumYY 0 ;
> set sumXY 0 ;
> set sumXXY 0 ;
> set sumYYY 0 ;
> set sumXXX 0 ;
> set sumXYY 0 ;
> set sumXXplusYY 0 ;
> set sumXXYplusYYY 0 ;
> set sumXXXplusXYY 0;

This is a good use for a foreach look setting variables by name:

foreach var {sumX sumY sumXX sumYY sumXY sumXXY sumYYY sumXXX sumXYY
sumXXplusYY sumXXYplusYYY sumXXXplusXYY} {
set $var 0

}

> set nsumXXplusYY -[expr {$sumXX + $sumYY}]
> set nsumXXYplusYYY -[expr {$sumXXY + $sumYYY}]
> set nsumXXXplusXYY -[expr {$sumXXX + $sumXYY}]

Don't place your negation outside of the expr expansion, that causes
your numbers to shimmer between numeric types and string types on each
loop cycle (i.e., your code wastes cycles needlessly converting
numerics to strings and back to numerics later). Perform the negation
inside the expression, then the variables will remain numeric types.

> -------------------
> domain error: argument not in valid range
> while executing
> "expr {double($A($k,$i)) / double($A($i,$i))}"
> (procedure "CircleFit" line 64)
> invoked from within
> "CircleFit {10 20 30 40} {10 20 30 40}"
> invoked from within
> "puts [CircleFit {10 20 30 40} {10 20 30 40}]"
> (file "testTCL.tcl" line 86)

Two ways to debug this are:

1) place a temporary puts statement (best to use stderr as the file
descriptor) just prior to this command that prints the values contained
in A($k,$i) and A($i,$i), then you can see what values are inside that
is causing trouble.

2) place a temporary catch around the expr, test the result of the
catch with if to detect an error condition, and puts the values of
A($k,$i) and A($i,$i) when an error condition was detected (as well as
halting at that point).

When you see what values were going in that caused the problem, you'll
be better able to trace back to how they were created.

[Vaibhav Kumar]

The code is failing for smaller values such as {10 20 30 40} but not failing for values such as {100 150 340 420}.

But the list returned by the function should be [xCenter yCenter Radius] but the values are not appropriate for the command :-

puts [CircleFit {100 150 340 420} {100 210 330 580}]

The result was :-
12.142117333727356 526.4655618850339 394.9843754642296

Can anyone help me in figuring out what's going wrong ?

[two...@gmail.com]

On Tuesday, February 26, 2019 at 9:23:26 AM UTC-8, Vaibhav Kumar wrote:

> Can anyone help me in figuring out what's going wrong ?

The error => domain error: argument not in valid range

can be seen when doing

() 4 % expr 0.0 / 0.0

domain error: argument not in valid range

As opposed to when you don't have both num and denom 0.0

() 5 % expr 5 / 0.0
Inf
() 6 % expr 5 / 0
divide by zero


This is what I believe your problem is at the point of the error.

[two...@gmail.com]

Also, for debugging, consider the [parray] command to
dump your entire array values.

So, placing a

parray A

at the point of the error should help.

[Arjen Markus]

Your first example consists of four collinear points - the article is explicit about that situation: all five methods will fail. If you think of what it means in a geometrical sense, then it is clear what is happening: the best "circle" will actually be a line through these points.

As for getting the wrong result: have you tried plotting the circle in relation to the points? What is the result you expect?

Regards,

Arjen

[Vaibhav Kumar]

I am using the code from the article
http://echosaint.blogspot.com/2008/04/best-fit-circle-using-least-squares.html?_sm_au_=iVVq7S4Rn6JsD373


It is failing for the points which are collinear.
But code is doinf OK for the non-collinear points.
What should i do ?

[Arjen Markus]

You can check if the divisor is zero and if so return an error code/message. You may also treat this as a special case - return a line instead of a circle.

It quite depends on what you want to use this for. If you want more considered advice, then tell us the purpose.

Regards,

Arjen

[Vaibhav Kumar]

If i get zero divisor means the points the collinear and i can draw circle accordingly. Thanks for the help

What if i had to draw ellipse using the same ?

[Arjen Markus]

On Wednesday, February 27, 2019 at 12:07:12 PM UTC+1, Vaibhav Kumar wrote:
> If i get zero divisor means the points the collinear and i can draw circle accordingly. Thanks for the help
>
> What if i had to draw ellipse using the same ?

Finding a general ellipse (that is oriented arbitrarily) is a completely different problem and I can just imagine that even an ellipse that has its axes parallel to the x- and y-axes will prove to be more complicated.

The first thing will be to define in what sense an ellipse fits the data points best. You can use the sum of the distances again, like for the circle, but you are probably stuck with a numerical method to approximate the optimum.

Ellipses aligned with the x- and y-axes can be drawn on the canvas using the oval items. More general ones will require a specific technique.

Regards,

Arjen

[keithv]

On Tuesday, February 26, 2019 at 1:07:27 AM UTC-8, Vaibhav Kumar wrote:
> Hi,
>
> Can someone guide me regarding how to create best fit circle using given set of 2D points ?

Check out my page "Smallest Enclosing Disc" (https://wiki.tcl-lang.org/page/Smallest+Enclosing+Disc), which I think does exactly what you're asking for.
Keith

[Vaibhav Kumar]

Please read the whole discussion carefully.
It consisits of link
https://pdfs.semanticscholar.org/faac/44067f04abf10af7dd583fca0c35c5937f95.pdf

related to algorithms for generating best fit circle.

[Kevin Kenny]

For the 'best fit to circle' problem, I've found

http://www.spaceroots.org/documents/circle/circle-fitting.pdf

to be useful. The Levenberg-Marquardt approach is reasonably robust, and where I've used it, I have physical constraints that make degenerate problems unlikely.