welcome

Welcome to GFA-Basic 32

Geometrics algorithms Help File !

This help file mainly includes procedures and functions written in GFA-Basic 32 concerning geometry but not only.
I wrote this help in order to compile a maximum of algorithms dedicated particularly to geometry. You will find each time an example of use. In all cases, the functions and procedures are not necessarily suitable for all forms of use. I currently use UDTs or CONSTANTs that are specific to me, the list of which you will find here : Specific data type and constants.
For general UDT manipulations look here : Manipulating specifics UDT.
I can only advise you to adapt them to your own way of coding and your specific UDTs.

I am far from having finished gathering all my codes, so I should publish a second more complete version during the year 2023.
In the meantime, dear friends, I wish you happy reading.

Thank's :
- Sjoule Hamstra, the one and only current programmer capable of developing and improving the GFA-Basic suite !
- Frank Ostrowki the creator of GFA-Basic.
- Jean-Marie Melanson for HSL spectrum correction, and all other many fondamental algorithm.
- Peter Heinzig for the fabulous and very effective first opengl demo and directsound demo.
- All contributors to the https://gb32.proboards.com/

Nicolas Rey Nicolas REY
France - 24/07/2023
Contact : scalion@free.fr

{Created by Nicolas Rey : 10/07/2023}


UDT and CONSTANTS

Specific data type and constants

UDT

Type LINE3D_STRUCT
  - XYZ_STRUCT p1, p2
EndType

Type XYZ_STRUCT
  - Double x, y, z
EndType

Type TRIANGLE3D_STRUCT
  - XYZ_STRUCT a, b, c
EndType

Type QUAD3D_STRUCT
  - XYZ_STRUCT a, b, c, d
EndType

Type SPHERE_STRUCT Packed 1
  - XYZ_STRUCT Center
  - Double Radius
EndType

CONSTANTS

Global Const REFRACTION_AIR = 1.000272
Global Const REFRACTION_EAU = 1.325
Global Const REFRACTION_DIAMANT = 2.4175


Manipulating specifics UDT

Proc Translate_QUAD3D(ByRef Q As QUAD3D_STRUCT, x#, y#, z#)

  Translate_XYZ Q.a, x, y, z
  Translate_XYZ Q.b, x, y, z
  Translate_XYZ Q.c, x, y, z
  Translate_XYZ Q.d, x, y, z

EndProc
Proc Translate_XYZ(ByRef P As XYZ_STRUCT, x#, y#, z#)

  P.x += x#
  P.y += y#
  P.z += z#

EndProc
Proc Turn_QUAD3D_Y(ByRef Q As QUAD3D_STRUCT, Angle)

  Turn_XYZ_Y Q.a, Angle
  Turn_XYZ_Y Q.b, Angle
  Turn_XYZ_Y Q.c, Angle
  Turn_XYZ_Y Q.d, Angle

EndProc
Proc Turn_XYZ_Y(ByRef P As XYZ_STRUCT, Angle)

  Turn P.x, P.z, Angle

EndProc
Proc Turn_QUAD3D_X(ByRef Q As QUAD3D_STRUCT, Angle)

  Turn_XYZ_X Q.a, Angle
  Turn_XYZ_X Q.b, Angle
  Turn_XYZ_X Q.c, Angle
  Turn_XYZ_X Q.d, Angle

EndProc
Proc Turn_XYZ_X(ByRef P As XYZ_STRUCT, Angle)

  Turn P.y, P.z, Angle

EndProc
Function XYZ_Sub(ByRef a As XYZ_STRUCT, ByRef b As XYZ_STRUCT) As XYZ_STRUCT

  XYZ_Sub.x = a.x - b.x
  XYZ_Sub.y = a.y - b.y
  XYZ_Sub.z = a.z - b.z

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc


Rotate a point

This procedure turn x,y arround 0,0 by radian angle a.

Proc Turn(ByRef x#, ByRef y#, a#)

  Static Double s, c, newx, newy

  s = Sin(a)
  c = Cos(a)
  newx = c * x + s * y
  newy = c * y - s * x
  x = newx
  y = newy

EndProc

Example :

Global xpoint As Double, Ypoint As Double, i As Long

OpenW 1
Do
  Sleep
Loop Until Me Is Nothing
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  Cls
  Circle x, y, 10
  xpoint = x - _X / 2
  Ypoint = y - _Y / 2

  ' Turning point of 0.4 radian 5 times
  For i = 1 To 5
    Turn xpoint, Ypoint, 0.4
    Line 0, _Y / 2, _X, _Y / 2
    Line _X / 2, 0, _X / 2, _Y
    PCircle _X / 2 + xpoint, _Y / 2 + Ypoint, 10
  Next i

EndSub


Transpose using a normal

Transpose point using a normal

This function performs the necessary 3D rotations so that a point defined in a space centered in (0,0,0) with a normal (0,1,0) corresponds to the otherwise normalized space.

Function XYZ_TranposeByNormale(ByRef Normale As XYZ_STRUCT, ByRef P As XYZ_STRUCT) As XYZ_STRUCT

  Static Double xz, s, c, z

  xz = Sqr(Normale.x * Normale.x + Normale.z * Normale.z)

  If xz <> 0
    s = Normale.x / xz
    c = Normale.z / xz
    z = Normale.y * P.z + xz * P.y
    XYZ_TranposeByNormale.y = Normale.y * P.y - xz * P.z
    XYZ_TranposeByNormale.x = c * P.x + s * z
    XYZ_TranposeByNormale.z = c * z - s * P.x
  Else
    XYZ_TranposeByNormale.z = Normale.y * P.z
    XYZ_TranposeByNormale.y = Normale.y * P.y
    XYZ_TranposeByNormale.x = P.x
  EndIf

EndFunc

Example

$Library "Direct2D"

Type XYZ_STRUCT
  - Double x, y, z
EndType

OpenW Center 1

Global Double x, y, cx, cy, cz, r, a, x1, y1, x2, y2, z2
Global Normale As XYZ_STRUCT
Global PA As XYZ_STRUCT, PB As XYZ_STRUCT
Global Cube(0 .. 7) As XYZ_STRUCT
Global CubeX(0 .. 7) As Double, CubeY(0 .. 7) As Double
Global Long i

Cube(0) = XYZ(-1, -1, -1)
Cube(1) = XYZ(-1, -1,  1)
Cube(2) = XYZ( 1, -1,  1)
Cube(3) = XYZ( 1, -1, -1)
Cube(4) = XYZ(-1,  1, -1)
Cube(5) = XYZ(-1,  1,  1)
Cube(6) = XYZ( 1,  1,  1)
Cube(7) = XYZ( 1,  1, -1)

Global Object Win1RT
Set Win1RT = D2GetRT(Win_1)

Ocx Timer FPS
FPS.Enabled = True
FPS.Interval = 20

Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer

  Static work = False

  If !work
    work = True
    DrawAll
    work = False
  EndIf

EndSub
Proc DrawAll

  ' Initialise a Direct2D drawing session with a dark blue background
  D2BeginDraw Win1RT, D2C_DarkBlue

  ' Compute 3D center and size (r) of figure.
  cx = _X / 2 : cy = _Y / 2 : cz = _Y / 2 + 300
  r = cz / 2

  ' Define a normal for transpose the space with the mouse position
  a = PI - ((MouseX - cx) / cx) * PI
  Normale = XYZ(Sin(a), (MouseY - cy) / 50,  Cos(a))
  XYZ_Normalise Normale

  ' Draw a representation of this normal with an arrow
  ' Pseudo3D is a simplist algorithm to simulate a 3D view
  Pseudo3D cx, cy, cz, x1, y1
  Pseudo3D cx + Normale.x * r, cy + Normale.y * r, x1 + Normale.z * r, x2, y2
  D2ForeColor = D2C_White
  D2DefLine , 2
  DrawArrow x1, y1, x2, y2

  ' Transpose the cube with the normale and draw it
  D2DefLine , 1
  D2ForeColor = D2C_Orange
  For i = 0 To 7
    PA = XYZ_TranposeByNormale(Normale, Cube(i))
    PA.x = cx + PA.x * r
    PA.y = cy + PA.y * r
    PA.z = cz + PA.z * r
    Pseudo3D PA.x, PA.y, PA.z, CubeX(i), CubeY(i)
  Next i
  D2Line CubeX(0), CubeY(0), CubeX(1), CubeY(1)
  D2Line CubeX(1), CubeY(1), CubeX(2), CubeY(2)
  D2Line CubeX(2), CubeY(2), CubeX(3), CubeY(3)
  D2Line CubeX(0), CubeY(0), CubeX(3), CubeY(3)
  D2Line CubeX(0), CubeY(0), CubeX(4), CubeY(4)
  D2Line CubeX(1), CubeY(1), CubeX(5), CubeY(5)
  D2Line CubeX(2), CubeY(2), CubeX(6), CubeY(6)
  D2Line CubeX(3), CubeY(3), CubeX(7), CubeY(7)
  D2Line CubeX(4), CubeY(4), CubeX(5), CubeY(5)
  D2Line CubeX(5), CubeY(5), CubeX(6), CubeY(6)
  D2Line CubeX(6), CubeY(6), CubeX(7), CubeY(7)
  D2Line CubeX(4), CubeY(4), CubeX(7), CubeY(7)

  ' Draw a few moving particles to understand the spatial transposition.
  ' Using internal GFABasic randomize is sometimes very useful.
  D2DefLine , 0
  Randomize 1.1
  x2 = Frac(Timer / 19.5)
  y2 = Frac(Timer / 20)
  z2 = Frac(Timer / 19)
  For i = 1 To 100
    PA = XYZ(-1 + Abs(2 - 4 * Frac(x2 + Rnd)), -1 + Abs(2 - 4 * Frac(y2 + Rnd)), -1 + Abs(2 - 4 * Frac(z2 + Rnd)))
    PB = XYZ_TranposeByNormale(Normale, PA)
    Pseudo3D PB.x * r + cx, PB.y * r + cy, PB.z * r + cz, x, y
    D2PCircle x, y, 4 - 2 * PB.z, , , , D2Brush(ARGB(100, 255, 255, 0))
  Next i

  ' Drawing !
  D2EndDraw

EndProc
Proc  Pseudo3D(x#, y#, z#, ByRef x3d#, ByRef y3d#)

  Static Double dx, dy, d

  d = 1 / (1 + Abs(z) / 500)
  dx = x - _X / 2
  dy = y - _Y / 2
  x3d = _X / 2 + dx * d
  y3d = _Y / 2 + dy * d

EndProc
Function XYZ_TranposeByNormale(ByRef Normale As XYZ_STRUCT, ByRef P As XYZ_STRUCT) As XYZ_STRUCT

  Static Double xz, s, c, z

  xz = Sqr(Normale.x * Normale.x + Normale.z * Normale.z)

  If xz <> 0
    s = Normale.x / xz
    c = Normale.z / xz
    z = Normale.y * P.z + xz * P.y
    XYZ_TranposeByNormale.y = Normale.y * P.y - xz * P.z
    XYZ_TranposeByNormale.x = c * P.x + s * z
    XYZ_TranposeByNormale.z = c * z - s * P.x
  Else
    XYZ_TranposeByNormale.z = Normale.y * P.z
    XYZ_TranposeByNormale.y = Normale.y * P.y
    XYZ_TranposeByNormale.x = P.x
  EndIf

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Proc DrawArrow(x1#, y1#, x2#, y2#, Optional size# = 15)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  size = Min(d, size)

  If d > 0

    lx = size * lx / d
    ly = size * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py
    D2Line x2, y2, px + ly , py - lx
    D2Line px + ly , py - lx, px - ly , py + lx
    D2Line px - ly , py + lx , x2, y2

  EndIf

EndProc


Square to disc transposition

Square to disc homogeneous transposition

This procedure give the homogeneous cordinates (in the sense of distributing) of the corresponding point in a disc.

Proc Square2Disc(sx#, sy#, ByRef dx#, ByRef dy#)

  Static Double x, y, h, s, d = Sqr(0.5)

  If sx = 0 And sy = 0
    dx = 0
    dy = 0
    Exit Proc
  EndIf
  h = Sqr(sx * sx + sy * sy)
  If Abs(sx) >= Abs(sy)
    sy = sy / (1 + (Abs(sx) / h - d))
    s = Iif(sx < 0, -1, 1)
    y = sy * s / sx
    x = s
  Else
    sx = sx / (1 + (Abs(sy) / h - d))
    s = Iif(sy < 0, -1, 1)
    x = sx * s / sy
    y = s
  EndIf
  h = Sqr(x * x + y * y)
  dx = sx / h
  dy = sy / h

EndProc

Example

Global Double x = -1, y = -1, dx, dy, p = 0.5
Global Long i

OpenW Center 1

Do

  For i = 1 To 1000

    If Abs(x) < 0.25 And Abs(y) < 0.25
      RGBColor &H20D0F0
    Else If x < 0 And y < 0 And ((x + 0.5) ^ 2 + (y + 0.5) ^ 2) < 0.0625
      RGBColor &hFF
    Else If x > 0 And y < 0
      RGBColor Iif(Frac(x * 4) - Frac(y * 4 ) < 1, &hFFFFFF, &HFF0F)
    Else If (Int(x * 10) + Int(y * 10)) % 2 = 0
      RGBColor RGB(128 + x * 127, 128 + y * 127, 255)
    Else
      RGBColor 0
    EndIf

    Plot _X / 2 - 250 + x * 200, _Y / 2 + y * 200
    Square2Disc x, y, dx, dy
    Plot _X / 2 + 250 + dx * 200, _Y / 2 + dy * 200

    x = x + p
    If x > 1
      x = -1 : y = y + p
      If y > 1  y = -1 : p = p / 2
    EndIf

  Next i

  Sleep

Loop Until Me Is Nothing
Proc Square2Disc(sx#, sy#, ByRef dx#, ByRef dy#)

  Static Double x, y, h, s, d = Sqr(0.5)

  If sx = 0 And sy = 0
    dx = 0
    dy = 0
    Exit Proc
  EndIf
  h = Sqr(sx * sx + sy * sy)
  If Abs(sx) >= Abs(sy)
    sy = sy / (1 + (Abs(sx) / h - d))
    s = Iif(sx < 0, -1, 1)
    y = sy * s / sx
    x = s
  Else
    sx = sx / (1 + (Abs(sy) / h - d))
    s = Iif(sy < 0, -1, 1)
    x = sx * s / sy
    y = s
  EndIf
  h = Sqr(x * x + y * y)
  dx = sx / h
  dy = sy / h

EndProc


Redefining Homogeneous Coordinates

Transform canonic cordinates with segment (x1,y1,x2,y2) for the given radius.

Proc Homonormalise(x1, y1, x2, y2, r, ByRef x As Double, ByRef y As Double)

  Local Double lx, ly

  y = y1 + (y - y1) * (y2 - y1) / r
  If y2 <> y1
    x = x + (y - y1) * (x2 - x1) / (y2 - y1)
  EndIf

EndProc

Example

$Library "Direct2D"
Global Double x, y, x1, y1, x2, y2, r, a
Global i As Long, win1RT As Object
OpenW Center 1
Set win1RT = D2GetRT(Me)
DrawAll
Do
  Sleep
Loop Until Me Is Nothing
Proc Homonormalise(x1, y1, x2, y2, r, ByRef x As Double, ByRef y As Double)

  Local Double lx, ly

  y = y1 + (y - y1) * (y2 - y1) / r
  If y2 <> y1
    x = x + (y - y1) * (x2 - x1) / (y2 - y1)
  EndIf

EndProc
Sub Win_1_Paint

  DrawAll

EndSub
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  DrawAll

EndSub
Proc DrawAll

  D2BeginDraw win1RT, D2C_DarkBlue
  r = _Y / 2.5
  x1 = _X / 2 : y1 = _Y / 2
  x2 = MouseX : y2 = MouseY
  D2ForeColor = RGB(0, 100, 100)
  D2Line x1, y1, x1, y1 + r
  D2Line x1, y1, x2, y2
  For i = 1 To 1000
    a = PI * i / 500
    Select Mod(i, 4)
    Case 0
      x = x1 + Sin(a) * r
      y = y1 + Cos(a) * r
    Case 1, 2
      x = x1 - r / 3 + r / 3 * Mod(i, 2) * 2 + Sin(a) * r * (0.02 + 0.25 * Abs(Mod(i, 30) > 1))
      y = y1 - r / 3 + Cos(a) * r * (0.02 + 0.25 * Abs(Mod(i, 30) > 1))
    Case 3
      x = x1 + Sin(a / 3 - 1) * r / 1.3
      y = y1 + Cos(a / 3 - 1) * r / 1.3
    EndSelect
    D2DefLine , 0
    D2Plot x, y, D2Brush(D2C_Yellow)
    ' r is the fixed radius reference
    Homonormalise x1, y1, x2, y2, r, x, y
    D2PCircle x, y, 3, , , , D2Brush(D2C_Green)
  Next i
  D2EndDraw

EndProc


Rotate Scale Move

Rotate, Scale and Move

This procedure rotate, scale and move the point (x,y).

Proc RotateScaleMove(xCenter, yCenter, xOffset, yOffset, Scaling, Rotate, ByRef x, ByRef y)

  Static Double s, c, newx, newy

  s = Sin(Rotate)
  c = Cos(Rotate)
  x -= xCenter
  y -= yCenter
  x *= Scaling
  y *= Scaling
  newx = c * x + s * y
  newy = c * y - s * x
  x = xCenter + xOffset + newx
  y = yCenter + yOffset + newy

EndProc

Example

$Library "Direct2D"
Global Win1Rt As Object
Global Dim xy() As Double
Array xy() = Mkd( _
  -1, -1, -1.5, -1.5, -2.5, -1.5, -3, -1, -3, 0.5, -2.5, 1, -1.5, 1, -1, 0.5, -1, 0, -1.5, 0, _
  -0.5, 1, -0.5, -1.5, 1, -1.5, -0.5, -0.5, 0.5, -0.5, _
  1.3, 1, 1.3, -0.5, 1.8, -1.5, 2.3, -1.5, 2.8, -0.5, 2.8, 1, 1.3, 0, 2.8, 0)
OpenW Center 1
Ocx Timer FPS
FPS.Interval = 20
FPS.Enabled = True
Set  Win1Rt = D2GetRT(Me)
Do
  Sleep
Loop Until Me Is Nothing
Proc RotateScaleMove(xCenter, yCenter, xOffset, yOffset, Scaling, Rotate, ByRef x, ByRef y)

  Static Double s, c, newx, newy

  s = Sin(Rotate)
  c = Cos(Rotate)
  x -= xCenter
  y -= yCenter
  x *= Scaling
  y *= Scaling
  newx = c * x + s * y
  newy = c * y - s * x
  x = xCenter + xOffset + newx
  y = yCenter + yOffset + newy

EndProc
Sub FPS_Timer
  Static Boolean Work = False
  If !Work
    Work = True
    DrawAll
    Work = False
  EndIf
EndSub
Proc DrawAll

  Local Long i
  Local Double x1, y1, x2, y2
  Local Double Rotation = Sin(Frac(Timer / 7) * 6.283)
  Local Double Scaling = 80 + Sin(Frac(Timer / 9) * 6.283) * 70

  D2BeginDraw Win1Rt, D2C_DarkBlue
  D2DefLine 4096, 5
  For i = 2 To Dim?(xy()) - 2 Step 2
    If i != 20 And i != 26 And i != 30 And i != 42
      x1 = xy(i - 2) , y1 = xy(i - 1), x2 = xy(i) , y2 = xy(i + 1)
      RotateScaleMove 0, 0, MouseX, MouseY, Scaling, Rotation, x1, y1
      RotateScaleMove 0, 0, MouseX, MouseY, Scaling, Rotation, x2, y2
      D2Line x1, y1, x2, y2, D2Brush(D2C_White)
    EndIf
  Next i
  D2EndDraw

EndProc


Converting HSL to RGB

This function return RGB color value from Hue Saturation and Luminosity (H,S,L).
The spectrum correction use the Jean-Marie Melanson's algorithm and the brightness curve is smoothed using a bezier formula.

Function HSL2RGB(H#, S#, L#) As Long

  Local Double r, g, b, m
  Local Double cst1 = 255
  Local Double cst2 = 0.29, cst3 = 0.58, cst4 = 0.13
  Local Double d1 = 1.999999999, d2 = 0.666666666, d3 = 1.333333333, d4 = 0.999999999

  H = Max(0, Min(1, H))
  S = Max(0, Min(1, S))
  L = Max(0, Min(1, L))

  r = Max(0, Min(1, Abs(Frac(d4 - H * d1)) * 3 - 1))
  g = Max(0, Min(1, Abs(Frac(d4 - Frac(H + d2) * d1)) * 3 - 1))
  b = Max(0, Min(1, Abs(Frac(d4 - Frac(H + d3) * d1)) * 3 - 1))

  ' Using Jean-Marie Melanson's algorithm for realistic spectrum :

  m = r * cst2 + g * cst3 + b * cst4

  r = m + S * (r - m)
  g = m + S * (g - m)
  b = m + S * (b - m)

  ' Bezier softing for luminosity :

  L = 2 * L * ((L - 1.5) * L + 1.5)

  If L > 1
    L = L - 1
    r = r + L * (1 - r)
    g = g + L * (1 - g)
    b = b + L * (1 - b)
  Else
    r = r * L
    g = g * L
    b = b * L
  EndIf

  Return RGB(r * cst1, g * cst1, b * cst1)

EndFunc

Example

$Library "Gdip.inc"
Global Bitmap As Long, AdrBitmap As Long, BIH As BITMAPINFOHEADER
BIH.biSize = SizeOf(BIH)
BIH.biWidth = 600
BIH.biHeight = 200
BIH.biPlanes = 1
BIH.biBitCount = 32
OpenW Center 1
Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)
Do
  If _Mess = WM_MOUSEMOVE DrawAll
  Sleep
Loop Until Me Is Nothing
FreeBmp Bitmap
Function HSL2RGB(H#, S#, L#) As Long

  Local Double r, g, b, m
  Local Double cst1 = 255
  Local Double cst2 = 0.29, cst3 = 0.58, cst4 = 0.13
  Local Double d1 = 1.999999999, d2 = 0.666666666, d3 = 1.333333333, d4 = 0.999999999

  H = Max(0, Min(1, H))
  S = Max(0, Min(1, S))
  L = Max(0, Min(1, L))

  r = Max(0, Min(1, Abs(Frac(d4 - H * d1)) * 3 - 1))
  g = Max(0, Min(1, Abs(Frac(d4 - Frac(H + d2) * d1)) * 3 - 1))
  b = Max(0, Min(1, Abs(Frac(d4 - Frac(H + d3) * d1)) * 3 - 1))

  ' Using Jean-Marie Melanson's algorithm for realistic spectrum :

  m = r * cst2 + g * cst3 + b * cst4

  r = m + S * (r - m)
  g = m + S * (g - m)
  b = m + S * (b - m)

  ' Bezier softing for luminosity :

  L = 2 * L * ((L - 1.5) * L + 1.5)

  If L > 1
    L = L - 1
    r = r + L * (1 - r)
    g = g + L * (1 - g)
    b = b + L * (1 - b)
  Else
    r = r * L
    g = g * L
    b = b * L
  EndIf

  Return RGB(r * cst1, g * cst1, b * cst1)

EndFunc
Proc DrawAll

  Local Double Variation
  Local Long Aw, X, Y

  Variation = MouseY / _Y
  Aw = AdrBitmap
  For Y = 0 To 199
    For X = 0 To 199
      LPoke Aw, Bswap3(HSL2RGB(Variation, X / 200, Y / 200))
      LPoke Aw + 800, Bswap3(HSL2RGB(X / 200, Variation, Y / 200))
      LPoke Aw + 1600, Bswap3(HSL2RGB(X / 200, Y / 200, Variation))
      Add Aw, 4
    Next X
    Add Aw, 1600
  Next Y

  X = _X / 2 - 300 : Y = _Y / 2 - 100
  Put X, Y, Bitmap
  Y = Y - TextHeight("?")
  Text X, Y, "HUE"
  Text X + 200, Y, "SATURATION"
  Text X + 400, Y, "LUMINOSITY"

EndProc


Computing ray reflection

Get incident vector (reflection)

This function compute the incident vector resulting of the incoming vector with a normal vector surface.

Function GetXYZReflect(ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double dp

  XYZ_Normalise Income
  XYZ_Normalise Normale

  dp = 2 * (Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)

  GetXYZReflect = XYZ(Income.x - Normale.x * dp, Income.y - Normale.y * dp, Income.z - Normale.z * dp)

EndFunc

Example

$Library "Direct2D"

Type XYZ_STRUCT
  - Double x, y, z
EndType
Global Dim XY() As D2D1_POINT_2F
Array XY() = Mks(-3, 0, -2, -1, -1.5, -2, 0.5, -1.5, 1.5, -1.5, 1.5, 0.5, 0.5, 1.5, 0, 3, -1, 1, -3, 0)
Global Long i, NCount = Dim?(XY()) - 1
Global Dim N(NCount) As XYZ_STRUCT
Global Object Geo, GeoSink, Win1RT
Global Normal As XYZ_STRUCT, Incoming As XYZ_STRUCT, Incident As XYZ_STRUCT
Global Double cx, cy

OpenW Center 1

cx = _X / 2 : cy = _Y / 2
For i = 0 To NCount - 1
  N(i) = XYZ(XY(i).y - XY(i + 1).y, XY(i + 1).x - XY(i).x, 0)
  XYZ_Normalise N(i)
  XY(i).x = cx + XY(i).x * 100
  XY(i).y = cy + XY(i).y * 100
Next i

Set Win1RT = D2GetRT(Me)
Set Geo = D2PathGeometry()
Set GeoSink = D2GeometrySink(Geo)
D2Sink_FillMode GeoSink, D2D1_FILL_MODE_WINDING
D2Sink_BeginFigure GeoSink, XY(NCount - 1).x, XY(NCount - 1).y
D2Sink_AddLines GeoSink, XY(), NCount - 1
D2Sink_EndFigure GeoSink, 1
D2Sink_Close GeoSink

Do
  DrawAll
  Sleep
Loop Until Me Is Nothing
Function GetXYZReflect(ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double dp

  XYZ_Normalise Income
  XYZ_Normalise Normale

  dp = 2 * (Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)

  GetXYZReflect = XYZ(Income.x - Normale.x * dp, Income.y - Normale.y * dp, Income.z - Normale.z * dp)

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Proc DrawArrow(x1#, y1#, x2#, y2#, C%)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  If d > 0

    lx = 15 * lx / d
    ly = 15 * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py, D2Brush(C)
    D2Line x2, y2, px + ly , py - lx, D2Brush(C)
    D2Line px + ly , py - lx ,  px - ly , py + lx, D2Brush(C)
    D2Line px - ly , py + lx , x2, y2, D2Brush(C)

  EndIf

EndProc
Proc DrawAll

  Local Double ix, iy, x1, y1, x2, y2, mx, my, Radius = 100

  D2BeginDraw Win1RT, D2C_DarkBlue

  D2FillGeometry Geo, D2Brush(ARGB(255, 100, 100, 100))

  For i = 0 To NCount - 1
    x1 = XY(i).x, y1 = XY(i).y, x2 = XY((i + 1) % NCount).x, y2 = XY((i + 1) % NCount).y
    mx = (x1 + x2) / 2
    my = (y1 + y2) / 2
    Incoming = XYZ(MouseX - cx, MouseY - cy, 0)
    DrawArrow cx, cy, cx + Incoming.x, cy + Incoming.y, D2C_Yellow
    If IntersectionSegments(MouseX, MouseY, cx, cy, x1, y1, x2, y2, ix, iy)
      Incident = GetXYZReflect(Incoming, N(i))
      D2DefLine , 2
      DrawArrow ix, iy, ix + Incident.x * Radius, iy + Incident.y * Radius, ARGB(255, 100, 255, 100)
      DrawArrow ix, iy, ix + N(i).x * Radius, iy + N(i).y * Radius, D2C_White
    Else
      D2DefLine , 1
      DrawArrow mx, my, mx + N(i).x * Radius / 2, my + N(i).y * Radius / 2, ARGB(255, 200, 200, 200)
    EndIf
  Next i

  D2EndDraw

EndProc
Function IntersectionSegments(x1#, y1#, x2#, y2#, x3#, y3#, x4#, y4#, ByRef ix#, ByRef iy#) As Boolean

  Local Double Lx2, Lx1, Ly, Lx, Ly1, Ly2, M, Ua, Ub, M1, M2

  Lx2 = x4 - x3
  Lx1 = x2 - x1
  Ly1 = y2 - y1
  Ly2 = y4 - y3

  M1 = Lx2 * Ly1
  M2 = Ly2 * Lx1

  If M1 <> M2  ' If the lines are not parallel

    Ly = y1 - y3
    Lx = x1 - x3
    M = M2 - M1

    Ua = ( Lx2 * Ly - Ly2 * Lx ) / M
    Ub = ( Lx1 * Ly - Ly1 * Lx ) / M

    If Ua > 0
      If Ua < 1
        If Ub > 0
          If Ub < 1
            ix = x1 + Ua * Lx1
            iy = y1 + Ua * Ly1
            Return True
          EndIf
        EndIf
      EndIf
    EndIf

  EndIf

  Return False

EndFunc


Compute ray refraction

This function compute the refraction with the Income ray vector and the Normale vector.
n1 is the exterior refractive index (that mean it's where Income ray), n2 is the interior refractive index.
Refractive index = 3e10 (light speed in cosmos) / light speed of the material.

Function GetXYZRefract(n1#, n2#, ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double c1, n, c2, cc
  Static Outgoing As XYZ_STRUCT

  XYZ_Normalise Income

  n1 = Max(1, n1)
  n2 = Max(1, n2)

  c1 = -(Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)
  n = n1 / n2
  c2 = Sqr(Max(0, 1 - (n * n) * (1 - (c1 * c1))))
  cc = n * c1 - c2
  Outgoing.x = (n * Income.x) + cc * Normale.x
  Outgoing.y = (n * Income.y) + cc * Normale.y
  Outgoing.z = (n * Income.z) + cc * Normale.z
  XYZ_Normalise Outgoing

  Return Outgoing

EndFunc

Example

$Library "Gdip.inc"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Global Const REFRACTION_AIR = 1.000272
Global Const REFRACTION_WATER = 1.325
Global Const REFRACTION_CRISTAL = 1.545
Global Const REFRACTION_DIAMOND = 2.4175
Global Bitmap As Long, AdrBitmap As Long, BIH As BITMAPINFOHEADER
BIH.biSize = SizeOf(BIH)
BIH.biWidth = 600
BIH.biHeight = 200
BIH.biPlanes = 1
BIH.biBitCount = 32
Global Dim Normal(0 .. 199, 0 .. 199) As XYZ_STRUCT
Global Long X, Y, i
Global Double dx, dy, d
For Y = 0 To 199
  For X = 0 To 199
    dy = (Y - 100)
    dx = (X - 100)
    d = Sqr(dx * dx + dy * dy)
    If d >= 80
      Normal(X, Y) =  XYZ(0, 0, -1)
    Else
      d = (80 - d) / 80
      Normal(X, Y) = XYZ(dx / 80, dy / 80, -1 + Sqr(1 - d * d))
      XYZ_Normalise Normal(X, Y)
    EndIf
  Next X
Next Y
Global Incoming As XYZ_STRUCT
Global Refract As XYZ_STRUCT

OpenW Center 1
Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)
Do
  If _Mess = WM_MOUSEMOVE DrawAll
  Sleep
Loop Until Me Is Nothing
FreeBmp Bitmap
Function GetXYZRefract(n1#, n2#, ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double c1, n, c2, cc
  Static Outgoing As XYZ_STRUCT

  XYZ_Normalise Income

  n1 = Max(1, n1)
  n2 = Max(1, n2)

  c1 = -(Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)
  n = n1 / n2
  c2 = Sqr(Max(0, 1 - (n * n) * (1 - (c1 * c1))))
  cc = n * c1 - c2
  Outgoing.x = (n * Income.x) + cc * Normale.x
  Outgoing.y = (n * Income.y) + cc * Normale.y
  Outgoing.z = (n * Income.z) + cc * Normale.z
  XYZ_Normalise Outgoing

  Return Outgoing

EndFunc
Proc DrawAll

  Local Long X, Y, C, Aw
  Local Double ix, iy, ox, oy
  Static Boolean Work = False

  If !Work
    Work = True
    ox = -MouseX
    oy = MouseY
    Aw = AdrBitmap
    For Y = 0 To 199
      For X = 0 To 199

        Refract = GetXYZRefract(REFRACTION_AIR, REFRACTION_WATER, Incoming, Normal(X, Y))
        C = CInt(Abs((Abs(Normal(X, Y).z) ^ 0.2) * 254))
        C = RGB(C, C, Iif(Normal(X, Y).z = -1, 255, 0))
        If Refract.z > 0.001
          ix = X + 80 * Refract.x / Refract.z
          iy = Y + 80 * Refract.y / Refract.z
        Else
          ix = X
          iy = Y
        EndIf
        LPoke Aw, Damier((ox + ix) / 50, (oy + iy) / 50) * C

        Refract = GetXYZRefract(REFRACTION_AIR, REFRACTION_DIAMOND , Incoming, Normal(X, Y))
        C = CInt(Abs((Abs(Normal(X, Y).z) ^ 0.2) * 254.0))
        C = RGB(Iif(Normal(X, Y).z = -1, 255, 0), C, C)
        If Refract.z > 0.001
          ix = X + 200 + 80 * Refract.x / Refract.z
          iy = Y + 80 * Refract.y / Refract.z
        Else
          ix = X + 200
          iy = Y
        EndIf
        LPoke Aw + 800, Damier((ox + ix) / 50, (oy + iy) / 50) * C

        Refract = GetXYZRefract(REFRACTION_WATER, REFRACTION_AIR, Incoming, Normal(X, Y))
        C = CInt(Abs((Abs(Normal(X, Y).z) ^ 2) * 254.0))
        C = RGB(C, Iif(Normal(X, Y).z = -1, 128 + C / 2, 255), C / 2 )
        If Refract.z > 0.001
          ix = X + 400 + 80 * Refract.x / Refract.z
          iy = Y + 80 * Refract.y / Refract.z
        Else
          ix = X + 400
          iy = Y
        EndIf
        LPoke Aw + 1600, Damier((ox + ix) / 50, (oy + iy) / 50) * C

        Add Aw, 4

      Next X
      Add Aw, 1600
    Next Y

    X = _X / 2 - 300 : Y = _Y / 2 - 100
    Put X, Y, Bitmap

    Y = Y - TextHeight("?")
    Text X, Y, "AIR/WATER"
    Text X + 200, Y, "AIR/DIAMOND"
    Text X + 400, Y, "WATER/AIR"

    Work = False
  EndIf

EndProc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Function Damier(ix#, iz#) As Long

  Local Long alpha = 1
  If Frac(Abs(ix)) < 0.5 alpha = 1 - alpha
  If Frac(Abs(iz)) < 0.5 alpha = 1 - alpha
  If iz < 0 alpha = 1 - alpha
  If ix < 0 alpha = 1 - alpha
  Return alpha

EndFunc


Deviation from normal

This function return the vector deviation arround the Normale (on the plane of it) rotating by G turns.

Function GetNormaleDeviation(ByRef Normale As XYZ_STRUCT, G As Double) As XYZ_STRUCT

  Local Double ax, ay, x, y, z, s, c, d, newx, newz
  Local Double p2 = 6.283185307179586476925286766559
  Local Double p = 3.1415926535897932384626433832795

  G = G * p2

  x = Sin(G) : y = 0 : z = Cos(G)

  d = Sqr(Normale.x * Normale.x + Normale.z * Normale.z)
  If Normale.y <> 0
    ax = Atan(d / Normale.y)
    If Normale.y < 0 ax = ax + p
    s = Sin(ax) : c = Cos(ax)
  Else
    s = 1 : c = 0
  EndIf

  newz = c * z + s * y
  y = c * y - s * z
  z = newz

  If Normale.z <> 0
    ay = Atan(Normale.x / Normale.z)
    If Normale.z < 0 ay = ay + p
    s = Sin(ay) : c = Cos(ay)
  Else If Normale.x < 0
    s = -1 : c = 0
  Else
    s = 1 : c = 0
  EndIf

  newx = c * x + s * z
  z = c * z - s * x
  x = newx

  GetNormaleDeviation.x = x
  GetNormaleDeviation.y = y
  GetNormaleDeviation.z = z

EndFunc

Example

$Library "Direct2D"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Global Object Win1RT
Global Normal As XYZ_STRUCT, Deviation As XYZ_STRUCT
Global Double G

OpenW Center 1
Set Win1RT = D2GetRT(Me)
Ocx Timer FPS
FPS.Interval = 20
FPS.Enabled = True
Do
  Sleep
Loop Until Me Is Nothing
Function GetNormaleDeviation(ByRef Normale As XYZ_STRUCT, G As Double) As XYZ_STRUCT

  Local Double ax, ay, x, y, z, s, c, d, newx, newz
  Local Double p2 = 6.283185307179586476925286766559
  Local Double p = 3.1415926535897932384626433832795

  G = G * p2

  x = Sin(G) : y = 0 : z = Cos(G)

  d = Sqr(Normale.x * Normale.x + Normale.z * Normale.z)
  If Normale.y <> 0
    ax = Atan(d / Normale.y)
    If Normale.y < 0 ax = ax + p
    s = Sin(ax) : c = Cos(ax)
  Else
    s = 1 : c = 0
  EndIf

  newz = c * z + s * y
  y = c * y - s * z
  z = newz

  If Normale.z <> 0
    ay = Atan(Normale.x / Normale.z)
    If Normale.z < 0 ay = ay + p
    s = Sin(ay) : c = Cos(ay)
  Else If Normale.x < 0
    s = -1 : c = 0
  Else
    s = 1 : c = 0
  EndIf

  newx = c * x + s * z
  z = c * z - s * x
  x = newx

  GetNormaleDeviation.x = x
  GetNormaleDeviation.y = y
  GetNormaleDeviation.z = z

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Sub FPS_Timer
  DrawAll
EndSub
Proc DrawAll

  Static Boolean Work = False
  Local Double a, cx, cy, cz, x, y, x1, y1, x2, y2, radius
  Local Long i

  If !Work

    Work = True

    G = Frac(Timer / 10)

    cx = _X / 2 : cy = _Y / 2
    cz = _Y
    radius = cz / 1.5
    D2BeginDraw Win1RT, D2C_DarkBlue

    a = -1.5 + 3 * (MouseX / _X)
    Normal = XYZ(Sin(a), 10 * ((MouseY / _Y) - 0.5), -Cos(a))
    XYZ_Normalise Normal

    D2DefLine , 0

    For i = 1 To 500
      Deviation = GetNormaleDeviation(Normal, i / 500)
      Pseudo3D cx + Deviation.x * radius, cy + Deviation.y * radius, cz + Deviation.z * radius, x, y
      D2PCircle x, y, 10 - (1 + Deviation.z) * 4 , , , , D2Brush(ARGB(20, 10, 255, 100))
    Next i

    D2DefLine , 1
    Pseudo3D cx, cy, cz, x1, y1
    Pseudo3D cx + Normal.x * radius, cy + Normal.y * radius, cz + Normal.z * radius, x2, y2
    DrawArrow x1, y1, x2, y2, D2C_White

    D2DefLine 4096, 3
    Deviation = GetNormaleDeviation(Normal, G)
    Pseudo3D cx + Deviation.x * radius, cy + Deviation.y * radius, cz + Deviation.z * radius, x2, y2
    DrawArrow x1, y1, x2, y2, ARGB(255, 255, 255, 0)

    D2SetFont "Lucida Console", 40
    D2DrawText 0, 0, _X, 100, "G = " & Left(Trim(Round(G, 3)) & "     ", 6) & " tr.", D2T_CENTER | D2T_PARCENTER, , D2Brush(D2C_White)

    D2EndDraw

    Work = False
  EndIf

EndProc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Proc DrawArrow(x1#, y1#, x2#, y2#, C%)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  If d > 0

    lx = 15 * lx / d
    ly = 15 * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py, D2Brush(C)
    D2Line x2, y2, px + ly , py - lx, D2Brush(C)
    D2Line px + ly , py - lx ,  px - ly , py + lx, D2Brush(C)
    D2Line px - ly , py + lx , x2, y2, D2Brush(C)

  EndIf

EndProc
Proc  Pseudo3D(x#, y#, z#, ByRef x3d#, ByRef y3d#)

  Static Double dx, dy, d

  d = 1 / (1 + Abs(z) / 500)
  dx = x - _X / 2
  dy = y - _Y / 2
  x3d = _X / 2 + dx * d
  y3d = _Y / 2 + dy * d

EndProc


Mandelbrot Fractal

This function return the total iterations before exiting the Mandelbrot fractal figure from x0,y0.
Because this total can never be reached we must indicate Itermax (the maximum counting iterations) before considerate it's infinite (inside the fractal surface).
If this maximum is not indicated the default value is 150.

Function Mandelbrot(x0 As Double, y0 As Double, Optional Itermax%=150) As Long Naked

  Static Double x, y, xnew, ynew
  Static Long Compteur

  x = 0
  y = 0
  xnew = 0
  ynew = 0
  For Compteur = 0 To Itermax
    ynew = 2 * xnew * ynew  + y0
    xnew = x - y + x0
    x = xnew * xnew
    y = ynew * ynew
    If x + y > 4 Return Compteur
  Next Compteur

  Return Itermax

EndFunc

Example

$Library "Gdip.inc"
Global Bitmap As Long, AdrBitmap As Long, BIH As BITMAPINFOHEADER
BIH.biSize = SizeOf(BIH)
BIH.biWidth = 600
BIH.biHeight = 600
BIH.biPlanes = 1
BIH.biBitCount = 32
OpenW Center 1
BackColor = colBtnFace
Ocx Command C_Reset = "Reset", 10, 10, 100, 30
Ocx Command C_Itermax = "Itermax + 100", 10, 50, 100, 30
Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)
Global Double t, x = -2, y = -2, Fx = 0, Fy = 0, P = 1 / 150, Zoom = 2
Global Long i, Tx = 0, Ty = 0, M, IterMax = 150, Aw = AdrBitmap, C

Do

  FractalTracing
  Put _X / 2 - 300, _Y / 2 - 300, Bitmap

  If Abs(MouseX - _X / 2) < 300 And Abs(MouseY - _Y / 2) < 300
    Clip _X / 2 - 300, _Y / 2 - 300, 600, 600
    RGBColor 0
    Box MouseX - 150, MouseY - 150, MouseX + 150, MouseY + 150
    RGBColor &hffffff
    Box MouseX - 151, MouseY - 151, MouseX + 151, MouseY + 151
    Clip Off
  EndIf

  Sleep
Loop Until Me Is Nothing
FreeBmp Bitmap
Sub C_Reset_Click

  x = -2, y = -2, Fx = 0, Fy = 0, P = 1 / 150, Zoom = 2
  Tx = 0, Ty = 0, IterMax = 150, Aw = AdrBitmap
  Me.Caption = "Zoom : " & 2 / Zoom

EndSub
Sub C_Itermax_Click

  Add IterMax, 100

EndSub
Sub Win_1_MouseUp(Button&, Shift&, x!, y!)

  If Abs(MouseX - _X / 2) < 300 And Abs(MouseY - _Y / 2) < 300
    Fx = Fx + (MouseX - _X / 2 ) / 300 * Zoom
    Fy = Fy + (_Y / 2 - MouseY) / 300 * Zoom
    Zoom = Zoom / 2
    x = Fx - Zoom
    y = Fy - Zoom
    P = P / 2
    Tx = 0
    Ty = 0
    Aw = AdrBitmap
    Me.Caption = "Zoom : " & 2 / Zoom
  EndIf

EndSub
Proc FractalTracing

  t = Timer
  Do
    M = Mandelbrot(x, y, IterMax)
    If M < IterMax
      C = RGB((M % 4) * 63, ( M % 5) * 51, (M % 7) * 36)
    Else
      C = 0
    EndIf
    LPoke Aw, C
    Add Aw, 4
    Inc Tx
    x += P
    If Tx = 600
      Tx = 0
      x = Fx - Zoom
      Inc Ty
      y += P
      If Ty = 600
        y = Fy - Zoom
        Ty = 0
        Aw = AdrBitmap
      EndIf
    EndIf
  Loop Until Timer - t > 0.2

EndProc
Function Mandelbrot(x0 As Double, y0 As Double, Optional Itermax%=150) As Long Naked

  Static Double x, y, xnew, ynew
  Static Long Compteur

  x = 0
  y = 0
  xnew = 0
  ynew = 0
  For Compteur = 0 To Itermax
    ynew = 2 * xnew * ynew  + y0
    xnew = x - y + x0
    x = xnew * xnew
    y = ynew * ynew
    If x + y > 4 Return Compteur
  Next Compteur

  Return Itermax

EndFunc


Random point on a sphere

This function returns a random point so the statistical distribution is homogeneous.

Function SphericalRandom() As XYZ_STRUCT

  Static Double h, shh, o, p = 3.1415926535897932384626433832795

  h = -1 + 2 * Rnd
  o = Rnd * p * 2
  shh = Sqr(1 - h * h)
  SphericalRandom.x = shh * Cos(o)
  SphericalRandom.y = shh * Sin(o)
  SphericalRandom.z = h

EndFunc

Example

Type XYZ_STRUCT
  - Double x, y, z
EndType
Global S As XYZ_STRUCT
Global Double Radius
Global Long i
OpenW Center 1
Radius = _X / 8
Do
  Text _X / 5, _Y / 2 - Radius - 20, "Front"
  Text _X / 2, _Y / 2 - Radius - 20, "Top"
  Text _X * 0.8, _Y / 2 - Radius - 20, "Left"
  For i = 1 To 10000
    S = SphericalRandom()
    Plot _X / 5   + S.x * Radius, _Y / 2 + S.y * Radius
    Plot _X / 2   + S.x * Radius, _Y / 2 + S.z * Radius
    Plot _X * 0.8 + S.y * Radius, _Y / 2 + S.z * Radius
  Next i
  Sleep
Loop Until Me Is Nothing
Function SphericalRandom() As XYZ_STRUCT

  Static Double h, shh, o, p = 3.1415926535897932384626433832795

  h = -1 + 2 * Rnd
  o = Rnd * p * 2
  shh = Sqr(1 - h * h)
  SphericalRandom.x = shh * Cos(o)
  SphericalRandom.y = shh * Sin(o)
  SphericalRandom.z = h

EndFunc


Plane Equation By Triangle

Obtain a plane equation using triangle

This function compute the plane equation a,b,c,d from the given triangle T.
This method for obtaining the equation of a plan has a small drawback: on values that are too large, you can be confronted with an overflow during multiplications.
If no overflow error occur the function return True.

Function PlaneEquation(x1#, y1#, z1#, x2#, y2#, z2#, x3#, y3#, z3#, ByRef a#, ByRef b#, ByRef c#, ByRef d#) As Boolean

  Try

    a = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
    b = z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2)
    c = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    d = x1 * (y2 * z3 - y3 * z2) + x2 * (y3 * z1 - y1 * z3) + x3 * (y1 * z2 - y2 * z1)
    d = -d

    Return True

  Catch

    Return False

  EndCatch

EndFunc

Example

$Library "Direct2D"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Type TRIANGLE3D_STRUCT
  - XYZ_STRUCT a, b, c
EndType
Global T As TRIANGLE3D_STRUCT
Global Double a, b, c, d, cx, cy, cz, RotationView
Global Object Win1RT
OpenW Center 1
Set Win1RT = D2GetRT(Win_1)
' Define a 3D center to draw all
cx = _X / 2 : cy = _Y / 2 : cz = _Y
Ocx Timer FPS : .Interval = 20 : .Enabled = True
Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer

  DrawAll

EndSub
Proc DrawGround

  Local Long i
  Local Double r, u

  r = cz / 1.5
  For i = -10 To 10
    u = i / 10
    DrawLine3D(XYZ(cx + u * r, cy + 400, cz + r), XYZ(cx + u * r, cy + 400, cz - r), D2C_Orange)
    DrawLine3D(XYZ(cx + r, cy + 400, cz + u * r), XYZ(cx - r, cy + 400, cz + u * r), D2C_Orange)
  Next i

EndProc
Proc DrawPlane

  Local Long i, j
  Local P(0 .. 3) As XYZ_STRUCT
  Local Double r

  r = cz / 1.5

  For j = 0 To 1

    P(0) = XYZ(cx - r, 0, cz - r)
    P(1) = XYZ(cx + r, 0, cz - r)
    P(2) = XYZ(cx + r, 0, cz + r)
    P(3) = XYZ(cx - r, 0, cz + r)

    ' Compute all Y with plane equation applied to X,Z
    For i = 0 To 3
      P(i).y = (-a * P(i).x - c * P(i).z - d) / b
    Next i

    DrawLine3D(P(0), P(1), D2C_Yellow)
    DrawLine3D(P(1), P(2), D2C_Yellow)
    DrawLine3D(P(2), P(3), D2C_Yellow)
    DrawLine3D(P(3), P(0), D2C_Yellow)

    r /= 1.05

  Next j

EndProc
Proc DrawAll

  Local Double ua, ub

  Static Boolean Work = False

  If !Work

    Work = True

    RotationView = Frac(Timer / 30) * 6.283

    ' Define an arbitrary triangle T in space
    ua = (MouseX / _X) * 6.283
    ub = -1.5 + 3 * (MouseY / _Y)
    T.a = XYZ(cx + Sin(ua) * 250, cy + Sin(ub) * 100, cz + Cos(ua) * 250)
    T.b = XYZ(cx + Sin(ua + 2) * 250 , cy + Sin(ub + 1) * 200, cz + Cos(ua + 2) * 250)
    T.c = XYZ(cx + Sin(ua + 4) * 250, cy + Sin(ub + 2) * 20, cz + Cos(ua + 4) * 250)

    ' Obtain the plane equation ax+bx+cx+d=0 where is the triangle T
    PlaneEquation(T.a.x, T.a.y, T.a.z, T.b.x, T.b.y, T.b.z, T.c.x, T.c.y, T.c.z, a, b, c, d)

    D2BeginDraw Win1RT, D2C_DarkBlue

    ' A ground like a grid 20x20
    DrawGround

    ' The triangle
    D2DefLine 4096, 5
    DrawLine3D(T.a, T.b, D2C_White)
    DrawLine3D(T.b, T.c, D2C_White)
    DrawLine3D(T.c, T.a, D2C_White)

    ' The plane delimited by ground limit
    DrawPlane

    D2EndDraw
    Work = False
  EndIf

EndProc
Function DrawLine3D(ByRef A As XYZ_STRUCT, ByRef B As XYZ_STRUCT, C As Long)

  Local Double x1, y1, x2, y2

  Pseudo3D A.x, A.y, A.z, x1, y1
  Pseudo3D B.x, B.y, B.z, x2, y2

  D2Line x1, y1, x2, y2, D2Brush(C)

EndFunc
Function PlaneEquation(x1#, y1#, z1#, x2#, y2#, z2#, x3#, y3#, z3#, ByRef a#, ByRef b#, ByRef c#, ByRef d#) As Boolean

  Try

    a = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
    b = z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2)
    c = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    d = x1 * (y2 * z3 - y3 * z2) + x2 * (y3 * z1 - y1 * z3) + x3 * (y1 * z2 - y2 * z1)
    d = -d

    Return True

  Catch

    Return False

  EndCatch

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Proc  Pseudo3D(x#, y#, z#, ByRef x3d#, ByRef y3d#)

  Static Double dx, dy, d

  x -= cx
  z -= cz
  Turn x, z, RotationView
  x += cx
  z += cz

  d = 1 / (1 + Abs(z) / 500)
  dx = x - _X / 2
  dy = y - _Y / 2
  x3d = _X / 2 + dx * d
  y3d = _Y / 2 + dy * d

EndProc
Proc Turn(ByRef x#, ByRef y#, a#)

  Static Double s, c, newx, newy

  s = Sin(a)
  c = Cos(a)
  newx = c * x + s * y
  newy = c * y - s * x
  x = newx
  y = newy

EndProc


LCM and GCD

LCM (PPCM) and GCD (PGCD)

These 2 functions PGCD() and PPCM() calculate the greatest common divisor (GCD) and the least common multiple (LCM) of 2 integers.

Function PPCM(A As Long, B As Long) As Long

  Local Long AM, BM

  If A = 0 Return 0
  If B = 0 Return 0

  A = Abs(A)
  B = Abs(B)
  AM = A
  BM = B

  Do
    If AM = BM  Return AM
    If AM < BM
      Add AM, A
    Else
      Add BM, B
    EndIf
  Loop

  Return AM

EndFunc
Function PGCD(A As Long, B As Long) As Long

  Local Long T

  A = Abs(A)
  B = Abs(B)

  While B <> 0
    T = B
    B = A % B
    A = T
  Wend

  Return A

EndFunc

Example

OpenW Center 1
AutoRedraw = True
BackColor = colBtnFace
Mode StrSpace 0
Global Long A, B, C, D, i
Global Double w
Global String T
w = _X - 50
A = 15
B = 12
C = PPCM(A, B)
D = PGCD(A, B)
FontSize = 20
RGBColor 0
For i = 1 To A
  RGBColor Mod(i, 2) * &H404040
  PBox 25 + (i - 1) * w / A, 25, 25 + i * w / A - 2, 75
Next i
For i = 1 To B
  RGBColor Mod(i, 2) * &H404040
  PBox 25 + (i - 1) * w / B, 80, 25 + i * w / B - 2, 130
Next i
RGBColor 0
PBox 25, 135, 25 + w - 2 , 185
For i = 1 To C
  RGBColor RGB(255, 0, 0)
  PCircle 25 + (i - 0.5) * w / C, 50, 5
  RGBColor RGB(255, 200, 0)
  PCircle 25 + (i - 0.5) * w / C, 105, 5
  RGBColor RGB(100, 255, 100)
  PCircle 25 + (i - 0.5) * w / C, 160, 5
Next i
RGBColor 0
PCircle _X / 3, 300, 70
PCircle _X * 0.66, 300, 70
RGBColor &hFF
For i = 1 To A
  PCircle _X / 3 + SinQ((i + 0.5) * 360 / A) * 50, 300 + CosQ((i + 0.5) * 360 / A) * 50, 5
Next i
RGBColor RGB(255, 200, 0)
For i = 1 To B
  PCircle _X * 0.66 + SinQ((i + 0.5) * 360 / B) * 50, 300 + CosQ((i + 0.5) * 360 / B) * 50, 5
Next i
RGBColor RGB(100, 255, 100)
DefLine , 3
For i = 1 To D
  Line _X / 3, 300, _X / 3 + SinQ(i * 360 / D) * 70, 300 + CosQ(i * 360 / D) * 70
  Line _X * 0.66, 300, _X * 0.66 + SinQ(i * 360 / D) * 70, 300 + CosQ(i * 360 / D) * 70
Next i
RGBColor 0
T = "LCM of " & A & " and " & B & " = " & C & " ; " &  C \ A & " * " & A & " = " &  C \ B & " * " & B & " = " & C
Text _X / 2 - TextWidth(T) / 2, 400, T
T = "GCD of " & A & " and " & B & " = " & D
Text _X / 2 - TextWidth(T) / 2, 440, T
Do
  Sleep
Loop Until Me Is Nothing
Function PPCM(A As Long, B As Long) As Long

  Local Long AM, BM

  If A = 0 Return 0
  If B = 0 Return 0

  A = Abs(A)
  B = Abs(B)
  AM = A
  BM = B

  Do
    If AM = BM  Return AM
    If AM < BM
      Add AM, A
    Else
      Add BM, B
    EndIf
  Loop

  Return AM

EndFunc
Function PGCD(A As Long, B As Long) As Long

  Local Long T

  A = Abs(A)
  B = Abs(B)

  While B <> 0
    T = B
    B = A % B
    A = T
  Wend

  Return A

EndFunc


Angle Segment

Compute angle of a line segment with Asin()

This function return the angle in radian of the segment x1,y1,x2,y2. Considering vertical top-down segment (0,0,0,1) is angle 0. The returned value >=0 and <=2*pi.

Function GetSegmentAngle(x1#, y1#, x2#, y2#) As Double

  Local Double lx, a, h
  Local Double p = 3.1415926535897932384626433832795
  Local Double p2 = 6.283185307179586476925286766559

  lx = x2 - x1
  h = Sqr(lx * lx + (y2 - y1) ^ 2)
  If h > 0
    lx = lx / h
    If y2 > y1
      a = Asin(lx)
    Else
      a = p - Asin(lx)
    EndIf
  EndIf

  If a < 0
    a = a + p2
  EndIf

  Return a

EndFunc

Example

Global Double x1, y1, x2, y2, Angle
OpenW Center 1
Do
  Sleep
Loop Until Me Is Nothing
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  DrawAll

EndSub
Proc DrawAll

  Cls
  x1 = _X / 2 : y1 = _Y / 2 : x2 = MouseX : y2 = MouseY
  Angle = GetSegmentAngle(x1, y1, x2, y2)
  DefLine , 0
  RGBColor RGB(0, 100, 0)
  PCircle x1, y1, 200, -90, -90 + Angle * 57.29
  RGBColor RGB(255, 0, 0)
  DefLine , 10

  Line x1, y1, x2, y2

EndProc
Function GetSegmentAngle(x1#, y1#, x2#, y2#) As Double

  Local Double lx, a, h
  Local Double p = 3.1415926535897932384626433832795
  Local Double p2 = 6.283185307179586476925286766559

  lx = x2 - x1
  h = Sqr(lx * lx + (y2 - y1) ^ 2)
  If h > 0
    lx = lx / h
    If y2 > y1
      a = Asin(lx)
    Else
      a = p - Asin(lx)
    EndIf
  EndIf

  If a < 0
    a = a + p2
  EndIf

  Return a

EndFunc


Interpolation of Z with X and Y

Interpolation of Z with X,Y

This function return the z4 value of a 3D triangle (x1,y1,z1)(x2,y2,z2)(x3,y3,z3) for a given (x4,y4) value.

Function InterpolationZfromXY(x1#, y1#, z1#, x2#, y2#, z2#, x3#, y3#, z3#, x4#, y4#) As Double

  Local Double xn, yn, zn, z4

  xn = (y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1)
  yn = (x3 - x1) * (z2 - z1) - (x2 - x1) * (z3 - z1)
  zn = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)

  z4 = (x1 * xn + y1 * yn + z1 * zn - x4 * xn - y4 * yn) / zn

  Return z4

EndFunc

Example

$Library "OpenGL.inc"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Global Long x, z
Global Player As XYZ_STRUCT
Global PlayerView As XYZ_STRUCT
Global Double yGround

' This example use OpenGl, it's much more convenient to render terrain
InitialisationOpenGL

' Define a random ground to test interpolation Y from XZ
' The Y ground is computed in DrawGround Procedure
' By call to the function InterpolationZfromXY with permutation between Z and Y
Global Dim Ground(0 .. 19, 0 .. 19) As Double
For z = 0 To 19
  For x = 0 To 19
    Ground(x, z) = Max(0, 1 - (Sqr((x - 10) ^ 2 + (z - 10) ^ 2) / 10))
    Ground(x, z) += Rnd * Ground(x, z) * 2
  Next x
Next z

' A first person evolution need a minimum player structure
Player = XYZ(1.4, 1, 1.5)
PlayerView = XYZ(0.7, 0, 0.7)

' Create a marvelous texture to render the nicest rendering possible.
Global Dim RGBADATA(0 .. 255, 0 .. 255) As Long
For z = 0 To 255
  For x = 0 To 255
    RGBADATA(x, z) =  &HFF000000 | ((Int(x / 64) + Int(z / 64))  % 2) * &HFFF000
    RGBADATA(x, z) |=  Sqr((x - 128) ^ 2 + (z - 128) ^ 2) * 1.4
    RGBADATA(x, z) |= RGB(Rand(100), Rand(100), Rand(100))
  Next x
Next z
Global Long Texture = RGBAtoTexture(V:RGBADATA(0, 0), 256, 256)
Erase RGBADATA()

glClearColor(0, 0, 0, 0)
glEnable(GL_DEPTH_TEST)

Do

  glViewport(0, 0, Max(1, _X), Max(1, _Y))
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()
  gluPerspective(45, Max(1, _X) / Max(1, _Y), 0.1, 1000)

  PlayerView.y = -((MouseY / _Y) - 0.5) / 2
  Player.y = Player.y + (yGround - Player.y) / 20
  gluLookAt(Player.x, Player.y + 0.5, Player.z, Player.x + PlayerView.x, Player.y + PlayerView.y + 0.3, Player.z + PlayerView.z, 0, 1, 0)

  ' Moving with keyboard Interaction
  If Btst(GetAsyncKeyState(VK_RIGHT), 15) = True Turn PlayerView.x, PlayerView.z, -0.005
  If Btst(GetAsyncKeyState(VK_LEFT), 15) = True Turn PlayerView.x, PlayerView.z, 0.005
  If Btst(GetAsyncKeyState(VK_UP), 15) = True
    Player.x = Player.x + PlayerView.x / 100
    Player.z = Player.z + PlayerView.z / 100
  EndIf
  If Btst(GetAsyncKeyState(VK_DOWN), 15) = True
    Player.x = Player.x - PlayerView.x / 100
    Player.z = Player.z - PlayerView.z / 100
  EndIf

  Player.x = Max(0.01, Min(18.99, Player.x))
  Player.z = Max(0.01, Min(18.99, Player.z))

  DrawGround

  SwapBuffers(OpenGL_DC)

  ' Usually it's better to use Sleep instead PeekEvent and have an OCX timer to draw screen
  ' For this exemple it's not important
  PeekEvent
Loop Until Me Is Nothing
wglDeleteContext(OpenGL_RDC) ' Never forget this!
Function InterpolationZfromXY(x1#, y1#, z1#, x2#, y2#, z2#, x3#, y3#, z3#, x4#, y4#) As Double

  Local Double xn, yn, zn, z4

  xn = (y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1)
  yn = (x3 - x1) * (z2 - z1) - (x2 - x1) * (z3 - z1)
  zn = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)

  z4 = (x1 * xn + y1 * yn + z1 * zn - x4 * xn - y4 * yn) / zn

  Return z4

EndFunc
Proc Turn(ByRef x#, ByRef y#, a#)

  Static Double s, c, newx, newy

  s = Sin(a)
  c = Cos(a)
  newx = c * x + s * y
  newy = c * y - s * x
  x = newx
  y = newy

EndProc
Proc InitialisationOpenGL

  Global Long OpenGL_DC, OpenGL_RDC, PixelFormat
  Global PFD As PIXELFORMATDESCRIPTOR

  Form Full OpenGLWin

  PFD.nSize = SizeOf(PFD)
  PFD.nVersion = 1
  PFD.dwFlags = PFD_SUPPORT_OPENGL | PFD_DOUBLEBUFFER | PFD_DRAW_TO_WINDOW
  PFD.cDepthBits = 32

  OpenGL_DC = Me.hDC
  PixelFormat = ChoosePixelFormat(OpenGL_DC, PFD)
  Void SetPixelFormat(OpenGL_DC, PixelFormat, PFD)
  OpenGL_RDC = wglCreateContext(OpenGL_DC)

  If !wglMakeCurrent(OpenGL_DC, OpenGL_RDC)
    Message "OpenGL can not be Initialised."
    End
  EndIf

EndProc
Proc DrawGround

  Local Long X, z
  Local Double x1, y1, z1, x2, y2, z2, x3, y3, z3, a
  glBindTexture(GL_TEXTURE_2D, Texture)
  glBegin(GL_TRIANGLES)
  yGround = 0

  For z = 0 To 18
    For X = 0 To 18

      a = Sqr(((X - 9) / 9) ^ 2 + ((z - 9) / 9) ^ 2) * 4

      glColor3d(0.5 + Sin(a) / 2, 0.5 + Sin(a + 1) / 2, 0.5 + Sin(a + 2) / 2)
      x1 = X : y1 = Ground(X, z) : z1 = z
      x2 = X + 1 : y2 = Ground(X + 1, z) : z2 = z
      x3 = X + 1 : y3 = Ground(X + 1, z + 1) : z3 = z + 1

      If PointInTriangleLeftRight(Player.x, Player.z, x1, z1, x2, z2, x3, z3)
        yGround = InterpolationZfromXY(x1, z1, y1, x2, z2, y2, x3, z3, y3, Player.x, Player.z)
      EndIf

      glTexCoord2d(0, 0) : glVertex3d(x1, y1, z1)
      glTexCoord2d(1, 0) : glVertex3d(x2, y2, z2)
      glTexCoord2d(1, 1) : glVertex3d(x3, y3, z3)

      x2 = X  : y2 = Ground(X , z + 1) : z2 = z + 1

      If PointInTriangleLeftRight(Player.x, Player.z, x1, z1, x2, z2, x3, z3)
        yGround = InterpolationZfromXY(x1, z1, y1, x2, z2, y2, x3, z3, y3, Player.x, Player.z)
      EndIf

      glTexCoord2d(0, 0) : glVertex3d(x1, y1, z1)
      glTexCoord2d(0, 1) : glVertex3d(x2, y2, z2)
      glTexCoord2d(1, 1) : glVertex3d(x3, y3, z3)

    Next X
  Next z

  glEnd()

EndProc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Function PointInTriangleLeftRight(Px#, Py#, x1#, y1#, x2#, y2#, x3#, y3#) As Boolean Naked

  Static Long s

  s = Sgn((x2 - x1) * (Py - y1) - (y2 - y1) * (Px - x1))
  If s <> Sgn((x3 - x2) * (Py - y2) - (y3 - y2) * (Px - x2)) Return False
  If s <> Sgn((x1 - x3) * (Py - y3) - (y1 - y3) * (Px - x3)) Return False

  Return True

EndFunc
Function RGBAtoTexture(Adr As Long, Width As Long, Height As Long) As Long

  Local Long TextureIndex

  glEnable(GL_TEXTURE_2D)
  glGenTextures(1, V:TextureIndex)
  glBindTexture(GL_TEXTURE_2D, TextureIndex)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, Width, Height, 0, GL_RGBA, GL_UNSIGNED_BYTE, Adr)

  Return TextureIndex

EndFunc


Compute bissector

Compute a bissector

This procedure calculates the position (ix,iy) of the intersection point of the bissector starting from (x2,y2) on the segment (x1,y1)(x3,y3).

Proc GetBissector(x1, y1, x2, y2, x3, y3, ByRef ix As Double, ByRef iy As Double)

  Local Double lx, ly, a, b, abSum, c, u

  lx = x2 - x1 : ly = y2 - y1
  a = Sqr(lx * lx + ly * ly)

  lx = x3 - x2 : ly = y3 - y2
  b = Sqr(lx * lx + ly * ly)

  abSum = a + b

  If abSum = 0
    ix = x2
    iy = y2
    Exit Proc
  EndIf

  lx = x1 - x3 : ly = y1 - y3
  c = Sqr(lx * lx + ly * ly)

  u = b / abSum

  ix = x3 + u * lx
  iy = y3 + u * ly

EndProc

Example

' Calling GetBissector() function to compute the inscribed circle in a triangle
$Library "Direct2D"
Global Double x1, y1, x2, y2, x3, y3, bx1, by1, bx2, by2, bx3, by3, Radius, ix, iy, x, y
Global Double Lx2, Lx1, Ly1, Ly2, U, M1, M2, h
Global Long PMOVE
Global Object Win1RT
OpenW Center 1
Set Win1RT = D2GetRT(Me)
x1 = 100 : y1 = 100 : x2 = 400 : y2 = 50 : x3 = 300 : y3 = 400
Do

  D2BeginDraw Win1RT, D2C_DarkBlue

  ' Interaction for moving points
  If MouseK = 0 PMOVE = 0
  If PMOVE > 0
    If PMOVE = 1 x1 = MouseX : y1 = MouseY
    If PMOVE = 2 x2 = MouseX : y2 = MouseY
    If PMOVE = 3 x3 = MouseX : y3 = MouseY
  EndIf
  If Abs(MouseX - x1) < 10 And Abs(MouseY - y1) < 10 D2PCircle x1, y1, 9, , , , D2Brush(D2C_Yellow)
  If Abs(MouseX - x2) < 10 And Abs(MouseY - y2) < 10 D2PCircle x2, y2, 9, , , , D2Brush(D2C_Yellow)
  If Abs(MouseX - x3) < 10 And Abs(MouseY - y3) < 10 D2PCircle x3, y3, 9, , , , D2Brush(D2C_Yellow)

  D2ForeColor = D2C_White

  ' Draw the triangle
  D2Line x1, y1, x2, y2
  D2Line x2, y2, x3, y3
  D2Line x3, y3, x1, y1

  ' Trace 3 bissectors from 3 points
  GetBissector x2, y2, x1, y1, x3, y3, bx1, by1
  D2Line x1, y1, bx1, by1, D2Brush(ARGB(255, 0, 255, 0))

  GetBissector x1, y1, x2, y2, x3, y3, bx2, by2
  D2Line x2, y2, bx2, by2, D2Brush(ARGB(255, 0, 255, 0))

  GetBissector x2, y2, x3, y3, x1, y1, bx3, by3
  D2Line x3, y3, bx3, by3, D2Brush(ARGB(255, 0, 255, 0))

  ' Calculate Bissectors Intersection
  Lx2 = bx2 - x2 : Lx1 = bx1 - x1
  Ly1 = by1 - y1 : Ly2 = by2 - y2
  M1 = Lx2 * Ly1 : M2 = Ly2 * Lx1
  U = Iif(M1 <> M2, (Lx2 * (y1 - y2) - Ly2 * (x1 - x2)) / (M2 - M1), 0.5)
  ix = x1 + U * Lx1 : iy = y1 + U * Ly1
  D2Circle ix, iy, 5

  'Calculate Radius = Perpendicular
  Lx1 = x2 - x1 : Ly1 = y2 - y1
  h = Lx1 * Lx1 + Ly1 * Ly1
  If h = 0
    Radius = 0
  Else
    U = (((ix - x1) * Lx1) + ((iy - y1) * Ly1)) / h
    x = (x1 + U * Lx1) : y = (y1 + U * Ly1)
    D2Line ix, iy, x, y, D2Brush(D2C_Gray)
    x = ix - x : y = iy - y
    Radius = Sqr(x * x + y * y)
  EndIf
  D2Circle ix, iy, Radius
  D2EndDraw

  Sleep

Loop Until Me Is Nothing
Sub Win_1_MouseDown(Button&, Shift&, x!, y!)

  If Abs(MouseX - x1) < 10 And Abs(MouseY - y1) < 10 PMOVE = 1
  If Abs(MouseX - x2) < 10 And Abs(MouseY - y2) < 10 PMOVE = 2
  If Abs(MouseX - x3) < 10 And Abs(MouseY - y3) < 10 PMOVE = 3

EndSub

Proc GetBissector(x1, y1, x2, y2, x3, y3, ByRef ix As Double, ByRef iy As Double)

  Local Double lx, ly, a, b, abSum, c, u

  lx = x2 - x1 : ly = y2 - y1
  a = Sqr(lx * lx + ly * ly)

  lx = x3 - x2 : ly = y3 - y2
  b = Sqr(lx * lx + ly * ly)

  abSum = a + b

  If abSum = 0
    ix = x2
    iy = y2
    Exit Proc
  EndIf

  lx = x1 - x3 : ly = y1 - y3
  c = Sqr(lx * lx + ly * ly)

  u = b / abSum

  ix = x3 + u * lx
  iy = y3 + u * ly

EndProc


Compute vector angle

This function use Atan() to obtain the radian angle of the vector (x,y).

Function GetAngleAtan(x#, y#) As Double

  If y <> 0
    GetAngleAtan = Atan(x / y)
    If y < 0 GetAngleAtan +=  PI
  Else If x < 0
    GetAngleAtan = -PI / 2
  Else
    GetAngleAtan = PI / 2
  EndIf

EndFunc

Example

' Get angle using Atan to simulate magnetic field
$Library "Direct2D"
OpenW Center 1
Global Object win1RT, GeoMagnet
Set win1RT = D2GetRT(Me)
Global Dim vx(0 .. 20, 0 .. 20) As Double
Global Dim vy(0 .. 20, 0 .. 20) As Double
Global Long X, Y, XX, YY
Global Double sx, sy, ox, oy, dx, dy, d, a
Set GeoMagnet = CreateGeoMagnet(0.25)
For Y = 1 To 19
  For X = 1 To 19
    a = Rnd * 6.283 : vx(X, Y) = Sin(a) : vy(X, Y) = Cos(a)
  Next X
Next Y
Do

  D2BeginDraw win1RT, D2C_DarkBlue

  D2ForeColor = D2C_White
  sy = _Y / 2 - 270
  For Y = 1 To 19
    sx = _X / 2 - 270
    For X = 1 To 19

      ' Angle between mini-rotor and magnet
      dx = MouseX - sx : dy = MouseY - sy
      a = GetAngleAtan(dx, dy) * 2 // Muliply angle by 2 = bipolarity

      ' Progressive orientation to have fun animation
      vx(X, Y) += Cos(a) / 20
      vy(X, Y) -= Sin(a) / 20
      d = Sqr(vx(X, Y) ^ 2 + vy(X, Y) ^ 2)
      If d > 0
        vx(X, Y) /= d : vy(X, Y) /= d
      EndIf

      D2Line sx, sy, sx + vx(X, Y) * 15, sy + vy(X, Y) * 15
      D2PCircle sx, sy, 2, , , , D2Brush(D2C_DarkOliveGreen)
      sx += 30
    Next X
    sy += 30
  Next Y

  ' Draw a magnet at mouse position
  D2DefLine , 0
  D2Transform D2TF_MOVE, MouseX, MouseY - 75
  D2FillGeometry GeoMagnet, D2Brush(ARGB(200, 255, 0, 0))
  D2PBox -62.5, 37.5, -25, 75, D2Brush(ARGB(200, 200, 200, 255))
  D2PBox 25, 37.5, 62.5, 75, D2Brush(ARGB(200, 200, 200, 255))

  D2EndDraw

  PeekEvent

Loop Until Me Is Nothing
Function GetAngleAtan(x#, y#) As Double

  If y <> 0
    GetAngleAtan = Atan(x / y)
    If y < 0 GetAngleAtan +=  PI
  Else If x < 0
    GetAngleAtan = -PI / 2
  Else
    GetAngleAtan = PI / 2
  EndIf

EndFunc

Function CreateGeoMagnet(Zoom As Double) As Object

  Local Object GeoSink , CreateGeoResult
  Local P2F As D2D1_POINT_2F
  Local ArcSegment As D2D1_ARC_SEGMENT
  Local Bezier As D2D1_BEZIER_SEGMENT
  Local QBezier As D2D1_QUADRATIC_BEZIER_SEGMENT
  Set CreateGeoResult = D2PathGeometry()
  Set GeoSink = D2GeometrySink(CreateGeoResult)
  ' Figure #1
  D2Sink_BeginFigure GeoSink, -250 * Zoom, -100 * Zoom
  ArcSegment.point.x = 250 * Zoom : ArcSegment.point.y = -100 * Zoom
  ArcSegment.size.width = 250 * Zoom
  ArcSegment.size.height = 250 * Zoom
  ArcSegment.rotationAngle = 0
  ArcSegment.sweepDirection = 1
  ArcSegment.arcSize = 0
  D2Sink_AddArc GeoSink, ArcSegment
  P2F.x = 250 * Zoom : P2F.y = 150 * Zoom
  D2Sink_AddLine GeoSink, P2F
  P2F.x = 100 * Zoom : P2F.y = 150 * Zoom
  D2Sink_AddLine GeoSink, P2F
  P2F.x = 100 * Zoom : P2F.y = -100 * Zoom
  D2Sink_AddLine GeoSink, P2F
  ArcSegment.point.x = -100 * Zoom : ArcSegment.point.y = -100 * Zoom
  ArcSegment.size.width = 100 * Zoom
  ArcSegment.size.height = 100 * Zoom
  ArcSegment.rotationAngle = 0
  ArcSegment.sweepDirection = 0
  ArcSegment.arcSize = 0
  D2Sink_AddArc GeoSink, ArcSegment
  P2F.x = -100 * Zoom : P2F.y = 150 * Zoom
  D2Sink_AddLine GeoSink, P2F
  P2F.x = -250 * Zoom : P2F.y = 150 * Zoom
  D2Sink_AddLine GeoSink, P2F
  D2Sink_EndFigure GeoSink
  D2Sink_Close GeoSink
  Set CreateGeoMagnet = CreateGeoResult

EndFunc


Bi-Sector Cosinus

Compute Bi-Vectors Cosinus

This function compute the cosinus of angle bteween vector (v1x,v1y,v1z) and vector (v2x,v2y,v2z).
The returned value is between -1 and 1.
To compute the cosinus of 2D vectors, pass simply zero on z values v1z and v2z.
This is particularly useful for evaluating a degree of exposure on a surface as part of a ray-tracing program.
Note vector are normalised inside the function and it's not necessary to normalise the vectors before calling it.

Function CosVV(v1x#, v1y#, v1z#, v2x#, v2y#, v2z#) As Double

  Local Double d

  d = Sqr(v1x * v1x + v1y * v1y + v1z * v1z)
  If d > 0
    v1x /= d
    v1y /= d
    v1z /= d
  Else If v1x <> 0
    v1x = Sgn(v1x)
    v1y = 0
    v1z = 0
  Else If v1y <> 0
    v1x = 0
    v1y = Sgn(v1y)
    v1z = 0
  Else If v1z <> 0
    v1x = 0
    v1y = 0
    v1z = Sgn(v1y)
  Else
    v1x = 1
    v1y = 0
    v1z = 0
  EndIf

  d = Sqr(v2x * v2x + v2y * v2y + v2z * v2z)
  If d > 0
    v2x /= d
    v2y /= d
    v2z /= d
  Else If v2x <> 0
    v2x = Sgn(v2x)
    v2y = 0
    v2z = 0
  Else If v2y <> 0
    v2x = 0
    v2y = Sgn(v2y)
    v2z = 0
  Else If v2z <> 0
    v2x = 0
    v2y = 0
    v2z = Sgn(v2y)
  Else
    v2x = 1
    v2y = 0
    v2z = 0
  EndIf

  Return (v1x * v2x + v1y * v2y + v1z * v2z)

EndFunc

Example

$Library "Gdip.inc"

Type XYZ_STRUCT
  - Double x, y, z
EndType

Global Long Idx, Aw , X , Y, Mx, My, Bitmap , AdrBitmap, WH = 400
Global Long BSize = WH * WH * 4
Global Double cvv

Global BIH As BITMAPINFOHEADER
Global Ray As XYZ_STRUCT

OpenW 1
Me.Adjust WH, WH
Me.Center True

' Create a bitmap and the corresponding array of normals
CreateSurface

Do

  Idx = 0
  Aw = AdrBitmap + WH * WH * 4
  Mx = (MouseX - (_X / 2 - WH / 2))
  My = ((_Y / 2 + WH / 2) - MouseY)

  For Y = 0 To WH - 1
    For X = 0 To WH - 1

      ' A Ray incoming from mouse position to pixel position
      Ray.x = X - Mx
      Ray.y = Y - My
      Ray.z = -WH / 10
      XYZ_Normalise Ray

      ' Compute the cosinus between a Ray and a normal of surface
      cvv = CosVV(N(Idx).x, N(Idx).y, N(Idx).z, Ray.x, Ray.y, Ray.z)

      ' Use this Cosinus to apply lighting exposure
      cvv =  Max(0, Min(1, ( 1 - (cvv + 1)  )))
      Poke Aw, cvv * Peek(Aw - BSize)
      Poke Aw + 1, cvv * Peek(Aw + 1 - BSize)
      Poke Aw + 2, cvv * Peek(Aw + 2 - BSize)
      If cvv > 0.75
        cvv = ((cvv - 0.75) * 4) ^ 5
        Poke Aw, Peek(Aw) + cvv * (255 - Peek(Aw))
        Poke Aw + 1, Peek(Aw + 1) + cvv * (255 - Peek(Aw + 1))
        Poke Aw + 2, Peek(Aw + 2) + cvv * (255 - Peek(Aw + 2))
      EndIf

      Inc Idx
      Add Aw, 4
    Next X

  Next Y

  ' Show the result
  Put _X / 2 - WH / 2, _Y / 2 - WH / 2, Bitmap

  Sleep

Loop Until Me Is Nothing
FreeBmp Bitmap
Proc CreateSurface

  Local Long X, Y, G, hDc, hFont
  Local R As RECT

  BIH.biSize = SizeOf(BIH)
  BIH.biWidth = WH
  BIH.biHeight = WH * 2
  BIH.biPlanes = 1
  BIH.biBitCount = 32
  Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)
  hDc = CreateCompatibleDC(Win_1.hDC)
  Void SelectObject(hDc, Bitmap)

  Aw = AdrBitmap
  For Y = 0 To WH - 1
    For X = 0 To WH - 1
      If X > 15 And X < WH - 15 And Y > 15 And Y < WH - 15
        G = (1 - Rnd * Rnd) * 200 * (1 - (Max(Abs(X - WH / 2), Abs(Y - WH / 2)) / WH) ^ 0.75)
      Else
        G = (Int(X + Y) % 10) * 20
      EndIf
      LPoke Aw, Bswap3(_RGB(G + Rnd * (X / WH) * 50,  G + Rnd * 21, G + Rnd * (Y / WH) * 50))
      Add Aw, 4
    Next X
  Next Y

  Contrast32Bits(AdrBitmap, WH, WH, 1.4)
  GaussianBlur32bits(AdrBitmap, WH, WH, 4)
  hFont = CreateFont(48, 0, 0, 0, FW_BOLD, False, True, False, DEFAULT_CHARSET, OUT_OUTLINE_PRECIS, CLIP_DEFAULT_PRECIS, 0, VARIABLE_PITCH, Wide("Impact"))
  Void SelectObject(hDc, hFont)
  R.Left = 0 : R.Top = WH : R.Right = WH - 1 : R.Bottom = WH * 2 - 1
  Void SetBkMode(hDc, TRANSPARENT)
  Void SetTextColor(hDc, 0)
  Void DrawText(hDc, "GFA-BASIC 32", 12, R, DT_CENTER | DT_VCENTER | DT_SINGLELINE)
  GaussianBlur32bits(AdrBitmap, WH, WH, 5)
  AutoLevel32Bits(AdrBitmap, WH, WH )
  Contrast32Bits(AdrBitmap, WH, WH, 10)
  Void SetTextColor(hDc, &h4FFF00)
  Void DrawText(hDc, "GFA-BASIC 32", 12, R, DT_CENTER | DT_VCENTER | DT_SINGLELINE)

  Global Dim N(WH * WH) As XYZ_STRUCT
  NormalsFromHeights(AdrBitmap, WH, WH, 32, 100, 3, N())

  Void DeleteObject(hFont)
  Void DeleteDC(hDc)

EndProc
Proc GaussianBlur32bits(Adr As Long, w As Long, h As Long, blurlevel As Double)

  Local Long X, Y, XX, YY, Radius, Aw, Ar
  Local Double Distance, Alpha

  If blurlevel = 0 Exit Proc

  Local Dim R(w, h) As Double
  Local Dim G(w, h) As Double
  Local Dim B(w, h) As Double
  Local Dim Sum(w, h) As Double

  ArrayFill R(), 0
  ArrayFill G(), 0
  ArrayFill B(), 0
  ArrayFill Sum(), 0

  Radius = blurlevel + 2
  Aw = Adr
  For Y = 0 To h - 1
    For X = 0 To w - 1
      For YY = Y - Radius To Y + Radius
        If YY >= 0 And YY < h
          For XX = X - Radius To X + Radius
            If XX >= 0 And XX < w
              Distance = Sqr((XX - X) ^ 2 + (YY - Y) ^ 2)
              If Distance <= blurlevel
                Alpha = (blurlevel - Distance) / blurlevel
                Sum(X, Y) += Alpha
                Ar = Adr + (XX + YY * w) * 4
                R(X, Y) += Alpha * Peek(Ar)
                G(X, Y) += Alpha * Peek(Ar + 1)
                B(X, Y) += Alpha * Peek(Ar + 2)
              EndIf
            EndIf
          Next XX
        EndIf
      Next YY
      Add Aw, 4
    Next X
  Next Y

  Aw = Adr
  For Y = 0 To h - 1
    For X = 0 To w - 1
      Poke Aw, R(X, Y) / Sum(X, Y)
      Poke Aw + 1, G(X, Y) / Sum(X, Y)
      Poke Aw + 2, B(X, Y) / Sum(X, Y)
      Add Aw, 4
    Next X
  Next Y

  Erase R(), G(), B(), Sum()

EndProc
Function AutoLevel32Bits(Adr As Long, W As Long, h As Long) As Boolean

  Local Long Channel
  Local Long ValueMin, ValueMax
  Local Long X, Y, Aw
  Local Double ValueMul

  Aw = Adr

  ValueMin = 255
  ValueMax = 0

  For Channel = 0 To 2
    For Aw = Adr + Channel To Adr + Channel + W * h * 4 - 4 Step 4
      ValueMin = Min(ValueMin, Peek(Aw))
      ValueMax = Max(ValueMax, Peek(Aw))
    Next Aw
  Next Channel

  If ValueMax > ValueMin
    For Channel = 0 To 2
      ValueMul = 255 / (ValueMax - ValueMin)
      For Aw = Adr + Channel To Adr + Channel + W * h * 4 - 4 Step 4
        Poke Aw, CByte((Peek(Aw) - ValueMin) * ValueMul)
      Next Aw
    Next Channel

  EndIf

EndFunc
Function Contrast32Bits(Adr As Long, W As Long, h As Long, intensity As Double) As Boolean

  Local Long valueMin, ValueMax
  Local Double ValueMed
  Local Long X, Y, Aw

  Aw = Adr

  valueMin = 255
  ValueMax = 0

  For Aw = Adr To Adr + W * h * 4 - 4 Step 4
    valueMin = Min(valueMin, Peek(Aw), Peek(Aw + 1), Peek(Aw + 2))
    ValueMax = Max(ValueMax, Peek(Aw), Peek(Aw + 1), Peek(Aw + 2))
  Next Aw

  ValueMed = (valueMin + ValueMax) / 2

  For Aw = Adr To Adr + W * h * 4 - 4 Step 4
    Poke Aw, CByte(Max(0, Min(255, ValueMed + (Peek(Aw) - ValueMed) * intensity)))
    Poke Aw + 1, CByte(Max(0, Min(255, ValueMed + (Peek(Aw + 1) - ValueMed) * intensity)))
    Poke Aw + 2, CByte(Max(0, Min(255, ValueMed + (Peek(Aw + 2) - ValueMed) * intensity)))
  Next Aw

EndFunc
Proc NormalsFromHeights(AdrBitmap As Long, Width As Long, Height As Long, BitCount As Long, Intensity As Double, Softening As Long, ByRef N() As XYZ_STRUCT)

  Local Long X, Y, X1, Y1, X2, Y2, i, j, S, D, idx, Ar, ArStep
  Local Double Z(0 .. Width - 1, 0 .. Height - 1), z1, z2
  Local Double nx, ny, nz, nd, ndSum
  Local Long Px(0 .. 8), Py(0 .. 8)
  Local Ns As XYZ_STRUCT

  Px(0) = -1 : Py(0) = -1
  Px(1) =  0 : Py(1) = -1
  Px(2) =  1 : Py(2) = -1
  Px(3) =  1 : Py(3) =  0
  Px(4) =  1 : Py(4) =  1
  Px(5) =  0 : Py(5) =  1
  Px(6) = -1 : Py(6) =  1
  Px(7) = -1 : Py(7) =  0
  Px(8) = -1 : Py(8) = -1

  If Dim?(N()) < Width * Height + 1
    ReDim N(Width * Height + 1)
  EndIf

  Ar = AdrBitmap
  ArStep = BitCount \ 8

  For Y = 0 To Height - 1
    For X = 0 To Width - 1
      Select BitCount
      Case 8
        Z(X, Y) = Intensity * Peek(Ar)
      Case 16
        Z(X, Y) = ((CPeek(Ar) >> 10) + And(CPeek(Ar) >> 5, &x11111) + And(CPeek(Ar), &x11111)) / 93
      Case 24, 32
        Z(X, Y) = Intensity * Max(Peek(Ar) , Peek(Ar + 1), Peek(Ar + 2)) / 255
      EndSelect
      Add Ar, ArStep
    Next X
  Next Y

  For Y = 0 To Height - 1
    For X = 0 To Width - 1

      ndSum = 0
      For i = 0 To 7

        X1 = X + Px(i)     : Y1 = Y + Py(i)
        X2 = X + Px(i + 1) : Y2 = Y + Py(i + 1)
        z1 = Iif(X1 >= 0 And X1 < Width And Y1 >= 0 And Y1 < Height, Z(X1, Y1), Z(X, Y))
        z2 = Iif(X2 >= 0 And X2 < Width And Y2 >= 0 And Y2 < Height, Z(X2, Y2), Z(X, Y))

        nx = Y1 * (z2 - Z(X, Y)) + Y2 * (Z(X, Y) - z1) + Y * (z1 - z2)
        ny = z1 * (X2 - X) + z2 * (X - X1) + Z(X, Y) * (X1 - X2)
        nz = X1 * (Y2 - Y) + X2 * (Y - Y1) + X * (Y1 - Y2)

        nd = Sqr(nx * nx + ny * ny + nz * nz)

        N(idx).x += nx / nd
        N(idx).y += ny / nd
        N(idx).z += nz / nd

        ndSum += nd

      Next i

      N(idx).x /= ndSum
      N(idx).y /= ndSum
      N(idx).z /= ndSum

      Inc idx

    Next X
  Next Y

  If Softening > 0

    Local Dim Nb(idx) As XYZ_STRUCT

    For S = 1 To Softening

      idx = 0
      For Y = 0 To Height - 1
        For X = 0 To Width - 1
          D = 0
          Nb(idx) = N(idx)
          For i = 0 To 7
            X1 = X + Px(i)     : Y1 = Y + Py(i)
            If X1 >= 0 And X1 < Width And Y1 >= 0 And Y1 < Height
              j = X1 + Y1 * Width
              Inc D
              Nb(idx).x += N(j).x
              Nb(idx).y += N(j).y
              Nb(idx).z += N(j).z
            EndIf
          Next i
          Nb(idx).x /= D
          Nb(idx).y /= D
          Nb(idx).z /= D
          Inc idx
        Next X
      Next Y

      For i = 0 To idx - 1
        N(i) = Nb(i)
        XYZ_Normalise N(i)
      Next i

    Next S

    Erase Nb()

  EndIf

EndProc
Function CosVV(v1x#, v1y#, v1z#, v2x#, v2y#, v2z#) As Double

  Local Double d

  d = Sqr(v1x * v1x + v1y * v1y + v1z * v1z)
  If d > 0
    v1x /= d
    v1y /= d
    v1z /= d
  Else If v1x <> 0
    v1x = Sgn(v1x)
    v1y = 0
    v1z = 0
  Else If v1y <> 0
    v1x = 0
    v1y = Sgn(v1y)
    v1z = 0
  Else If v1z <> 0
    v1x = 0
    v1y = 0
    v1z = Sgn(v1z)
  Else
    v1x = 1
    v1y = 0
    v1z = 0
  EndIf

  d = Sqr(v2x * v2x + v2y * v2y + v2z * v2z)
  If d > 0
    v2x /= d
    v2y /= d
    v2z /= d
  Else If v2x <> 0
    v2x = Sgn(v2x)
    v2y = 0
    v2z = 0
  Else If v2y <> 0
    v2x = 0
    v2y = Sgn(v2y)
    v2z = 0
  Else If v2z <> 0
    v2x = 0
    v2y = 0
    v2z = Sgn(v2z)
  Else
    v2x = 1
    v2y = 0
    v2z = 0
  EndIf

  Return (v1x * v2x + v1y * v2y + v1z * v2z)

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc


Computing Triangle Area

This function return the area of the given triangle (x1,y1,x2,y2,x3,y3) using the Heron's formula.

Function TriangleArea(x1#, y1#, x2#, y2#, x3#, y3#) As Double

  Static Double a, b, c, p

  a = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
  b = Sqr((x3 - x2) ^ 2 + (y3 - y2) ^ 2)
  c = Sqr((x1 - x3) ^ 2 + (y1 - y3) ^ 2)
  p = (a + b + c) / 2

  Return Sqr(p * (p - a) * (p - b) * (p - c))

EndFunc

Example

$Library "Direct2D"
Global Vertice(0 .. 3) As D2D1_POINT_2F
Global Long VerticeMoving = 0, VerticeHover = 0
Global Double Area
OpenW Center 1
Global Object Win1RT, Geo, GeoSink
Set Win1RT = D2GetRT(Me)
Array Vertice() = Mks(0, 0, _X / 4, _Y / 10, _X / 1.5, _Y / 8, _X / 2.5, _Y / 1.2)
Do

  Area = TriangleArea(Vertice(1).x, Vertice(1).y, Vertice(2).x, Vertice(2).y, Vertice(3).x, Vertice(3).y)

  DrawAll

  Sleep

Loop Until Me Is Nothing
Function TriangleArea(x1#, y1#, x2#, y2#, x3#, y3#) As Double

  Static Double a, b, c, p

  a = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
  b = Sqr((x3 - x2) ^ 2 + (y3 - y2) ^ 2)
  c = Sqr((x1 - x3) ^ 2 + (y1 - y3) ^ 2)
  p = (a + b + c) / 2

  Return Sqr(p * (p - a) * (p - b) * (p - c))

EndFunc
Proc DrawAll

  Local Long i

  D2BeginDraw Win1RT, D2C_DarkBlue

  For i = 1 To 3
    D2Line Vertice(i).x, Vertice(i).y, Vertice(1 + (i + 1) % 3).x, Vertice(1 + (i + 1) % 3).y, D2Brush(D2C_White)
  Next i
  For i = 1 To 3
    D2PCircle Vertice(i).x, Vertice(i).y, 5, , , , D2Brush(Iif(i = Max(VerticeMoving, VerticeHover), D2C_Yellow, D2C_Gray))
  Next i
  Set Geo = D2PathGeometry()
  Set GeoSink = D2GeometrySink(Geo)
  D2Sink_BeginFigure GeoSink, Vertice(1).x, Vertice(1).y
  D2Sink_AddLine GeoSink, Vertice(2)
  D2Sink_AddLine GeoSink, Vertice(3)
  D2Sink_EndFigure GeoSink
  D2Sink_Close GeoSink
  D2FillGeometry Geo, D2Brush(ARGB(100, 100, 200, 200))
  D2SetFont "Arial", 40
  D2DrawText 5, _Y / 2, _X - 5, _Y - 5, "Triangle Area = " & Round(Area, 2) & " px", D2T_LEFT | D2T_PARBOTTOM , , D2Brush(D2C_White)
  D2EndDraw

EndProc
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  Local Long i

  If VerticeMoving = 0
    VerticeHover = 0
    For i = 1 To 3
      If (MouseX - Vertice(i).x) ^ 2 + (MouseY - Vertice(i).y) ^ 2 < 100
        VerticeHover = i
      EndIf
    Next i
  Else
    Vertice(VerticeMoving).x = MouseX
    Vertice(VerticeMoving).y = MouseY
  EndIf

EndSub
Sub Win_1_MouseDown(Button&, Shift&, x!, y!)

  Local Long i
  For i = 1 To 3
    If (MouseX - Vertice(i).x) ^ 2 + (MouseY - Vertice(i).y) ^ 2 < 100
      VerticeMoving = i
    EndIf
  Next i

EndSub
Sub Win_1_MouseUp(Button&, Shift&, x!, y!)

  VerticeMoving = 0

EndSub


Point In Triangle

Determine if a point is in a triangle

This function return True if the point Px,Py is in the triangle x1,y1,x2,y2,x3,y3. The algorithm is based on the segment left/right side test.

Function PointInTriangleLeftRight(Px#, Py#, x1#, y1#, x2#, y2#, x3#, y3#) As Boolean Naked

  Static Long s

  s = Sgn((x2 - x1) * (Py - y1) - (y2 - y1) * (Px - x1))
  If s <> Sgn((x3 - x2) * (Py - y2) - (y3 - y2) * (Px - x2)) Return False
  If s <> Sgn((x1 - x3) * (Py - y3) - (y1 - y3) * (Px - x3)) Return False

  Return True

EndFunc

Example

OpenW Center 1
Global V(0 .. 2) As POINT, i As Long
Global Double x, y

OpenW Center 1
Array V() = Mkl(_X / 4, _Y / 10, _X / 1.5, _Y / 8, _X / 2.5, _Y / 1.2)
Do

  Draw V(1).x, V(1).y To V(2).x, V(2).y To V(0).x, V(0).y To V(1).x, V(1).y
  For i = 1 To 1000
    x = Rnd * 100 - 50
    y = MouseY + 2 * (Rnd - 0.5) * Sqr(2500 - x * x)
    x = MouseX + x
    If PointInTriangleLeftRight(x, y, V(1).x, V(1).y, V(2).x, V(2).y, V(0).x, V(0).y)
      Plot x, y, RGB(255, 0, 0)
    Else
      Plot x, y, RGB(0, 255, 0)
    EndIf
  Next i

  Sleep

Loop Until Me Is Nothing
Function PointInTriangleLeftRight(Px#, Py#, x1#, y1#, x2#, y2#, x3#, y3#) As Boolean Naked

  Static Long s

  s = Sgn((x2 - x1) * (Py - y1) - (y2 - y1) * (Px - x1))
  If s <> Sgn((x3 - x2) * (Py - y2) - (y3 - y2) * (Px - x2)) Return False
  If s <> Sgn((x1 - x3) * (Py - y3) - (y1 - y3) * (Px - x3)) Return False

  Return True

EndFunc


Intersection Circles

Intersection between 2 circles

This function return True if intersection found between circles x1,y1,r1 and x2,y2,r2.
Then 2 intersecting points are returned by ix1,iy1 and ix2,iy2.

Function IntersectionsCircles(x1#, y1#, r1#, x2#, y2#, r2#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Local Double a, b, d, h, dx, dy, x3, y3

  dx = x2 - x1
  dy = y2 - y1
  d = Sqr(dx * dx + dy * dy)

  If d > r1 + r2  Return False
  If d < Abs(r2 - r1) Return False
  If d = 0 Return False
  r1 = r1 * r1
  a = (r1 - r2 * r2 + d * d) / (2 * d)
  h = r1 - a * a
  If h < 0 Return False

  h = Sqr(h)
  dx /= d
  dy /= d
  x3 = x1 + a * dx
  y3 = y1 + a * dy
  dx *= h
  dy *= h
  ix1 = x3 + dy
  iy1 = y3 - dx
  ix2 = x3 - dy
  iy2 = y3 + dx

  Return True

EndFunc

Example

$Library "Direct2D"
Type Bubble_Struct
  - Double x, y, r, dx, dy
EndType
Global Object Win1RT
Global Long i, Bubbles = 20
Global Dim Bubble(1 .. Bubbles) As Bubble_Struct
OpenW Center 1
For i = 1 To Bubbles
  Bubble(i).x = Rnd * _X
  Bubble(i).y = Rnd * _Y
  Bubble(i).dx = Rnd - 0.5
  Bubble(i).dy = Rnd - 0.5
  Bubble(i).r = 50 + Rnd * 100
Next i
Set Win1RT = D2GetRT(Me)
Ocx Timer FPS : .Interval = 20 : .Enabled = True
Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer

  DrawAll

EndSub
Function IntersectionsCircles(x1#, y1#, r1#, x2#, y2#, r2#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Local Double a, b, d, h, dx, dy, x3, y3

  dx = x2 - x1
  dy = y2 - y1
  d = Sqr(dx * dx + dy * dy)

  If d > r1 + r2  Return False
  If d < Abs(r2 - r1) Return False
  If d = 0 Return False
  r1 = r1 * r1
  a = (r1 - r2 * r2 + d * d) / (2 * d)
  h = r1 - a * a
  If h < 0 Return False

  h = Sqr(h)
  dx /= d
  dy /= d
  x3 = x1 + a * dx
  y3 = y1 + a * dy
  dx *= h
  dy *= h
  ix1 = x3 + dy
  iy1 = y3 - dx
  ix2 = x3 - dy
  iy2 = y3 + dx

  Return True

EndFunc
Proc DrawAll

  Local Long i, j
  Local Double ix1, iy1, ix2, iy2, d
  Static Boolean Work = False
  If Work Exit Proc
  Work = True

  D2BeginDraw Win1RT, D2C_DarkBlue
  D2ForeColor = D2C_White
  For i = 1 To Bubbles
    D2DefLine , 0
    D2PCircle Bubble(i).x, Bubble(i).y, Bubble(i).r, , , , D2Brush(ARGB(100, 255, (i * 7) % 200, (i * 13) % 200))
    Bubble(i).x += Bubble(i).dx
    If Bubble(i).x < 0 Or Bubble(i).x > _X
      Bubble(i).dx = -Bubble(i).dx
      Bubble(i).x = Max(0, Min(_X, Bubble(i).x))
    EndIf
    Bubble(i).y += Bubble(i).dy
    If Bubble(i).y < 0 Or Bubble(i).y > _Y
      If Bubble(i).y > _Y Bubble(i).dy = -Bubble(i).dy
      Bubble(i).y = Max(0, Min(_Y, Bubble(i).y))
    EndIf
    For j = i + 1 To Bubbles
      If IntersectionsCircles(Bubble(i).x, Bubble(i).y, Bubble(i).r, Bubble(j).x, Bubble(j).y, Bubble(j).r, ix1, iy1, ix2, iy2)
        D2DefLine , 1
        D2Line ix1, iy1, ix2, iy2, D2Brush(D2C_White)
        d = Sqr((ix2 - ix1) ^ 2 + (iy2 - iy1) ^ 2) - Min(Bubble(i).r, Bubble(j).r)
        Bubble(i).dx += Sgn(d) * (Bubble(i).x - Bubble(j).x) / 1000
        Bubble(j).dx -= Sgn(d) * (Bubble(i).x - Bubble(j).x) / 1000
        Bubble(i).dy += Sgn(d) * (Bubble(i).y - Bubble(j).y) / 1000
        Bubble(j).dy -= Sgn(d) * (Bubble(i).y - Bubble(j).y) / 1000
      EndIf
    Next j
    Bubble(i).dx /= 1.01
    Bubble(i).dy /= 1.01
  Next i
  D2EndDraw

  Work = False

EndProc


Box intersection

Rectangles Intersection

This function return True is 2 rectangles (x1,y1,x2,y2) (x3,y3,x4,y4) have a commun area.The common area is returned by (x5,y5,x6,y6).
It does the same thing as the IntersectRect() function except that it works on doubles.

Function IntersectionBoxBox(x1#, y1#, x2#, y2#, x3#, y3#, x4#, y4#, ByRef x5 As Double, ByRef y5 As Double, ByRef x6 As Double, ByRef y6 As Double)  As Boolean

  If Abs(x1 + x2 - x3 - x4) - Abs(x2 - x1) - Abs(x4 - x3) < 0
    If Abs(y1 + y2 - y3 - y4) - Abs(y2 - y1) - Abs(y4 - y3) < 0
      x5 = Max(Min(x1, x2), Min(x3, x4))
      y5 = Max(Min(y1, y2), Min(y3, y4))
      x6 = Min(Max(x1, x2), Max(x3, x4))
      y6 = Min(Max(y1, y2), Max(y3, y4))
      Return True
    EndIf
  EndIf

  Return False

EndFunc

Example

$Library "Direct2D"
OpenW Center 1
Global Object Win1RT
Set Win1RT = D2GetRT(Me)
Ocx Timer FPS
FPS.Enabled = True
FPS.Interval = 30
Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer

  DrawAll

EndSub
Proc DrawAll

  Local Double x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6, u
  Local R1 As RECT, R2 As RECT, R3 As RECT
  D2BeginDraw Win1RT, D2C_DarkBlue

  x1 = _X / 2 - 200
  y1 = _Y / 2 - 100
  x2 = x1 + 400
  y2 = y1 + 200
  x3 = MouseX - 80
  y3 = MouseY - 180
  x4 = x3 + 160
  y4 = y3 + 360
  D2DefLine , 2
  D2ForeColor = D2C_White
  D2Box x1, y1, x2, y2
  D2Box x3, y3, x4, y4
  If IntersectionBoxBox(x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6)
    u = Frac(Timer / 2)
    D2DefLine , 0
    D2PBox x5, y5, x6, y6, D2Brush(ARGB(150, 0, 0, 0))
    D2PBox x5 + u * 10, y5 - u * 10, x6 + u * 10, y6 - u * 10, D2Brush(ARGB(255 - u * 250, 255, 200, 50))
  EndIf

  ' The Winuser function IntersectRect() do the same thing working with integers :
  R1.Left = x1 : R1.Top = y1 : R1.Right = x2 : R1.Bottom = y2
  R2.Left = x3 : R2.Top = y3 : R2.Right = x4 : R2.Bottom = y4
  If IntersectRect(*R3, R1, R2)
    D2DefLine 2, 2
    D2Box R3.Left, R3.Top, R3.Right, R3.Bottom, D2Brush(D2C_Red)
  EndIf

  D2EndDraw

EndProc
Function IntersectionBoxBox(x1#, y1#, x2#, y2#, x3#, y3#, x4#, y4#, ByRef x5 As Double, ByRef y5 As Double, ByRef x6 As Double, ByRef y6 As Double)  As Boolean

  If Abs(x1 + x2 - x3 - x4) - Abs(x2 - x1) - Abs(x4 - x3) < 0
    If Abs(y1 + y2 - y3 - y4) - Abs(y2 - y1) - Abs(y4 - y3) < 0
      x5 = Max(Min(x1, x2), Min(x3, x4))
      y5 = Max(Min(y1, y2), Min(y3, y4))
      x6 = Min(Max(x1, x2), Max(x3, x4))
      y6 = Min(Max(y1, y2), Max(y3, y4))
      Return True
    EndIf
  EndIf

  Return False

EndFunc


Intersection 2 segments

Intersection between 2 lines segments

This function return True if exist intersection between segments (x1,y1,x2,y2) (x3,y3,x4,y4). Then ix,iy is the intersecting point.

Function IntersectionSegments(x1#, y1#, x2#, y2#, x3#, y3#, x4#, y4#, ByRef ix#, ByRef iy#) As Boolean

  Local Double Lx2, Lx1, Ly, Lx, Ly1, Ly2, M, Ua, Ub, M1, M2

  Lx2 = x4 - x3
  Lx1 = x2 - x1
  Ly1 = y2 - y1
  Ly2 = y4 - y3

  M1 = Lx2 * Ly1
  M2 = Ly2 * Lx1

  If M1 <> M2  ' If the lines are not parallel

    Ly = y1 - y3
    Lx = x1 - x3
    M = M2 - M1

    Ua = ( Lx2 * Ly - Ly2 * Lx ) / M
    Ub = ( Lx1 * Ly - Ly1 * Lx ) / M

    If Ua > 0
      If Ua < 1
        If Ub > 0
          If Ub < 1
            ix = x1 + Ua * Lx1
            iy = y1 + Ua * Ly1
            Return True
          EndIf
        EndIf
      EndIf
    EndIf

  EndIf

  Return False

EndFunc

Example

$Library "Direct2D"
OpenW Center 1
Global Object Win1RT
Set Win1RT = D2GetRT(Me)
Do
  DrawAll
  Sleep
Loop Until Me Is Nothing
Proc DrawAll

  Local Double cx, cy, x1, y1, x2, y2, x3, y3, x4, y4, ix, iy, a, pa, Radius,  decal
  Local Double s, c, ox1, oy1, ox2, oy2, NewX
  Static Long N = 16
  Local Long i, j

  pa = 2 * PI / N
  s = Sin(pa)
  c = Cos(pa)
  D2BeginDraw Win1RT, D2C_DarkBlue
  cx = _X / 2 : cy = _Y / 2
  D2DefLine , 1
  ox1 = 0 : oy1 = _Y / 2.1
  ox2 = MouseX - cx : oy2 = MouseY - cy

  ' Draw lines segments
  x1 = ox1 : y1 = oy1 : x2 = ox2 : y2 = oy2
  For i = 1 To N
    D2Line cx + x1, cy + y1, cx + x2, cy + y2, D2Brush(D2C_White)
    NewX = c * x1 + s * y1 : y1 = c * y1 - s * x1 : x1 = NewX
    NewX = c * x2 + s * y2 : y2 = c * y2 - s * x2 : x2 = NewX
  Next i

  ' Draw intersections between all segments
  D2DefLine , 3
  x1 = ox1 : y1 = oy1 : x2 = ox2 : y2 = oy2
  For i = 1 To N
    x3 = x1 : y3 = y1 : x4 = x2 : y4 = y2
    For j = i - 1 DownTo 1
      NewX = c * x3 - s * y3 : y3 = c * y3 + s * x3 : x3 = NewX
      NewX = c * x4 - s * y4 : y4 = c * y4 + s * x4 : x4 = NewX
      If IntersectionSegments(x1, y1, x2, y2, x3, y3, x4, y4, ix, iy)
        D2Circle cx + ix, cy + iy, 5, , , , D2Brush(ARGB(200, 255, 255, 100))
      EndIf
    Next j
    NewX = c * x1 + s * y1 : y1 = c * y1 - s * x1 : x1 = NewX
    NewX = c * x2 + s * y2 : y2 = c * y2 - s * x2 : x2 = NewX
  Next i

  D2EndDraw

EndProc
Function IntersectionSegments(x1#, y1#, x2#, y2#, x3#, y3#, x4#, y4#, ByRef ix#, ByRef iy#) As Boolean

  Local Double Lx2, Lx1, Ly, Lx, Ly1, Ly2, M, Ua, Ub, M1, M2

  Lx2 = x4 - x3
  Lx1 = x2 - x1
  Ly1 = y2 - y1
  Ly2 = y4 - y3

  M1 = Lx2 * Ly1
  M2 = Ly2 * Lx1

  If M1 <> M2  ' If the lines are not parallel

    Ly = y1 - y3
    Lx = x1 - x3
    M = M2 - M1

    Ua = ( Lx2 * Ly - Ly2 * Lx ) / M
    Ub = ( Lx1 * Ly - Ly1 * Lx ) / M

    If Ua > 0
      If Ua < 1
        If Ub > 0
          If Ub < 1
            ix = x1 + Ua * Lx1
            iy = y1 + Ua * Ly1
            Return True
          EndIf
        EndIf
      EndIf
    EndIf

  EndIf

  Return False

EndFunc


Intersection Line Ellipse

Intersection between a line and an ellipse

This function return True if line defined by segment (x1,y1,x2,y2) intersect the ellipse centered to (cx,cy) with radius rx and ry.
Ao (Angle Orientation) is angle in Degree that indicate the rotation of ellipse arround the center.
This ellipse definition corresponding to the D2D1_ARC_SEGMENT specification.
Then (ix1,iy1) and (ix2,iy2) are the coordinates of intersections.

Function IntersectionLineEllipse(x1#, y1#, x2#, y2#, cx#, cy#, rx#, ry#, Ao#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Local Double a, b, c, v, U, lx, ly, r, m, cao, sao, xt, yt

  m = rx / ry
  r = rx

  Ao = Ao * PI / 180
  cao = Cos(Ao)
  sao = Sin(Ao)

  xt = cx + sao * (x1 - cx) - cao * (y1 - cy)
  yt = cy + sao * (y1 - cy) + cao * (x1 - cx)
  x1 = xt
  y1 = yt

  xt = cx + sao * (x2 - cx) - cao * (y2 - cy)
  yt = cy + sao * (y2 - cy) + cao * (x2 - cx)
  x2 = xt
  y2 = yt

  x1 = cx + (x1 - cx) * m
  x2 = cx + (x2 - cx) * m

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    a = a * 2
    v = Sqr(v)
    U = (-b - v) / a
    ix1 = cx + ((x1 + U * lx) - cx) / m
    iy1 = y1 + U * ly
    xt = cx + sao * (ix1 - cx) + cao * (iy1 - cy)
    yt = cy + sao * (iy1 - cy) - cao * (ix1 - cx)
    ix1 = xt
    iy1 = yt

    U = (-b + v) / a
    ix2 = cx + ((x1 + U * lx) - cx) / m
    iy2 = y1 + U * ly
    xt = cx + sao * (ix2 - cx) + cao * (iy2 - cy)
    yt = cy + sao * (iy2 - cy) - cao * (ix2 - cx)
    ix2 = xt
    iy2 = yt

    Return True

  EndIf

  Return False

EndFunc

Example

$Library "Direct2D"

Global Double ECenterX, ECenterY, ERadiusX, ERadiusY, EOrientation
Global Double RayX1, RayY1, RayX2, RayY2
Global Double ix1, iy1, ix2, iy2
Global Long PMoving = 0, PHover = 0
Global Object Win1RT

OpenW Center 1
Set Win1RT = D2GetRT(Me)

' Define the ellipse coordinates and radius
ECenterX = _X / 2 : ECenterY = _Y / 2
ERadiusX = _Y / 2.1 : ERadiusY = ERadiusX / 2

' The line segment who can be moved by user
RayX1 = _X / 8 : RayY1 = _Y / 2 : RayX2 = _X / 1.25 : RayY2 = _Y / 2

Ocx Timer FPS
FPS.Interval = 30 : FPS.Enabled = True
Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer

  DrawAll

EndSub
Function IntersectionLineEllipse(x1#, y1#, x2#, y2#, cx#, cy#, rx#, ry#, Ao#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Local Double a, b, c, v, U, lx, ly, r, m, cao, sao, xt, yt

  m = rx / ry
  r = rx

  Ao = Ao * PI / 180
  cao = Cos(Ao)
  sao = Sin(Ao)

  xt = cx + sao * (x1 - cx) - cao * (y1 - cy)
  yt = cy + sao * (y1 - cy) + cao * (x1 - cx)
  x1 = xt
  y1 = yt

  xt = cx + sao * (x2 - cx) - cao * (y2 - cy)
  yt = cy + sao * (y2 - cy) + cao * (x2 - cx)
  x2 = xt
  y2 = yt

  x1 = cx + (x1 - cx) * m
  x2 = cx + (x2 - cx) * m

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    a = a * 2
    v = Sqr(v)
    U = (-b - v) / a
    ix1 = cx + ((x1 + U * lx) - cx) / m
    iy1 = y1 + U * ly
    xt = cx + sao * (ix1 - cx) + cao * (iy1 - cy)
    yt = cy + sao * (iy1 - cy) - cao * (ix1 - cx)
    ix1 = xt
    iy1 = yt

    U = (-b + v) / a
    ix2 = cx + ((x1 + U * lx) - cx) / m
    iy2 = y1 + U * ly
    xt = cx + sao * (ix2 - cx) + cao * (iy2 - cy)
    yt = cy + sao * (iy2 - cy) - cao * (ix2 - cx)
    ix2 = xt
    iy2 = yt

    Return True

  EndIf

  Return False

EndFunc
Proc DrawAll

  Local Long i

  EOrientation = Frac(Timer / 60) * 360

  D2BeginDraw Win1RT, D2C_DarkBlue

  ' Draw the segment of the line
  DrawArrow RayX1, RayY1, RayX2, RayY2, D2C_Yellow

  ' Draw the ellipse
  D2ForeColor = D2C_White
  D2TransformSave
  D2Transform D2TF_ROTATE, EOrientation, ECenterX, ECenterY
  D2PEllipse ECenterX, ECenterY, ERadiusX, ERadiusY, , , , D2Brush(ARGB(100, 50, 50, 255))
  D2TransformRestore

  ' Draw intersections if found
  If IntersectionLineEllipse(RayX1, RayY1, RayX2, RayY2, ECenterX, ECenterY, ERadiusX, ERadiusY, EOrientation, ix1, iy1, ix2, iy2)
    D2PCircle ix1, iy1, 10, , , , D2Brush(D2C_Blue)
    D2PCircle ix2, iy2, 10, , , , D2Brush(D2C_Red)
    D2Text ix1 + 14, iy1 - 5, "First intersection"
    D2Text ix2 + 14, iy2 - 5, "Second intersection"
  EndIf

  ' Highlighting hovered points of line segment
  Select PHover
  Case 1 : D2PCircle RayX1, RayY1, 5, , , , D2Brush(D2C_Yellow)
  Case 2 : D2PCircle RayX2, RayY2, 5, , , , D2Brush(D2C_Yellow)
  EndSelect

  D2EndDraw

EndProc
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  Local Double d1, d2

  If PMoving = 0
    d1 = Sqr((MouseX - RayX1) ^ 2 + (MouseY - RayY1) ^ 2)
    d2 = Sqr((MouseX - RayX2) ^ 2 + (MouseY - RayY2) ^ 2)
    If d1 < d2 And d1 < 15
      PHover = 1
    Else If d2 < d1 And d2 < 15
      PHover = 2
    Else
      PHover = 0
    EndIf
  Else
    Select PMoving
    Case 1 : RayX1 = MouseX : RayY1 = MouseY
    Case 2 : RayX2 = MouseX : RayY2 = MouseY
    EndSelect
  EndIf

EndSub
Sub Win_1_MouseDown(Button&, Shift&, x!, y!)

  PMoving = PHover

EndSub
Sub Win_1_MouseUp(Button&, Shift&, x!, y!)

  PMoving = 0

EndSub
Proc DrawArrow(x1#, y1#, x2#, y2#, C%)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  If d > 0

    lx = 15 * lx / d
    ly = 15 * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py, D2Brush(C)
    D2Line x2, y2, px + ly , py - lx, D2Brush(C)
    D2Line px + ly , py - lx ,  px - ly , py + lx, D2Brush(C)
    D2Line px - ly , py + lx , x2, y2, D2Brush(C)

  EndIf

EndProc


Intersection line circle

Intersection between a line and a circle

This function return True if line defined by segment (x1,y1,x2,y2) intersect circle (cx,cy,r) then (ix1,iy1) and (ix2,iy2) are the coordinates of intersections.

Function IntersectionLineCircle(x1#, y1#, x2#, y2#, cx#, cy#, r#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Static Double a, b, c, v, U, lx, ly

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    a *= 2

    U = (-b - Sqr(v)) / a
    ix1 = x1 + U * lx
    iy1 = y1 + U * ly

    U = (-b + Sqr(v)) / a
    ix2 = x1 + U * lx
    iy2 = y1 + U * ly

    Return True

  EndIf

  Return False

EndFunc

Example

$Library "Direct2D"

Global Double CenterX, CenterY, Radius
Global Double RayX1, RayY1, RayX2, RayY2
Global Double ix1, iy1, ix2, iy2
Global Long PMoving = 0, PHover = 0
Global Object Win1RT

OpenW Center 1
Set Win1RT = D2GetRT(Me)

' Define the circle coordinates and radius
CenterX = _X / 2 : CenterY = _Y / 2 : Radius = _Y / 2.5

' The line segment who can be moved by user
RayX1 = _X / 8 : RayY1 = _Y / 2 : RayX2 = _X / 1.25 : RayY2 = _Y / 2

Ocx Timer FPS
FPS.Interval = 30 : FPS.Enabled = True
Do
  Sleep
Loop Until Me Is Nothing
Sub FPS_Timer
  DrawAll
EndSub
Function IntersectionLineCircle(x1#, y1#, x2#, y2#, cx#, cy#, r#, ByRef ix1#, ByRef iy1#, ByRef ix2#, ByRef iy2#) As Boolean

  Static Double a, b, c, v, U, lx, ly

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    a *= 2

    U = (-b - Sqr(v)) / a
    ix1 = x1 + U * lx
    iy1 = y1 + U * ly

    U = (-b + Sqr(v)) / a
    ix2 = x1 + U * lx
    iy2 = y1 + U * ly

    Return True

  EndIf

  Return False

EndFunc
Proc DrawAll

  Local Long i

  D2BeginDraw Win1RT, D2C_DarkBlue

  ' Draw the segment of the line
  DrawArrow RayX1, RayY1, RayX2, RayY2, D2C_Yellow

  ' Draw the circle
  D2ForeColor = D2C_White
  D2PCircle CenterX, CenterY, Radius, , , , D2Brush(ARGB(100, 50, 50, 255))

  ' Draw intersections if found
  If IntersectionLineCircle(RayX1, RayY1, RayX2, RayY2, CenterX, CenterY, Radius, ix1, iy1, ix2, iy2)
    D2PCircle ix1, iy1, 10, , , , D2Brush(D2C_Blue)
    D2PCircle ix2, iy2, 10, , , , D2Brush(D2C_Red)
    D2Text ix1 + 14, iy1 - 5, "First intersection"
    D2Text ix2 + 14, iy2 - 5, "Second intersection"
  EndIf

  ' Highlighting hovered points of line segment
  Select PHover
  Case 1 : D2PCircle RayX1, RayY1, 5, , , , D2Brush(D2C_Yellow)
  Case 2 : D2PCircle RayX2, RayY2, 5, , , , D2Brush(D2C_Yellow)
  EndSelect

  D2EndDraw

EndProc
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  Local Double d1, d2

  If PMoving = 0
    d1 = Sqr((MouseX - RayX1) ^ 2 + (MouseY - RayY1) ^ 2)
    d2 = Sqr((MouseX - RayX2) ^ 2 + (MouseY - RayY2) ^ 2)
    If d1 < d2 And d1 < 15
      PHover = 1
    Else If d2 < d1 And d2 < 15
      PHover = 2
    Else
      PHover = 0
    EndIf
  Else
    Select PMoving
    Case 1 : RayX1 = MouseX : RayY1 = MouseY
    Case 2 : RayX2 = MouseX : RayY2 = MouseY
    EndSelect
  EndIf

EndSub
Sub Win_1_MouseDown(Button&, Shift&, x!, y!)

  PMoving = PHover

EndSub
Sub Win_1_MouseUp(Button&, Shift&, x!, y!)

  PMoving = 0

EndSub
Proc DrawArrow(x1#, y1#, x2#, y2#, C%)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  If d > 0

    lx = 15 * lx / d
    ly = 15 * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py, D2Brush(C)
    D2Line x2, y2, px + ly , py - lx, D2Brush(C)
    D2Line px + ly , py - lx ,  px - ly , py + lx, D2Brush(C)
    D2Line px - ly , py + lx , x2, y2, D2Brush(C)

  EndIf

EndProc


Intersection Line Segment Circle

Intersection between line segment and circle

This function return True if segment (x1,y1,x2,y2) intersect the circle (cx,cy,cr) then (ix,iy) is the coordinates of intersection.

Function IntersectionSegmentCircle(x1#, y1#, x2#, y2#, cx#, cy#, r#, ByRef ix#, ByRef iy#) As Boolean

  Static Double a, b, c, v, U, lx, ly

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    v = Sqr(v)
    a = a * 2

    U = (-b - v)  / a

    If U >= 0 And U <= 1

      ix = x1 + U * lx
      iy = y1 + U * ly

      Return True

    EndIf

    U = (-b + v)  / a

    If U >= 0 And U <= 1

      ix = x1 + U * lx
      iy = y1 + U * ly

      Return True

    EndIf

  EndIf

  Return False

EndFunc

Example

$Library "Direct2D"

Global Double CenterX, CenterY, Radius
Global Double RayX1, RayY1, RayX2, RayY2
Global Double ix, iy
Global Long PMoving = 0, PHover = 0
Global Object Win1RT

OpenW Center 1
Set Win1RT = D2GetRT(Me)

' Define the circle coordinates and radius
CenterX = _X / 2 : CenterY = _Y / 2 : Radius = _Y / 2.5

' The segment who can be moved by user
RayX1 = _X / 8 : RayY1 = _Y / 2 : RayX2 = _X / 1.25 : RayY2 = _Y / 2
Do
  DrawAll
  Sleep
Loop Until Me Is Nothing
Function IntersectionSegmentCircle(x1#, y1#, x2#, y2#, cx#, cy#, r#, ByRef ix#, ByRef iy#) As Boolean

  Static Double a, b, c, v, U, lx, ly

  lx = x2 - x1
  ly = y2 - y1

  a = lx * lx + ly * ly

  If a = 0 Return False

  b = 2 * ( lx * (x1 - cx) + ly * (y1 - cy) )

  c = cx * cx + cy * cy + x1 * x1 + y1 * y1 - 2 * (cx * x1 + cy * y1 ) - r * r

  v = b * b - 4 * a * c

  If v >= 0

    v = Sqr(v)
    a = a * 2

    U = (-b - v)  / a

    If U >= 0 And U <= 1

      ix = x1 + U * lx
      iy = y1 + U * ly

      Return True

    EndIf

    U = (-b + v)  / a

    If U >= 0 And U <= 1

      ix = x1 + U * lx
      iy = y1 + U * ly

      Return True

    EndIf

  EndIf

  Return False

EndFunc

Proc DrawAll

  Local Long i

  D2BeginDraw Win1RT, D2C_DarkBlue

  ' Draw the segment
  DrawArrow RayX1, RayY1, RayX2, RayY2, D2C_Yellow

  ' Draw the circle
  D2ForeColor = D2C_White
  D2PCircle CenterX, CenterY, Radius, , , , D2Brush(ARGB(100, 50, 50, 255))

  ' Draw intersection if found
  If IntersectionSegmentCircle(RayX1, RayY1, RayX2, RayY2, CenterX, CenterY, Radius, ix, iy)
    D2PCircle ix, iy, 10, , , , D2Brush(D2C_Blue)
    D2Text ix + 14, iy - 5, "Intersection"
  EndIf

  ' Highlighting hovered points of segment
  Select PHover
  Case 1 : D2PCircle RayX1, RayY1, 5, , , , D2Brush(D2C_Yellow)
  Case 2 : D2PCircle RayX2, RayY2, 5, , , , D2Brush(D2C_Yellow)
  EndSelect

  D2EndDraw

EndProc
Sub Win_1_MouseMove(Button&, Shift&, x!, y!)

  Local Double d1, d2

  If PMoving = 0
    d1 = Sqr((MouseX - RayX1) ^ 2 + (MouseY - RayY1) ^ 2)
    d2 = Sqr((MouseX - RayX2) ^ 2 + (MouseY - RayY2) ^ 2)
    If d1 < d2 And d1 < 15
      PHover = 1
    Else If d2 < d1 And d2 < 15
      PHover = 2
    Else
      PHover = 0
    EndIf
  Else
    Select PMoving
    Case 1 : RayX1 = MouseX : RayY1 = MouseY
    Case 2 : RayX2 = MouseX : RayY2 = MouseY
    EndSelect
  EndIf

EndSub
Sub Win_1_MouseDown(Button&, Shift&, x!, y!)

  PMoving = PHover

EndSub
Sub Win_1_MouseUp(Button&, Shift&, x!, y!)

  PMoving = 0

EndSub
Proc DrawArrow(x1#, y1#, x2#, y2#, C%)

  Local Double lx, ly, px, py, d

  lx = x2 - x1
  ly = y2 - y1
  d = Sqr(lx * lx + ly * ly)

  If d > 0

    lx = 15 * lx / d
    ly = 15 * ly / d
    px = x2 - lx
    py = y2 - ly
    lx = lx / 2.5
    ly = ly / 2.5

    D2Line x1, y1, px, py, D2Brush(C)
    D2Line x2, y2, px + ly , py - lx, D2Brush(C)
    D2Line px + ly , py - lx ,  px - ly , py + lx, D2Brush(C)
    D2Line px - ly , py + lx , x2, y2, D2Brush(C)

  EndIf

EndProc


Intersection Line Sphere

This function return True if line defined by segment (x1,y1,z1)(x2,y2,z2) intersecting sphere (x3,y3,r).
Then (ix,iy,iz) is the intersection and (nx,ny,nz) the normal from surface.

Function IntersectionLineSphere(x1#, y1#, z1#, x2#, y2#, z2#, x3#, y3#, z3#, r#, ByRef ix#, ByRef iy#, ByRef iz#, ByRef nx#, ByRef ny#, ByRef nz#) As Boolean

  Static Double b, c, d, v, U
  Static Double lx, ly, lz
  Static Double iox, ioy, ioz
  Static Double rr
  Static Long InSide

  iox = x1 - x3
  ioy = y1 - y3
  ioz = z1 - z3

  rr = r * r

  lx = x2 - x1
  ly = y2 - y1
  lz = z2 - z1

  d = Sqr(lx * lx + ly * ly + lz * lz)
  If d > 0
    lx /= d
    ly /= d
    lz /= d
  Else If lx <> 0
    lx = Sgn(lx)
    ly = 0
    lz = 0
  Else If ly <> 0
    lx = 0
    ly = Sgn(ly)
    lz = 0
  Else If lz <> 0
    lx = 0
    ly = 0
    lz = Sgn(ly)
  Else
    lx = 1
    ly = 0
    lz = 0
  EndIf

  b = 2 * (lx * iox + ly * ioy + lz * ioz)
  c = x3 ^ 2 + y3 ^ 2 + z3 ^ 2 + x1 ^ 2 + y1 ^ 2 + z1 ^ 2 - 2 * (x3 * x1 + y3 * y1 + z3 * z1 ) - rr
  v = b * b - 4 * c

  Clr ix, iy, iz
  If v >= 0
    v = Sqr(v)
    If iox * iox + ioy * ioy + ioz * ioz < rr
      U = Max(-b - v, -b + v) / 2
      InSide = True
    Else
      U = Min(-b - v, -b + v) / 2
      InSide = False
    EndIf
    If U >= 0
      ix = x1 + U * lx
      iy = y1 + U * ly
      iz = z1 + U * lz
      If InSide = True
        nx = x3 - ix
        ny = y3 - iy
        nz = z3 - iz
      Else
        nx = ix - x3
        ny = iy - y3
        nz = iz - z3
      EndIf
      Normalise3D nx, ny, nz
      Return True
    EndIf
  EndIf

  Return False

EndFunc

Example

$Library "Gdip.inc"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Type LINE3D_STRUCT
  - XYZ_STRUCT p1, p2
EndType
Type SPHERE_STRUCT Packed 1
  - XYZ_STRUCT Center
  - Double Radius
EndType
Global Long X, Y, Aw, Bitmap, AdrBitmap, WBitmap = 320 , HBitmap = 240
Global BIH As BITMAPINFOHEADER
Global Sphere As SPHERE_STRUCT
Global Intersection As XYZ_STRUCT
Global Normal As XYZ_STRUCT
Global SpotLight As XYZ_STRUCT
Global Reflect As XYZ_STRUCT
Global Double wd2 = WBitmap / 2, hd2 = HBitmap / 2, cvv, u, xp, yp, zp, RotationTexture, GroundScrollingX, GroundScrollingZ
Global Double RotationTextureS, RotationTextureC
BIH.biSize = SizeOf(BIH)
BIH.biWidth = WBitmap
BIH.biHeight = HBitmap
BIH.biPlanes = 1
BIH.biBitCount = 32
Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)
Sphere.Center = XYZ(0, -50, 300)
Sphere.Radius = 100
SpotLight = XYZ(500, 4000, -4000)
OpenW Center 1
BackColor = RGB(0, 0, 64)
Do

  RotationTexture = Frac(Timer / 20) * 6.283
  RotationTextureS = Sin(RotationTexture) : RotationTextureC = Cos(RotationTexture)
  GroundScrollingX = Frac(Timer / 5)
  GroundScrollingZ = Frac(Timer / 5)
  Sphere.Center = XYZ(500 * (MouseX - _X / 2) / _X, -50 + Abs(Sin(3.1416 * Frac(Timer / 2)) * 80), 300 + Max(0, _Y - MouseY))

  Aw = AdrBitmap
  For Y = 0 To HBitmap - 1
    For X = 0 To WBitmap - 1
      If IntersectionLineSphere(0, 0, 0, X - wd2, Y - hd2, 300, Sphere, Intersection, Normal)

        cvv = CosVV(Normal.x, Normal.y, Normal.z, SpotLight.x, SpotLight.y, SpotLight.z)
        LPoke Aw, TextureSphere(Intersection.x, Intersection.y, Intersection.z)
        u = ((1 + cvv) / 2) ^ 2
        Poke Aw, u * Peek(Aw)
        Poke Aw + 1, u * Peek(Aw + 1)
        Poke Aw + 2, u * Peek(Aw + 2)
        If cvv > 0.5
          u = 2 * (cvv - 0.5)
          Poke Aw, Min(255, Peek(Aw) + u * 100)
          Poke Aw + 1, Min(255, Peek(Aw + 1) + u * 80)
          Poke Aw + 2, Peek(Aw + 2) + (255 - Peek(Aw + 2)) * u
        EndIf

        Reflect = GetXYZReflect(XYZ(X - wd2, Y - hd2, 300), Normal)

      Else

        LPoke Aw, &h4444ff

        Reflect = XYZ(X - WBitmap / 2, Y - HBitmap / 2, 300)
        Intersection = XYZ(0, 0, 0)

      EndIf

      If Reflect.y < 0.01
        If Intersection.y + Reflect.y <> 0// -149.99
          yp = (-150 - Intersection.y) / (Reflect.y)
          xp = Intersection.x + Reflect.x * yp
          zp = Intersection.z + Reflect.z * yp
          yp = -150
          If Damier(GroundScrollingX + xp / 100, -150, GroundScrollingZ + zp / 100)
            u = 1 - Max(0, 1 - Sqr(zp * zp + xp * xp) / 2000) / 2
            Poke Aw, Peek(Aw) * u
            Poke Aw + 1, Peek(Aw + 1) * u
            Poke Aw + 2, Peek(Aw + 2) * u
          Else
            u = Max(0, 1 - Sqr(zp * zp + xp * xp) / 2000) * 50
            Poke Aw, Min(255, Peek(Aw) + u)
            Poke Aw + 1, Min(255, Peek(Aw + 1) + u)
            Poke Aw + 2, Min(255, Peek(Aw + 2) + u)
          EndIf
          If IntersectionLineSphere(xp, yp, zp,  SpotLight.x ,  SpotLight.y ,  SpotLight.z , Sphere, Intersection, Normal)
            Poke Aw, Peek(Aw) / 2
            Poke Aw + 1, Peek(Aw + 1) / 2
            Poke Aw + 2, Peek(Aw + 2) / 2
          EndIf
        EndIf
      EndIf

      Add Aw, 4
    Next X
  Next Y
  Put _X / 2 - WBitmap / 2, _Y / 2 - HBitmap / 2, Bitmap
  PeekEvent
Loop Until Me Is Nothing
FreeBmp Bitmap
Function TextureSphere(x, y, z) As Long

  x -= Sphere.Center.x
  y -= Sphere.Center.y
  z -= Sphere.Center.z
  x /= 100
  y /= 100
  z /= 100

  If Damier(RotationTextureS * x + RotationTextureC * z, y, RotationTextureS * z - RotationTextureC * x) = 0
    Return &HFF0000
  Else
    Return &HFFFF44
  EndIf

EndFunc
Function IntersectionLineSphere(x1#, y1#, z1#, x2#, y2#, z2#, ByRef S As SPHERE_STRUCT, ByRef Intersection As XYZ_STRUCT, ByRef Normal As XYZ_STRUCT) As Boolean

  Static Double b, c, d, v, U
  Static Double lx, ly, lz
  Static Double iox, ioy, ioz
  Static Double rr
  Static Long InSide

  iox = x1 - S.Center.x
  ioy = y1 - S.Center.y
  ioz = z1 - S.Center.z

  rr = S.Radius * S.Radius

  lx = x2 - x1
  ly = y2 - y1
  lz = z2 - z1

  d = Sqr(lx * lx + ly * ly + lz * lz)
  If d > 0
    lx /= d
    ly /= d
    lz /= d
  Else If lx <> 0
    lx = Sgn(lx)
    ly = 0
    lz = 0
  Else If ly <> 0
    lx = 0
    ly = Sgn(ly)
    lz = 0
  Else If lz <> 0
    lx = 0
    ly = 0
    lz = Sgn(ly)
  Else
    lx = 1
    ly = 0
    lz = 0
  EndIf

  b = 2 * (lx * iox + ly * ioy + lz * ioz)
  c = S.Center.x ^ 2 + S.Center.y ^ 2 + S.Center.z ^ 2 + x1 ^ 2 + y1 ^ 2 + z1 ^ 2 - 2 * (S.Center.x * x1 + S.Center.y * y1 + S.Center.z * z1 ) - rr
  v = b * b - 4 * c

  If v >= 0
    v = Sqr(v)
    If iox * iox + ioy * ioy + ioz * ioz < rr
      U = Max(-b - v, -b + v) / 2
      InSide = True
    Else
      U = Min(-b - v, -b + v) / 2
      InSide = False
    EndIf
    If U >= 0
      Intersection.x = x1 + U * lx
      Intersection.y = y1 + U * ly
      Intersection.z = z1 + U * lz
      If InSide = True
        Normal.x = S.Center.x - Intersection.x
        Normal.y = S.Center.y - Intersection.y
        Normal.z = S.Center.z - Intersection.z
      Else
        Normal.x = Intersection.x - S.Center.x
        Normal.y = Intersection.y - S.Center.y
        Normal.z = Intersection.z - S.Center.z
      EndIf
      XYZ_Normalise Normal
      Return True
    EndIf
  EndIf

  Return False

EndFunc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Function CosVV(v1x#, v1y#, v1z#, v2x#, v2y#, v2z#) As Double

  Local Double d

  d = Sqr(v1x * v1x + v1y * v1y + v1z * v1z)
  If d > 0
    v1x /= d
    v1y /= d
    v1z /= d
  Else If v1x <> 0
    v1x = Sgn(v1x)
    v1y = 0
    v1z = 0
  Else If v1y <> 0
    v1x = 0
    v1y = Sgn(v1y)
    v1z = 0
  Else If v1z <> 0
    v1x = 0
    v1y = 0
    v1z = Sgn(v1y)
  Else
    v1x = 1
    v1y = 0
    v1z = 0
  EndIf

  d = Sqr(v2x * v2x + v2y * v2y + v2z * v2z)
  If d > 0
    v2x /= d
    v2y /= d
    v2z /= d
  Else If v2x <> 0
    v2x = Sgn(v2x)
    v2y = 0
    v2z = 0
  Else If v2y <> 0
    v2x = 0
    v2y = Sgn(v2y)
    v2z = 0
  Else If v2z <> 0
    v2x = 0
    v2y = 0
    v2z = Sgn(v2y)
  Else
    v2x = 1
    v2y = 0
    v2z = 0
  EndIf

  Return (v1x * v2x + v1y * v2y + v1z * v2z)

EndFunc
Function Damier(ix#, iy, iz#) As Long

  Local Long alpha = 1

  If Frac(Abs(ix)) < 0.5 alpha = 1 - alpha
  If Frac(Abs(iy)) < 0.5 alpha = 1 - alpha
  If Frac(Abs(iz)) < 0.5 alpha = 1 - alpha
  If iz < 0 alpha = 1 - alpha
  If iy < 0 alpha = 1 - alpha
  If ix < 0 alpha = 1 - alpha

  Return alpha

EndFunc
Function GetXYZReflect(ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double dp

  XYZ_Normalise Income
  XYZ_Normalise Normale

  dp = 2 * (Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)

  GetXYZReflect = XYZ(Income.x - Normale.x * dp, Income.y - Normale.y * dp, Income.z - Normale.z * dp)

EndFunc


Intersection of a plane and a line

This function return True when L intersect plane equation P then Intesercetion is the instersection and Side = -1 match the negative face or 1 for the positive face.

Function IntersectionLinePlane(ByRef L As LINE3D_STRUCT, ByRef P As PLANE_STRUCT, ByRef Intersection As XYZ_STRUCT, ByRef Side%) As Boolean

  Local Double U, d

  U = (P.a * L.p1.x) + (P.b * L.p1.y) + (P.c * L.p1.z) + P.d

  d = (P.a * (L.p1.x - L.p2.x) + P.b * (L.p1.y - L.p2.y) + P.c * (L.p1.z - L.p2.z))

  If U <= 0
    Side = 1
    If d >= 0 Return False
  Else
    Side = -1
    If d < 0 Return False
  EndIf

  If Abs(d) < 0.00001 Return False
  U = U / d

  Intersection.x = L.p1.x + U * (L.p2.x - L.p1.x)
  Intersection.y = L.p1.y + U * (L.p2.y - L.p1.y)
  Intersection.z = L.p1.z + U * (L.p2.z - L.p1.z)

  Return True

EndFunc

Example

$Library "Gdip.inc"
Type XYZ_STRUCT
  - Double x, y, z
EndType
Type LINE3D_STRUCT
  - XYZ_STRUCT p1, p2
EndType
Type PLANE_STRUCT
  - Double a, b, c, d
EndType
Global Long Faces = 12, i, j, Side(0 .. Faces), X, Y, Aw, Bitmap, AdrBitmap, WBitmap = 320, HBitmap = 240
Global BIH As BITMAPINFOHEADER
Global Plane(0 .. Faces) As PLANE_STRUCT
Global Intersection(0 .. Faces) As XYZ_STRUCT
Global IntersectionOk As XYZ_STRUCT
Global Normal As XYZ_STRUCT
Global SpotLight As XYZ_STRUCT
Global Reflect As XYZ_STRUCT
Global Ray As LINE3D_STRUCT, RefRay As LINE3D_STRUCT
Global Double wd2 = WBitmap / 2, hd2 = HBitmap / 2, cvv, u, xp, yp, zp, Rotation
Global Boolean Ok

BIH.biSize = SizeOf(BIH)
BIH.biWidth = WBitmap
BIH.biHeight = HBitmap
BIH.biPlanes = 1
BIH.biBitCount = 32
Bitmap =  CreateDIBSection(Null, *BIH, DIB_RGB_COLORS, V:AdrBitmap, Null, 0)

OpenW Center 1
BackColor = RGB(0, 0, 64)

' A spot light source
SpotLight = XYZ(500, 3000, -2000)

' Position of eye camera
Ray.p1 = XYZ(0, 0, -190)

' Define the plane of Ground
Plane(0).a = 0, Plane(0).b = 1, Plane(0).c = 0, Plane(0).d = 30

' Define the face planes of a Dodecahedron
Global DodecSize = 20
Plane(1).a = 0, Plane(1).b = -1, Plane(1).c = 0, Plane(1).d = DodecSize
Plane(12).a = 0, Plane(12).b = 1, Plane(12).c = 0, Plane(12).d = DodecSize
For i = 2 To 11
  u = i * 0.6283
  Plane(i).a = Sin(u) * Sqr(0.8)
  Plane(i).b = (-1 + (i % 2) * 2) * Sqr(0.2)
  Plane(i).c = Cos(u) * Sqr(0.8)
  Plane(i).d = DodecSize
Next i

Do

  For i = 1 To Faces
    Turn Plane(i).a, Plane(i).b, 0.0052
    Turn Plane(i).a, Plane(i).c, 0.006
    Turn Plane(i).b, Plane(i).c, 0.0045
  Next i

  Aw = AdrBitmap
  For Y = 0 To HBitmap - 1
    For X = 0 To WBitmap - 1

      Ray.p2 = XYZ(X - wd2, Y - hd2, 300)

      ' Recensing all intersections
      For i = 1 To Faces
        ~IntersectionLinePlane(Ray, Plane(i), Intersection(i), Side(i))
      Next i

      LPoke Aw, 0

      ' Determine if a plane intersection point is on a dodecahedron face
      ' the way is to look if all others faces are in the negative space
      Ok = False
      For i = 1 To Faces
        If Side(i) = 1 ' positive space
          Ok = True
          For j = 1 To Faces
            If j <> i
              If SidePlane(Intersection(i), Plane(j)) = 1 ' if an other is positive then fail
                Ok = False : j = Faces
              EndIf
            EndIf
          Next j
          If Ok
            Normal = XYZ(Plane(i).a, Plane(i).b, Plane(i).c)
            IntersectionOk = Intersection(i)
            ' Colors of the dodecahedrond faces
            LPoke Aw, Choose(1 + (i % 7), &HFF0000, &HFFFF00, &H00FF00, &H0000FF, &HFF00FF, &H00FFFF, &HFFFFFF)
            i = Faces
          EndIf
        EndIf
      Next i

      If Ok ' One face intercepted

        ' Compute a lighting
        cvv = CosVV(Normal.x, Normal.y, Normal.z, SpotLight.x, SpotLight.y, SpotLight.z)
        u = ((1 - cvv) / 2) ^ 2
        Poke Aw, u * Peek(Aw)
        Poke Aw + 1, u * Peek(Aw + 1)
        Poke Aw + 2, u * Peek(Aw + 2)
        If u > 0.8
          u = (u - 0.8) / 0.2
          Poke Aw, Min(255, Peek(Aw) + u * 200)
          Poke Aw + 1, Min(255, Peek(Aw + 1) + u * 200)
          Poke Aw + 2, Min(255, Peek(Aw + 2) + u * 200)
        EndIf

        ' Add little reflection (ground or sky)
        Reflect = GetXYZReflect(XYZ(X - wd2, Y - hd2, 300), Normal)
        RefRay = LINE3D(IntersectionOk, XYZ(IntersectionOk.x + Reflect.x, IntersectionOk.y + Reflect.y, IntersectionOk.z + Reflect.z))
        If IntersectionLinePlane(RefRay, Plane(0), Intersection(0), Side(0))
          If Damier(Intersection(0).x / 100, 0, Intersection(0).z / 100) = 1
            Poke Aw, Peek(Aw) / 2
            Poke Aw + 1, Peek(Aw + 1) / 2
            Poke Aw + 2, Peek(Aw + 2) / 2
          Else
            Poke Aw, Min(255, Peek(Aw + 1) + 50)
            Poke Aw + 1, Min(255, Peek(Aw + 1) + 50)
            Poke Aw + 2, Min(255, Peek(Aw + 2) + 50)
          EndIf
        Else
          Poke Aw , (Peek(Aw) + (255 - (Max(0, Reflect.y ) ^ 0.5) * 255)) / 2
          Poke Aw + 1, Peek(Aw + 1) / 2
          Poke Aw + 2, Peek(Aw + 2) / 2
        EndIf

      Else ' The ground (plane(0)) or the sky

        If IntersectionLinePlane(Ray, Plane(0), Intersection(0), Side(0))
          If Damier(Intersection(0).x / 100, 0, Intersection(0).z / 100) = 1
            LPoke Aw, &h101010
          Else
            LPoke Aw, &h959595
          EndIf
        Else
          Poke Aw , Max(0, Min(255, 255 - Abs(Ray.p1.y - Ray.p2.y)))
        EndIf

      EndIf

      Add Aw, 4

    Next X
  Next Y

  Put _X / 2 - WBitmap / 2, _Y / 2 - HBitmap / 2, Bitmap

  PeekEvent

Loop Until Me Is Nothing

FreeBmp Bitmap
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc
Function LINE3D(ByRef P1 As XYZ_STRUCT, ByRef P2 As XYZ_STRUCT) As LINE3D_STRUCT

  LINE3D.p1 = P1
  LINE3D.p2 = P2

EndFunc
Function XYZ(x As Double, y As Double, z As Double) As XYZ_STRUCT

  XYZ.x = x
  XYZ.y = y
  XYZ.z = z

EndFunc
Function CosVV(v1x#, v1y#, v1z#, v2x#, v2y#, v2z#) As Double

  Local Double d

  d = Sqr(v1x * v1x + v1y * v1y + v1z * v1z)
  If d > 0
    v1x /= d
    v1y /= d
    v1z /= d
  Else If v1x <> 0
    v1x = Sgn(v1x)
    v1y = 0
    v1z = 0
  Else If v1y <> 0
    v1x = 0
    v1y = Sgn(v1y)
    v1z = 0
  Else If v1z <> 0
    v1x = 0
    v1y = 0
    v1z = Sgn(v1y)
  Else
    v1x = 1
    v1y = 0
    v1z = 0
  EndIf

  d = Sqr(v2x * v2x + v2y * v2y + v2z * v2z)
  If d > 0
    v2x /= d
    v2y /= d
    v2z /= d
  Else If v2x <> 0
    v2x = Sgn(v2x)
    v2y = 0
    v2z = 0
  Else If v2y <> 0
    v2x = 0
    v2y = Sgn(v2y)
    v2z = 0
  Else If v2z <> 0
    v2x = 0
    v2y = 0
    v2z = Sgn(v2y)
  Else
    v2x = 1
    v2y = 0
    v2z = 0
  EndIf

  Return (v1x * v2x + v1y * v2y + v1z * v2z)

EndFunc
Function Damier(ix#, iy, iz#) As Long

  Local Long alpha = 1

  ix = ix + iz : iz = iz - (ix - iz) ' 1/4 turn
  If Frac(Abs(ix)) < 0.5 alpha = 1 - alpha
  If Frac(Abs(iy)) < 0.5 alpha = 1 - alpha
  If Frac(Abs(iz)) < 0.5 alpha = 1 - alpha
  If iz < 0 alpha = 1 - alpha
  If iy < 0 alpha = 1 - alpha
  If ix < 0 alpha = 1 - alpha

  Return alpha

EndFunc
Function GetXYZReflect(ByRef Income As XYZ_STRUCT, ByRef Normale As XYZ_STRUCT) As XYZ_STRUCT

  Static Double dp

  XYZ_Normalise Income
  XYZ_Normalise Normale

  dp = 2 * (Normale.x * Income.x + Normale.y * Income.y + Normale.z * Income.z)

  GetXYZReflect = XYZ(Income.x - Normale.x * dp, Income.y - Normale.y * dp, Income.z - Normale.z * dp)

EndFunc
Function IntersectionLinePlane(ByRef L As LINE3D_STRUCT, ByRef P As PLANE_STRUCT, ByRef Intersection As XYZ_STRUCT, ByRef Side%) As Boolean

  Local Double U, d

  U = (P.a * L.p1.x) + (P.b * L.p1.y) + (P.c * L.p1.z) + P.d

  d = (P.a * (L.p1.x - L.p2.x) + P.b * (L.p1.y - L.p2.y) + P.c * (L.p1.z - L.p2.z))

  If U <= 0
    Side = 1
    If d >= 0 Return False
  Else
    Side = -1
    If d < 0 Return False
  EndIf

  If Abs(d) < 0.00001 Return False
  U = U / d

  Intersection.x = L.p1.x + U * (L.p2.x - L.p1.x)
  Intersection.y = L.p1.y + U * (L.p2.y - L.p1.y)
  Intersection.z = L.p1.z + U * (L.p2.z - L.p1.z)

  Return True

EndFunc
Function SidePlane(ByRef Pt As XYZ_STRUCT, ByRef P As PLANE_STRUCT) As Long Naked
  If (P.a * Pt.x) + (P.b * Pt.y) + (P.c * Pt.z) + P.d < 0
    Return 1
  EndIf
  Return -1
EndFunc
Proc Turn(ByRef x#, ByRef y#, a#)

  Static Double s, c, newx, newy

  s = Sin(a)
  c = Cos(a)
  newx = c * x + s * y
  newy = c * y - s * x
  x = newx
  y = newy

EndProc


Distance Point Box

Distance between a point and a rectangle

This function return the distance from a point x,y to a rectangle x1,y1,x2,y2. The point ix,iy is the nearest point.

Function DistanceFromBox(x#, y#, x1#, y1#, x2#, y2#, ByRef ix#, ByRef iy#) As Double

  Local Double dx, dy

  If x1 > x2 Swap x1, x2
  If y1 > y2 Swap y1, y2

  If x > x1
    If x < x2
      If y > y1
        If y < y2

          ' In the box
          dy = y - y1
          iy = y1
          If dy > y2 - y
            dy = y2 - y
            iy = y2
          EndIf

          dx = x - x1
          ix = x1
          If dx > x2 - x
            dx = x2 - x
            ix = x2
          EndIf

          If dx < dy
            iy = y
            Return -dx
          Else
            ix = x
            Return -dy
          EndIf

        EndIf
      EndIf
    EndIf
  EndIf

  If x < x1 And y > y1 And y < y2
    ix = x1
    iy = y
    Return x1 - x
  EndIf

  If x > x2 And y > y1 And y < y2
    ix = x2
    iy = y
    Return x - x2
  EndIf

  If y < y1 And x > x1 And x < x2
    ix = x
    iy = y1
    Return y1 - y
  EndIf

  If y > y2 And x > x1 And x < x2
    ix = x
    iy = y2
    Return y - y2
  EndIf

  If x < x1
    ix = x1
  Else
    ix = x2
  EndIf

  If y < y1
    iy = y1
  Else
    iy = y2
  EndIf

  dx = ix - x
  dy = iy - y

  Return Sqr(dx * dx + dy * dy)

EndFunc


Geodesic distance 2 points on sphere

Geodesic distance between 2 points on a sphere

This function return the geodesic distance of 2 points on the surface of a sphere with radius =1.
It is not necessary to normalise A and B.
Note if A and B are from a sphere with a radius <> 1 you must multiply the result by the radius.
The sphere center is implicitly (0,0,0).

Function SphericalDistance(ByRef A As XYZ_STRUCT, ByRef B As XYZ_STRUCT) As Double

  Local Double C

  C = CosVV(A.x, A.y, A.z, B.x, B.y, B.z)
  C = Max(-1, Min(1, C))

  Return Acos(C) / 2

EndFunc


Distance Point Segment

Distance Point Line Segment

This function return the shortest distance between point (px,py) and the line segment (x1,y1,x2,y2) then (ix,iy) is the nearest point of the line segment.

Function DistancePointSegment(x1#, y1#, x2#, y2#, px#, py#, ByRef ix#, ByRef iy#) As Double

  Local Double  h, U

  h = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
  If h NEAR 0
    ix = (x1 + x2) / 2 : iy = (y1 + y2) / 2
    Return Sqr((ix - px) ^ 2 + (iy - py) ^ 2)
  EndIf
  U = (((px - x1) * (x2 - x1)) + ((py - y1) * (y2 - y1)))
  U = U / (h ^ 2)
  U = Max(0, Min(1, U))
  ix = x1 + U * (x2 - x1)
  iy = y1 + U * (y2 - y1)

  Return Sqr((ix - px) ^ 2 + (iy - py) ^ 2)

EndFunc


Distance Point Circle

Distance between a point and a circle

This function compute the distance from the given point x,y to the circle defined at center cx,cy of radius cr. The nearest point on circle is returned by ix,iy.

Function DistanceFromCircle(x, y, cx, cy, cr, ByRef ix As Double, ByRef iy As Double) As Double

  Local Double lx, ly, d

  lx = x - cx
  ly = y - cy
  d = Sqr(lx * lx + ly * ly)

  If d = 0
    ix = cx + cr
    iy = cy
  Else
    ix = cx + cr * lx / d
    iy = cy + cr * ly / d
  EndIf

  Return d - cr

EndFunc


Distance between point and ellipse

This function return the distance from a given point (Px,Py) to the ellipse of center (Cx,Cy) and radius (Rx,Ry) turned by a degree angle Ao.
Then the nearest point is (ix,iy).
This function uses dichotomous reduction.

Function DistancePointEllipse(Px, Py, Cx, Cy, Rx, Ry, Ao, ByRef ix As Double, ByRef iy As Double) As Double

  Local Long i, q, N
  Local Double a, pa, sao, cao, aStart, aBest, lx, ly, d, dmin, dbest, xt, yt, ie , corde, f, fx, fy, d1, d2

  Px -= Cx
  Py -= Cy

  corde = 2 * Max(Rx, Ry)
  Ao = Ao * PI / 180
  sao = Sin(Ao)
  cao = Cos(Ao)

  If Rx > Ry
    f = Sqr(Rx * Rx - Ry * Ry)
    fx = cao * f
    fy = sao * f
  Else
    f = Sqr(Ry * Ry - Rx * Rx)
    fx = -sao * f
    fy = cao * f
  EndIf

  lx = Px - fx  : ly = Py - fy  : d1 = Sqr(lx * lx + ly * ly)
  lx = Px + fx  : ly = Py + fy  : d2 = Sqr(lx * lx + ly * ly)

  If d1 + d2 < corde
    ie = -1
  Else
    ie = 1
  EndIf

  xt =  cao * Px + sao * Py
  yt =  cao * Py - sao * Px
  Px = xt
  Py = yt

  dbest = _maxDbl
  a = 0

  For q = 0 To 1

    pa =  PI / 3
    dmin = _maxDbl
    For N = 0 To 50
      For i = 0 To 2
        a = a + pa
        lx = Sin(a) * Rx - Px
        ly = Cos(a) * Ry - Py
        d = lx * lx + ly * ly
        If d < dmin
          dmin = d
          aBest = a
        EndIf
      Next i

      pa = pa / 2
      a = aBest - pa * 1.5
    Next N

    a = PI

    If dmin < dbest
      ix = Sin(aBest) * Rx
      iy = Cos(aBest) * Ry
      dbest = dmin
    EndIf

  Next q

  xt =  cao * ix - sao * iy
  yt =  cao * iy + sao * ix
  ix = Cx + xt
  iy = Cy + yt

  Return ie * Sqr(dbest)

EndFunc


Example for all distances functions

Example with all distance functions

$Library "Direct2D"
Type CIRCLE_STRUCT Packed 1
  - Single Cx, Cy, Cr
EndType
Type ELLIPSE_STRUCT Packed 1
  - Single Cx, Cy, Rx, Ry, Ao
EndType
Type RECTANGLE_STRUCT Packed 1
  - Single x1, y1, x2, y2
EndType
Type XYZ_STRUCT Packed 1
  - Double x, y, z
EndType
Type SPHERE_STRUCT Packed 1
  - XYZ_STRUCT Center, PointA, PointB
  - Double Radius
EndType
Global Object Win1RT, Geo, GeoSink
Global Dim DemoRectangle(0) As RECTANGLE_STRUCT
Global Dim DemoCircle(0) As CIRCLE_STRUCT
Global Dim DemoEllipse(0) As ELLIPSE_STRUCT
Global Dim DemoSphere(0) As SPHERE_STRUCT
OpenW Center 1
Set Win1RT = D2GetRT(Me)
Array DemoRectangle() = Mks(_X / 8 - _X / 16, _Y / 2 - 130, _X / 8 + _X / 16, _Y / 2 + 130)
Array DemoEllipse() = Mks(_X / 2.66, _Y / 2, 100, 50, 20)
Array DemoCircle() = Mks(_X / 1.6, _Y / 2, 100)
Array DemoSphere() = Mkd(_X / 1.14, _Y / 2, 0, 0, 0, 0,  - Sqr(1 / 3) * 120, Sqr(1 / 3) * 120, Sqr(1 / 3) * 120, 120)
Do
  DrawAll
  Sleep
Loop Until Me Is Nothing
Proc DrawAll

  Static A As D2D1_ARC_SEGMENT
  Local Double Distance, ix, iy, RotationY, RotationXZ
  Local Long Mx, My

  Mx = MouseX : My = MouseY

  D2BeginDraw Win1RT, D2C_DarkBlue

  D2ForeColor = D2C_White
  D2TransformSave
  D2Transform D2TF_ROTATE, DemoEllipse(0).Ao, DemoEllipse(0).Cx, DemoEllipse(0).Cy
  D2Ellipse DemoEllipse(0).Cx, DemoEllipse(0).Cy, DemoEllipse(0).Rx, DemoEllipse(0).Ry
  D2TransformRestore

  D2Box DemoRectangle(0).x1, DemoRectangle(0).y1, DemoRectangle(0).x2, DemoRectangle(0).y2
  D2Circle DemoCircle(0).Cx, DemoCircle(0).Cy, DemoCircle(0).Cr

  DemoSphere(0).PointA.x = Mx - DemoSphere(0).Center.x
  DemoSphere(0).PointA.y = My - DemoSphere(0).Center.y
  DemoSphere(0).PointA.z = DemoSphere(0).Radius
  XYZ_Normalise DemoSphere(0).PointA
  DemoSphere(0).PointA.x *= DemoSphere(0).Radius
  DemoSphere(0).PointA.y *= DemoSphere(0).Radius
  DemoSphere(0).PointA.z *= DemoSphere(0).Radius
  DrawSphere DemoSphere(0)

  D2DefLine 2, 2
  D2SetFont "Arial", 20
  Distance = DistanceFromBox(Mx, My, DemoRectangle(0).x1, DemoRectangle(0).y1, DemoRectangle(0).x2, DemoRectangle(0).y2, ix, iy)
  D2Line Mx, My, ix, iy, D2Brush(ARGB(240, 255, 220, 50))
  D2DrawText DemoRectangle(0).x1, DemoRectangle(0).y1, DemoRectangle(0).x2, DemoRectangle(0).y2, Round(Distance, 2), D2T_CENTER | D2T_PARCENTER

  Distance = DistancePointEllipse(Mx, My, DemoEllipse(0).Cx, DemoEllipse(0).Cy, DemoEllipse(0).Rx, DemoEllipse(0).Ry, DemoEllipse(0).Ao, ix, iy)
  D2Line Mx, My, ix, iy, D2Brush(ARGB(240, 255, 220, 50))
  D2DrawText DemoEllipse(0).Cx - 50, DemoEllipse(0).Cy - 50, DemoEllipse(0).Cx + 50, DemoEllipse(0).Cy + 50, Round(Distance, 2) , D2T_CENTER | D2T_PARCENTER

  Distance = DistanceFromCircle(Mx, My, DemoCircle(0).Cx, DemoCircle(0).Cy, DemoCircle(0).Cr, ix, iy)
  D2Line Mx, My, ix, iy, D2Brush(ARGB(240, 255, 220, 50))
  D2DrawText DemoCircle(0).Cx - 50, DemoCircle(0).Cy - 50, DemoCircle(0).Cx + 50, DemoCircle(0).Cy + 50, Round(Distance, 2) , D2T_CENTER | D2T_PARCENTER

  Distance = SphericalDistance(DemoSphere(0).PointA, DemoSphere(0).PointB) * DemoSphere(0).Radius
  D2DrawText DemoSphere(0).Center.x - 50, DemoSphere(0).Center.y - 50, DemoSphere(0).Center.x + 50, DemoSphere(0).Center.y + 50, Round(Distance, 2) , D2T_CENTER | D2T_PARCENTER

  D2EndDraw

EndProc
Function DistanceFromBox(x#, y#, x1#, y1#, x2#, y2#, ByRef ix#, ByRef iy#) As Double

  Local Double dx, dy

  If x1 > x2 Swap x1, x2
  If y1 > y2 Swap y1, y2

  If x > x1
    If x < x2
      If y > y1
        If y < y2

          ' In the box
          dy = y - y1
          iy = y1
          If dy > y2 - y
            dy = y2 - y
            iy = y2
          EndIf

          dx = x - x1
          ix = x1
          If dx > x2 - x
            dx = x2 - x
            ix = x2
          EndIf

          If dx < dy
            iy = y
            Return -dx
          Else
            ix = x
            Return -dy
          EndIf

        EndIf
      EndIf
    EndIf
  EndIf

  If x <= x1 And y > y1 And y < y2
    ix = x1
    iy = y
    Return x1 - x
  EndIf

  If x >= x2 And y > y1 And y < y2
    ix = x2
    iy = y
    Return x - x2
  EndIf

  If y <= y1 And x > x1 And x < x2
    ix = x
    iy = y1
    Return y1 - y
  EndIf

  If y >= y2 And x > x1 And x < x2
    ix = x
    iy = y2
    Return y - y2
  EndIf

  If x < x1
    ix = x1
  Else
    ix = x2
  EndIf

  If y < y1
    iy = y1
  Else
    iy = y2
  EndIf

  dx = ix - x
  dy = iy - y

  Return Sqr(dx * dx + dy * dy)

EndFunc
Function SphericalDistance(ByRef A As XYZ_STRUCT, ByRef B As XYZ_STRUCT) As Double

  Local Double C

  C = CosVV(A.x, A.y, A.z, B.x, B.y, B.z)
  C = Max(-1, Min(1, C))

  Return Acos(C) / 2

EndFunc
Function DistancePointSegment(x1#, y1#, x2#, y2#, px#, py#, ByRef ix#, ByRef iy#) As Double

  Local Double  h, U

  h = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
  If h NEAR 0
    ix = (x1 + x2) / 2 : iy = (y1 + y2) / 2
    Return Sqr((ix - px) ^ 2 + (iy - py) ^ 2)
  EndIf
  U = (((px - x1) * (x2 - x1)) + ((py - y1) * (y2 - y1)))
  U = U / (h ^ 2)
  U = Max(0, Min(1, U))
  ix = x1 + U * (x2 - x1)
  iy = y1 + U * (y2 - y1)

  Return Sqr((ix - px) ^ 2 + (iy - py) ^ 2)

EndFunc
Function DistanceFromCircle(x, y, cx, cy, cr, ByRef ix As Double, ByRef iy As Double) As Double

  Local Double lx, ly, d

  lx = x - cx
  ly = y - cy
  d = Sqr(lx * lx + ly * ly)

  If d = 0
    ix = cx + cr
    iy = cy
  Else
    ix = cx + cr * lx / d
    iy = cy + cr * ly / d
  EndIf

  Return d - cr

EndFunc
Function DistancePointEllipse(Px, Py, Cx, Cy, Rx, Ry, Ao, ByRef ix As Double, ByRef iy As Double) As Double

  Local Long i, q, N
  Local Double a, pa, sao, cao, aStart, aBest, lx, ly, d, dmin, dbest, xt, yt, ie , corde, f, fx, fy, d1, d2

  Px -= Cx
  Py -= Cy

  corde = 2 * Max(Rx, Ry)
  Ao = Ao * PI / 180
  sao = Sin(Ao)
  cao = Cos(Ao)

  If Rx > Ry
    f = Sqr(Rx * Rx - Ry * Ry)
    fx = cao * f
    fy = sao * f
  Else
    f = Sqr(Ry * Ry - Rx * Rx)
    fx = -sao * f
    fy = cao * f
  EndIf

  lx = Px - fx  : ly = Py - fy  : d1 = Sqr(lx * lx + ly * ly)
  lx = Px + fx  : ly = Py + fy  : d2 = Sqr(lx * lx + ly * ly)

  If d1 + d2 < corde
    ie = -1
  Else
    ie = 1
  EndIf

  xt =  cao * Px + sao * Py
  yt =  cao * Py - sao * Px
  Px = xt
  Py = yt

  dbest = _maxDbl
  a = 0

  For q = 0 To 1

    pa =  PI / 3
    dmin = _maxDbl
    For N = 0 To 50
      For i = 0 To 2
        a = a + pa
        lx = Sin(a) * Rx - Px
        ly = Cos(a) * Ry - Py
        d = lx * lx + ly * ly
        If d < dmin
          dmin = d
          aBest = a
        EndIf
      Next i

      pa = pa / 2
      a = aBest - pa * 1.5
    Next N

    a = PI

    If dmin < dbest
      ix = Sin(aBest) * Rx
      iy = Cos(aBest) * Ry
      dbest = dmin
    EndIf

  Next q

  xt =  cao * ix - sao * iy
  yt =  cao * iy + sao * ix
  ix = Cx + xt
  iy = Cy + yt

  Return ie * Sqr(dbest)

EndFunc
Function CosVV(v1x#, v1y#, v1z#, v2x#, v2y#, v2z#) As Double

  Local Double d

  d = Sqr(v1x * v1x + v1y * v1y + v1z * v1z)
  If d > 0
    v1x /= d
    v1y /= d
    v1z /= d
  Else If v1x <> 0
    v1x = Sgn(v1x)
    v1y = 0
    v1z = 0
  Else If v1y <> 0
    v1x = 0
    v1y = Sgn(v1y)
    v1z = 0
  Else If v1z <> 0
    v1x = 0
    v1y = 0
    v1z = Sgn(v1y)
  Else
    v1x = 1
    v1y = 0
    v1z = 0
  EndIf

  d = Sqr(v2x * v2x + v2y * v2y + v2z * v2z)
  If d > 0
    v2x /= d
    v2y /= d
    v2z /= d
  Else If v2x <> 0
    v2x = Sgn(v2x)
    v2y = 0
    v2z = 0
  Else If v2y <> 0
    v2x = 0
    v2y = Sgn(v2y)
    v2z = 0
  Else If v2z <> 0
    v2x = 0
    v2y = 0
    v2z = Sgn(v2y)
  Else
    v2x = 1
    v2y = 0
    v2z = 0
  EndIf

  Return (v1x * v2x + v1y * v2y + v1z * v2z)

EndFunc
Proc DrawSphere(ByRef S As SPHERE_STRUCT)

  Local Double m, x, y, z, px, py, pz, nx, ny, nz, nd
  Local Long i
  Local GS(1 .. 3) As D2Gradient

  GS(1).position = 0.0
  GS(1).argbcolor = D2C_White
  GS(2).position = 0.2
  GS(2).argbcolor = D2C_Green
  GS(3).position = 1.0
  GS(3).argbcolor = D2C_Black

  D2PCircle S.Center.x, S.Center.y, S.Radius, , , , D2BrushGradient(1, GS(), 3, S.Center.x, S.Center.y, S.Radius / 2, -S.Radius / 2, 0.75, S.Radius * 1.5, S.Radius * 1.5)

  x = S.PointA.x : y = S.PointA.y : z = S.PointA.z
  px = (S.PointB.x - x) / 200 : py = (S.PointB.y - y) / 200 : pz = (S.PointB.z - z) / 200
  For i = 0 To 200
    x += px : y += py : z += pz
    nx = x : ny = y : nz = nz
    nd = Sqr(x * x + y * y + z * z)
    If nd > 0
      nx /= nd
      ny /= nd
      nz /= nd
    EndIf
    D2Plot S.Center.x + nx * S.Radius, S.Center.y + ny * S.Radius, D2Brush(ARGB(240, 255, 220, 50))
  Next i

  m = 1 + S.PointA.z / S.Radius
  D2PCircle S.Center.x + S.PointA.x, S.Center.y + S.PointA.y, 5 + m * 3, , , , D2Brush(ARGB(200, 255, 255, 0))
  m = 1 + S.PointB.z / S.Radius
  D2PCircle S.Center.x + S.PointB.x, S.Center.y + S.PointB.y, 5 + m * 3, , , , D2Brush(ARGB(200, 0, 255, 0))

EndProc
Proc Turn(ByRef x#, ByRef y#, a#)

  Static Double s, c, newx, newy

  s = Sin(a)
  c = Cos(a)
  newx = c * x + s * y
  newy = c * y - s * x
  x = newx
  y = newy

EndProc
Proc XYZ_Normalise(ByRef a As XYZ_STRUCT)

  Static Double d

  d = Sqr(a.x * a.x + a.y * a.y + a.z * a.z)

  If d > 0
    a.x /= d
    a.y /= d
    a.z /= d
  Else If a.x < 0
    a.x = -1
    a.y = 0
    a.z = 0
  Else If a.x > 0
    a.x = 1
    a.y = 0
    a.z = 0
  Else If a.y < 0
    a.x = 0
    a.y = -1
    a.z = 0
  Else If a.y > 0
    a.x = 0
    a.y = 1
    a.z = 0
  Else If a.z < 0
    a.x = 0
    a.y = 0
    a.z = -1
  Else
    a.x = 0
    a.y = 0
    a.z = 1
  EndIf

EndProc