ympyrän sovitus 2D

tsuppa kirjoitti:

Löysin sattumoisin koodin joka sovittaa annetuista pisteistä ympyrän. Tällaista kaipaisinkin, vaan eipä toimi, ei. Osaisiko joku neuvoa mikä mättää? (Epäilen jotain yksinkertaista virhettä jota tollo ei tajua).
Tökkää jo riville:"Private Type XYPair"

Alla koodi:

Option Explicit

Const PRECSN = 0.00000001 ' precision for solving circles
Const INCR_STEP = 0.1 ' multiplier for decrementing search increment for solving circles

Private Type XYPair
X1 As Double
Y1 As Double
X2 As Double
Y2 As Double
D As Double
End Type



Private Function Radius_Err(dblXY() As Double, dblX As Double, dblY As Double) As Double()
' given an array of XY points and X and Y center coordinates, calculate the average
' radius and average error (as a positive number)
Dim i As Integer
Dim j As Integer
Dim dblDist() As Double ' array to hold individual distances
Dim dblRad As Double ' intermediate to calculate radius
Dim dblSquares As Double ' intermediate to calculate error
Dim dblRadErr(1) As Double ' array to return radius and error

j = UBound(dblXY, 2) ' how many XY pairs were passed
ReDim dblDist(j)
' loop through the pairs and calculate the individual distances and
' the sum of those distances
For i = 0 To j
dblDist(i) = Sqr((dblXY(0, i) - dblX) ^ 2 (dblXY(1, i) - dblY) ^ 2)
dblRad = dblRad dblDist(i)
Next i
dblRadErr(0) = dblRad / CDbl(j) ' the average distance = radius
' loop through and sum the individual errors as squares
For i = 0 To j
dblSquares = dblSquares (dblDist(i) - dblRadErr(0)) ^ 2
Next i
dblRadErr(1) = Sqr(dblSquares) / CDbl(j) ' the average error
Radius_Err = dblRadErr ' pass the array
End Function

Private Function Solve_Circle(dblXY() As Double, dblPrec As Double) As Double()
'given an array of XY points and a precision for iteration, calculate the XY center
' point and the radius
Dim dblX3 As Double
Dim dblY3 As Double
Dim dblMA As Double
Dim dblMB As Double
Dim dblR As Double
Dim dblNewD As Double
Dim dblLastD As Double
Dim dblNewX As Double
Dim dblNewY As Double
Dim dblLastX As Double
Dim dblLastY As Double
Dim dblXYR(2) As Double
Dim dblIncr As Double
Dim dblErrRad() As Double
Dim dblLastErr As Double
Dim intConvFlg As Integer
Dim Pairs() As XYPair
Dim ThePair As XYPair
Dim i As Long
Dim j As Long
Dim k As Long

' find the two points that have the geatest distance between them
k = -1
ReDim Pairs(0)
' go through all the permutations and calculate the distance
For i = 0 To UBound(dblXY, 2)
For j = 0 To UBound(dblXY, 2)
If i j Then
k = k 1
ReDim Preserve Pairs(k)
With Pairs(k)
.X1 = dblXY(0, i)
.Y1 = dblXY(1, i)
.X2 = dblXY(0, j)
.Y2 = dblXY(1, j)
.D = Sqr((.X1 - .X2) ^ 2 (.Y1 - .Y2) ^ 2)
End With
End If
Next j
Next i
' find the the pair largest distance
For i = 0 To UBound(Pairs)
If Pairs(i).D > ThePair.D Then
With Pairs(i)
ThePair.X1 = .X1
ThePair.Y1 = .Y1
ThePair.X2 = .X2
ThePair.Y2 = .Y2
ThePair.D = .D
End With
End If
Next i
' find the point that is farthest from a line between the previous two points
With ThePair
For i = 0 To UBound(dblXY, 2)
dblNewD = DistToSegment(dblXY(0, i), dblXY(1, i), .X1, .Y1, .X2, .Y2)
If dblNewD > dblLastD Then
dblX3 = dblXY(0, i)
dblY3 = dblXY(1, i)
dblLastD = dblNewD
End If
Next i
' calculate the center and radius given the three points described above
If .X2 - .X1 0# Then
dblMA = (.Y2 - .Y1) / (.X2 - .X1)
Else
dblMA = (.Y2 - .Y1) / 1E-16
End If
If dblMA = 0 Then dblMA = 1E-16
If dblX3 - .X2 0# Then
dblMB = (dblY3 - .Y2) / (dblX3 - .X2)
Else
dblMB = (dblY3 - .Y2) / 1E-16
End If
If dblMB = 0 Then dblMB = 1E-16
dblLastX = (dblMA * dblMB * (.Y1 - dblY3) dblMB * (.X1 .X2) - dblMA * (.X2 dblX3)) / 2# / (dblMB - dblMA)
dblLastY = (dblLastX - (.X1 .X2) / 2#) / dblMA (.Y1 .Y2) / 2
dblR = Sqr((dblLastX - .X1) ^ 2 (dblLastY - .Y1) ^ 2)
End With
dblIncr = dblR / 2#
' get starting value for error
dblErrRad = Radius_Err(dblXY, dblLastX, dblLastY)
dblLastErr = dblErrRad(1)
' loop so long as the increment is larger than the required precision
' increment X and test, if improved, continue incrementing and testing
' if made worse, increment and test the other direction, if improved, continue
' when changing X no longer helps, do the same thing with Y
' when Y no longer helps, check to see if an improvement was made, if improved,
' loop but don't change the casting increment
' once the loop can be completed with no improvement, reduce the increment (usually
' by a factor of 10)
While dblIncr > dblPrec
intConvFlg = 0 ' set flag to 0 = no improvement
dblNewX = dblLastX dblIncr ' new X to be tested
dblErrRad = Radius_Err(dblXY, dblNewX, dblLastY) ' get the new error
If dblErrRad(1) < dblLastErr Then ' if improved
' the improvement must be greater than the precision
While dblErrRad(1) < dblLastErr - dblPrec ' while continuing to improve
intConvFlg = 1 ' set flag for improved, if we get here, we've improved at least once
dblLastX = dblNewX ' last X becomes the X just tried
dblNewX = dblLastX dblIncr ' calculate a new X
dblLastErr = dblErrRad(1) ' last error becomes the error just returned
dblErrRad = Radius_Err(dblXY, dblNewX, dblLastY) ' get a new error value
Wend
ElseIf dblErrRad(1) > dblLastErr Then ' if made worse
dblNewX = dblLastX - dblIncr
dblErrRad = Radius_Err(dblXY, dblNewX, dblLastY)
While dblErrRad(1) < dblLastErr - dblPrec
intConvFlg = 1
dblLastX = dblNewX
dblNewX = dblLastX - dblIncr
dblLastErr = dblErrRad(1)
dblErrRad = Radius_Err(dblXY, dblNewX, dblLastY)
Wend
End If
dblNewY = dblLastY dblIncr
dblErrRad = Radius_Err(dblXY, dblLastX, dblNewY)
If dblErrRad(1) < dblLastErr Then
While dblErrRad(1) < dblLastErr - dblPrec
intConvFlg = 1
dblLastY = dblNewY
dblNewY = dblLastY dblIncr
dblLastErr = dblErrRad(1)
dblErrRad = Radius_Err(dblXY, dblLastX, dblNewY)
Wend
ElseIf dblErrRad(1) > dblLastErr Then
dblNewY = dblLastY - dblIncr
dblErrRad = Radius_Err(dblXY, dblLastX, dblNewY)
While dblErrRad(1) < dblLastErr - dblPrec
intConvFlg = 1
dblLastY = dblNewY
dblNewY = dblLastY - dblIncr
dblLastErr = dblErrRad(1)
dblErrRad = Radius_Err(dblXY, dblLastX, dblNewY)
Wend
End If
' if improved, don't change the casting value
If intConvFlg = 0 Then dblIncr = dblIncr * INCR_STEP
Wend
' load the array
dblXYR(0) = dblLastX
dblXYR(1) = dblLastY
dblXYR(2) = dblErrRad(0)
' pass the array
Solve_Circle = dblXYR
End Function

Public Function DistToSegment(ByVal px As Double, ByVal py _
As Double, ByVal X1 As Double, ByVal Y1 As Double, _
ByVal X2 As Double, ByVal Y2 As Double) As Double
' distance between a line segment and a point
Dim dx As Double
Dim dy As Double
Dim T As Double

dx = X2 - X1
dy = Y2 - Y1
If dx = 0 And dy = 0 Then
' It's a point not a line segment.
dx = px - X1
dy = py - Y1
DistToSegment = Sqr(dx * dx dy * dy)
Exit Function
End If

If dx dy 0 Then
T = (px py - X1 - Y1) / (dx dy)
Else
T = (px py - X1 - Y1) / 1E-16
End If

If T < 0 Then
dx = px - X1
dy = py - Y1
ElseIf T > 1 Then
dx = px - X2
dy = py - Y2
Else
X2 = X1 T * dx
Y2 = Y1 T * dy
dx = px - X2
dy = py - Y2
End If
DistToSegment = Sqr(dx * dx dy * dy)
End Function

kirjoitti:

Anonyymi

Kirjoita kommentti

5000

Lähetä