All code posted here is released to the Public Domain.
Code:
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
' BSurface.bas                    by Jordi Vallès      version 1g     16/09/2008
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
'  A Bézier surface is formed as the cartesian product of the blending functions
'  of two orthogonal Bézier curves.
'
'  The theory about splines, patches and surfaces created from Bézier curves I
'  have learned in Alan Watt book "3D Computer Graphics" 3rd edition, basically
'  the 4 first chapters.
'
'  This sample shows one or two Bézier patches. Each patch is defined as a matrix
'  of 4 x 4 points, and all operations are based over the points of this matrix.
'
'  Best explanation about Bézier Surfaces is in:
'  [URL="http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/AV0405/DONAVANIK/bezier.html"]Bézier Surfaces[/URL]
'  and Links referenced in this page, specially the written by Paul Bourke.
'
'  Some routines are adapted to PB from ones of P. Bourke, D. H. Frost and
'  others, mainly from educational sites.
'
'  Trackbars are created with PB Graphics. Code adapted from models provided by
'  Raffaello Bervini.
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
'  - All code posted here is released to the Public Domain.
'  - Program developed and tested with PowerBASIC for Windows (PB/Win 8.04) on
'    PC HP Pavilion m7760 1.80 GHz with Windows Vista Home Premium.
'  - Only PB Graphic package is used to draw information generated by program.
'  - Four Graphic windows are created and attached in this program.
'  - 16/09/2008 Recompiled and tested with PB/Win 9.0 without problems.
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
' To be solved:
'  - Continuity across the two patches.
'  - Some irregular behaviour on axis rotation.
'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
' SED_PBWIN

#Compile Exe "BSurface.exe"
#Dim All

#Include "Win32Api.inc"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%ID_GRAPHIC    = 1100
%ID_DRAWCP     = 1111
%ID_DRAWGR     = 1112
%ID_DRAWSH     = 1113
%ID_RESET      = 1114
%ID_EXIT       = 1115
%ID_PATCH1     = 1116
%ID_PATCH2     = 1117
%ID_INFO1      = 1118
%ID_INFO2      = 1119

%ID_XSTP       = 3001      'my controls starts always on 3000
%ID_YSTP       = 3002
%ID_ZSTP       = 3003

%LTRED         = &h5555FF???
%SLIDER        = &hEEFAFA???
%CANVAS        = &h400000???
$TITLE         = "Bèzier surface  -  1g"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Macro Pi       =  3.141592653589793#   'some comment?
Macro Pi180    =  0.017453292519943#   'Pi/180

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Type POLYARRAY
   cnt As Long
   x1 As Single  :  y1 As Single
   x2 As Single  :  y2 As Single
   x3 As Single  :  y3 As Single
   x4 As Single  :  y4 As Single
End Type

Type JCONTROL
   handler  As Dword
   idc      As Long
   vtot     As Long
   divis    As Long
   txt      As String * 4
End Type

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Global OldGraphicProc           As Dword        'old process, used in subclass
Global hDlg, hGraphic           As Dword
Global drawCP, drawGR, drawSH   As Long         'main switches for draw
Global rl, rh                   As Long         'graphic space dimension
Global X(), Y(), Z()            As Single       'control points or major nodes
Global pX(), pY()               As Long         'projected control points
Global bX(), bY(), bZ()         As Single       'compute patch vertices
Global pbX(), pbY()             As Long         'projected patch vertices
Global Trnf()                   As Single       'from 3D to 2D
Global numPatches               As Long         'number of patches, 1 or 2
Global steps, maxSteps, mouseSw As Long
Global eyeX, eyeY, eyeZ         As Single       'used to pass from 3D to 2D
Global factor                   As Single
Global currentI, currentJ       As Long         'current cursor position on graphic space
Global centerX, centerY         As Long         'start position of patch
Global svX, svY, svZ            As Long
Global jcount                   As Long
Global jCtl()                   As JCONTROL

Global mu00, mu01, mu02, mu10, mu11, mu12, mu20, mu21, mu22 As Single
Global sv00, sv01, sv02, sv10, sv11, sv12, sv20, sv21, sv22 As Single

'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Declare Function ComputeNodePoints(u As Single, v As Single, P() As Single, n As Long) As Single
Declare Function CreateTrackbarHControl(hWnd As Dword, id As Long, posx As Long, posy As Long, wide As Long, high As Long, _
                                        mn As Long, mx As Long, per As Long, divis As Long, txt As String) As Dword
Declare Function TrackbarHControl(hWnd As Dword, id As Long, posmx As Long, wantMouse As Long) As Long

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function PBMain()
   Local signature As Asciiz * 90

  '--- only one instance of this program is allowed ---
   signature = $TITLE + $TITLE + $TITLE
   If CreateMutex(ByVal 0, 1, signature) Then
      If GetLastError = %ERROR_ALREADY_EXISTS Then Exit Function
   End If

   Call InitialParms

   Dialog New Pixels, 0, $TITLE, , , 506, 604, %WS_CAPTION OR %WS_SYSMENU To hDlg

   Control Add Graphic,  hDlg, %ID_GRAPHIC,"", 3, 3, 500, 500, %SS_NOTIFY
   Control Get Client    hDlg, %ID_GRAPHIC To rl, rh    'get client (canvas) size
   Graphic Scale (0, 0) - (rl, rh)                      'scale to pixel coordinate system

   Control Add Button,   hDlg, %ID_RESET,   "Reset",                 436, 546,  66,  26
   Control Add Button,   hDlg, %ID_EXIT,    "Quit",                  436, 574,  66,  26

   Control Add Option,   hDlg, %ID_PATCH1,  "One Patch",             288, 506, 100,  16, %WS_GROUP
   Control Add Option,   hDlg, %ID_PATCH2,  "Two Patches",           288, 522, 100,  16
   Control Add CheckBox, hDlg, %ID_DRAWCP,  "Draw Control Points",   288, 546, 138,  16
   Control Add CheckBox, hDlg, %ID_DRAWGR,  "Draw Grid Detail",      288, 565, 138,  16
   Control Add CheckBox, hDlg, %ID_DRAWSH,  "Draw Gradient Shade",   288, 584, 138,  16

   Control Set Option    hDlg, %ID_PATCH1, %ID_PATCH1, %ID_PATCH2

   Control Add Label,    hDlg, %ID_INFO1,   "Pos:",                   20, 504,  80,  14
   Control Add Label,    hDlg, %ID_INFO2,   "",                      150, 504, 120,  14

   Control Set Check     hDlg, %ID_DRAWCP, 1
   Control Set Check     hDlg, %ID_DRAWGR, 0
   Control Disable       hDlg, %ID_DRAWSH

   '--- tracbars creation ---
   CreateTrackbarHControl(hDlg, %ID_XSTP,  4, 522, 260, 20, -100, 100, 0, 20, "X = ")
   CreateTrackbarHControl(hDlg, %ID_YSTP,  4, 550, 260, 20, -100, 100, 0, 20, "Y = ")
   CreateTrackbarHControl(hDlg, %ID_ZSTP,  4, 578, 260, 20, -100, 100, 0, 20, "Z = ")

   '--- main graphic activation ---
   Graphic Attach hDlg, %ID_GRAPHIC, Redraw                          'ATTACH
   Graphic Color %WHITE, %CANVAS

   Graphic Box (0, 0) - (rl, rh), 6, %WHITE, %CANVAS                 'my GRAPHIC CLEAR

   Dialog Show Modal hDlg Call DlgProc
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
CallBack Function DlgProc() As Long
   Local CtrlID As Long

   Select Case As Long CbMsg
      Case %WM_INITDIALOG
         Control Handle    CbHndl, %ID_GRAPHIC To hGraphic
         OldGraphicProc = SetWindowLong(hGraphic, %GWL_WNDPROC, CodePtr(GraphicProc))
         Control Get Check CbHndl, %ID_DRAWCP To drawCP
         Control Get Check CbHndl, %ID_DRAWGR To drawGR
         Sleep 100
         TrackbarHControl  CbHndl, %ID_XSTP, 30, %FALSE
         TrackbarHControl  CbHndl, %ID_YSTP,  0, %FALSE
         TrackbarHControl  CbHndl, %ID_ZSTP,  0, %FALSE
         Sleep 100
         Call InitControlPoints
         Call RotateAxisX(30)
         Call PaintPatch(%TRUE)

      Case %WM_SETCURSOR
         If CbWParam <> hGraphic Then
            Control Set Text CbHndl, %ID_INFO1, "Pos:"
            CtrlID = GetDlgCtrlID(CbWParam) : GetAsyncKeyState(%VK_LBUTTON)
            Select Case CtrlID
               Case %ID_XSTP To %ID_ZSTP
                  If GetAsyncKeyState(%VK_LBUTTON) <> 0 Then
                     Select Case CtrlId
                        Case %ID_XSTP
                           Call RotateAxisX(-svX)
                           Call RotateAxisX(TrackbarHControl(CbHndl, CtrlID, 0, %TRUE))
                        Case %ID_YSTP
                           Call RotateAxisY(-svY)
                           Call RotateAxisY(TrackbarHControl(CbHndl, CtrlID, 0, %TRUE))
                        Case %ID_ZSTP
                           Call RotateAxisZ(-svZ)
                           Call RotateAxisZ(TrackbarHControl(CbHndl, CtrlID, 0, %TRUE))
                     End Select
                  End If
                  Call PaintPatch(%TRUE)
            End Select
         End If

      Case %WM_COMMAND
         If CbCtlMsg <> %BN_CLICKED Then Exit Select
         Select Case CbCtl
            Case %ID_PATCH1
               numPatches = 1
               Control Send CbHndl, %ID_RESET, %BM_CLICK, 0, 0
            Case %ID_PATCH2
               numPatches = 2
               Control Send CbHndl, %ID_RESET, %BM_CLICK, 0, 0
            Case %ID_DRAWCP
               Control Get Check CbHndl, %ID_DRAWCP To drawCP
               Call PaintPatch(%TRUE)
            Case %ID_DRAWGR
               Control Get Check CbHndl, %ID_DRAWGR To drawGR
               If drawGR Then
                  Control Enable hDlg, %ID_DRAWSH
               Else
                  Control Disable hDlg, %ID_DRAWSH
               End If
               Call PaintPatch(%TRUE)
            Case %ID_DRAWSH
               Control Get Check CbHndl, %ID_DRAWSH To drawSH
               Call PaintPatch(%TRUE)
            Case %IDCANCEL, %ID_EXIT
               Dialog End CbHndl, 0
               Function = 1 : Exit Function
            Case %ID_RESET
               Control Set Check CbHndl, %ID_DRAWCP, 1  :  drawCP = 1
               Control Set Check CbHndl, %ID_DRAWGR, 0  :  drawGR = 0
               Control Set Check CbHndl, %ID_DRAWSH, 0  :  drawSH = 0
               Call InitControlPoints
               Call RotateAxisX(30)
               Call PaintPatch(%TRUE)
               TrackbarHControl CbHndl, %ID_XSTP, 30, %FALSE
               TrackbarHControl CbHndl, %ID_YSTP,  0, %FALSE
               TrackbarHControl CbHndl, %ID_ZSTP,  0, %FALSE
         End Select

      Case %WM_SYSCOMMAND
         If CbWParam = %SC_CLOSE Then
            Dialog End CbHndl, 0
            Function = 1
            Exit Function
         End If

      Case %WM_NCACTIVATE
         Static hWndsvFocus As Dword
         If IsFalse CbWParam Then
            hWndsvFocus = GetFocus()
         ElseIf hWndsvFocus Then
            SetFocus(hWndsvFocus)
            hWndsvFocus = 0
         End If

      Case %WM_DESTROY
         If OldGraphicProc Then SetWindowLong hGraphic, %GWL_WNDPROC, OldGraphicProc
   End Select
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function GraphicProc(Byval hWnd As Dword, Byval wMsg As Dword, Byval wParam As Dword, Byval lParam As Long) As Long
   Select Case wMsg
      Case %WM_LBUTTONDOWN  :  Call LeftButtonDown (hWnd, wParam, lParam)
      Case %WM_MOUSEMOVE    :  Call MouseMove      (hWnd, wParam, lParam)
      Case %WM_LBUTTONUP    :  Call LeftButtonUp   (hWnd, wParam, lParam)
   End Select
   Function = CallWindowProc(OldGraphicProc, hWnd, wMsg, wParam, lParam)
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub LeftButtonDown(Byval hWnd As Dword, Byval wParam As Dword, Byval lParam As Long)
   Register i As Long, j As Long, t AS Long
   Local pt As POINTAPI

   GetCursorPos pt
   ScreenToClient hWnd, pt
   mouseSw = %TRUE

   'look for a match
   For t = 0 To 3
      For i = 0 To numPatches*3
         For j = 0 To 3
            If ((pt.x-t <= pX(i,j)) And (pX(i,j) <= pt.x+t) And (pt.y-t <= pY(i,j)) And (pY(i,j) <= pt.y+t)) Then
               currentI = i  :  currentJ = j
               Control Set Text hDlg, %ID_INFO2, "Node" + Str$(i) + Str$(j) + " caught"
               Sleep 1
               Exit Sub
            End If
         Next j
      Next i
   Next k
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub LeftButtonUp(Byval hWnd As Dword, Byval wParam As Dword, Byval lParam As Long)
   currentI = -1
   Call PaintPatch(%TRUE)
   Control Set Text hDlg, %ID_INFO2, Space$(20)
   Control Set Text hDlg, %ID_INFO1, "Pos:"
   mouseSw = %FALSE
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub MouseMove(Byval hWnd As Dword, Byval wParam As Dword, Byval lParam As Long)
   Local pt As POINTAPI

   If mouseSw Then Control Set Text hDlg, %ID_INFO1, "Pos:" + Str$(Lo(Integer, lParam)) + "," + Str$(Hi(Integer, lParam))
   If currentI = -1 Then Exit Sub
   GetCursorPos pt
   ScreenToClient hWnd, pt

    X(currentI,currentJ) = eyeX + ((pt.x-centerX) * (Z(currentI,currentJ)-eyeZ) / -factor)
   Y(currentI,currentJ) = eyeY + ((centerY-pt.y) * (Z(currentI,currentJ)-eyeZ) / -factor)
   Call PaintPatch(%TRUE)
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub PaintPatch(wantErase As Long)
   Register i As Long, j As Long, n As Long
   Local u, v, step1 As Single
   Local u1, v1 As Long
   Local poly As POLYARRAY

   Graphic Attach hDlg, %ID_GRAPHIC, Redraw        'get graphic control, ATTACH

   If wantErase Then Graphic Box (0, 0) - (rl, rh), 6, %WHITE, %CANVAS      'my GRAPHIC CLEAR

   Call DrawControlPoints

   If currentI = -1 Then
      For n = 0 to numPatches-1
         step1 = 1.0 / steps
         u1 = 0  :  v1 = 0

         'create space subdivisions from Control Points
         'the result is the green grid (4.2.2)
         For u = 0 To 1.0001 Step step1
            For v = 0 To 1.0001 Step step1
               'compute node points
               bX(u1,v1) = ComputeNodePoints(u, v, X(), n*3)
               bY(u1,v1) = ComputeNodePoints(u, v, Y(), n*3)
               bZ(u1,v1) = ComputeNodePoints(u, v, Z(), n*3)
               Incr v1
            Next v
            v1 = 0
            Incr u1
         Next u

         'transform from 3D to 2D
         For i = 0 To steps
            For j = 0 To steps
               Call Interpolation(bX(i,j), bY(i,j), bZ(i,j))
               pbX(i,j) = ((Trnf(0)-eyeX) / (Trnf(2)-eyeZ) * -factor + 0.5) + centerX
               pbY(i,j) = centerY - ((Trnf(1)-eyeY) / (Trnf(2)-eyeZ) * -factor + 0.5)
            Next j
         Next i

         'draw wireframe, if selected
         If drawGR Then
            For i = 0 to steps
               For j = 0 To steps
                  If i < steps Then
                     Graphic Line (pbX(i,j), pbY(i,j)) - (pbX(i+1,j), pbY(i+1,j)), %GREEN
                  End If
                  If j < steps Then
                     Graphic Line (pbX(i,j), pbY(i,j)) - (pbX(i,j+1), pbY(i,j+1)), %GREEN
                  End If
               Next j
            Next i

            'draw gradient shade, if selected
            'needs drawGR switch on (wireframe)
            If drawSH Then
               poly.cnt = 4
               For i = 0 to steps
                  For j = 0 To steps
                     If i < steps And j < steps Then
                        poly.x1 = pbX(i,j)     :  poly.y1 = pbY(i,j)
                        poly.x2 = pbX(i+1,j)   :  poly.y2 = pbY(i+1,j)
                        poly.x4 = pbX(i,j+1)   :  poly.y4 = pbY(i,j+1)
                        poly.x3 = pbX(i+1,j+1) :  poly.y3 = pbY(i+1,j+1)
                        u = 140 + bz(i,j)*10
                        Graphic Polygon poly, RGB(u/2,u,20), RGB(u/2,u,0)
                     End If
                  Next j
               Next i
            End If
         End If

      Next n
   End If
   Graphic Redraw

End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub DrawControlPoints()
   Register i As Long, j As Long

   If drawCP Then
      For i = 0 To numPatches * 3
         For j = 0 To 3
            'calculate control points positions
            'transform from 3D to 2D and after display it
            Call Interpolation(X(i,j), Y(i,j), Z(i,j))
            pX(i,j) = ((Trnf(0) - eyeX) / (Trnf(2) - eyeZ) * -factor + 0.5) + centerX
            pY(i,j) = centerY - ((Trnf(1) - eyeY) / (Trnf(2) - eyeZ) * -factor + 0.5)
         Next j
      Next i
      For i = 0 to numPatches * 3
         For j = 0 to 3
            If (i < numPatches*3) Then
               Graphic Line (pX(i,j), pY(i,j)) - (pX(i+1,j), pY(i+1,j)), %GRAY
            End If
            If (j < 3) Then
               Graphic Line (pX(i,j),   pY(i,j)) - (pX(i,j+1), pY(i,j+1)), %GRAY
            End If
            Graphic Box (pX(i,j)-2, pY(i,j)-2) - (pX(i,j)+2, pY(i,j)+2) , , -1, -1
         Next j
      Next i
   End If
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function ComputeNodePoints(u0 As Single, v0 As Single, P() As Single, n As Long) As Single
   'u0 and v0 must be values between 0 and 1

   'Calculate the heights of Control Points of a Bézier patch using the Bernstein cubic
   'polynomials method, based on descriptions (3.1 and 3.4) of book referenced. Also
   'some samples and comments of Paul Bourke and Daniel H. Frost, help me a lot.
   'Not is the best but is faster.

   'This routine is called 1452 times per patch for each change on Control Point matrix,
   'just when the left button is released. See PaintPatch subroutine.
   '           1452 = 3 * 484 times for each axis
   '            484 = 22 * 22 steps to form the grid

   'To investigate: Cubic Bézier curve used in Postscript, developed by Adobe,
   '                according to some authors is the best and fast method.

   Local u1, v1, v1squared, v1cubed, v0squared, v0cubed As Single

   u1        = 1 - u0
   v1        = 1 - v0
   v1squared = v1*v1
   v1cubed   = v1*v1squared
   v0squared = v0*v0
   v0cubed   = v0*v0squared

   Function = _
      u1*u1*u1 * _
         (P(n+0,0) * v1cubed + 3*P(n+0,1) * v1squared * v0 + 3*P(n+0,2) * v1 * v0squared + P(n+0,3) * v0cubed) + _
      3*u1*u1*u0 * _
         (P(n+1,0) * v1cubed + 3*P(n+1,1) * v1squared * v0 + 3*P(n+1,2) * v1 * v0squared + P(n+1,3) * v0cubed) + _
      3*u1*u0*u0 * _
         (P(n+2,0) * v1cubed + 3*P(n+2,1) * v1squared * v0 + 3*P(n+2,2) * v1 * v0squared + P(n+2,3) * v0cubed) + _
      u0*u0*u0 * _
         (P(n+3,0) * v1cubed + 3*P(n+3,1) * v1squared * v0 + 3*P(n+3,2) * v1 * v0squared + P(n+3,3) * v0cubed)
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub InitialParms()
   numPatches = 1
   centerX = 250  :  centerY = 200
   maxSteps = 50
   currentI = -1    '-1 means no control point is current

   Redim X(0 To 6, 0 To 3) As Single
   Redim Y(0 To 6, 0 To 3) As Single
   Redim Z(0 To 6, 0 To 3) As Single
   Redim pX(0 To 6, 0 To 3) As Long
   Redim pY(0 To 6, 0 To 3) As Long
   Redim Trnf(0 to 2) As Single
   Redim bx(0 To maxSteps, 0 To maxSteps)
   Redim by(0 To maxSteps, 0 To maxSteps)
   Redim bz(0 To maxSteps, 0 To maxSteps)
   Redim pbx(0 To maxSteps, 0 To maxSteps)
   Redim pby(0 To maxSteps, 0 To maxSteps)
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub InitControlPoints()
   Register j As Long, i As Long

   steps      = 21
   factor     = 600
   eyeX = 0  :  eyeY = 5  :  eyeZ = 50

   'set mu ranges
   mu00 = 1.0  :  mu11 = 1.0  :  mu22 = 1.0
   mu01 = 0.0  :  mu02 = 0.0  :  mu10 = 0.0
   mu12 = 0.0  :  mu20 = 0.0  :  mu21 = 0.0

   'control points initial values
   X(0,0) = -30  :  Y(0,0) = 0  :  Z(0,0) =   15
   X(0,1) = -30  :  Y(0,1) = 0  :  Z(0,1) =    5
   X(0,2) = -30  :  Y(0,2) = 0  :  Z(0,2) =   -5
   X(0,3) = -30  :  Y(0,3) = 0  :  Z(0,3) =  -15
   X(1,0) = -20  :  Y(1,0) = 0  :  Z(1,0) =   15
   X(1,1) = -20  :  Y(1,1) = 0  :  Z(1,1) =    5
   X(1,2) = -20  :  Y(1,2) = 0  :  Z(1,2) =   -5
   X(1,3) = -20  :  Y(1,3) = 0  :  Z(1,3) =  -15
   X(2,0) = -10  :  Y(2,0) = 0  :  Z(2,0) =   15
   X(2,1) = -10  :  Y(2,1) = 0  :  Z(2,1) =    5
   X(2,2) = -10  :  Y(2,2) = 0  :  Z(2,2) =   -5
   X(2,3) = -10  :  Y(2,3) = 0  :  Z(2,3) =  -15
   X(3,0) =   0  :  Y(3,0) = 0  :  Z(3,0) =   15
   X(3,1) =   0  :  Y(3,1) = 0  :  Z(3,1) =    5
   X(3,2) =   0  :  Y(3,2) = 0  :  Z(3,2) =   -5
   X(3,3) =   0  :  Y(3,3) = 0  :  Z(3,3) =  -15
   X(4,0) =  10  :  Y(4,0) = 0  :  Z(4,0) =   15
   X(4,1) =  10  :  Y(4,1) = 0  :  Z(4,1) =    5
   X(4,2) =  10  :  Y(4,2) = 0  :  Z(4,2) =   -5
   X(4,3) =  10  :  Y(4,3) = 0  :  Z(4,3) =  -15
   X(5,0) =  20  :  Y(5,0) = 0  :  Z(5,0) =   15
   X(5,1) =  20  :  Y(5,1) = 0  :  Z(5,1) =    5
   X(5,2) =  20  :  Y(5,2) = 0  :  Z(5,2) =   -5
   X(5,3) =  20  :  Y(5,3) = 0  :  Z(5,3) =  -15
   X(6,0) =  30  :  Y(6,0) = 0  :  Z(6,0) =   15
   X(6,1) =  30  :  Y(6,1) = 0  :  Z(6,1) =    5
   X(6,2) =  30  :  Y(6,2) = 0  :  Z(6,2) =   -5
   X(6,3) =  30  :  Y(6,3) = 0  :  Z(6,3) =  -15

   If numPatches = 1 Then   'if only 1 patch the first half is needed must be centered
      For i = 0 to 3
         For j = 0 to 3
            X(i,j) = X(i,j) +15
         Next j
      Next i
   End If

   For i = 0 to 6           'if only 1 patch only the half is useful
      For j = 0 to 3
         X(i,j) = X(i,j) / 2
         Z(i,j) = Z(i,j) / 2
      Next j
   Next i
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub Interpolation(xx As Single, yy As Single, zz As Single)
   'three control point Bezier interpolation
    'mu ranges from 0 to 1, start to end of the curve
   'adapted from Paul Bourke
   Trnf(0) = xx * mu00 + yy * mu10 + zz * mu20
   Trnf(1) = xx * mu01 + yy * mu11 + zz * mu21
   Trnf(2) = xx * mu02 + yy * mu12 + zz * mu22
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub RotateAxisX(degrees As Long)
   Local radians, cn, sn As Single
   svX = degrees
   radians = degrees * Pi180
   sn = Sin(radians)
   cn = Cos(radians)
   Call MuCopyToSave
   mu01 = sv01 * cn + sv02 * -sn
   mu02 = sv01 * sn + sv02 *  cn
   mu11 = sv11 * cn + sv12 * -sn
   mu12 = sv11 * sn + sv12 *  cn
   mu21 = sv21 * cn + sv22 * -sn
   mu22 = sv21 * sn + sv22 *  cn
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub RotateAxisY(degrees As Long)
   Local radians, cn, sn As Single
   svY = degrees
   radians = degrees * Pi180
   sn = Sin(radians)
   cn = Cos(radians)
   Call MuCopyToSave
   mu00 = sv00 *  cn + sv02 * sn
   mu02 = sv00 * -sn + sv02 * cn
   mu10 = sv10 *  cn + sv12 * sn
   mu12 = sv10 * -sn + sv12 * cn
   mu20 = sv20 *  cn + sv22 * sn
   mu22 = sv20 * -sn + sv22 * cn
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub RotateAxisZ(degrees As Long)
   Local radians, cn, sn As Single
   svZ = degrees
   radians = degrees * Pi180
   sn = Sin(radians)
   cn = Cos(radians)
   Call MuCopyToSave
   mu00 = sv00 * cn + sv01 * -sn
   mu01 = sv00 * sn + sv01 *  cn
   mu10 = sv10 * cn + sv11 * -sn
   mu11 = sv10 * sn + sv11 *  cn
   mu20 = sv20 * cn + sv21 * -sn
   mu21 = sv20 * sn + sv21 *  cn
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Sub MuCopyToSave()
   sv00 = mu00  :  sv01 = mu01  :  sv02 = mu02
   sv10 = mu10  :  sv11 = mu11  :  sv12 = mu12
   sv20 = mu20  :  sv21 = mu21  :  sv22 = mu22
End Sub

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function CreateTrackbarHControl(hWnd As Dword, id As Long, posx As Long, posy As Long, wide As Long, high As Long, _
                             mn As Long, mx As Long, per As Long, divis As Long, txt As String) As Dword
   Local rll, rhh As Long

   Incr jcount
   Redim Preserve jCtl(1 To jcount) As JCONTROL
   jCtl(jcount).idc   = id -3000                          'my controls starts always on 3000
   jCtl(jcount).vtot  = Abs(mx-mn)                        'range
   jCtl(jcount).divis = divis                             'number of divisions
   jCtl(jcount).txt   = txt                               'text to display inside trackbar
   Control Add Graphic, hWnd, id, "", posx, posy, wide, high, %SS_NOTIFY
   Control Get Client   hWnd, id To rll, rhh              'get client (canvas) size
   Graphic Scale (0, 0) - (rll, rhh)                      'scale to pixel coordinate
   Graphic Attach hWnd, id, ReDraw                        'ATTACH
   Graphic Font "Verdana", 8, 0
   Graphic Width 3
   Graphic Box (0,0) - (wide, high), ,%GRAY, %SLIDER
   Graphic ReDraw
   Control Handle hWnd, id To jCtl(jcount).handler

   FUNCTION = jCtl(jcount).handler
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
Function TrackbarHControl(hWnd As Dword, id As Long, posmx As Long, wantMouse As Long) As Long
   Local pt As POINTAPI
   Local handler As Dword
   Local d As Single
   Local j, k, v, mx, my, wide, high, posx, posy, idc, sx, sy As Long
   Local txt As String

   idc = id -3000                                         'my controls starts always on 3000
   handler = jCtl(idc).handler                            'get appropriate trackbar handler
   v = jCtl(idc).vtot                                     'range to be covered by slider

   Graphic Attach handler, id, ReDraw                     'get graphic control, ATTACH
   Graphic Get Client To wide, high
   Control Get Loc jCtl(idc).handler, id To posx, posy

   If wantMouse Then                                      'mouse position requested?
      GetCursorPos pt
      ScreenToClient handler, pt
      mx = pt.x-posx  :  my = pt.y-posy
      Control Set Text hDlg, %ID_INFO1, "Pos:" + Str$(mx) + "," + Str$(my)
      txt = jCtl(idc).txt + Ltrim$(Str$(Round (Min(v, v*(pt.x-posx) / wide),0) -v/2))
      FUNCTION = Round(Min(v, v*(pt.x-posx) / wide),0) -v/2
   Else                                                   'position is passed as parameter
      d = wide / v
      mx = (wide/2) + (posmx*d)
      txt = jCtl(idc).txt + Ltrim$(Str$(posmx))
      FUNCTION = posmx
   End If

   Graphic Width 3                                            'draw frame and
   Graphic Box (0,0) - (wide, high), ,%GRAY, %SLIDER          'clear previous info

   Graphic Width 1
   Graphic Box (mx-7, posy+3) - (mx+8, posy+high-3), ,%GRAY   'slider box
   Graphic Line (mx, posy+3) - (mx, posy+high-3), %LTRED      'slider line

   Graphic Text Size txt To sx, sy                            'set text
   Graphic Set Pos (wide/2-sx/2, high/2-sy/2-2)               '
   Graphic Color %BLUE, -2                                    '
   Graphic Print txt                                          '

   If jCtl(idc).divis > 0 Then                                'divisions required?
      k = (wide-6) / jCtl(idc).divis
      For j = 1 To jCtl(idc).divis -1
         Graphic Line (posx+k*j, posy+high-6) - (posx+k*j, posy+high-2)
      Next j
   End If

   Graphic Redraw                                             'draw all
End Function

'¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
'eof