2008年4月11日 星期五

2008年4月10日 星期四

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] }

A New Face of this Blog

Finally, I decided to switch this blog space to my experimental "nerd base". The first reason is that I think it can be a good review to track all the works I've done. The other reason is I hope it can help people who may meet the same problems and is looking for solutions. I'll only post articles related to my research, study, and work. Please feel free to leave a message if you have any idea or better solution, or, just make a comment^^. Enjoy~~~~