Procedure to define a best-fit Plane through a set of points

Everything related to 3D programming
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Procedure to define a best-fit Plane through a set of points

Post by IdeasVacuum »

I am getting sets of point coordinates (between 100 and 500 points) that define a plane - they are not all exactly on the plane, there are some outliers. The plane they define is typically at a slight compound angle to one of the Origin planes.

Does anyone have a PB Procedure to define a best-fit plane (orthogonal distance regression plane)?

I have been searching for solutions in plain C, but most examples are in C++, Python, Rust etc which use there own function libraries. The best example I could find: http://mathforum.org/library/drmath/view/63765.html
However, I do not see how to turn that into something that works :?

The final description of the plane I'm looking for is the Plane Normal and point V0

A small sample of points (40) that are representative:

Code: Select all

-3.0,196.0,4.0
3.0,197.0,5.0
4.0,197.0,2.0
4.0,197.0,-2.0
1.0,197.0,-3.0
-4.0,197.0,-3.0
-4.0,197.0,0.0
0.0,197.0,0.0
1.0,197.0,3.0
-6.0,196.0,4.0
-1.14,197.22,1.38
-1.28,197.53,-2.66
-3.26,197.68,-5.63
1.99,197.89,-5.41
3.53,198.19,-8.33
-0.15,198.15,-9.88
-0.35,197.89,-6.74
-3.97,198.02,-10.3
-5.44,197.66,-6.62
-6.09,197.21,-1.28
0.04,196.86,6.63
-3.36,196.6,8.05
-6.91,196.52,7.12
2.26,197.48,0.0
3.9,197.96,-5.2
1.31,198.1,-8.44
-2.16,197.94,-8.35
-5.54,197.49,-4.44
-0.46,196.99,4.66
-0.56,197.73,-4.7
-6.43,197.0,1.3
-4.45,196.98,2.53
-5.13,196.66,6.28
-1.82,196.8,6.32
4.24,197.57,-0.05
-7.56,196.7,4.43
-1.38,196.7,7.89
2.26,196.84,8.14
4.31,197.05,6.58
5.57,197.35,3.49
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
applePi
Addict
Addict
Posts: 1404
Joined: Sun Jun 25, 2006 7:28 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by applePi »

Hi IdeasVacuum
if someone ever is able to design a PB procedure to find the best fit plane for your 3d data, i suggest to use this PB Ogre approach to display the data and the plane, since it is much easier (too much) than OpenGL and more friendly
i will use spheres to draw the data points since the data count is small, and will use for the best fit plane this equation:
1842.08+0.374578*x-9.33871*y (more about this equation at the end)

Code: Select all

#CameraSpeed = 1
Declare DrawPlane()

Structure vertices
  x.f
  y.f
  z.f
EndStructure

#texture = 1
Global plane.i, planeMesh.i
Global sphere.i
Global sphereMesh.i

Define.f KeyX, KeyY
Global Dim vertex.vertices(39) 

InitEngine3D()
InitSprite()
InitKeyboard()
InitMouse()

OpenWindow(0, 0, 0, 800, 600, "use the Mouse with the keyboard to navigate the scene ..... Esc to exit", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, 800, 600)

Add3DArchive(".", #PB_3DArchive_FileSystem)
Add3DArchive(#PB_Compiler_Home + "Examples/3D/Data/Textures", #PB_3DArchive_FileSystem)
Add3DArchive(#PB_Compiler_Home + "Examples/3D/Data/Scripts", #PB_3DArchive_FileSystem)
Parse3DScripts()

CreateMaterial(#texture, LoadTexture(#texture, "wood.jpg"))
MaterialCullingMode(#texture, #PB_Material_NoCulling)

CreateCamera(0, 0, 0, 100, 100)
MoveCamera(0, 28, 199, 0)
CameraLookAt(0, 0,196,0)

;CreateLight(#Light, Color [, x, y, z [, Flags]])
CreateLight(1, RGB(220,20,120), 0,200,30)

DrawPlane()

sphereMesh = CreateSphere(#PB_Any, 0.2, 8,8)
sphere = CreateEntity(#PB_Any, MeshID(sphereMesh), #PB_Material_None, 0, 0, 0)


For i=0 To 39 ; create 40 spheres to represent Data
    
  Read.f vertex(i)\x
  Read.f vertex(i)\y
  Read.f vertex(i)\z
   
    ;Debug vertex(i)\y
  CreateEntity(i, MeshID(sphereMesh), #PB_Material_None, vertex(i)\x, vertex(i)\y, vertex(i)\z)
  AttachEntityObject(sphere, "", EntityID(i))
Next
  
CreateLine3D(#PB_Any , -200, 0, 0, RGB(255,   0,   0), 200,  0,  0, RGB(255,   0,   0))  ; Axis X
CreateLine3D(#PB_Any , 0, -200, 0, RGB(  0, 255,   0),  0, 200,  0, RGB(  0, 255,   0))  ; Axis Y
CreateLine3D(#PB_Any , 0, 0, -200, RGB(  0,   0, 255),  0,  0, 200, RGB(  0,   0, 255))  ; Axis Z
    
    
Repeat
    Repeat
    event = WindowEvent()
  Until event = 0
  If ExamineMouse()
        MouseX = -MouseDeltaX() * #CameraSpeed * 0.2
        MouseY = -MouseDeltaY() * #CameraSpeed * 0.2
      EndIf
  ShowCursor_(0)
    
; Use arrow keys and mouse to rotate and fly in/out
  ExamineKeyboard()
      
        If KeyboardPushed(#PB_Key_Left)
          KeyX = -#CameraSpeed
        ElseIf KeyboardPushed(#PB_Key_Right)
          KeyX = #CameraSpeed
        Else
          KeyX = 0
        EndIf
        
        If KeyboardPushed(#PB_Key_Up)
          KeyY = -#CameraSpeed
        ElseIf KeyboardPushed(#PB_Key_Down)
          KeyY = #CameraSpeed
        Else
          KeyY = 0
        EndIf

        ;RotateEntity(sphere, 0,0.3,0, #PB_Relative)
        ;RotateEntity(plane, 0,0.3,0, #PB_Relative)

  RotateCamera(0, MouseY, MouseX, 0, #PB_Relative)
  MoveCamera  (0, KeyX, 0, KeyY)
      
  
  RenderWorld()
  
  FlipBuffers()
  
  ExamineKeyboard()
  
Until KeyboardPushed(#PB_Key_Escape)

Procedure DrawPlane()
  Protected.l a, b, Nb
  Protected.w P1, P2, P3, P4
  Protected.f x, z, y
  NbX=600
  ;NbZ=600
  NbY=600
  planeMesh = CreateMesh(#PB_Any, #PB_Mesh_LineStrip, #PB_Mesh_Static)
  
  ;xMin.f = -7.56 :xMax.f = 5.57: yMin.f = 196: yMax = 200
  xMin.f = -10 : xMax.f = 10: yMin.f = 196 : yMax = 200
  rangeX.f = xMax - xMin
  stepX.f = rangeX / NbX
  rangeY.f = yMax - yMin
  stepY.f = rangeY / NbY
  ;Debug StrF(stepY)
  x.f = xMin:  y.f = yMin
  
  For b=0 To NbY
    
    For a=0 To NbX
      
        z.f = 0.374578*x -9.33871*y +1842.08
                      
        MeshVertexPosition(x, y, z) 
        MeshVertexTextureCoordinate(a/NbX, b/NbY)
        x.f + stepX
    Next a
    
    x = xMin
    y.f + stepY
  Next b
  
  FinishMesh(#True)
  SetMeshMaterial(planeMesh, MaterialID(#texture))
  plane = CreateEntity(#PB_Any, MeshID(planeMesh), MaterialID(#texture))
  
EndProcedure

DataSection
Data.f -3.0,196.0,4.0, 3.0,197.0,5.0, 4.0,197.0,2.0, 4.0,197.0,-2.0, 1.0,197.0,-3.0, -4.0,197.0,-3.0, -4.0,197.0,0.0, 0.0,197.0,0.0, 1.0,197.0,3.0, -6.0,196.0,4.0, -1.14,197.22,1.38, -1.28,197.53,-2.66, -3.26,197.68,-5.63, 1.99,197.89,-5.41, 3.53,198.19,-8.33, -0.15,198.15,-9.88, -0.35,197.89,-6.74, -3.97,198.02,-10.3, -5.44,197.66,-6.62, -6.09,197.21,-1.28, 0.04,196.86,6.63, -3.36,196.6,8.05, -6.91,196.52,7.12, 2.26,197.48,0.0, 3.9,197.96,-5.2, 1.31,198.1,-8.44, -2.16,197.94,-8.35, -5.54,197.49,-4.44, -0.46,196.99,4.66, -0.56,197.73,-4.7, -6.43,197.0,1.3, -4.45,196.98,2.53, -5.13,196.66,6.28, -1.82,196.8,6.32, 4.24,197.57,-0.05, -7.56,196.7,4.43, -1.38,196.7,7.89, 2.26,196.84,8.14, 4.31,197.05,6.58, 5.57,197.35,3.49
EndDataSection

download this tutorial about how to find the equation of the best fit plane with mathematica v4.1
http://s000.tinyupload.com/index.php?fi ... 7702916022
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

Hi ApplePi

Many many thanks for your efforts, you never fail to amaze 8)

I can see that the function used is not exposed (potentially an iterative least-squares approach). However, it has to be part of my program, so I can't put it to use.

Concerning OpenGL v Ogre, it's easy to display the points in OpenGL, as per your examples on the forum (if required, usually just a just means to an end). Like CAD programs, it's also easy to display a representative plane, using a single square quad.

I have been experimenting with the removal of outliers (for the real world larger point sets). Works quite well by sorting the points for each axis and "trimming" off a percentage of them, as the worst cases are around the point cloud "edges". So, better data, but I still don't know how to define a best fit plane with them :(
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
applePi
Addict
Addict
Posts: 1404
Joined: Sun Jun 25, 2006 7:28 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by applePi »

Like CAD programs, it's also easy to display a representative plane, using a single square quad.
you are right, there is no need to draw thousands of point to draw a plane, just 4 points to draw it. since we know previously x and y of the points at the edges (it is our choice) we need only to find the z using the Equation
the site http://al-roomi.org/3DPlot/index.html can draw such equation to show the plane but the coordinates is like mathematica, why the people don't adhere to one opinion !!! i lost easily in the space, so any change in the coordinates directions makes my brain flip flop.
note that we need always to direct the camera to the objects correctly look lines 36-37 else we will see nothing
here is the code again :

Code: Select all

#CameraSpeed = 1

Declare DrawBetterPlane()

Structure vertices
  x.f
  y.f
  z.f
EndStructure

#texture = 1
Global plane.i, planeMesh.i
Global sphere.i
Global sphereMesh.i

Define.f KeyX, KeyY
Global Dim vertex.vertices(39) 

InitEngine3D()
InitSprite()
InitKeyboard()
InitMouse()

OpenWindow(0, 0, 0, 800, 600, "use the Mouse with the keyboard to navigate the scene ..... Esc to exit", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
OpenWindowedScreen(WindowID(0), 0, 0, 800, 600)

Add3DArchive(".", #PB_3DArchive_FileSystem)
Add3DArchive(#PB_Compiler_Home + "Examples/3D/Data/Textures", #PB_3DArchive_FileSystem)
Add3DArchive(#PB_Compiler_Home + "Examples/3D/Data/Scripts", #PB_3DArchive_FileSystem)
Parse3DScripts()

CreateMaterial(#texture, LoadTexture(#texture, "Geebee2.bmp"))
MaterialCullingMode(#texture, #PB_Material_NoCulling)

CreateCamera(0, 0, 0, 100, 100)
MoveCamera(0, 28, 199, 0)
CameraLookAt(0, 0,196,0)

;CreateLight(#Light, Color [, x, y, z [, Flags]])
CreateLight(1, RGB(220,20,120), 0,200,30)

DrawBetterPlane()

sphereMesh = CreateSphere(#PB_Any, 0.2, 8,8)
sphere = CreateEntity(#PB_Any, MeshID(sphereMesh), #PB_Material_None, 0, 0, 0)


For i=0 To 39 ; create 40 spheres to represent Data
    
  Read.f vertex(i)\x
  Read.f vertex(i)\y
  Read.f vertex(i)\z
   
    ;Debug vertex(i)\y
  CreateEntity(i, MeshID(sphereMesh), #PB_Material_None, vertex(i)\x, vertex(i)\y, vertex(i)\z)
  AttachEntityObject(sphere, "", EntityID(i))
Next
  
CreateLine3D(#PB_Any , -200, 0, 0, RGB(255,   0,   0), 200,  0,  0, RGB(255,   0,   0))  ; Axis X
CreateLine3D(#PB_Any , 0, -200, 0, RGB(  0, 255,   0),  0, 200,  0, RGB(  0, 255,   0))  ; Axis Y
CreateLine3D(#PB_Any , 0, 0, -200, RGB(  0,   0, 255),  0,  0, 200, RGB(  0,   0, 255))  ; Axis Z
    
    
Repeat
    Repeat
    event = WindowEvent()
  Until event = 0
  If ExamineMouse()
        MouseX = -MouseDeltaX() * #CameraSpeed * 0.2
        MouseY = -MouseDeltaY() * #CameraSpeed * 0.2
      EndIf
  ShowCursor_(0)
    
; Use arrow keys and mouse to rotate and fly in/out
  ExamineKeyboard()
      
        If KeyboardPushed(#PB_Key_Left)
          KeyX = -#CameraSpeed
        ElseIf KeyboardPushed(#PB_Key_Right)
          KeyX = #CameraSpeed
        Else
          KeyX = 0
        EndIf
        
        If KeyboardPushed(#PB_Key_Up)
          KeyY = -#CameraSpeed
        ElseIf KeyboardPushed(#PB_Key_Down)
          KeyY = #CameraSpeed
        Else
          KeyY = 0
        EndIf

        ;RotateEntity(sphere, 0,0.3,0, #PB_Relative)
        ;RotateEntity(plane, 0,0.3,0, #PB_Relative)

  RotateCamera(0, MouseY, MouseX, 0, #PB_Relative)
  MoveCamera  (0, KeyX, 0, KeyY)
      
  
  RenderWorld()
  
  FlipBuffers()
  
  ExamineKeyboard()
  
Until KeyboardPushed(#PB_Key_Escape)

Procedure DrawBetterPlane()
  Dim planeZ.f(3)
  Protected.f x, z, y
  x=-10: y=196
  planeZ(0)= 0.374578*x -9.33871*y +1842.08
  x=10: y=196
  planeZ(1)= 0.374578*x -9.33871*y +1842.08
  
  x=-10: y=200
  planeZ(2)= 0.374578*x -9.33871*y +1842.08
  x=10: y=200
  planeZ(3)= 0.374578*x -9.33871*y +1842.08
  
  planeMesh = CreateMesh(#PB_Any, #PB_Mesh_TriangleList, #PB_Mesh_Static) 
  MeshVertexPosition(-10, 196, planeZ(0)) ;index 0
  MeshVertexNormal(0,1,0) 
  MeshVertexTextureCoordinate(0, 1)
  
  MeshVertexPosition(10, 196, planeZ(1))  ; index 1
  MeshVertexNormal(0,1,0)
  MeshVertexTextureCoordinate(1, 1)
  
  MeshVertexPosition(-10, 200, planeZ(2)) ; index 2
  MeshVertexNormal(0,1,0)
  MeshVertexTextureCoordinate(0, 0)
  
  MeshVertexPosition(10, 200, planeZ(3))  ; index 3
  MeshVertexNormal(0,1,0)
  MeshVertexTextureCoordinate(1, 0)
  
  MeshFace(0, 1, 3)
  MeshFace(3, 2, 0)
  
  
  FinishMesh(#True)
  SetMeshMaterial(planeMesh, MaterialID(#texture))
  plane = CreateEntity(#PB_Any, MeshID(planeMesh), MaterialID(#texture))
EndProcedure

DataSection
Data.f -3.0,196.0,4.0, 3.0,197.0,5.0, 4.0,197.0,2.0, 4.0,197.0,-2.0, 1.0,197.0,-3.0, -4.0,197.0,-3.0, -4.0,197.0,0.0, 0.0,197.0,0.0, 1.0,197.0,3.0, -6.0,196.0,4.0, -1.14,197.22,1.38, -1.28,197.53,-2.66, -3.26,197.68,-5.63, 1.99,197.89,-5.41, 3.53,198.19,-8.33, -0.15,198.15,-9.88, -0.35,197.89,-6.74, -3.97,198.02,-10.3, -5.44,197.66,-6.62, -6.09,197.21,-1.28, 0.04,196.86,6.63, -3.36,196.6,8.05, -6.91,196.52,7.12, 2.26,197.48,0.0, 3.9,197.96,-5.2, 1.31,198.1,-8.44, -2.16,197.94,-8.35, -5.54,197.49,-4.44, -0.46,196.99,4.66, -0.56,197.73,-4.7, -6.43,197.0,1.3, -4.45,196.98,2.53, -5.13,196.66,6.28, -1.82,196.8,6.32, 4.24,197.57,-0.05, -7.56,196.7,4.43, -1.38,196.7,7.89, 2.26,196.84,8.14, 4.31,197.05,6.58, 5.57,197.35,3.49
EndDataSection
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

Hi ApplePi

That's a very handy website, you do have a knack of finding treasure!

Edit: The equation is not accurate if the plane is at a compound angle though ("Z" plane at angle in X and in Y), which is where an orthogonal calculation wins-out.....
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
alter Mann
User
User
Posts: 39
Joined: Fri Oct 17, 2014 8:52 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by alter Mann »

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.

For your example the code works.

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.
Last edited by alter Mann on Thu Jan 16, 2020 9:32 am, edited 4 times in total.
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

Hi alter Mann

That's a huge amount of code, thank you for your hard work. Straight away you can recognise the C heritage :)

The demo points set is solved very well. I'm going to try it with more point sets.

Thanks again.
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

Hi alter Mann

Fallen at the first hurdle - not sure if the code fails or the result is simply off-screen. Mind you, it's very late here, my brain is asleep :mrgreen:

POINTS.DAT
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
alter Mann
User
User
Posts: 39
Joined: Fri Oct 17, 2014 8:52 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by alter Mann »

Yes, a lot of code for such a "simple" problem :wink: .

I've changed the code above. It was a simple glOrtho problem. The result was correct but you couldn't see it. And I forgot a CloseFile.
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

That's brilliant alter Mann, thanks. I probably would not have spotted that, still in zombie mode at the moment.
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
alter Mann
User
User
Posts: 39
Joined: Fri Oct 17, 2014 8:52 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by alter Mann »

I played a little bit and added zoom,center and identification functions (code above).
Hope you enjoy :)
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

:mrgreen: It's going to take me a while to understand it.
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
alter Mann
User
User
Posts: 39
Joined: Fri Oct 17, 2014 8:52 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by alter Mann »

I've added some comments (not for calculation of eigenvectors - this is numerical math). Maybe it makes something clearer, maybe it confuse you more :mrgreen:. (I'm not native in this language)
IdeasVacuum
Always Here
Always Here
Posts: 6425
Joined: Fri Oct 23, 2009 2:33 am
Location: Wales, UK
Contact:

Re: Procedure to define a best-fit Plane through a set of po

Post by IdeasVacuum »

...that's much appreciated. Actually your code solves more problems than best-fit planes (working with huge sets of points now), the zooming is better than it is in my project and the ability to identify individual points is very useful too. 8)
IdeasVacuum
If it sounds simple, you have not grasped the complexity.
alter Mann
User
User
Posts: 39
Joined: Fri Oct 17, 2014 8:52 pm

Re: Procedure to define a best-fit Plane through a set of po

Post by alter Mann »

If you have questions about triangles, lines, points, etc., don't hesitate to ask. I have been working on such problems for years. :)
Post Reply