@James D Jarvis here is another version but has a problem with the seam, maybe if the image were lengthened with a mirror image?
Code: (Select All)
_Title "Contour Plot 3S: spacebar for New World, escape to quit" 'b+ trans from SdlBasic to QB64, original from QB64.net
' 2021-11-17 Contour Plot 3: planets
' 2021-11-18 3S try adding shade
' Contour plot using Data Points by Mrwhy
'ref: SdlBasic forum, AndyA, 9-13-2016 http://sdlbasic.epizy.com/showthread.php?tid=302
' AndyA ref:
'http://www.qb64.net/forum/index.php?topic=3714.msg37019#msg37019
'https://en.wikipedia.org/wiki/Richard_V._Southwell
'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
Dim Shared sw, sh, nC, nP
sw = 800: sh = 400: nC = 32: nP = 12 ' screen Width and Height, number of Colors, number of Points
Screen _NewImage(sw, sh, 32)
_ScreenMove 200, 60 ' at least 60 on each side and to and bottom
_PrintMode _KeepBackground
ReDim Shared pal(nC) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
'color mixing green blue
For i = 1 To nC
pal(i) = _RGB32(0, i * 255 / nC, 128 / nC * (nC - i))
Next
restart:
ReDim h(sw, sh)
Cls
For i = 1 To nP ' new data points and color assignments
x(i) = Int((sw - 30) * Rnd + 15): y(i) = Int((sh - 30) * Rnd + 15): z(i) = Int(Rnd * (nC - 2) + 1)
Next
ztot = 0
'display height or potential data points on screen
For i = 1 To nP
fcirc x(i), sh - y(i), 3, pal(z(i))
ztot = ztot + z(i)
Next
Color &HFFFFFF40
_Font 8
s$ = "Calculating Contour Map"
x = (_Width - _PrintWidth(s$)) / 2
_PrintString (x, _Height - 14), s$
Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
'initialize some variables
zmean = ztot / nP
wo = sw * sh / nP
w = wo / 2
'generate initial Error estimates
For i = 1 To nP
E(i) = z(i) - zmean
Next
'Legend
'begin Relaxation (iterative method)
For jj = 1 To 9 * nP 'find max error point
emax = 0
For i = 1 To nP
If Abs(E(i)) > emax Then
emax = Abs(E(i))
ii = i
End If
k = E(ii)
Next
'fixit
For i = 1 To nP
dx = x(i) - x(ii)
dy = y(i) - y(ii)
dsq = dx * dx + dy * dy
E(i) = E(i) - k * Exp(-(dsq / w))
If i = ii Then
'update map with revised height or potential estimates for each pixel
For fy = 1 To (sh - 1)
For fx = 1 To (sw - 1)
dx = fx - x(ii)
dy = fy - y(ii)
dsq2 = dx * dx + dy * dy
dy = sh - fy
h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
Next
Next
End If
Next
Next
'Draw calculated contour map
Color pal(h(1, 1) + zmean)
Line (0, 0)-(sw - 1, 0)
Line (0, 0)-(0, sh - 1)
For fy = 1 To sh - 1
For fx = 1 To sw - 1
If Int(h(fx, fy) + zmean) > nC Then
c~& = pal(nC)
ElseIf Int(h(fx, fy) + zmean) < 1 Then
c~& = pal(1)
Else
c~& = pal(Int(h(fx, fy) + zmean))
End If
PSet (fx, fy), c~&
Next
Next
'display height or potential data points on contour map
'For i = 1 To nP
' 'fcirc x(i), sh - y(i), 3, &HFF000000
' fcirc x(i), sh - y(i), 2, pal(z(i))
'Next
'Legend
surface& = _NewImage(_Width, _Height, 32)
mid = Int(sw / 2)
_PutImage (mid - 2, 0)-(sw - 1, sh), 0, surface&
_PutImage (0, 0)-(mid + 2, sh), 0, surface&, (sw - 1, 0)-(0, sh)
_PutImage , surface&, 0
For y = 0 To sh ' fix seam !
PSet (mid + 2, y), Point(mid - 1, y)
PSet (mid + 1, y), Point(mid - 1, y)
PSet (mid, y), Point(mid - 1, y)
Next
_PutImage , 0, surface&
Cls
'_PutImage , surface&, 0 'check
Color &HFF000000 ' be rid of white points at poles
r = sh / 2.5: xc = sw / 2: yc = sh / 2: xo = 0: sr = 200
Do
For y = -r To r
x1 = Sqr(r * r - y * y)
tv = (_Asin(y / r) + 1.5) / 3
For x = -x1 + 1 To x1
tu = (_Asin(x / x1) + 1.5) / 6
_Source surface&
pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
_Dest 0
dist = (x * x + y * y) / (r * r) ' for shading
PSet (x + xc, y + yc), pc~&
PSet (x + xc, y + yc), _RGBA32(0, 0, 0, 180 * (dist ^ 2)) ' for shading
Next x
Next y
xo = xo - 1 + sw
xo = xo Mod sw
If (nC + 1) * 12 <= sh Then Legend
_Display
_Limit 30
Loop Until _KeyDown(32) Or _KeyDown(27)
_FreeImage surface&
If _KeyDown(32) Then GoTo restart
'show height or potential value for each color in map
Sub Legend
posy = 12
' just draw a balck box where legend is going
Line (0, 0)-(37, (nC + 1) * 12), &HFFFF0000, BF
Line (1, 1)-(36, (nC + 1) * 12 - 2), &HFFFFFF00, B
Color &HFFFFFF00
_Font 8
_PrintString (7, 3), "Key"
For nn = 1 To nC
_PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
posy = posy + 12
Next
End Sub
Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
Dim Radius As Long, RadiusError As Long
Dim X As Long, Y As Long
Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
If Radius = 0 Then PSet (CX, CY), C: Exit Sub
Line (CX - X, CY)-(CX + X, CY), C, BF
While X > Y
RadiusError = RadiusError + Y * 2 + 1
If RadiusError >= 0 Then
If X <> Y + 1 Then
Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
End If
X = X - 1
RadiusError = RadiusError - X * 2
End If
Y = Y + 1
Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
Wend
End Sub
b = b + ...