LU Decomposition

LU decomposition is a method to solve matrix questions. By finding the "pivot", the original matrix can be transfer to an upper triangular matrix. The unknown matrix then can be easily solved by simple algebra calculation. Take the following matrix for example, the c3Z=C, b2Y+b3Z=B, a1X+a2Y+a3Z=A. Put Z=C/c3 into above equations, the X and Y can be easily solved.

Tcl/Tk Code:
## LU Decomposition. Find X if AX=B by LU decomposition. After calculation X=B.
for {set i 0} {$i <= 1} {incr i 1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} {
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)}]
}

Best Fit Circle Using Least Squares Method

For best fitting a circle with a set of data points, the least squares method has been considered as one of the most useful method. The least squares method is based on minimizing the sum of mean square distance from the circle perimeter to data points. By given n points (xi, yi), 1<=i<=n, the objective function can be defined as

where di is the distance from the point (xi, yi) to the perimeter. The equation of circle is

where (a, b) is the circle center and R is the radius. Then the di is

The equation of the least squares method can be rewritten as

where B=-2a, C=-2b, and D=a^2+b^2-R^2. Next, differentiating F with respect to B, C, D yields a system of linear equations

B, C, and D can be solved by solving matrix

Let

The equation of matrix calculation can be rewritten as

The X can be solved by calculating the equation

where A^-1 is the inverse matrix of A. The matrix X can be solved using LU decomposition to get B, C, and D. The x coordinate of the center (a) equals to -B/2, the y coordinate of the center (b) equals to -C/2, and the radius (R) equals to

Reference:

N. Chernov and C. Lesort, Least squares fitting of circles, J. Math. Imaging & Vision 23, 2005, 239-252.

Code:

Tcl/Tk

## Circular regression (least squares method) by given the coordinates of a set of data points.

## Return the center points (XC , YC), and radius R. => CircleFit X Y returns XC YC R

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 1} { 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 ## Gaussian elimination. Find X if AX=B by LU decomposition. After calculation X=B. for {set i 0} {$i <= 1} {incr i 1} { for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} { 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] }