2008年4月11日 星期五
2008年4月10日 星期四
LU Decomposition
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
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
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
LetThe 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 toset 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] }