It took a lot of time to convert old C-code (for calculating eigenvalues and eigenvectors) to PB.
I've put your point data into a file.
Code: Select all
EnableExplicit
#PB_DBL_MAX = 1.7976931348623158e+308 ; max double value
#PB_DBL_EPS = 2.2204460492503131e-016 ; smallest such that 1.0+#PB_DBL_EPS != 1.0
#MAXIT = 50 ; Max. Iterationszahl pro EW
Macro QUAD(wert)
((wert)*(wert))
EndMacro
Procedure.l Num_ComDiv (ar.d,ai.d,br.d,bi.d,*cr.double,*ci.double)
Protected tmp.d
If br = 0.0 And bi = 0.0
ProcedureReturn 1
EndIf
If Abs(br) > Abs(bi)
tmp = bi / br
br = tmp * bi + br
*cr\d = (ar + tmp * ai) / br
*ci\d = (ai - tmp * ar) / br
Else
tmp = br / bi
bi = tmp * br + bi
*cr\d = (tmp * ar + ai) / bi
*ci\d = (tmp * ai - ar) / bi
EndIf
ProcedureReturn 0
EndProcedure
Procedure.d Num_ComAbs (ar.d,ai.d)
If ar = 0.0 And ai = 0.0
ProcedureReturn 0.0
ElseIf ai = 0.0
ProcedureReturn Abs(ar)
ElseIf ar = 0.0
ProcedureReturn Abs(ai)
EndIf
ar = Abs(ar)
ai = Abs(ai)
If ai > ar
Swap ai, ar
EndIf
ProcedureReturn ar * Sqr(1.0 + ai / ar * ai / ar)
EndProcedure
Procedure Num_Balance (n.i,Array mat.d(2),Array scal.d(1),*low.integer,*high.integer,basis.i)
Protected i.i, j.i, iter.i, k.i, m.i
Protected b2.d, r.d, c.d, f.d, g.d, s.d
b2 = (basis * basis)
m = 0
k = n - 1
Repeat
iter = #False
For j = k To 0 Step -1
r = 0.0
For i = 0 To k Step 1
If i <> j
r + Abs(mat(j,i))
EndIf
Next i
If r = 0.0
scal(k) = j
If j <> k
For i = 0 To k Step 1
Swap mat(i,j) , mat(i,k)
Next i
For i = m To n-1 Step 1
Swap mat(j,i), mat(k,i)
Next i
EndIf
k - 1
iter = #True
EndIf
Next j
Until iter = #False
Repeat
iter = #False
For j = m To k Step 1
c = 0.0
For i = m To k Step 1
If i <> j
c + Abs(mat(i,j))
EndIf
Next i
If c = 0.0
scal(m) = j
If j <> m
For i = 0 To k Step 1
Swap mat(i,j) , mat(i,m)
Next i
For i = m To n-1 Step 1
Swap mat(j,i), mat(m,i)
Next i
EndIf
m + 1
iter = #True
EndIf
Next j
Until iter = #False
*low\i = m
*high\i = k
For i = m To k Step 1
scal(i) = 1.0
Next i
Repeat
iter = #False
For i = m To k Step 1
c = 0.0
r = 0.0
For j = m To k Step 1
If j <> i
c + Abs(mat(j,i))
r + Abs(mat(i,j))
EndIf
Next j
g = r / basis
f = 1.0
s = c + r
While c < g
f * basis
c * b2
Wend
g = r * basis
While c >= g
f / basis
c / b2
Wend
If (c + r) / f < 0.95 * s
g = 1.0 / f
scal(i) * f
iter = #True
For j = m To n-1 Step 1
mat(i,j) * g
Next j
For j = 0 To k Step 1
mat(j,i) * f
Next j
EndIf
Next i
Until iter = #False
ProcedureReturn 0
EndProcedure
Procedure Num_BalBack (n.i,low.i,high.i,Array scal.d(1),Array eivec.d(2))
Protected i.i, j.i, k.i
Protected s.d
For i = low To high Step 1
s = scal(i)
For j = 0 To n-1 Step 1
eivec(i,j) * s
Next j
Next i
For i = low - 1 To 0 Step -1
k = Int(scal(i))
If k <> i
For j = 0 To n-1 Step 1
Swap eivec(i,j), eivec(k,j)
Next j
EndIf
Next i
For i = high + 1 To n-1 Step 1
k = Int(scal(i))
If k <> i
For j = 0 To n-1 Step 1
Swap eivec(i,j), eivec(k,j)
Next j
EndIf
Next i
ProcedureReturn 0
EndProcedure
Procedure.i Num_ElmHes (n.i,low.i,high.i,Array mat.d(2),Array perm.i(1))
Protected i.i, j.i, m.i
Protected x.d, y.d
For m = low + 1 To high-1 Step 1
i = m
x = 0.0
For j = m To high Step 1
If Abs(mat(j,m-1)) > Abs (x)
x = mat(j,m-1)
i = j
EndIf
Next j
perm(m) = i
If i <> m
For j = m - 1 To n-1 Step 1
Swap mat(i,j), mat(m,j)
Next j
For j = 0 To high Step 1
Swap mat(j,i), mat(j,m)
Next j
EndIf
If x <> 0.0
For i = m + 1 To high Step 1
y = mat(i,m-1)
If y <> 0.0
mat(i,m-1) = y / x
y = mat(i,m-1)
For j = m To n-1 Step 1
mat(i,j) - y * mat(m,j)
Next j
For j = 0 To high Step 1
mat(j,m) + y * mat(j,i)
Next j
EndIf
Next i
EndIf
Next m
ProcedureReturn 0
EndProcedure
Procedure.i Num_ElmTrans (n.i,low.i,high.i,Array mat.d(2),Array perm.i(1),Array h.d(2))
Protected k.i, i.i, j.i
For i = 0 To n-1 Step 1
For k = 0 To n-1 Step 1
h(i,k) = 0.0
Next k
h(i,i) = 1.0
Next i
For i = high - 1 To low+1 Step -1
j = perm(i)
For k = i + 1 To high Step 1
h(k,i) = mat(k,i-1)
Next k
If i <> j
For k = i To high Step 1
h(i,k) = h(j,k)
h(j,k) = 0.0
Next k
h(j,i) = 1.0
EndIf
Next i
ProcedureReturn 0
EndProcedure
Procedure.i Num_OrtHes (n.i,low.i,high.i,Array mat.d(2),Array d.d(1))
Protected i.i, j.i, m.i
Protected s.d,x.d = 0.0,y.d,eps.d
eps = 128.0 * #PB_DBL_EPS
For m = low + 1 To high-1 Step 1
y = 0.0
For i = high To m Step -1
x = mat(i,m - 1)
d(i) = x
y = y + x * x
Next i
If y <= eps
s = 0.0
Else
If x >= 0.0
s = -Sqr(y)
Else
s = Sqr(y)
EndIf
y - x * s
d(m) = x - s
For j = m To n-1 Step 1
x = 0.0
For i = high To m Step -1
x + d(i) * mat(i,j)
Next i
x / y
For i = m To high Step 1
mat(i,j) - x * d(i)
Next i
Next j
For i = 0 To high Step 1
x = 0.0
For j = high To m Step -1
x + d(j) * mat(i,j)
Next j
x / y
For j = m To high Step 1
mat(i,j) - x * d(j)
Next j
Next i
EndIf
mat(m,m - 1) = s
Next m
ProcedureReturn 0
EndProcedure
Procedure.l Num_OrtTrans (n.i,low.i,high.i,Array mat.d(2),Array d.d(1),Array v.d(2))
Protected i.i, j.i, m.i
Protected x.d,y.d
For i = 0 To n-1 Step 1
For j = 0 To n-1 Step 1
v(i,j) = 0.0
Next j
v(i,i) = 1.0
Next i
For m = high - 1 To low+1 Step -1
y = mat(m,m - 1)
If y <> 0.0
y * d(m)
For i = m + 1 To high Step 1
d(i) = mat(i,m - 1)
Next i
For j = m To high Step 1
x = 0.0
For i = m To high Step 1
x + d(i) * v(i,j)
Next i
x / y
For i = m To high Step 1
v(i,j) + x * d(i)
Next i
Next j
EndIf
Next m
ProcedureReturn 0
EndProcedure
Procedure Num_HqrVec (n.i,low.i,high.i,Array h.d(2),Array wr.d(1),Array wi.d(1),Array eivec.d(2))
Protected i.i, j.i, k.i, l.i, m.i, en.i, na.i
Protected p.d, q.d, r.d = 0.0, s.d = 0.0, t.d, w.d, x.d, y.d, z.d = 0.0, ra.d, sa.d, vr.d, vi.d, norm.d
norm = 0.0
For i = 0 To n-1 Step 1
For j = i To n-1 Step 1
norm + Abs(h(i,j))
Next j
Next i
If norm = 0.0
ProcedureReturn 1
EndIf
For en = n - 1 To 0 Step -1
p = wr(en)
q = wi(en)
na = en - 1
If q = 0.0
m = en
h(en,en) = 1.0
For i = na To 0 Step -1
w = h(i,i) - p
r = h(i,en)
For j = m To na Step 1
r + h(i,j) * h(j,en)
Next j
If wi(i) < 0.0
z = w
s = r
Else
m = i
If wi(i) = 0.0
If w = 0.0
w = #PB_DBL_EPS * norm
EndIf
h(i,en) = -r / w
Else
x = h(i,i+1)
y = h(i+1,i)
q = QUAD(wr(i) - p) + QUAD(wi(i))
h(i,en) = (x * s - z * r) / q
If Abs(x) > Abs(z)
h(i+1,en) = (-r -w * t) / x
Else
h(i+1,en) = (-s -y * t) / z
EndIf
EndIf
EndIf
Next i
ElseIf q < 0.0
m = na
If Abs(h(en,na)) > Abs(h(na,en))
h(na,na) = - (h(en,en) - p) / h(en,na)
h(na,en) = - q / h(en,na)
Else
Num_ComDiv(-h(na,en), 0.0, h(na,na)-p, q, @h(na,na), @h(na,en))
EndIf
h(en,na) = 1.0
h(en,en) = 0.0
For i = na - 1 To 0 Step -1
w = h(i,i) - p
ra = h(i,en)
sa = 0.0
For j = m To na Step 1
ra + h(i,j) * h(j,na)
sa + h(i,j) * h(j,en)
Next j
If wi(i) < 0.0
z = w
r = ra
s = sa
Else
m = i
If wi(i) = 0.0
Num_ComDiv (-ra, -sa, w, q, @h(i,na), @h(i,en))
Else
x = h(i,i+1)
y = h(i+1,i)
vr = QUAD(wr(i) - p) + QUAD (wi(i)) - QUAD (q)
vi = 2.0 * q * (wr(i) - p);
If vr = 0.0 And vi = 0.0
vr = #PB_DBL_EPS * norm * (Abs(w) + Abs(q) + Abs(x) + Abs(y) + Abs(z))
EndIf
Num_ComDiv (x * r - z * ra + q * sa, x * s - z * sa -q * ra, vr, vi, @h(i,na), @h(i,en))
If Abs(x) > Abs(z) + Abs(q)
h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x
h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x
Else
Num_ComDiv (-r - y * h(i,na), -s - y * h(i,en), z, q, @h(i+1,na), @h(i+1,en))
EndIf
EndIf
EndIf
Next i
EndIf
Next en
For i = 0 To n-1 Step 1
If i < low Or i > high
For k = i + 1 To n-1 Step 1
eivec(i,k) = h(i,k)
Next k
EndIf
Next i
For j = n - 1 To low Step -1
If j <= high
m = j
Else
m = high
EndIf
If wi(j) < 0.0
l = j - 1
For i = low To high Step 1
z = 0.0
y = 0.0
For k = low To m Step 1
y + eivec(i,k) * h(k,l)
z + eivec(i,k) * h(k,j)
Next k
eivec(i,l) = y
eivec(i,j) = z
Next i
Else
If wi(j) = 0.0
For i = low To high Step 1
z = 0.0
For k = low To m Step 1
z + eivec(i,k) * h(k,j)
Next k
eivec(i,j) = z
Next i
EndIf
EndIf
Next j
ProcedureReturn 0
EndProcedure
Procedure.i Num_Hqr2 (vec.i,n.i,low.i,high.i,Array h.d(2),Array wr.d(1),Array wi.d(1),Array eivec.d(2),Array cnt.i(1))
Protected i.i, j.i, na.i, en.i, iter.i, k.i, l.i, m.i
Protected p.d = 0.0, q.d = 0.0, r = 0.0, s.d, t.d, w.d, x.d, y.d, z.d
For i = 0 To n-1 Step 1
If i < low Or i > high
wr(i) = h(i,i)
wi(i) = 0.0
cnt(i) = 0
EndIf
Next i
en = high
t = 0.0
While en >= low
iter = 0
na = en - 1
Repeat
For l = en To low+1 Step -1
If Abs(h(l,l-1)) <= #PB_DBL_EPS * (Abs(h(l-1,l-1)) + Abs(h(l,l)))
Break
EndIf
Next l
x = h(en,en)
If l = en
h(en,en) = x + t
wr(en) = h(en,en)
wi(en) = 0.0
cnt(en) = iter
en - 1
Break
EndIf
y = h(na,na)
w = h(en,na) * h(na,en)
If l = na
p = (y - x) * 0.5
q = p * p + w
z = Sqr(Abs(q))
x + t
h(en,en) = X
h(na,na) = y + t
cnt(en) = -iter
cnt(na) = iter
If q >= 0.0
If p < 0.0
z = p - z
Else
z = p + z
EndIf
wr(na) = x + z
s = x - w / z
wr(en) = s
wi(na) = 0.0
wi(en) = 0.0
x = h(en,na)
r = Sqr(x * x + z * z)
If vec > 0
p = x / r
q = z / r
For j = na To n - 1 Step 1
z = h(na,j)
h(na,j) = q * z + p * h(en,j)
h(en,j) = q * h(en,j) - p * z
Next j
For i = 0 To en Step 1
z = h(i,na)
h(i,na) = q * z + p * h(i,en)
h(i,en) = q * h(i,en) - p * z
Next i
For i = low To high Step 1
z = eivec(i,na)
eivec(i,na) = q * z + p * eivec(i,en)
eivec(i,en) = q * eivec(i,en) - p * z
Next i
EndIf
Else
wr(en) = x + p
wr(na) = wr(en)
wi(na) = z
wi(en) = - z
EndIf
en - 2
Break
EndIf
If iter >= #MAXIT
cnt(en) = #MAXIT + 1
ProcedureReturn en
EndIf
If iter <> 0 And (iter%10) = 0
t + x
For i = low To en Step 1
h(i,i) - x
Next i
s = Abs(h(en,na)) + Abs(h(na,en-2))
y = 0.75 * s
x = y
w = -0.4375 * s * s
EndIf
iter + 1
For m = en - 2 To l Step -1
z = h(m,m)
r = x - z
s = y - z
p = ( r * s - w ) / h(m+1,m) + h(m,m+1)
q = h(m + 1,m + 1) - z - r - s
r = h(m + 2,m + 1)
s = Abs(p) + Abs(q) + Abs(r)
p / s
q / s
r / s
If m = l
Break
EndIf
If Abs(h(m,m-1)) * (Abs(q) + Abs(r)) <= #PB_DBL_EPS * Abs(p) * ( Abs(h(m-1,m-1)) + Abs(z) + Abs(h(m+1,m+1)))
Break
EndIf
Next m
For i = m + 2 To en Step 1
h(i,i-2) = 0.0
Next i
For i = m + 3 To en Step 1
h(i,i-3) = 0.0
Next i
For k = m To na Step 1
If k <> m
p = h(k,k-1)
q = h(k+1,k-1)
If k<> na
r = h(k+2,k-1)
Else
r = 0.0
EndIf
x = Abs(p) + Abs(q) + Abs(r)
If x = 0.0
Continue
EndIf
p / x
q / x
r / x
EndIf
s = Sqr(p * p + q * q + r * r)
If p < 0.0
s = -s
EndIf
If k <> m
h(k,k-1) = -s * x
ElseIf l <> m
h(k,k-1) = -h(k,k-1)
EndIf
p + s
x = p / s
y = q / s
z = r / s
q / p
r / p
For j = k To n-1 Step 1
p = h(k,j) + q * h(k+1,j)
If k <> na
p + r * h(k+2,j)
h(k+2,j) - p * z
EndIf
h(k+1,j) - p * y
h(k,j) - p * x
Next j
If k + 3 < en
j = k + 3
Else
j = en
EndIf
For i = 0 To j Step 1
p = x * h(i,k) + y * h(i,k+1)
If k <> na
p + z * h(i,k+2)
h(i,k+2) - p * r
EndIf
h(i,k+1) - p * q
h(i,k) - p
Next i
If vec <> 0
For i = low To high Step 1
p = x * eivec(i,k) + y * eivec(i,k+1)
If k <> na
p + z * eivec(i,k+2)
eivec(i,k+2) - p * r
EndIf
eivec(i,k+1) - p * q
eivec(i,k) - p
Next i
EndIf
Next k
ForEver
Wend
If vec <> 0
If Num_HqrVec (n, low, high, h(), wr(), wi(), eivec())
ProcedureReturn 99
EndIf
EndIf
ProcedureReturn 0
EndProcedure
Procedure.i Num_Norm1 (n.i,Array v.d(2),Array wi.d(1))
Protected i.l, j.l
Protected maxi.d, tr.d, ti.d
If n < 1
ProcedureReturn 1
EndIf
For j = 0 To n-1 Step 1
If wi(j) = 0.0
maxi = v(0,j)
For i = 1 To n-1 Step 1
If Abs(v(i,j)) > Abs(maxi)
maxi = v(i,j)
EndIf
Next i
If maxi <> 0.0
maxi = 1.0 / maxi
For i = 0 To n-1 Step 1
v(i,j) * maxi
Next i
EndIf
Else
tr = v(0,j)
ti = v(0,j+1)
For i = 1 To n-1 Step 1
If Num_ComAbs(v(i,j), v(i,j+1)) > Num_ComAbs(tr, ti)
tr = v(i,j)
ti = v(i,j+1)
EndIf
Next i
If tr <> 0.0 Or ti <> 0.0
For i = 0 To n-1 Step 1
Num_ComDiv(v(i,j), v(i,j+1), tr, ti, @v(i,j), @v(i,j+1))
Next i
EndIf
j + 1
EndIf
Next j
ProcedureReturn 0
EndProcedure
Procedure.i Num_Eigen (vec.i,ortho.i,ev_norm.i,n.i,Array mat.d(2),Array eivec.d(2),Array valre.d(1),Array valim.d(1),Array cnt.i(1))
Protected i.i, low.i, high.i, rc.i=SizeOf(Integer)
Protected Dim scale.d(n)
Protected Dim d.d(n)
If n < 1
ProcedureReturn 1
EndIf
For i = 0 To n-1 Step 1
cnt(i) = 0
Next i
If n = 1
eivec(0,0) = 1.0
valre(0) = mat(0,0)
valim(0) = 0.0
ProcedureReturn 0
EndIf
rc = Num_Balance (n, mat(), scale(),@low, @high, rc)
If rc
ProcedureReturn 100 + rc
EndIf
If ortho
rc = Num_OrtHes(n, low, high, mat(), d())
Else
rc = Num_ElmHes(n, low, high, mat(), cnt())
EndIf
If rc
ProcedureReturn 200 + rc
EndIf
If vec
If ortho
rc = Num_OrtTrans(n, low, high, mat(), d(), eivec())
Else
rc = Num_ElmTrans(n, low, high, mat(), cnt(), eivec())
EndIf
If rc
ProcedureReturn 300 + rc
EndIf
EndIf
rc = Num_Hqr2 (vec, n, low, high, mat(), valre(), valim(), eivec(), cnt())
If rc
ProcedureReturn 400 + rc
EndIf
If vec
rc = Num_BalBack(n, low, high, scale(), eivec())
If rc
ProcedureReturn 500 + rc
EndIf
If ev_norm
rc = Num_Norm1 (n, eivec(), valim())
If rc
ProcedureReturn 600 + rc
EndIf
EndIf
EndIf
ProcedureReturn 0
EndProcedure
;=========================================================================================
Procedure ReadDATFile(zFile.s,Array dPoints.d(1))
Protected iFile.i, iFormat.i, iAnz.i=0
Protected zZeile.s, zX.s, zY.s, zZ.s
iFile = OpenFile(#PB_Any,zFile)
If iFile
iFormat = ReadStringFormat(iFile)
If iFile
While Not Eof(iFile)
zZeile = ReadString(iFile,iFormat)
zX = StringField(zZeile,1,",")
zY = StringField(zZeile,2,",")
zZ = StringField(zZeile,3,",")
If zX <> "" And zY <> "" And zZ <> ""
ReDim dPoints(iAnz*3+2)
dPoints(iAnz*3+0) = ValD(zX)
dPoints(iAnz*3+1) = ValD(zY)
dPoints(iAnz*3+2) = ValD(zZ)
iAnz + 1
EndIf
Wend
CloseFile(iFile)
ProcedureReturn 0
EndIf
EndIf
ProcedureReturn 1
EndProcedure
Procedure DATcalcSize (Array adPoints.d(1),Array adMinMax.d(1),Array adSP.d(1))
Protected i.i, iAnz.i=ArraySize(adPoints())
adSP(0) = adPoints(0) : adSP(1) = adPoints(1) : adSP(2) = adPoints(2)
adMinMax(0) = adPoints(0) : adMinMax(1) = adPoints(1) : adMinMax(2) = adPoints(2)
adMinMax(3) = adPoints(0) : adMinMax(4) = adPoints(1) : adMinMax(5) = adPoints(2)
For i=3 To iAnz Step 1
adSP(i%3) + adPoints(i)
If adMinMax(i%3 ) > adPoints(i) : adMinMax(i%3 ) = adPoints(i) : EndIf
If adMinMax(i%3+3) < adPoints(i) : adMinMax(i%3+3) = adPoints(i) : EndIf
Next i
iAnz + 1
iAnz / 3
adSP(0) / iAnz
adSP(1) / iAnz
adSP(2) / iAnz
EndProcedure
Procedure DATdrawPoints(Array adPoints.d(1))
Protected i.i, iAnz.i=(ArraySize(adPoints())+1)/3
Protected Dim lafCol.f(3)
lafCol(0) = 1.0
lafCol(1) = 0.6
lafCol(2) = 0.2
lafCol(3) = 1.0
glMaterialfv_(#GL_FRONT_AND_BACK,#GL_AMBIENT_AND_DIFFUSE,lafCol())
glPointSize_(4)
glBegin_(#GL_POINTS)
For i=0 To iAnz-1 Step 1
glVertex3dv_(@adPoints(i*3))
Next i
glEnd_()
EndProcedure
Procedure CalcMat(Array adSP.d(1), Array adPoints.d(1), Array adMat.d(1))
Protected i.i, iAnz.i=(ArraySize(adPoints())+1)/3-1
adMat(0) = 0.0 : adMat(1) = 0.0 : adMat(2) = 0.0 : adMat(3) = 0.0 : adMat(4) = 0.0 : adMat(5) = 0.0 : adMat(6) = 0.0 : adMat(7) = 0.0 : adMat(8) = 0.0
For i=0 To iAnz Step 1
adMat(0) + (adPoints(i*3+0)-adSP(0))*(adPoints(i*3+0)-adSP(0))
adMat(1) + (adPoints(i*3+0)-adSP(0))*(adPoints(i*3+1)-adSP(1))
adMat(2) + (adPoints(i*3+0)-adSP(0))*(adPoints(i*3+2)-adSP(2))
adMat(3) + (adPoints(i*3+1)-adSP(1))*(adPoints(i*3+0)-adSP(0))
adMat(4) + (adPoints(i*3+1)-adSP(1))*(adPoints(i*3+1)-adSP(1))
adMat(5) + (adPoints(i*3+1)-adSP(1))*(adPoints(i*3+2)-adSP(2))
adMat(6) + (adPoints(i*3+2)-adSP(2))*(adPoints(i*3+0)-adSP(0))
adMat(7) + (adPoints(i*3+2)-adSP(2))*(adPoints(i*3+1)-adSP(1))
adMat(8) + (adPoints(i*3+2)-adSP(2))*(adPoints(i*3+2)-adSP(2))
Next i
iAnz + 1
adMat(0) / iAnz
adMat(1) / iAnz
adMat(2) / iAnz
adMat(3) / iAnz
adMat(4) / iAnz
adMat(5) / iAnz
adMat(6) / iAnz
adMat(7) / iAnz
adMat(8) / iAnz
EndProcedure
Procedure VekProd (Array adV1.d(1), Array adV2.d(1), Array adV3.d(1))
adV3(0) = adV1(1)*adV2(2)-adV1(2)*adV2(1)
adV3(1) = adV1(2)*adV2(0)-adV1(0)*adV2(2)
adV3(2) = adV1(0)*adV2(1)-adV1(1)*adV2(0)
EndProcedure
Procedure DrawPlane( Array adP.d(1), Array adV.d(1), dR.d)
Protected Dim adV1.d(2)
Protected Dim adV2.d(2)
Protected Dim lafCol.f(3)
Protected iH.i
If Abs(adV(0)) <= Abs(adV(1)) And Abs(adV(0)) <= Abs(adV(2))
iH = 0
ElseIf Abs(adV(1)) <= Abs(adV(0)) And Abs(adV(1)) <= Abs(adV(2))
iH = 1
ElseIf Abs(adV(2)) <= Abs(adV(0)) And Abs(adV(2)) <= Abs(adV(1))
iH = 2
EndIf
adV2(0) = 0.0 : adV2(2) = 0.0 : adV2(2) = 0.0
adV2(iH) = 1.0
VekProd (adV(),adV2(),adV1())
VekProd (adV1(),adV(),adV2())
lafCol(0) = 0.0
lafCol(1) = 0.4
lafCol(2) = 0.4
lafCol(3) = 0.1
dR * 0.5
glMaterialfv_(#GL_FRONT_AND_BACK,#GL_AMBIENT_AND_DIFFUSE,lafCol())
glBegin_(#GL_QUADS)
glVertex3d_(adP(0)-adV1(0)*dR-adV2(0)*dR,adP(1)-adV1(1)*dR-adV2(1)*dR,adP(2)-adV1(2)*dR-adV2(2)*dR)
glVertex3d_(adP(0)+adV1(0)*dR-adV2(0)*dR,adP(1)+adV1(1)*dR-adV2(1)*dR,adP(2)+adV1(2)*dR-adV2(2)*dR)
glVertex3d_(adP(0)+adV1(0)*dR+adV2(0)*dR,adP(1)+adV1(1)*dR+adV2(1)*dR,adP(2)+adV1(2)*dR+adV2(2)*dR)
glVertex3d_(adP(0)-adV1(0)*dR+adV2(0)*dR,adP(1)-adV1(1)*dR+adV2(1)*dR,adP(2)-adV1(2)*dR+adV2(2)*dR)
glEnd_()
EndProcedure
Procedure DistancePointLine ( Array adQ.d(1), Array adP.d(1), Array adV.d(1), *dU.double, *dA.double)
Protected Dim adU.d(2)
Protected Dim adN.d(2)
Protected dL
dL = adV(0)*adV(0) + adV(1)*adV(1) + adV(2)*adV(2)
*dA\d = 1.0e20
*dU\d = 0.0
adU(0) = adQ(0) - adP(0)
adU(1) = adQ(1) - adP(1)
adU(2) = adQ(2) - adP(2)
If dL > 0.0
adN(0) = adU(1)*adV(2) - adU(2)*adV(1)
adN(1) = adU(2)*adV(0) - adU(0)*adV(2)
adN(2) = adU(0)*adV(1) - adU(1)*adV(0)
*dA\d = (adN(0)*adN(0) + adN(1)*adN(1) + adN(2)*adN(2)) / dL
*dU\d = (adU(0)*adV(0) + adU(1)*adV(1) + adU(2)*adV(2)) / dL
Else
*dA\d = adU(0)*adU(0) + adU(1)*adU(1) + adU(2)*adU(2)
EndIf
*dA\d = Sqr(*dA\d)
EndProcedure
Procedure MatInverse (Array adMatIn.d(1), Array adMatout.d(1))
Protected dD.d
dD = adMatIn(0)*(adMatIn(5)*(adMatIn(10)*adMatIn(15)-adMatIn(11)*adMatIn(14))+adMatIn(6)*(adMatIn(11)*adMatIn(13)-adMatIn( 9)*adMatIn(15))+adMatIn(7)*(adMatIn( 9)*adMatIn(14)-adMatIn(10)*adMatIn(13)))
dD -adMatIn(1)*(adMatIn(4)*(adMatIn(10)*adMatIn(15)-adMatIn(11)*adMatIn(14))+adMatIn(6)*(adMatIn(11)*adMatIn(12)-adMatIn( 8)*adMatIn(15))+adMatIn(7)*(adMatIn( 8)*adMatIn(14)-adMatIn(10)*adMatIn(12)))
dD +adMatIn(2)*(adMatIn(4)*(adMatIn( 9)*adMatIn(15)-adMatIn(11)*adMatIn(13))+adMatIn(5)*(adMatIn(11)*adMatIn(12)-adMatIn( 8)*adMatIn(15))+adMatIn(7)*(adMatIn( 8)*adMatIn(13)-adMatIn( 9)*adMatIn(12)))
dD -adMatIn(3)*(adMatIn(4)*(adMatIn( 9)*adMatIn(14)-adMatIn(10)*adMatIn(13))+adMatIn(5)*(adMatIn(10)*adMatIn(12)-adMatIn( 8)*adMatIn(14))+adMatIn(6)*(adMatIn( 8)*adMatIn(13)-adMatIn( 9)*adMatIn(12)))
If Abs(dD)<1.0e-20
adMatOut( 0) = 1.0 : adMatOut( 5) = 1.0 : adMatOut(10) = 1.0 : adMatOut(15) = 1.0
adMatOut( 1) = 0.0 : adMatOut( 2) = 0.0 : adMatOut( 3) = 0.0 : adMatOut( 4) = 0.0
adMatOut( 6) = 0.0 : adMatOut( 7) = 0.0 : adMatOut( 8) = 0.0 : adMatOut( 9) = 0.0
adMatOut(11) = 0.0 : adMatOut(12) = 0.0 : adMatOut(13) = 0.0 : adMatOut(14) = 0.0
Else
dD = 1.0/dD
adMatOut( 0) = (adMatIn(5)*(adMatIn(10)*adMatIn(15)-adMatIn(11)*adMatIn(14))+adMatIn(6)*(adMatIn(11)*adMatIn(13)-adMatIn( 9)*adMatIn(15))+adMatIn(7)*(adMatIn( 9)*adMatIn(14)-adMatIn(10)*adMatIn(13)))*dD
adMatOut( 1) = (adMatIn(1)*(adMatIn(11)*adMatIn(14)-adMatIn(10)*adMatIn(15))+adMatIn(2)*(adMatIn( 9)*adMatIn(15)-adMatIn(11)*adMatIn(13))+adMatIn(3)*(adMatIn(10)*adMatIn(13)-adMatIn( 9)*adMatIn(14)))*dD
adMatOut( 2) = (adMatIn(1)*(adMatIn( 6)*adMatIn(15)-adMatIn( 7)*adMatIn(14))+adMatIn(2)*(adMatIn( 7)*adMatIn(13)-adMatIn( 5)*adMatIn(15))+adMatIn(3)*(adMatIn( 5)*adMatIn(14)-adMatIn( 6)*adMatIn(13)))*dD
adMatOut( 3) = (adMatIn(1)*(adMatIn( 7)*adMatIn(10)-adMatIn( 6)*adMatIn(11))+adMatIn(2)*(adMatIn( 5)*adMatIn(11)-adMatIn( 7)*adMatIn( 9))+adMatIn(3)*(adMatIn( 6)*adMatIn( 9)-adMatIn( 5)*adMatIn(10)))*dD
adMatOut( 4) = (adMatIn(4)*(adMatIn(11)*adMatIn(14)-adMatIn(10)*adMatIn(15))+adMatIn(6)*(adMatIn( 8)*adMatIn(15)-adMatIn(11)*adMatIn(12))+adMatIn(7)*(adMatIn(10)*adMatIn(12)-adMatIn( 8)*adMatIn(14)))*dD
adMatOut( 5) = (adMatIn(0)*(adMatIn(10)*adMatIn(15)-adMatIn(11)*adMatIn(14))+adMatIn(2)*(adMatIn(11)*adMatIn(12)-adMatIn( 8)*adMatIn(15))+adMatIn(3)*(adMatIn( 8)*adMatIn(14)-adMatIn(10)*adMatIn(12)))*dD
adMatOut( 6) = (adMatIn(0)*(adMatIn( 7)*adMatIn(14)-adMatIn( 6)*adMatIn(15))+adMatIn(2)*(adMatIn( 4)*adMatIn(15)-adMatIn( 7)*adMatIn(12))+adMatIn(3)*(adMatIn( 6)*adMatIn(12)-adMatIn( 4)*adMatIn(14)))*dD
adMatOut( 7) = (adMatIn(0)*(adMatIn( 6)*adMatIn(11)-adMatIn( 7)*adMatIn(10))+adMatIn(2)*(adMatIn( 7)*adMatIn( 8)-adMatIn( 4)*adMatIn(11))+adMatIn(3)*(adMatIn( 4)*adMatIn(10)-adMatIn( 6)*adMatIn( 8)))*dD
adMatOut( 8) = (adMatIn(4)*(adMatIn( 9)*adMatIn(15)-adMatIn(11)*adMatIn(13))+adMatIn(5)*(adMatIn(11)*adMatIn(12)-adMatIn( 8)*adMatIn(15))+adMatIn(7)*(adMatIn( 8)*adMatIn(13)-adMatIn( 9)*adMatIn(12)))*dD
adMatOut( 9) = (adMatIn(0)*(adMatIn(11)*adMatIn(13)-adMatIn( 9)*adMatIn(15))+adMatIn(1)*(adMatIn( 8)*adMatIn(15)-adMatIn(11)*adMatIn(12))+adMatIn(3)*(adMatIn( 9)*adMatIn(12)-adMatIn( 8)*adMatIn(13)))*dD
adMatOut(10) = (adMatIn(0)*(adMatIn( 5)*adMatIn(15)-adMatIn( 7)*adMatIn(13))+adMatIn(1)*(adMatIn( 7)*adMatIn(12)-adMatIn( 4)*adMatIn(15))+adMatIn(3)*(adMatIn( 4)*adMatIn(13)-adMatIn( 5)*adMatIn(12)))*dD
adMatOut(11) = (adMatIn(0)*(adMatIn( 7)*adMatIn( 9)-adMatIn( 5)*adMatIn(11))+adMatIn(1)*(adMatIn( 4)*adMatIn(11)-adMatIn( 7)*adMatIn( 8))+adMatIn(3)*(adMatIn( 5)*adMatIn( 8)-adMatIn( 4)*adMatIn( 9)))*dD
adMatOut(12) = (adMatIn(4)*(adMatIn(10)*adMatIn(13)-adMatIn( 9)*adMatIn(14))+adMatIn(5)*(adMatIn( 8)*adMatIn(14)-adMatIn(10)*adMatIn(12))+adMatIn(6)*(adMatIn( 9)*adMatIn(12)-adMatIn( 8)*adMatIn(13)))*dD
adMatOut(13) = (adMatIn(0)*(adMatIn( 9)*adMatIn(14)-adMatIn(10)*adMatIn(13))+adMatIn(1)*(adMatIn(10)*adMatIn(12)-adMatIn( 8)*adMatIn(14))+adMatIn(2)*(adMatIn( 8)*adMatIn(13)-adMatIn( 9)*adMatIn(12)))*dD
adMatOut(14) = (adMatIn(0)*(adMatIn( 6)*adMatIn(13)-adMatIn( 5)*adMatIn(14))+adMatIn(1)*(adMatIn( 4)*adMatIn(14)-adMatIn( 6)*adMatIn(12))+adMatIn(2)*(adMatIn( 5)*adMatIn(12)-adMatIn( 4)*adMatIn(13)))*dD
adMatOut(15) = (adMatIn(0)*(adMatIn( 5)*adMatIn(10)-adMatIn( 6)*adMatIn( 9))+adMatIn(1)*(adMatIn( 6)*adMatIn( 8)-adMatIn( 4)*adMatIn(10))+adMatIn(2)*(adMatIn( 4)*adMatIn( 9)-adMatIn( 5)*adMatIn( 8)))*dD
EndIf
EndProcedure
Procedure VecMultMat (Array adVecIn.d(1), Array adMat.d(1), Array adVecOut.d(1))
Protected i.i, j.i, k.i
adVecOut(0) = 0.0 : adVecOut(1) = 0.0 : adVecOut(2) = 0.0 : adVecOut(3) = 0.0
For i.i=0 To 15 Step 1
j = i % 4
k = i / 4
adVecOut(j) + adVecIn(k) * adMat(i)
Next i
EndProcedure
Procedure Main()
Protected zFile.s="C:\POINTS.DAT"
Protected Dim adPoints.d(1)
Protected Dim adMinMax.d(6)
Protected Dim adSP.d(3)
Protected Dim adMP.d(3)
Protected Dim adP.d(3)
Protected Dim adQ.d(3)
Protected Dim adU.d(3)
Protected Dim adV.d(3)
Protected Dim adW.d(3)
Protected Dim adMat.d(9)
Protected Dim adMat1.d(16)
Protected Dim adMat2.d(16)
Protected Dim adMat3.d(16)
Protected Dim adMatE.d(2,2)
Protected Dim adEigen.d(2,2)
Protected Dim adRe.d(2)
Protected Dim adIm.d(2)
Protected Dim aiIt.i(2)
Protected Dim lafCol.f(3)
Protected Dim lafPos.f(3)
Protected iWin.i,iOGL.i,iEvent.i,iR.i=3,iW.i,i.i,j.i
Protected dR.d,dX.d,dY.d,dF.d,dXmin.d,dYmin.d,dXmax.d,dYmax.d,dZmin.d,dZmax.d,dFz.d,dM.d
iWin = OpenWindow(#PB_Any, 0, 0, 820, 820, "Point Viewer")
zFile = OpenFileRequester("DAT File",GetHomeDirectory(),"DAT-File|*.DAT|Alle|*.*",0)
If ReadDATFile(zFile.s,adPoints()) <> 0
ProcedureReturn 1
EndIf
If ArraySize(adPoints()) < 8
ProcedureReturn 2
EndIf
DATcalcSize (adPoints(),adMinMax(),adSP())
adMP(0) = (adMinMax(0)+adMinMax(3)) * 0.5
adMP(1) = (adMinMax(1)+adMinMax(4)) * 0.5
adMP(2) = (adMinMax(2)+adMinMax(5)) * 0.5
adV(0) = adMinMax(3)-adMP(0)
adV(1) = adMinMax(4)-adMP(1)
adV(2) = adMinMax(5)-adMP(2)
dR = Sqr(adV(0)*adV(0)+adV(1)*adV(1)+adV(2)*adV(2))*1.4
dFz = 2*dR
If Abs(adMinMax(0)) > dFz : dFz = Abs(adMinMax(0)) : EndIf
If Abs(adMinMax(1)) > dFz : dFz = Abs(adMinMax(1)) : EndIf
If Abs(adMinMax(2)) > dFz : dFz = Abs(adMinMax(2)) : EndIf
If Abs(adMinMax(3)) > dFz : dFz = Abs(adMinMax(3)) : EndIf
If Abs(adMinMax(4)) > dFz : dFz = Abs(adMinMax(4)) : EndIf
If Abs(adMinMax(5)) > dFz : dFz = Abs(adMinMax(5)) : EndIf
CalcMat(adSP(),adPoints(),adMat())
adMatE(0,0) = adMat(0) : adMatE(0,1) = adMat(1) : adMatE(0,2) = adMat(2)
adMatE(1,0) = adMat(3) : adMatE(1,1) = adMat(4) : adMatE(1,2) = adMat(5)
adMatE(2,0) = adMat(6) : adMatE(2,1) = adMat(7) : adMatE(2,2) = adMat(8)
If Num_Eigen(1,1,0,3,adMatE(),adEigen(),adRe(),adIm(),aiIt()) ; calculate eigenvectors and eigenvalues
End
EndIf
dM = #PB_DBL_MAX
j = -1
For i=0 To 2
If adIm(i) = 0.0 And adRe(i) < dM
dM = adRe(i)
j = i
EndIf
Next i
If j = -1
End
EndIf
adV(0) = adEigen(0,j) ; the solution is the eigenvector of the smalest eigenvalue
adV(1) = adEigen(1,j)
adV(2) = adEigen(2,j)
dM = Sqr(adV(0)*adV(0)+adV(1)*adV(1)+adV(2)*adV(2))
adV(0) / dM
adV(1) / dM
adV(2) / dM
adU(0) = 0.0 : adU(1) = 0.0 : adU(2) = 1.0 : adU(3) = 0.0 ; initialize vector for ident
If iWin
SetWindowColor(iWin, RGB(220,220,220))
iOGL = OpenGLGadget(#PB_Any, 10, 10, 800 , 800 , #PB_OpenGL_Keyboard)
SetWindowTitle(iWin, "Calc Plane from "+Str((ArraySize(adPoints())+1)/3)+" points" + " Zoom")
SetGadgetAttribute(iOGL,#PB_OpenGL_Cursor,#PB_Cursor_Cross)
glMatrixMode_(#GL_MODELVIEW)
glLoadIdentity_()
glMatrixMode_(#GL_PROJECTION)
glLoadIdentity_()
glViewport_(0, 0, GadgetWidth(iOGL), GadgetHeight(iOGL))
glOrtho_(adMP(0)-dR,adMP(0)+dR,adMP(1)-dR,adMP(1)+dR,adMP(2)-2*dFz,adMP(2)+2*dFz)
glEnable_(#GL_DEPTH_TEST)
glPolygonMode_(#GL_FRONT_AND_BACK, #GL_FILL )
glDepthFunc_(#GL_LEQUAL)
glClearDepth_(1.0)
glEnable_(#GL_DEPTH_TEST)
glLightModeli_(#GL_LIGHT_MODEL_TWO_SIDE,#True)
glEnable_(#GL_LIGHTING) ;Enable Lighting
glEnable_(#GL_LIGHT0)
glEnable_(#GL_BLEND)
glShadeModel_(#GL_SMOOTH)
lafPos(0) = 10*dR
lafPos(1) = 10*dR
lafPos(2) = 10*dR
lafPos(3) = 0
glLightfv_(#GL_LIGHT0,#GL_POSITION,lafPos())
lafPos(0) = -10*dR
lafPos(1) = 10*dR
lafPos(2) = 10*dR
lafPos(3) = 0
glLightfv_(#GL_LIGHT1,#GL_POSITION,lafPos())
lafCol(0) = 1.0
lafCol(1) = 1.0
lafCol(2) = 1.0
lafCol(3) = 1.0
glColorMaterial_(#GL_FRONT_AND_BACK,#GL_AMBIENT_AND_DIFFUSE)
glLightfv_(#GL_LIGHT0,#GL_AMBIENT,@lafCol())
glLightfv_(#GL_LIGHT0,#GL_DIFFUSE,@lafCol())
glEnable_(#GL_NORMALIZE)
SetGadgetAttribute(iOGL, #PB_OpenGL_SetContext, #True)
glClear_ (#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)
glClearColor_ (0.8, 0.8, 0.8, 0.0)
DATdrawPoints(adPoints())
DrawPlane (adSP(),adV(),dR)
SetGadgetAttribute(iOGL, #PB_OpenGL_FlipBuffers, #True)
Repeat
iEvent = WaitWindowEvent()
If EventWindow() = iWin
If EventGadget() = iOGL
Select EventType()
Case #PB_EventType_LeftButtonUp ; ident
dX = GetGadgetAttribute(iOGL,#PB_OpenGL_MouseX)
dY = GetGadgetAttribute(iOGL,#PB_OpenGL_MouseY)
adP(0) = -1.0 + 2*(0 + dX / GadgetWidth (iOGL)) ; calculate mouse position in range[-1,1]
adP(1) = -1.0 + 2*(1 - dY / GadgetHeight(iOGL))
adP(2) = 0.0
adP(3) = 1.0
glMatrixMode_(#GL_PROJECTION)
glGetDoublev_(#GL_PROJECTION_MATRIX,@adMat1(0))
glMatrixMode_(#GL_MODELVIEW)
glGetDoublev_(#GL_MODELVIEW_MATRIX,@adMat2(0))
MatInverse(adMat1(),adMat3()) ;calculate the inverse of projection matrix
VecMultMat(adP(),adMat3(),adQ()) ;multiply mouse position with inverse matrix
MatInverse(adMat2(),adMat3()) ;calculate the inverse of modelview matrix
VecMultMat(adQ(),adMat3(),adP()) ; calculate world coordiantes of mouse position
VecMultMat(adU(),adMat3(),adW()) ; calculate ident direction
; gluUnProject_(dX,GadgetHeight(iOGL)-dY,0.0,adMat2(),adMat1(),alView(),@adQ(0),@adQ(1),@adQ(2)) , also a posibillity for calculation of world coordinates from mouse position
j = (ArraySize(adPoints())+1)/3-1
dF = #PB_DBL_MAX
For i=0 To j ; find the nearest point to the ident direction
adQ(0) = adPoints(i*3+0)
adQ(1) = adPoints(i*3+1)
adQ(2) = adPoints(i*3+2)
DistancePointLine (adQ(),adP(),adW(),@dX,@dY)
If dY < dF
dF = dY
iW = i
EndIf
Next i
MessageRequester ("Point","Nr. "+Str(iW) + "(X="+StrD(adPoints(iW*3+0),3)+",Y="+StrD(adPoints(iW*3+1),3)+",Z="+StrD(adPoints(iW*3+1),3)+")",#PB_MessageRequester_Ok|#PB_MessageRequester_Info)
Case #PB_EventType_MouseWheel ; rotate
dF = GetGadgetAttribute(iOGL,#PB_OpenGL_WheelDelta)
glMatrixMode_(#GL_MODELVIEW)
glTranslated_(adMP(0),adMP(1),adMP(2))
If iR=2 : glRotated_(dF,0,0,1) : EndIf
If iR=1 : glRotated_(dF,0,1,0) : EndIf
If iR=0 : glRotated_(dF,1,0,0) : EndIf
glTranslated_(-adMP(0),-adMP(1),-adMP(2))
If iR=3 ; zoom at mouse position
glMatrixMode_(#GL_MODELVIEW)
glGetDoublev_(#GL_MODELVIEW_MATRIX,@adMat2(0))
glMatrixMode_(#GL_PROJECTION)
glGetDoublev_(#GL_PROJECTION_MATRIX,@adMat1(0))
MatInverse(adMat1(),adMat3())
adP(0) = -1.0
adP(1) = -1.0
adP(2) = -1.0
adP(3) = 1.0
VecMultMat(adP(),adMat3(),adW()) ; calculate the actual min values for glOrtho
dXmin = adW(0)
dYmin = adW(1)
dZmin = -adW(2)
adP(0) = 1.0
adP(1) = 1.0
adP(2) = 1.0
adP(3) = 1.0
VecMultMat(adP(),adMat3(),adW()) ; calculate the actual max values for glOrtho
dXmax = adW(0)
dYmax = adW(1)
dZmax = -adW(2)
dF = 1.0+dF*0.15 ; define the zoom factor (15%)
dX = GetGadgetAttribute(iOGL,#PB_OpenGL_MouseX)
dY = GetGadgetAttribute(iOGL,#PB_OpenGL_MouseY)
adP(0) = -1.0 + 2*(0 + dX / GadgetWidth (iOGL))
adP(1) = -1.0 + 2*(1 - dY / GadgetHeight(iOGL))
adP(2) = 0.0
adP(3) = 1.0
VecMultMat(adP(),adMat3(),adW()) ; calculate the actual mouse position in projected coordinates
dXmin = adW(0)-dF*(adW(0)-dXmin) ; calculate the new min/max values for glOrtho (differenz from mouse position to actual values)
dXmax = adW(0)+dF*(dXmax-adW(0))
dYmin = adW(1)-dF*(adW(1)-dYmin)
dYmax = adW(1)+dF*(dYmax-adW(1))
glLoadIdentity_()
glOrtho_ (dXmin,dXmax,dYmin,dYmax,dZmin,dZmax)
glMatrixMode_(#GL_MODELVIEW)
glLoadMatrixd_(adMat2())
dF = 1.0
EndIf
Case #PB_EventType_MiddleButtonUp ;center
glMatrixMode_(#GL_MODELVIEW)
glGetDoublev_(#GL_MODELVIEW_MATRIX,@adMat2(0)) ; save the modelview matrix
glMatrixMode_(#GL_PROJECTION)
glLoadIdentity_()
glOrtho_(adMP(0)-dR,adMP(0)+dR,adMP(1)-dR,adMP(1)+dR,adMP(2)-2*dFz,adMP(2)+2*dFz) ; restore the glOrtho from global min/max
glMatrixMode_(#GL_MODELVIEW)
glLoadMatrixd_(adMat2()) ; restore the modelview matrix
Case #PB_EventType_RightButtonUp ; change rotation direction
iR + 1
iR % 4
If iR=2 : SetWindowTitle(iWin, "Calc Plane from "+Str((ArraySize(adPoints())+1)/3)+" points" +" Rotate (z)") : EndIf
If iR=1 : SetWindowTitle(iWin, "Calc Plane from "+Str((ArraySize(adPoints())+1)/3)+" points" +" Rotate (y)") : EndIf
If iR=0 : SetWindowTitle(iWin, "Calc Plane from "+Str((ArraySize(adPoints())+1)/3)+" points" +" Rotate (x)") : EndIf
If iR=3 : SetWindowTitle(iWin, "Calc Plane from "+Str((ArraySize(adPoints())+1)/3)+" points" +" Zoom") : EndIf
EndSelect
EndIf
EndIf
SetGadgetAttribute(iOGL, #PB_OpenGL_SetContext, #True)
glClear_ (#GL_COLOR_BUFFER_BIT | #GL_DEPTH_BUFFER_BIT)
glClearColor_ (0.8, 0.8, 0.8, 0.0)
DATdrawPoints(adPoints())
DrawPlane (adSP(),adV(),dR)
SetGadgetAttribute(iOGL, #PB_OpenGL_FlipBuffers, #True)
Until iEvent = #PB_Event_CloseWindow
EndIf
EndProcedure
Main()
It should be more then 3 or 4 points. Points and the resulting plane are shown in a OpenGL-Gadget.
Use right mouse button to change rotation direction and mouse wheel to rotate.