Author Topic: Contour Plot  (Read 6490 times)

0 Members and 1 Guest are viewing this topic.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Contour Plot
« on: November 16, 2021, 08:52:56 pm »
Ah, found some Contour plot code from SdlBasic which turns out started at Q B 6 4 .net [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there] by Mrwhy.

Here's how I found it SdlBasic Forum:
Code: [Select]
        ' Contour plot using Data Points by Mrwhy
       'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302

'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
'https://en.wikipedia.org/wiki/Richard_V._Southwell
'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)

Option QBASIC
Common sw, sh
'============================
' change size of output screen here
'============================
smallPlot = 1
'============================
If smallPlot = 1 Then
    'take about 25 secs to calc
    sw = 240 : sh = 240
Else
    'large plot takes over 2 minutes to calc
    sw = 640 : sh = 480
End If

SetDisplay(sw, sh, 32, 1)
SetCaption("Contour Map")

Dim pal[23], x[50], y[50], z[100], E[50], h[sw, sh]

'red to yellow palette
pal[0]  = 0        : pal[1]  = 0x800000
pal[2]  = 0xA00000 : pal[3]  = 0xC00000
pal[4]  = 0xD00000 : pal[5]  = 0xFF0000
pal[6]  = 0xFF2000 : pal[7]  = 0xFF4000
pal[8]  = 0xFF6000 : pal[9]  = 0xFF8000
pal[10] = 0xFFA000 : pal[11] = 0xFFC000
pal[12] = 0xFFE000 : pal[13] = 0xFFFF00
pal[14] = 0xFFFF20 : pal[15] = 0xFFFF40

'============================================
' number of height/potential points in data
'============================================
numPoints = 16
'============================================

If smallPlot = 1 Then
    x[1]  = 161 : y[1]  = 120 : z[1]  = 11
    x[2]  =  80 : y[2]  =  80 : z[2]  =  9
    x[3]  =  40 : y[3]  =  40 : z[3]  =  6
    x[4]  = 122 : y[4]  =  20 : z[4]  =  3
    x[5]  = 117 : y[5]  =  72 : z[5]  =  4
    x[6]  = 178 : y[6]  = 124 : z[6]  =  7
    x[7]  = 140 : y[7]  = 173 : z[7]  =  9
    x[8]  = 194 : y[8]  =  38 : z[8]  = 13
    x[9]  =  65 : y[9]  =  65 : z[9]  = 12
    x[10] = 135 : y[10] =  60 : z[10] =  9
    x[11] = 180 : y[11] = 160 : z[11] =  8
    x[12] =  70 : y[12] = 135 : z[12] = 11
    x[13] = 165 : y[13] = 150 : z[13] =  6
    x[14] = 225 : y[14] =  90 : z[14] =  2
    x[15] = 190 : y[15] = 178 : z[15] = 12
    x[16] = 135 : y[16] = 135 : z[16] = 10
Else
    x[1]  = 279: y[1]  = 220: z[1]  = 11    '1
    x[2]  = 160: y[2]  = 160: z[2]  =  9    '2
    x[3]  =  80: y[3]  =  80: z[3]  =  6    '3
    x[4]  = 144: y[4]  =  40: z[4]  =  7    '4
    x[5]  = 350: y[5]  = 158: z[5]  =  4    '5
    x[6]  = 356: y[6]  = 248: z[6]  =  7    '6
    x[7]  = 280: y[7]  = 347: z[7]  =  9    '7
    x[8]  = 370: y[8]  =  39: z[8]  = 13    '8
    x[9]  = 130: y[9]  = 230: z[9]  = 12    '9
    x[10] = 270: y[10] = 120: z[10] =  9    '10
    x[11] = 360: y[11] = 320: z[11] =  8    '11
    x[12] = 140: y[12] = 270: z[12] = 11    '12
    x[13] = 530: y[13] = 300: z[13] =  6    '13
    x[14] = 458: y[14] = 143: z[14] =  2    '14
    x[15] = 380: y[15] = 370: z[15] = 12    '15
    x[16] = 270: y[16] = 270: z[16] = 10    '16
End If

ztot = 0
'display height or potential data points on screen
For i = 1 To numPoints
    Ink (pal[z[i]])
   FillCircle (x[i], sh - y[i], 3)
   ztot = ztot + z[i]
Next

Ink(0xFFFF40)
Text(48,2,12,"Calculating contour map")
Box(40, 1, 212, 16)
WaitVBL

'initialize some variables
zmean = ztot / numPoints
wo = sw * sh / numPoints
w = wo / 2.0

'generate initial Error estimates
For i = 1 To numPoints
   E[i] = z[i] - zmean
Next
Legend()
WaitVBL

'begin Relaxation (iterative method)
For jj = 1 To 9 * numPoints 'find max error point
   emax = 0
   For i = 1 To numPoints
       If Abs(E[i]) > emax Then
            emax = Abs(E[i])
            ii = i
        End If
       k = E[ii]
   Next
   'fixit
   For i = 1 To numPoints
       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
Ink(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
       clr = pal[h[fx, fy] + zmean]
       Plot (fx, fy, clr)
   Next
Next

'display height or potential data points on contour map
For i = 1 To numPoints
    clr = pal[z[i]]
    Ink(0)
   FillCircle (x[i], sh - y[i], 3)
    Ink(clr)
    FillCircle (x[i], sh - y[i], 2)
Next
Legend()

Waitkey(27)
End

'show height or potential value for each color in map
Sub Legend()
Ink(0xFFFF40)
Text(5, 45, 10, "Key")
posy = 58
For nn = 1 To 15
    Ink(0xFFFF40)
    Text(0, posy, 10, Str$(nn))
    Ink(0)
   Bar(14, posy + 3, 30, posy + 9)
    Ink(pal[nn])
   Bar (16, posy + 5, 28, posy + 7)
   posy = posy + 12
Next
End Sub


And my QB64 translation, mostly changing the Legend:
Code: QB64: [Select]
  1. _Title "Contour Map" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2.  
  3. ' Contour plot using Data Points by Mrwhy
  4. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  5. ' AndyA ref:
  6. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  7. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  8. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  9.  
  10. Dim Shared sw, sh
  11. '============================
  12. ' change size of output screen here
  13. '============================
  14. smallPlot = 0
  15. '============================
  16. If smallPlot = 1 Then
  17.     'take about 25 secs to calc
  18.     sw = 240: sh = 240
  19.     'large plot takes over 2 minutes to calc
  20.     sw = 640: sh = 480
  21.  
  22. Screen _NewImage(sw, sh, 32)
  23. _ScreenMove 300, 100
  24.  
  25. Dim Shared pal(23) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  26.  
  27. 'red to yellow palette
  28. pal(0) = &HFF000000: pal(1) = &HFF800000
  29. pal(2) = &HFFA00000: pal(3) = &HFFC00000
  30. pal(4) = &HFFD00000: pal(5) = &HFFFF0000
  31. pal(6) = &HFFFF2000: pal(7) = &HFFFF4000
  32. pal(8) = &HFFFF6000: pal(9) = &HFFFF8000
  33. pal(10) = &HFFFFA000: pal(11) = &HFFFFC000
  34. pal(12) = &HFFFFE000: pal(13) = &HFFFFFF00
  35. pal(14) = &HFFFFFF20: pal(15) = &HFFFFFF40
  36.  
  37. '============================================
  38. ' number of height/potential points in data
  39. '============================================
  40. numPoints = 16
  41. '============================================
  42.  
  43. If smallPlot = 1 Then
  44.     x(1) = 161: y(1) = 120: z(1) = 11
  45.     x(2) = 80: y(2) = 80: z(2) = 9
  46.     x(3) = 40: y(3) = 40: z(3) = 6
  47.     x(4) = 122: y(4) = 20: z(4) = 3
  48.     x(5) = 117: y(5) = 72: z(5) = 4
  49.     x(6) = 178: y(6) = 124: z(6) = 7
  50.     x(7) = 140: y(7) = 173: z(7) = 9
  51.     x(8) = 194: y(8) = 38: z(8) = 13
  52.     x(9) = 65: y(9) = 65: z(9) = 12
  53.     x(10) = 135: y(10) = 60: z(10) = 9
  54.     x(11) = 180: y(11) = 160: z(11) = 8
  55.     x(12) = 70: y(12) = 135: z(12) = 11
  56.     x(13) = 165: y(13) = 150: z(13) = 6
  57.     x(14) = 225: y(14) = 90: z(14) = 2
  58.     x(15) = 190: y(15) = 178: z(15) = 12
  59.     x(16) = 135: y(16) = 135: z(16) = 10
  60.     x(1) = 279: y(1) = 220: z(1) = 11 '1
  61.     x(2) = 160: y(2) = 160: z(2) = 9 '2
  62.     x(3) = 80: y(3) = 80: z(3) = 6 '3
  63.     x(4) = 144: y(4) = 40: z(4) = 7 '4
  64.     x(5) = 350: y(5) = 158: z(5) = 4 '5
  65.     x(6) = 356: y(6) = 248: z(6) = 7 '6
  66.     x(7) = 280: y(7) = 347: z(7) = 9 '7
  67.     x(8) = 370: y(8) = 39: z(8) = 13 '8
  68.     x(9) = 130: y(9) = 230: z(9) = 12 '9
  69.     x(10) = 270: y(10) = 120: z(10) = 9 '10
  70.     x(11) = 360: y(11) = 320: z(11) = 8 '11
  71.     x(12) = 140: y(12) = 270: z(12) = 11 '12
  72.     x(13) = 530: y(13) = 300: z(13) = 6 '13
  73.     x(14) = 458: y(14) = 143: z(14) = 2 '14
  74.     x(15) = 380: y(15) = 370: z(15) = 12 '15
  75.     x(16) = 270: y(16) = 270: z(16) = 10 '16
  76.  
  77. ztot = 0
  78. 'display height or potential data points on screen
  79. For i = 1 To numPoints
  80.     fcirc x(i), sh - y(i), 3, pal(z(i))
  81.     ztot = ztot + z(i)
  82.  
  83. Color &HFFFFFF40
  84.  
  85. s$ = "Calculating Contour Map"
  86. x = (_Width - _PrintWidth(s$)) / 2
  87. _PrintString (x, _Height - 14), s$
  88. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  89.  
  90. 'initialize some variables
  91. zmean = ztot / numPoints
  92. wo = sw * sh / numPoints
  93. w = wo / 2.0
  94.  
  95. 'generate initial Error estimates
  96. For i = 1 To numPoints
  97.     E(i) = z(i) - zmean
  98. Legend
  99.  
  100. 'begin Relaxation (iterative method)
  101. For jj = 1 To 9 * numPoints 'find max error point
  102.     emax = 0
  103.     For i = 1 To numPoints
  104.         If Abs(E(i)) > emax Then
  105.             emax = Abs(E(i))
  106.             ii = i
  107.         End If
  108.         k = E(ii)
  109.     Next
  110.     'fixit
  111.     For i = 1 To numPoints
  112.         dx = x(i) - x(ii)
  113.         dy = y(i) - y(ii)
  114.         dsq = dx * dx + dy * dy
  115.         E(i) = E(i) - k * Exp(-(dsq / w))
  116.         If i = ii Then
  117.             'update map with revised height or potential estimates for each pixel
  118.             For fy = 1 To (sh - 1)
  119.                 For fx = 1 To (sw - 1)
  120.                     dx = fx - x(ii)
  121.                     dy = fy - y(ii)
  122.                     dsq2 = dx * dx + dy * dy
  123.                     dy = sh - fy
  124.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  125.                 Next
  126.             Next
  127.         End If
  128.     Next
  129.  
  130. 'Draw calculated contour map
  131. Color pal(h(1, 1) + zmean)
  132. Line (0, 0)-(sw - 1, 0)
  133. Line (0, 0)-(0, sh - 1)
  134. For fy = 1 To sh - 1
  135.     For fx = 1 To sw - 1
  136.         PSet (fx, fy), pal(h(fx, fy) + zmean)
  137.     Next
  138.  
  139. 'display height or potential data points on contour map
  140. For i = 1 To numPoints
  141.     fcirc x(i), sh - y(i), 3, pal(0)
  142.     fcirc x(i), sh - y(i), 2, pal(z(i))
  143. Legend
  144.  
  145. 'show height or potential value for each color in map
  146. Sub Legend
  147.     posy = 20
  148.     ' just draw a balck box where legend is going
  149.     Line (0, 0)-(37, 200), pal(0), BF
  150.     Line (1, 1)-(36, 199), &HFFFFFF40, B
  151.     Color &HFFFFFF40
  152.     _Font 8
  153.     _PrintString (7, 3), "Key"
  154.     For nn = 1 To 15
  155.         Color &HFFFFFF40
  156.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  157.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  158.         posy = posy + 12
  159.     Next
  160.  
  161. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  162.     Dim Radius As Long, RadiusError As Long
  163.     Dim X As Long, Y As Long
  164.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  165.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  166.     Line (CX - X, CY)-(CX + X, CY), C, BF
  167.     While X > Y
  168.         RadiusError = RadiusError + Y * 2 + 1
  169.         If RadiusError >= 0 Then
  170.             If X <> Y + 1 Then
  171.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  172.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  173.             End If
  174.             X = X - 1
  175.             RadiusError = RadiusError - X * 2
  176.         End If
  177.         Y = Y + 1
  178.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  179.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  180.     Wend
  181.  
  182.  

 
Contour Map.PNG

Offline _vince

  • Seasoned Forum Regular
  • Posts: 422
    • View Profile
Re: Contour Plot
« Reply #1 on: November 16, 2021, 10:02:28 pm »
Nice, I would love to score some of zom-B's (if I remember correctly, couldve been someone else) software rendered 3D work from the original dot net. I remember there was a rotating yellow cup and also bouncing oranges and probably more.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #2 on: November 16, 2021, 10:43:08 pm »
Working towards landscape generation:
Code: QB64: [Select]
  1. _Title "Contour Plot 2" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2. ' Contour plot using Data Points by Mrwhy
  3. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  4. ' AndyA ref:
  5. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  6. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  7. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  8.  
  9. Dim Shared sw, sh
  10. sw = 640: sh = 480
  11. Screen _NewImage(sw, sh, 32)
  12. _ScreenMove 300, 100
  13. ReDim Shared pal(23) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  14.  
  15. 'blue green
  16. For i = 1 To 16
  17.     If i > 5 Then r = i * 15 Else r = 0
  18.     pal(i) = _RGB32(r, i * 10, 5 * (14 - i))
  19. '============================================
  20. ' number of height/potential points in data
  21. '============================================
  22. numPoints = 16
  23. '============================================
  24. restart:
  25. ReDim h(sw, sh)
  26. For i = 1 To 16
  27.     x(i) = Int((sw - 10) * Rnd + 5): y(i) = Int((sh - 10) * Rnd + 5): z(i) = Int(Rnd * 12 + 2)
  28. ztot = 0
  29. 'display height or potential data points on screen
  30. For i = 1 To numPoints
  31.     fcirc x(i), sh - y(i), 3, pal(z(i))
  32.     ztot = ztot + z(i)
  33.  
  34. Color &HFFFFFF40
  35.  
  36. s$ = "Calculating Contour Map"
  37. x = (_Width - _PrintWidth(s$)) / 2
  38. _PrintString (x, _Height - 14), s$
  39. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  40.  
  41. 'initialize some variables
  42. zmean = ztot / numPoints
  43. wo = sw * sh / numPoints
  44. w = wo / 2
  45.  
  46. 'generate initial Error estimates
  47. For i = 1 To numPoints
  48.     E(i) = z(i) - zmean
  49. Legend
  50.  
  51. 'begin Relaxation (iterative method)
  52. For jj = 1 To 9 * numPoints 'find max error point
  53.     emax = 0
  54.     For i = 1 To numPoints
  55.         If Abs(E(i)) > emax Then
  56.             emax = Abs(E(i))
  57.             ii = i
  58.         End If
  59.         k = E(ii)
  60.     Next
  61.     'fixit
  62.     For i = 1 To numPoints
  63.         dx = x(i) - x(ii)
  64.         dy = y(i) - y(ii)
  65.         dsq = dx * dx + dy * dy
  66.         E(i) = E(i) - k * Exp(-(dsq / w))
  67.         If i = ii Then
  68.             'update map with revised height or potential estimates for each pixel
  69.             For fy = 1 To (sh - 1)
  70.                 For fx = 1 To (sw - 1)
  71.                     dx = fx - x(ii)
  72.                     dy = fy - y(ii)
  73.                     dsq2 = dx * dx + dy * dy
  74.                     dy = sh - fy
  75.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  76.                 Next
  77.             Next
  78.         End If
  79.     Next
  80.  
  81. 'Draw calculated contour map
  82. Color pal(h(1, 1) + zmean)
  83. Line (0, 0)-(sw - 1, 0)
  84. Line (0, 0)-(0, sh - 1)
  85. For fy = 1 To sh - 1
  86.     For fx = 1 To sw - 1
  87.         If Int(h(fx, fy) + zmean) > 15 Then
  88.             c~& = pal(15)
  89.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  90.             c~& = pal(1)
  91.         Else
  92.             c~& = pal(Int(h(fx, fy) + zmean))
  93.         End If
  94.         PSet (fx, fy), c~&
  95.     Next
  96.  
  97. 'display height or potential data points on contour map
  98. For i = 1 To numPoints
  99.     'fcirc x(i), sh - y(i), 3, &HFF000000
  100.     fcirc x(i), sh - y(i), 2, pal(z(i))
  101. Legend
  102. GoTo restart
  103.  
  104. 'show height or potential value for each color in map
  105. Sub Legend
  106.     posy = 20
  107.     ' just draw a balck box where legend is going
  108.     Line (0, 0)-(37, 200), &HFFFF0000, BF
  109.     Line (1, 1)-(36, 199), &HFFFFFF00, B
  110.     Color &HFFFFFF00
  111.     _Font 8
  112.     _PrintString (7, 3), "Key"
  113.     For nn = 1 To 15
  114.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  115.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  116.         posy = posy + 12
  117.     Next
  118.  
  119. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  120.     Dim Radius As Long, RadiusError As Long
  121.     Dim X As Long, Y As Long
  122.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  123.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  124.     Line (CX - X, CY)-(CX + X, CY), C, BF
  125.     While X > Y
  126.         RadiusError = RadiusError + Y * 2 + 1
  127.         If RadiusError >= 0 Then
  128.             If X <> Y + 1 Then
  129.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  130.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  131.             End If
  132.             X = X - 1
  133.             RadiusError = RadiusError - X * 2
  134.         End If
  135.         Y = Y + 1
  136.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  137.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  138.     Wend
  139.  
  140.  

Contour Plot 2.PNG
* Contour Plot 2.PNG (Filesize: 37.33 KB, Dimensions: 639x507, Views: 222)

Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 1091
  • he lives
    • View Profile
Re: Contour Plot
« Reply #3 on: November 16, 2021, 10:47:56 pm »
Hat tip to you bplus! Everyone makes terrain using Perlin noise, I thought I was the only one doing relaxation methods. Did you know you're actually solving physics problems that way? You are 100% allowed to regard your high and low spots as "fixed charges", and then the relaxed field in between will be the true field.

Vince - you're damn right, at least I think. I'm pretty sure it was zom-b. I think these days, all we have is his raytracer.

You're not done when it works, you're done when it's right.

Offline _vince

  • Seasoned Forum Regular
  • Posts: 422
    • View Profile
Re: Contour Plot
« Reply #4 on: November 17, 2021, 04:37:00 am »
I now remember that it was harixxx who was doing the 3D stuff.  I can't believe I didnt save all of those demos but here is one i found

You need to have the qb64icon bitmap in the same directory to run this:
Code: QB64: [Select]
  1. _TITLE "Demo of the Versatile _PUTIMAGE by harixxx"
  2. SCREEN _NEWIMAGE(640, 480, 32)
  3.  
  4. TYPE point2d
  5. TYPE point3d
  6. TYPE face4points
  7.  
  8. REDIM texture&(0, 0)
  9. 'loadTexture "picture1.jpg", texture&(), imgWidth, imgHeight
  10. loadTexture "qb64icon.bmp", texture&(), imgWidth, imgHeight
  11. DIM grid(imgWidth, imgWidth) AS point2d
  12.  
  13. READ points
  14. DIM vec2d(points) AS point2d, vec3d(points) AS point3d
  15. FOR i = 1 TO points
  16.    READ vec3d(i).x, vec3d(i).y, vec3d(i).z
  17. READ faces
  18. DIM cubeside(faces) AS face4points
  19. FOR i = 1 TO faces
  20.    READ cubeside(i).a, cubeside(i).b, cubeside(i).c, cubeside(i).d
  21.  
  22. 'points data
  23. '3d points for cube (8 points x,y,z)
  24. data 100,100,100,  -100,100,100,  -100,100,-100,  100,100,-100
  25. data 100,-100,100, -100,-100,100, -100,-100,-100, 100,-100,-100
  26. '4 points face connection for cube (6 sides a,b,c,d)
  27. data 1,3,2,4,  1,8,4,5,  8,6,7,5,  6,3,7,2,  3,8,7,4,  1,6,5,2
  28.  
  29. centerx = _WIDTH / 2
  30. centery = _HEIGHT / 2
  31.    'update rotation
  32.    xr = (xr + 2) MOD 360
  33.    yr = (yr + 4) MOD 360
  34.    zr = (zr + 1) MOD 360
  35.    zoom = (zoom + 1) MOD 360
  36.    '3d to 2d rotation
  37.    FOR i = 1 TO points
  38.       x = vec3d(i).x
  39.       y = vec3d(i).y
  40.       z = vec3d(i).z
  41.       rot3dto2d x, y, z, xr, yr, zr, nx, ny, centerx, centery, sinT(zoom)
  42.       vec2d(i).x = nx
  43.       vec2d(i).y = ny
  44.    NEXT
  45.    CLS
  46.    '3d cube setup
  47.    FOR i = 1 TO faces
  48.       x1 = vec2d(cubeside(i).a).x: y1 = vec2d(cubeside(i).a).y
  49.       x2 = vec2d(cubeside(i).b).x: y2 = vec2d(cubeside(i).b).y
  50.       x3 = vec2d(cubeside(i).c).x: y3 = vec2d(cubeside(i).c).y
  51.       x4 = vec2d(cubeside(i).d).x: y4 = vec2d(cubeside(i).d).y
  52.       'find the middle of convex 3d vectors
  53.       IF convex3d(x1, y1, x2, y2, x3, y3) >= 0 THEN
  54.          sx1 = (x2 - x3) / imgWidth
  55.          sy1 = (y2 - y3) / imgHeight
  56.          sx2 = (x4 - x1) / imgWidth
  57.          sy2 = (y4 - y1) / imgHeight
  58.          'grids setup
  59.          FOR x = 0 TO imgWidth
  60.             'grid steps size
  61.             tx1 = sx1 * x + vec2d(cubeside(i).c).x
  62.             ty1 = sy1 * x + vec2d(cubeside(i).c).y
  63.             tx2 = (sx2 * x + vec2d(cubeside(i).a).x - tx1) / imgWidth
  64.             ty2 = (sy2 * x + vec2d(cubeside(i).a).y - ty1) / imgHeight
  65.             FOR y = 0 TO imgHeight
  66.                grid(x, y).x = tx2 * y + tx1
  67.                grid(x, y).y = ty2 * y + ty1
  68.          NEXT y, x
  69.          'draw 3d cube
  70.          FOR y = 0 TO imgHeight - 1
  71.             FOR x = 0 TO imgWidth - 1
  72.                x1 = grid(x, y).x: y1 = grid(x, y).y
  73.                x2 = grid(x, y + 1).x: y2 = grid(x, y + 1).y
  74.                x3 = grid(x + 1, y + 1).x: y3 = grid(x + 1, y + 1).y
  75.                x4 = grid(x + 1, y).x: y4 = grid(x + 1, y).y
  76.                'draw 1 quadrangle polygon by connecting 2 triangles polygon
  77.                IF x1 > 0 OR x2 > 0 AND x3 > 0 OR x4 > 0 THEN
  78.                   triangle x1, y1, x2, y2, x3, y3, texture&(x, y)
  79.                   triangle x1, y1, x4, y4, x3, y3, texture&(x, y)
  80.                END IF
  81.          NEXT x, y
  82.       END IF
  83.    NEXT
  84.  
  85.    _LIMIT 30
  86.  
  87. '================= sines & cosines table
  88. FUNCTION sinT (n)
  89. sinT = SIN(3.141593 / 180 * n)
  90. FUNCTION cosT (n)
  91. cosT = COS(3.141593 / 180 * n)
  92.  
  93. '======================== convex 3d
  94. 'routines to find middle of convex 3d vectors
  95. '(http://en.wikipedia.org/wiki/Convex_function)
  96. FUNCTION convex3d (x1, y1, x2, y2, x3, y3)
  97. convex3d = (x1 - x2) * (y3 - y2) - (x3 - x2) * (y1 - y2)
  98.  
  99. '======================== 3d to 2d rotation
  100. 'thank's to Sami Kyostila 12-24-1997 for 3d to 2d rotation
  101. 'but i don't know who the first creator for this code :D
  102. SUB rot3dto2d (x, y, z, rx, ry, rz, nx, ny, cx, cy, sz)
  103. x1 = x: y1 = y: z1 = z
  104. xx = x1
  105. yy = y1 * cosT(rx) + z1 * sinT(rx)
  106. zz = z1 * cosT(rx) - y1 * sinT(rx)
  107. y1 = yy
  108. x1 = xx * cosT(ry) - zz * sinT(ry)
  109. z1 = xx * sinT(ry) + zz * cosT(ry)
  110. zz = z1
  111. xx = x1 * cosT(rz) - y1 * sinT(rz)
  112. yy = x1 * sinT(rz) + y1 * cosT(rz)
  113. fc = zz / 500 + 1
  114. nx = xx / fc * sz * 1.2 + cx
  115. ny = yy / fc * sz * 1.2 + cy
  116.  
  117. '============================= make texture from file
  118. 'make REDIM array as long before calling this routine
  119. 'because array updated and contain of 32bits color
  120. 'example: 'REDIM texture&(0, 0)' or 'REDIM texture(0, 0) as long'
  121. SUB loadTexture (imageFile$, textureName&(), imgWidth, imgHeight)
  122. myImg& = _LOADIMAGE(imageFile$)
  123. imgWidth = _WIDTH(myImg&)
  124. imgHeight = _HEIGHT(myImg&)
  125. REDIM textureName&(imgWidth, imgHeight)
  126. _SOURCE myImg&
  127. FOR y = 0 TO imgHeight
  128.    FOR x = 0 TO imgWidth
  129.       textureName&(x, y) = POINT(x, y)
  130. NEXT x, y
  131.  
  132. '=============================== solid triangle
  133. '/------------ my solid triangle v1.1b
  134. SUB triangle (x1%, y1%, x2%, y2%, x3%, y3%, c&)
  135. px1 = x1%: px2 = x2%: px3 = x3%
  136. py1 = y1%: py2 = y2%: py3 = y3%
  137. '(thank's to relminator for 'sort y' correction
  138. ' i see your comment on the mono&disco module)
  139. 'sort y
  140. IF py1 > py2 THEN SWAP py1, py2: SWAP px1, px2
  141. IF py1 > py3 THEN SWAP py1, py3: SWAP px1, px3
  142. IF py2 > py3 THEN SWAP py2, py3: SWAP px2, px3
  143. 'alpha
  144. A = px1: B = px1: C = px2
  145. 'delta
  146. dA = (px1 - px3) / (py1 - py3)
  147. dB = (px1 - px2) / (py1 - py2)
  148. dC = (px2 - px3) / (py2 - py3)
  149. 'wow! your triangle shading close to mine
  150. '(thank's to wiki for trigonometry
  151. ' http://en.wikipedia.org/wiki/Triangle#Using_trigonometry)
  152. 'triangle shading
  153. FOR y% = py1 TO py3
  154.    xx1% = A
  155.    IF y% < py2 THEN xx2% = B: B = B + dB ELSE xx2% = C: C = C + dC
  156.    LINE (xx1%, y%)-(xx2%, y%), c&
  157.    A = A + dA
  158.  
  159.  

* qb64icon.bmp (Filesize: 2.05 KB, Dimensions: 32x32, Views: 474)

Offline DANILIN

  • Forum Regular
  • Posts: 128
    • View Profile
    • Danilin youtube
Re: Contour Plot
« Reply #5 on: November 17, 2021, 06:51:23 am »
use pseudo 3d relief qb64 my milli program:

on: March 29, 2019, 07:01:25 AM

https://qb64.org/forum/index.php?topic=702.msg103948#msg103948

 
reliefqb.gif


New: Relief 3d multivariate parametric
https://qb64.org/forum/index.php?topic=4398.0
« Last Edit: November 17, 2021, 04:50:37 pm by DANILIN »
Russia looks world from future. big data is peace data.
https://youtube.com/playlist?list=PLBBTP9oVY7IagpH0g9FNUQ8JqmHwxDDDB
i never recommend anything to anyone and always write only about myself

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #6 on: November 17, 2021, 02:35:40 pm »
Plot #3 Earthlike planets
Code: QB64: [Select]
  1. _Title "Contour Plot 3: spacebar for New World, escape to quit" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2.  
  3. ' Contour plot using Data Points by Mrwhy
  4. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  5. ' AndyA ref:
  6. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  7. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  8. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  9.  
  10. Dim Shared sw, sh, nC, nP
  11. sw = 800: sh = 400: nC = 32: nP = 12 ' screen Width and Height, number of Colors, number of Points
  12. Screen _NewImage(sw, sh, 32)
  13. _ScreenMove 200, 60 ' at least 60 on each side and to and bottom
  14. ReDim Shared pal(nC) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  15.  
  16. 'color mixing green blue
  17. For i = 1 To nC
  18.     pal(i) = _RGB32(0, i * 255 / nC, 128 / nC * (nC - i))
  19.  
  20. restart:
  21. ReDim h(sw, sh)
  22. For i = 1 To nP ' new data points and color assignments
  23.     x(i) = Int((sw - 30) * Rnd + 15): y(i) = Int((sh - 30) * Rnd + 15): z(i) = Int(Rnd * (nC - 2) + 1)
  24. ztot = 0
  25. 'display height or potential data points on screen
  26. For i = 1 To nP
  27.     fcirc x(i), sh - y(i), 3, pal(z(i))
  28.     ztot = ztot + z(i)
  29.  
  30. Color &HFFFFFF40
  31.  
  32. s$ = "Calculating Contour Map"
  33. x = (_Width - _PrintWidth(s$)) / 2
  34. _PrintString (x, _Height - 14), s$
  35. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  36.  
  37. 'initialize some variables
  38. zmean = ztot / nP
  39. wo = sw * sh / nP
  40. w = wo / 2
  41.  
  42. 'generate initial Error estimates
  43. For i = 1 To nP
  44.     E(i) = z(i) - zmean
  45. 'Legend
  46.  
  47. 'begin Relaxation (iterative method)
  48. For jj = 1 To 9 * nP 'find max error point
  49.     emax = 0
  50.     For i = 1 To nP
  51.         If Abs(E(i)) > emax Then
  52.             emax = Abs(E(i))
  53.             ii = i
  54.         End If
  55.         k = E(ii)
  56.     Next
  57.     'fixit
  58.     For i = 1 To nP
  59.         dx = x(i) - x(ii)
  60.         dy = y(i) - y(ii)
  61.         dsq = dx * dx + dy * dy
  62.         E(i) = E(i) - k * Exp(-(dsq / w))
  63.         If i = ii Then
  64.             'update map with revised height or potential estimates for each pixel
  65.             For fy = 1 To (sh - 1)
  66.                 For fx = 1 To (sw - 1)
  67.                     dx = fx - x(ii)
  68.                     dy = fy - y(ii)
  69.                     dsq2 = dx * dx + dy * dy
  70.                     dy = sh - fy
  71.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  72.                 Next
  73.             Next
  74.         End If
  75.     Next
  76.  
  77. 'Draw calculated contour map
  78. Color pal(h(1, 1) + zmean)
  79. Line (0, 0)-(sw - 1, 0)
  80. Line (0, 0)-(0, sh - 1)
  81. For fy = 1 To sh - 1
  82.     For fx = 1 To sw - 1
  83.         If Int(h(fx, fy) + zmean) > nC Then
  84.             c~& = pal(nC)
  85.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  86.             c~& = pal(1)
  87.         Else
  88.             c~& = pal(Int(h(fx, fy) + zmean))
  89.         End If
  90.         PSet (fx, fy), c~&
  91.     Next
  92.  
  93. 'display height or potential data points on contour map
  94. 'For i = 1 To nP
  95. '    'fcirc x(i), sh - y(i), 3, &HFF000000
  96. '    fcirc x(i), sh - y(i), 2, pal(z(i))
  97. 'Next
  98. 'Legend
  99.  
  100. surface& = _NewImage(_Width, _Height, 32)
  101. mid = Int(sw / 2)
  102. _PutImage (mid - 2, 0)-(sw - 1, sh), 0, surface&
  103. _PutImage (0, 0)-(mid + 2, sh), 0, surface&, (sw - 1, 0)-(0, sh)
  104. _PutImage , surface&, 0
  105. For y = 0 To sh ' fix seam !
  106.     PSet (mid + 2, y), Point(mid - 1, y)
  107.     PSet (mid + 1, y), Point(mid - 1, y)
  108.     PSet (mid, y), Point(mid - 1, y)
  109. _PutImage , 0, surface&
  110. '_PutImage , surface&, 0  'check
  111.  
  112. Color &HFF000000 ' be rid of white points at poles
  113. r = sh / 2.5: xc = sw / 2: yc = sh / 2: xo = 0: sr = 200
  114.     For y = -r To r
  115.         x1 = Sqr(r * r - y * y)
  116.         tv = (_Asin(y / r) + 1.5) / 3
  117.         For x = -x1 To x1
  118.             tu = (_Asin(x / x1) + 1.5) / 6
  119.             _Source surface&
  120.             pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  121.             _Dest 0
  122.             PSet (x + xc, y + yc), pc~&
  123.         Next x
  124.     Next y
  125.     xo = xo - 1 + sw
  126.     xo = xo Mod sw
  127.     If (nC + 1) * 12 <= sh Then Legend
  128.     _Display
  129.     _Limit 30
  130. _FreeImage surface&
  131. If _KeyDown(32) Then GoTo restart
  132.  
  133.  
  134. 'show height or potential value for each color in map
  135. Sub Legend
  136.     posy = 12
  137.     ' just draw a balck box where legend is going
  138.     Line (0, 0)-(37, (nC + 1) * 12), &HFFFF0000, BF
  139.     Line (1, 1)-(36, (nC + 1) * 12 - 2), &HFFFFFF00, B
  140.     Color &HFFFFFF00
  141.     _Font 8
  142.     _PrintString (7, 3), "Key"
  143.     For nn = 1 To nC
  144.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  145.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  146.         posy = posy + 12
  147.     Next
  148.  
  149. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  150.     Dim Radius As Long, RadiusError As Long
  151.     Dim X As Long, Y As Long
  152.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  153.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  154.     Line (CX - X, CY)-(CX + X, CY), C, BF
  155.     While X > Y
  156.         RadiusError = RadiusError + Y * 2 + 1
  157.         If RadiusError >= 0 Then
  158.             If X <> Y + 1 Then
  159.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  160.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  161.             End If
  162.             X = X - 1
  163.             RadiusError = RadiusError - X * 2
  164.         End If
  165.         Y = Y + 1
  166.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  167.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  168.     Wend
  169.  
  170.  

 
Plot #3.PNG

Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 1091
  • he lives
    • View Profile
Re: Contour Plot
« Reply #7 on: November 17, 2021, 08:33:01 pm »
bplus, my maaaaaaan....

nice work. say, give this a run, and cross your eyes so that the two planets coincide. stereoscopic 3d baby!

Code: QB64: [Select]
  1. _Title "Contour Plot 3: spacebar for New World, escape to quit" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2.  
  3. ' Contour plot using Data Points by Mrwhy
  4. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  5. ' AndyA ref:
  6. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  7. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  8. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  9.  
  10. Dim Shared sw, sh, nC, nP
  11. sw = 800: sh = 400: nC = 32: nP = 12 ' screen Width and Height, number of Colors, number of Points
  12. Screen _NewImage(sw, sh, 32)
  13. _ScreenMove 200, 60 ' at least 60 on each side and to and bottom
  14. ReDim Shared pal(nC) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  15.  
  16. 'color mixing green blue
  17. For i = 1 To nC
  18.     pal(i) = _RGB32(0, i * 255 / nC, 128 / nC * (nC - i))
  19.  
  20. restart:
  21. ReDim h(sw, sh)
  22. For i = 1 To nP ' new data points and color assignments
  23.     x(i) = Int((sw - 30) * Rnd + 15): y(i) = Int((sh - 30) * Rnd + 15): z(i) = Int(Rnd * (nC - 2) + 1)
  24. ztot = 0
  25. 'display height or potential data points on screen
  26. For i = 1 To nP
  27.     fcirc x(i), sh - y(i), 3, pal(z(i))
  28.     ztot = ztot + z(i)
  29.  
  30. Color &HFFFFFF40
  31.  
  32. s$ = "Calculating Contour Map"
  33. x = (_Width - _PrintWidth(s$)) / 2
  34. _PrintString (x, _Height - 14), s$
  35. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  36.  
  37. 'initialize some variables
  38. zmean = ztot / nP
  39. wo = sw * sh / nP
  40. w = wo / 2
  41.  
  42. 'generate initial Error estimates
  43. For i = 1 To nP
  44.     E(i) = z(i) - zmean
  45. 'Legend
  46.  
  47. 'begin Relaxation (iterative method)
  48. For jj = 1 To 9 * nP 'find max error point
  49.     emax = 0
  50.     For i = 1 To nP
  51.         If Abs(E(i)) > emax Then
  52.             emax = Abs(E(i))
  53.             ii = i
  54.         End If
  55.         k = E(ii)
  56.     Next
  57.     'fixit
  58.     For i = 1 To nP
  59.         dx = x(i) - x(ii)
  60.         dy = y(i) - y(ii)
  61.         dsq = dx * dx + dy * dy
  62.         E(i) = E(i) - k * Exp(-(dsq / w))
  63.         If i = ii Then
  64.             'update map with revised height or potential estimates for each pixel
  65.             For fy = 1 To (sh - 1)
  66.                 For fx = 1 To (sw - 1)
  67.                     dx = fx - x(ii)
  68.                     dy = fy - y(ii)
  69.                     dsq2 = dx * dx + dy * dy
  70.                     dy = sh - fy
  71.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  72.                 Next
  73.             Next
  74.         End If
  75.     Next
  76.  
  77. 'Draw calculated contour map
  78. Color pal(h(1, 1) + zmean)
  79. Line (0, 0)-(sw - 1, 0)
  80. Line (0, 0)-(0, sh - 1)
  81. For fy = 1 To sh - 1
  82.     For fx = 1 To sw - 1
  83.         If Int(h(fx, fy) + zmean) > nC Then
  84.             c~& = pal(nC)
  85.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  86.             c~& = pal(1)
  87.         Else
  88.             c~& = pal(Int(h(fx, fy) + zmean))
  89.         End If
  90.         PSet (fx, fy), c~&
  91.     Next
  92.  
  93. 'display height or potential data points on contour map
  94. 'For i = 1 To nP
  95. '    'fcirc x(i), sh - y(i), 3, &HFF000000
  96. '    fcirc x(i), sh - y(i), 2, pal(z(i))
  97. 'Next
  98. 'Legend
  99.  
  100. surface& = _NewImage(_Width, _Height, 32)
  101. mid = Int(sw / 2)
  102. _PutImage (mid - 2, 0)-(sw - 1, sh), 0, surface&
  103. _PutImage (0, 0)-(mid + 2, sh), 0, surface&, (sw - 1, 0)-(0, sh)
  104. _PutImage , surface&, 0
  105. For y = 0 To sh ' fix seam !
  106.     PSet (mid + 2, y), Point(mid - 1, y)
  107.     PSet (mid + 1, y), Point(mid - 1, y)
  108.     PSet (mid, y), Point(mid - 1, y)
  109. _PutImage , 0, surface&
  110. '_PutImage , surface&, 0  'check
  111.  
  112. Color &HFF000000 ' be rid of white points at poles
  113. r = sh / 2.5: xc = sw / 2: yc = sh / 2: xo = 0: sr = 200
  114.     For y = -r To r
  115.         x1 = Sqr(r * r - y * y)
  116.         tv = (_Asin(y / r) + 1.5) / 3
  117.         For x = -x1 To x1
  118.             tu = (_Asin(x / x1) + 1.5) / 6 - .05
  119.             _Source surface&
  120.             pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  121.             _Dest 0
  122.             PSet (180 + x + xc, y + yc), pc~&
  123.         Next x
  124.     Next y
  125.     xo = xo - 1 + sw
  126.     xo = xo Mod sw
  127.     If (nC + 1) * 12 <= sh Then Legend
  128.  
  129.     For y = -r To r
  130.         x1 = Sqr(r * r - y * y)
  131.         tv = (_Asin(y / r) + 1.5) / 3
  132.         For x = -x1 To x1
  133.             tu = (_Asin(x / x1) + 1.5) / 6
  134.             _Source surface&
  135.             pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  136.             _Dest 0
  137.             PSet (-180 + x + xc, y + yc), pc~&
  138.         Next x
  139.     Next y
  140.     xo = xo - 1 + sw
  141.     xo = xo Mod sw
  142.     If (nC + 1) * 12 <= sh Then Legend
  143.  
  144.  
  145.  
  146.  
  147.  
  148.     _Display
  149.     _Limit 30
  150. _FreeImage surface&
  151. If _KeyDown(32) Then GoTo restart
  152.  
  153.  
  154. 'show height or potential value for each color in map
  155. Sub Legend
  156.     posy = 12
  157.     ' just draw a balck box where legend is going
  158.     Line (0, 0)-(37, (nC + 1) * 12), &HFFFF0000, BF
  159.     Line (1, 1)-(36, (nC + 1) * 12 - 2), &HFFFFFF00, B
  160.     Color &HFFFFFF00
  161.     _Font 8
  162.     _PrintString (7, 3), "Key"
  163.     For nn = 1 To nC
  164.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  165.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  166.         posy = posy + 12
  167.     Next
  168.  
  169. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  170.     Dim Radius As Long, RadiusError As Long
  171.     Dim X As Long, Y As Long
  172.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  173.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  174.     Line (CX - X, CY)-(CX + X, CY), C, BF
  175.     While X > Y
  176.         RadiusError = RadiusError + Y * 2 + 1
  177.         If RadiusError >= 0 Then
  178.             If X <> Y + 1 Then
  179.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  180.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  181.             End If
  182.             X = X - 1
  183.             RadiusError = RadiusError - X * 2
  184.         End If
  185.         Y = Y + 1
  186.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  187.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  188.     Wend
  189.  
  190.  
  191.  

 
ss1.png
You're not done when it works, you're done when it's right.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #8 on: November 17, 2021, 09:37:05 pm »
Could not see it, so I started looking for very old 3D glasses I found the other day... nah that's not going to work either... How about super imposing the images advance then for a moment flick back then advance further... try to send the eye and brain 2 different messages.

Eh? use spacebar to toggle, jiggle mode. Watch Title bar if you can't tell the mode:

Code: QB64: [Select]
  1. _Title "Contour Plot 3: spacebar for New World, escape to quit" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2.  
  3. ' Contour plot using Data Points by Mrwhy
  4. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  5. ' AndyA ref:
  6. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  7. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  8. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  9.  
  10. Dim Shared sw, sh, nC, nP
  11. sw = 800: sh = 400: nC = 32: nP = 12 ' screen Width and Height, number of Colors, number of Points
  12. Screen _NewImage(sw, sh, 32)
  13. _ScreenMove 200, 60 ' at least 60 on each side and to and bottom
  14. ReDim Shared pal(nC) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  15.  
  16. 'color mixing green blue
  17. For i = 1 To nC
  18.     pal(i) = _RGB32(0, i * 255 / nC, 128 / nC * (nC - i))
  19.  
  20. restart:
  21. ReDim h(sw, sh)
  22. For i = 1 To nP ' new data points and color assignments
  23.     x(i) = Int((sw - 30) * Rnd + 15): y(i) = Int((sh - 30) * Rnd + 15): z(i) = Int(Rnd * (nC - 2) + 1)
  24. ztot = 0
  25. 'display height or potential data points on screen
  26. For i = 1 To nP
  27.     fcirc x(i), sh - y(i), 3, pal(z(i))
  28.     ztot = ztot + z(i)
  29.  
  30. Color &HFFFFFF40
  31.  
  32. s$ = "Calculating Contour Map"
  33. x = (_Width - _PrintWidth(s$)) / 2
  34. _PrintString (x, _Height - 14), s$
  35. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  36.  
  37. 'initialize some variables
  38. zmean = ztot / nP
  39. wo = sw * sh / nP
  40. w = wo / 2
  41.  
  42. 'generate initial Error estimates
  43. For i = 1 To nP
  44.     E(i) = z(i) - zmean
  45. 'Legend
  46.  
  47. 'begin Relaxation (iterative method)
  48. For jj = 1 To 9 * nP 'find max error point
  49.     emax = 0
  50.     For i = 1 To nP
  51.         If Abs(E(i)) > emax Then
  52.             emax = Abs(E(i))
  53.             ii = i
  54.         End If
  55.         k = E(ii)
  56.     Next
  57.     'fixit
  58.     For i = 1 To nP
  59.         dx = x(i) - x(ii)
  60.         dy = y(i) - y(ii)
  61.         dsq = dx * dx + dy * dy
  62.         E(i) = E(i) - k * Exp(-(dsq / w))
  63.         If i = ii Then
  64.             'update map with revised height or potential estimates for each pixel
  65.             For fy = 1 To (sh - 1)
  66.                 For fx = 1 To (sw - 1)
  67.                     dx = fx - x(ii)
  68.                     dy = fy - y(ii)
  69.                     dsq2 = dx * dx + dy * dy
  70.                     dy = sh - fy
  71.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  72.                 Next
  73.             Next
  74.         End If
  75.     Next
  76.  
  77. 'Draw calculated contour map
  78. Color pal(h(1, 1) + zmean)
  79. Line (0, 0)-(sw - 1, 0)
  80. Line (0, 0)-(0, sh - 1)
  81. For fy = 1 To sh - 1
  82.     For fx = 1 To sw - 1
  83.         If Int(h(fx, fy) + zmean) > nC Then
  84.             c~& = pal(nC)
  85.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  86.             c~& = pal(1)
  87.         Else
  88.             c~& = pal(Int(h(fx, fy) + zmean))
  89.         End If
  90.         PSet (fx, fy), c~&
  91.     Next
  92.  
  93. 'display height or potential data points on contour map
  94. 'For i = 1 To nP
  95. '    'fcirc x(i), sh - y(i), 3, &HFF000000
  96. '    fcirc x(i), sh - y(i), 2, pal(z(i))
  97. 'Next
  98. 'Legend
  99.  
  100. surface& = _NewImage(_Width, _Height, 32)
  101. mid = Int(sw / 2)
  102. _PutImage (mid - 2, 0)-(sw - 1, sh), 0, surface&
  103. _PutImage (0, 0)-(mid + 2, sh), 0, surface&, (sw - 1, 0)-(0, sh)
  104. _PutImage , surface&, 0
  105. For y = 0 To sh ' fix seam !
  106.     PSet (mid + 2, y), Point(mid - 1, y)
  107.     PSet (mid + 1, y), Point(mid - 1, y)
  108.     PSet (mid, y), Point(mid - 1, y)
  109. _PutImage , 0, surface&
  110. '_PutImage , surface&, 0  'check
  111.  
  112. Color &HFFFFFF00 ' be rid of white points at poles , eh! don't know where the white cresents are coming from!!!
  113. r = sh / 2.5: xc = sw / 2: yc = sh / 2: xo = 0: sr = 200
  114.     'Cls
  115.     If InKey$ = " " Then
  116.         If jiggle Then
  117.             jiggle = 0
  118.             _Title "Contour Plot 3: spacebar for New World, escape to quit"
  119.         Else
  120.             jiggle = -1
  121.             _Title "Contour Plot 3: spacebar for New World, escape to quit, jiggle: ON"
  122.         End If
  123.     End If
  124.     For y = -r To r
  125.         x1 = Sqr(r * r - y * y)
  126.         tv = (_Asin(y / r) + 1.5) / 3
  127.         For x = -x1 To x1
  128.             tu = (_Asin(x / x1) + 1.5) / 6
  129.             _Source surface&
  130.             pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  131.             _Dest 0
  132.             PSet (x + xc, y + yc), pc~&
  133.         Next x
  134.     Next y
  135.     xo = xo - 1 + sw
  136.     xo = xo Mod sw
  137.     If (nC + 1) * 12 <= sh Then Legend
  138.     _Display ' this advance display should dwell a little longer and next will fall back
  139.  
  140.     ' ++++++++++++++++++++ press spacebar to toggle jiggle and no jiggle
  141.     If jiggle Then
  142.  
  143.         For y = -r To r
  144.             x1 = Sqr(r * r - y * y)
  145.             tv = (_Asin(y / r) + 1.5) / 3
  146.             For x = -x1 To x1
  147.                 tu = (_Asin(x / x1) + 1.5) / 6 - .005 ' <<<<< jiggle to send 2 signals to the eyes and brain
  148.                 _Source surface&
  149.                 pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  150.                 _Dest 0
  151.                 PSet (x + xc, y + yc), pc~& ' <<<<< jiggle to send 2 signals to the eyes and brain
  152.             Next x
  153.         Next y
  154.         xo = xo - 1 + sw
  155.         xo = xo Mod sw
  156.         'If (nC + 1) * 12 <= sh Then Legend  'one legend per frame is enough
  157.         _Display
  158.         '_Limit 30  ' go as fast as you can QB64!!!
  159.     End If
  160. _FreeImage surface&
  161. If _KeyDown(32) Then GoTo restart
  162.  
  163.  
  164. 'show height or potential value for each color in map
  165. Sub Legend
  166.     posy = 12
  167.     ' just draw a balck box where legend is going
  168.     Line (0, 0)-(37, (nC + 1) * 12), &HFFFF0000, BF
  169.     Line (1, 1)-(36, (nC + 1) * 12 - 2), &HFFFFFF00, B
  170.     Color &HFFFFFF00
  171.     _Font 8
  172.     _PrintString (7, 3), "Key"
  173.     For nn = 1 To nC
  174.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  175.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  176.         posy = posy + 12
  177.     Next
  178.  
  179. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  180.     Dim Radius As Long, RadiusError As Long
  181.     Dim X As Long, Y As Long
  182.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  183.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  184.     Line (CX - X, CY)-(CX + X, CY), C, BF
  185.     While X > Y
  186.         RadiusError = RadiusError + Y * 2 + 1
  187.         If RadiusError >= 0 Then
  188.             If X <> Y + 1 Then
  189.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  190.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  191.             End If
  192.             X = X - 1
  193.             RadiusError = RadiusError - X * 2
  194.         End If
  195.         Y = Y + 1
  196.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  197.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  198.     Wend
  199.  
  200.  
  201.  

Mostly just get car sick :(


Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 1091
  • he lives
    • View Profile
Re: Contour Plot
« Reply #9 on: November 18, 2021, 08:13:36 am »
bplus you're scaring me!

The other day you see a splotch of yellow in a program that nobody else could see, and three days later you aren't seeing in 3d? How close do you sit to your monitor? Do you view it from straight on? It's a flatscreen right? You also have... ahem... have two working eyes I presume?

Nah in all seriousness, I know how to make it work with color-filtered 3D glasses. (If your glasses look ordinary but still claim to be 3D glasses, they are probably the wrong type. We want different colors for each eye, not a different polarization for each eye.) Anyway, plot your planet for one eye with all red removed, and for the other eye, plot with all blue removed and red restored. Give that a try maybe? Along those lines, you're like 2 steps away from creating your own magic eye images with this program.

You're not done when it works, you're done when it's right.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #10 on: November 18, 2021, 12:08:40 pm »
bplus you're scaring me!

The other day you see a splotch of yellow in a program that nobody else could see, and three days later you aren't seeing in 3d? How close do you sit to your monitor? Do you view it from straight on? It's a flatscreen right? You also have... ahem... have two working eyes I presume?

Nah in all seriousness, I know how to make it work with color-filtered 3D glasses. (If your glasses look ordinary but still claim to be 3D glasses, they are probably the wrong type. We want different colors for each eye, not a different polarization for each eye.) Anyway, plot your planet for one eye with all red removed, and for the other eye, plot with all blue removed and red restored. Give that a try maybe? Along those lines, you're like 2 steps away from creating your own magic eye images with this program.

Two working eyes yes, working together, no. One does close up work like reading and the other use to do distance stuff but getting bad.

Yellow thing I think is optical illusion but you have to run it at just the right speed, it's like seeing dots where a white line grid lays over a dark background. You look at any intersection directly and nothing just intersect of 2 white lines but from a distance dots that disappear and reappear as you move your eyes over the image. Not going to argue that one because the better effect is captured in my avatar.

I don't know how I am supposed to see 3D Stereo when 2 very separate globes are offset in spinning positions even if one is red and one is not. I think they have to mostly overlap and appear to each eye as from a slightly different angle, offsetting the rotation won't do it, IMHO.

I think shading would do a better job of 3D work, oh an idea...!

Does anyone else see 3D from those 2 very separated globes in STx version?


Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #11 on: November 18, 2021, 12:57:55 pm »
Here is better 3D effect with shading:
Code: QB64: [Select]
  1. _Title "Contour Plot 3S: spacebar for New World, escape to quit" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  2. ' 2021-11-17 Contour Plot 3: planets
  3. ' 2021-11-18 3S try adding shade
  4.  
  5. ' Contour plot using Data Points by Mrwhy
  6. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  7. ' AndyA ref:
  8. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  9. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  10. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  11.  
  12. Dim Shared sw, sh, nC, nP
  13. sw = 800: sh = 400: nC = 32: nP = 12 ' screen Width and Height, number of Colors, number of Points
  14. Screen _NewImage(sw, sh, 32)
  15. _ScreenMove 200, 60 ' at least 60 on each side and to and bottom
  16. ReDim Shared pal(nC) As _Unsigned Long, x(50), y(50), z(100), E(50), h(sw, sh)
  17.  
  18. 'color mixing green blue
  19. For i = 1 To nC
  20.     pal(i) = _RGB32(0, i * 255 / nC, 128 / nC * (nC - i))
  21.  
  22. restart:
  23. ReDim h(sw, sh)
  24. For i = 1 To nP ' new data points and color assignments
  25.     x(i) = Int((sw - 30) * Rnd + 15): y(i) = Int((sh - 30) * Rnd + 15): z(i) = Int(Rnd * (nC - 2) + 1)
  26. ztot = 0
  27. 'display height or potential data points on screen
  28. For i = 1 To nP
  29.     fcirc x(i), sh - y(i), 3, pal(z(i))
  30.     ztot = ztot + z(i)
  31.  
  32. Color &HFFFFFF40
  33.  
  34. s$ = "Calculating Contour Map"
  35. x = (_Width - _PrintWidth(s$)) / 2
  36. _PrintString (x, _Height - 14), s$
  37. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  38.  
  39. 'initialize some variables
  40. zmean = ztot / nP
  41. wo = sw * sh / nP
  42. w = wo / 2
  43.  
  44. 'generate initial Error estimates
  45. For i = 1 To nP
  46.     E(i) = z(i) - zmean
  47. 'Legend
  48.  
  49. 'begin Relaxation (iterative method)
  50. For jj = 1 To 9 * nP 'find max error point
  51.     emax = 0
  52.     For i = 1 To nP
  53.         If Abs(E(i)) > emax Then
  54.             emax = Abs(E(i))
  55.             ii = i
  56.         End If
  57.         k = E(ii)
  58.     Next
  59.     'fixit
  60.     For i = 1 To nP
  61.         dx = x(i) - x(ii)
  62.         dy = y(i) - y(ii)
  63.         dsq = dx * dx + dy * dy
  64.         E(i) = E(i) - k * Exp(-(dsq / w))
  65.         If i = ii Then
  66.             'update map with revised height or potential estimates for each pixel
  67.             For fy = 1 To (sh - 1)
  68.                 For fx = 1 To (sw - 1)
  69.                     dx = fx - x(ii)
  70.                     dy = fy - y(ii)
  71.                     dsq2 = dx * dx + dy * dy
  72.                     dy = sh - fy
  73.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  74.                 Next
  75.             Next
  76.         End If
  77.     Next
  78.  
  79. 'Draw calculated contour map
  80. Color pal(h(1, 1) + zmean)
  81. Line (0, 0)-(sw - 1, 0)
  82. Line (0, 0)-(0, sh - 1)
  83. For fy = 1 To sh - 1
  84.     For fx = 1 To sw - 1
  85.         If Int(h(fx, fy) + zmean) > nC Then
  86.             c~& = pal(nC)
  87.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  88.             c~& = pal(1)
  89.         Else
  90.             c~& = pal(Int(h(fx, fy) + zmean))
  91.         End If
  92.         PSet (fx, fy), c~&
  93.     Next
  94.  
  95. 'display height or potential data points on contour map
  96. 'For i = 1 To nP
  97. '    'fcirc x(i), sh - y(i), 3, &HFF000000
  98. '    fcirc x(i), sh - y(i), 2, pal(z(i))
  99. 'Next
  100. 'Legend
  101.  
  102. surface& = _NewImage(_Width, _Height, 32)
  103. mid = Int(sw / 2)
  104. _PutImage (mid - 2, 0)-(sw - 1, sh), 0, surface&
  105. _PutImage (0, 0)-(mid + 2, sh), 0, surface&, (sw - 1, 0)-(0, sh)
  106. _PutImage , surface&, 0
  107. For y = 0 To sh ' fix seam !
  108.     PSet (mid + 2, y), Point(mid - 1, y)
  109.     PSet (mid + 1, y), Point(mid - 1, y)
  110.     PSet (mid, y), Point(mid - 1, y)
  111. _PutImage , 0, surface&
  112. '_PutImage , surface&, 0  'check
  113.  
  114. Color &HFF000000 ' be rid of white points at poles
  115. r = sh / 2.5: xc = sw / 2: yc = sh / 2: xo = 0: sr = 200
  116.     For y = -r To r
  117.         x1 = Sqr(r * r - y * y)
  118.         tv = (_Asin(y / r) + 1.5) / 3
  119.         For x = -x1 + 1 To x1
  120.             tu = (_Asin(x / x1) + 1.5) / 6
  121.             _Source surface&
  122.             pc~& = Point((xo + tu * (sw - 1)) Mod sw, tv * sh)
  123.             _Dest 0
  124.             dist = (x * x + y * y) / (r * r) ' for shading
  125.             PSet (x + xc, y + yc), pc~&
  126.             PSet (x + xc, y + yc), _RGBA32(0, 0, 0, 180 * (dist ^ 2)) ' for shading
  127.         Next x
  128.     Next y
  129.     xo = xo - 1 + sw
  130.     xo = xo Mod sw
  131.     If (nC + 1) * 12 <= sh Then Legend
  132.     _Display
  133.     _Limit 30
  134. _FreeImage surface&
  135. If _KeyDown(32) Then GoTo restart
  136.  
  137.  
  138. 'show height or potential value for each color in map
  139. Sub Legend
  140.     posy = 12
  141.     ' just draw a balck box where legend is going
  142.     Line (0, 0)-(37, (nC + 1) * 12), &HFFFF0000, BF
  143.     Line (1, 1)-(36, (nC + 1) * 12 - 2), &HFFFFFF00, B
  144.     Color &HFFFFFF00
  145.     _Font 8
  146.     _PrintString (7, 3), "Key"
  147.     For nn = 1 To nC
  148.         _PrintString (4, posy), Right$(" " + _Trim$(Str$(nn)), 2)
  149.         Line (21, posy + 1)-(32, posy + 5), pal(nn), BF
  150.         posy = posy + 12
  151.     Next
  152.  
  153. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  154.     Dim Radius As Long, RadiusError As Long
  155.     Dim X As Long, Y As Long
  156.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  157.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  158.     Line (CX - X, CY)-(CX + X, CY), C, BF
  159.     While X > Y
  160.         RadiusError = RadiusError + Y * 2 + 1
  161.         If RadiusError >= 0 Then
  162.             If X <> Y + 1 Then
  163.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  164.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  165.             End If
  166.             X = X - 1
  167.             RadiusError = RadiusError - X * 2
  168.         End If
  169.         Y = Y + 1
  170.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  171.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  172.     Wend
  173.  
  174.  

Contour Plot 3S.PNG
* Contour Plot 3S.PNG (Filesize: 53.34 KB, Dimensions: 803x426, Views: 93)

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #12 on: November 19, 2021, 09:38:05 pm »
Back to Contour #2 plotting, here is 3D Charting of Height Plane for contour map generated:
Code: QB64: [Select]
  1. _Title "3D Chart Height Plane: press any for New Chart" ' b+ 2021-11-19
  2. ' merge cube it with
  3. ' _Title "Contour Plot 2" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  4. ' Contour plot using Data Points by Mrwhy
  5. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  6. ' AndyA ref:
  7. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  8. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  9. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  10.  
  11. Const Xmax = 800, Ymax = 600
  12. ' normal x, y axis plane starts back and workes forward as 45 degrees with z height
  13. Const fSide = 15.73
  14. Const lSide = fSide * Sqr(2) / 3.9
  15. Const yOff = .67 * Ymax
  16. Dim Shared As _Unsigned Long fFace, sFace, tFace
  17. fFace = _RGB32(80, 0, 0)
  18. sFace = _RGB32(255, 40, 0)
  19.  
  20.  
  21. Screen _NewImage(Xmax, Ymax, 32)
  22. _ScreenMove 200, 60
  23.  
  24. 'blue green
  25. For i = 1 To 16
  26.     If i > 5 Then r = (i - 5) * 20 Else r = 0
  27.     pal(i) = _RGB32(r, i * 16, 5 * (14 - i))
  28. '============================================
  29. ' number of height/potential points in data
  30. '============================================
  31. stepper = 50 ' 40, 80, 160    20 is going to
  32. nP = (.5 * Xmax / stepper) * (.5 * Ymax / stepper)
  33. '============================================
  34. ReDim x(nP), y(nP), z(nP), E(nP), h(.5 * Xmax, .5 * Ymax)
  35. xoo = .25 * Xmax: yoo = 50 * Ymax
  36. restart:
  37. ReDim h(.5 * Xmax, .5 * Ymax)
  38.  
  39. i = 0
  40. For yy = stepper / 2 To .5 * Ymax Step stepper ' spread points evenly
  41.     For xx = stepper / 2 To .5 * Xmax Step stepper
  42.         i = i + 1
  43.         x(i) = xx: y(i) = yy: z(i) = Int(Rnd * 12 + 2)
  44.     Next
  45. 'Print nP, i     'check
  46. 'End
  47.  
  48. ztot = 0
  49. 'display height or potential data points on screen
  50. For i = 1 To nP
  51.     'fcirc x(i), .5*Ymax - y(i) , 3, pal(z(i))
  52.     ztot = ztot + z(i)
  53.  
  54. Color &HFFFFFF40
  55.  
  56. s$ = "Calculating Next Contour Map"
  57. x = (_Width - _PrintWidth(s$)) / 2
  58. _PrintString (x, _Height - 14), s$
  59. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  60.  
  61.  
  62. 'initialize some variables
  63. zmean = ztot / nP
  64. wo = .5 * Xmax * .5 * Ymax / nP
  65. w = wo / 2
  66.  
  67. 'generate initial Error estimates
  68. For i = 1 To nP
  69.     E(i) = z(i) - zmean
  70. Legend
  71.  
  72. 'begin Relaxation (iterative method)
  73. For jj = 1 To 9 * nP 'find max error point  >> 9 to 16?
  74.     emax = 0
  75.     For i = 1 To nP
  76.         If Abs(E(i)) > emax Then
  77.             emax = Abs(E(i))
  78.             ii = i
  79.         End If
  80.         k = E(ii)
  81.     Next
  82.     'fixit
  83.     For i = 1 To nP
  84.         dx = x(i) - x(ii)
  85.         dy = y(i) - y(ii)
  86.         dsq = dx * dx + dy * dy
  87.         E(i) = E(i) - k * Exp(-(dsq / w))
  88.         If i = ii Then
  89.             'update map with revised height or potential estimates for each pixel
  90.             For fy = 1 To (.5 * Ymax - 1)
  91.                 For fx = 1 To (.5 * Xmax - 1)
  92.                     dx = fx - x(ii)
  93.                     dy = fy - y(ii)
  94.                     dsq2 = dx * dx + dy * dy
  95.                     dy = .5 * Ymax - fy
  96.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  97.                 Next
  98.             Next
  99.         End If
  100.     Next
  101. 'Draw calculated contour map
  102. Color pal(h(1, 1) + zmean)
  103. 'Line (200, 10)-(.5 * Xmax - 1, .5 * Ymax - 1), , B
  104. 'Line (0, 0)-(0, .5 * Ymax - 1)
  105.  
  106. For fy = 1 To .5 * Ymax - 1
  107.     For fx = 1 To .5 * Xmax - 1
  108.         If Int(h(fx, fy) + zmean) > 15 Then
  109.             c~& = pal(15)
  110.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  111.             c~& = pal(1)
  112.         Else
  113.             c~& = pal(Int(h(fx, fy) + zmean))
  114.         End If
  115.         PSet (fx + 200, fy + 10), c~&
  116.     Next
  117.  
  118. 'display height or potential data points on contour map
  119. For i = 1 To nP
  120.     'fcirc x(i), Ymax - y(i), 3, &HFF000000
  121.     fcirc x(i) + 200, .5 * Ymax - y(i) + 10, 2, pal(z(i))
  122. Legend
  123. For r = 0 To 29
  124.     For c = 39 To 0 Step -1
  125.         fx = c * 10 + 5: fy = r * 10 + 5
  126.         If Int(h(fx, fy) + zmean) > 15 Then
  127.             ch = 15: c~& = pal(15)
  128.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  129.             ch = 1: c~& = pal(1)
  130.         Else
  131.             ch = Int(h(fx, fy) + zmean): c~& = pal(Int(h(fx, fy) + zmean))
  132.         End If
  133.         cube c, r, 5 * ch, c~&
  134.     Next
  135. GoTo restart
  136.  
  137. Sub cube (col, row, height, tFace As _Unsigned Long)
  138.     x0 = (row + 1) * lSide + col * fSide ' bottom left corner of front face
  139.     y0 = yOff + (row + 1) * lSide
  140.     x1 = x0 '                               left corner above
  141.     y1 = y0 - height
  142.     x2 = x0 + fSide '                    right bottom corner
  143.     y2 = y0
  144.     x4 = x2 '                      right top corner of front face
  145.     y4 = y1
  146.     Line (x1, y1)-(x2, y2), fFace, BF
  147.  
  148.     ' side face
  149.     x1 = x0 - lSide
  150.     y1 = y0 - lSide
  151.     x2 = x1
  152.     y2 = y1 - height
  153.     x3 = x0
  154.     y3 = y0 - height
  155.     fquad x0, y0, x1, y1, x2, y2, x3, y3, sFace
  156.  
  157.     ' top face
  158.     x0 = x3
  159.     y0 = y3
  160.     x1 = x2
  161.     y1 = y2
  162.     x2 = x1 + fSide
  163.     y2 = y1
  164.     fquad x0, y0, x1, y1, x2, y2, x4, y4, tFace
  165.  
  166. Sub ftri (x1, y1, x2, y2, x3, y3, K As _Unsigned Long)
  167.     Dim D As Long
  168.     Static a&
  169.     D = _Dest
  170.     If a& = 0 Then a& = _NewImage(1, 1, 32)
  171.     _Dest a&
  172.     PSet (0, 0), K
  173.     _Dest D
  174.     _MapTriangle _Seamless(0, 0)-(0, 0)-(0, 0), a& To(x1, y1)-(x2, y2)-(x3, y3)
  175.  
  176. 'need 4 non linear points (not all on 1 line) list then clockwise so x2, y2 is opposite of x4, y4
  177. Sub fquad (x1, y1, x2, y2, x3, y3, x4, y4, K As _Unsigned Long)
  178.     ftri x1, y1, x2, y2, x4, y4, K
  179.     ftri x3, y3, x2, y2, x4, y4, K
  180.  
  181. 'show height or potential value for each color in map
  182. Sub Legend
  183.     posy = 20: xoff = _Width - 100: yoo = 50
  184.     ' just draw a balck box where legend is going
  185.     Line (xoff, yoo)-(xoff + 37, yoo + 200), &HFFFF0000, BF
  186.     Line (xoff + 1, yoo + 1)-(xoff + 36, yoo + 199), &HFFFFFF00, B
  187.     Color &HFFFFFF00
  188.     _Font 8
  189.     _PrintString (xoff + 7, yoo + 3), "Key"
  190.     For nn = 1 To 15
  191.         _PrintString (xoff + 4, yoo + posy), Right$(" " + _Trim$(Str$(nn)), 2)
  192.         Line (xoff + 21, yoo + posy + 1)-(xoff + 32, yoo + posy + 5), pal(nn), BF
  193.         posy = posy + 12
  194.     Next
  195.  
  196. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  197.     Dim Radius As Long, RadiusError As Long
  198.     Dim X As Long, Y As Long
  199.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  200.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  201.     Line (CX - X, CY)-(CX + X, CY), C, BF
  202.     While X > Y
  203.         RadiusError = RadiusError + Y * 2 + 1
  204.         If RadiusError >= 0 Then
  205.             If X <> Y + 1 Then
  206.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  207.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  208.             End If
  209.             X = X - 1
  210.             RadiusError = RadiusError - X * 2
  211.         End If
  212.         Y = Y + 1
  213.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  214.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  215.     Wend
  216.  
  217.  

 
3D Chart Height Plane.PNG



Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Contour Plot
« Reply #13 on: November 20, 2021, 12:02:50 pm »
OK I halved the cube size a couple of times:
Code: QB64: [Select]
  1. _Title "3D Chart Height Plane #2: press any for New Chart" ' b+ 2021-11-19
  2. ' merge cube it with
  3. ' _Title "Contour Plot 2" 'b+ trans from SdlBasic to QB64, original from [abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]
  4. ' Contour plot using Data Points by Mrwhy
  5. 'ref: SdlBasic forum, AndyA, 9-13-2016  http://sdlbasic.epizy.com/showthread.php?tid=302
  6. ' AndyA ref:
  7. 'http://www.[abandoned, outdated and now likely malicious qb64 dot net website - don’t go there]/forum/index.php?topic=3714.msg37019#msg37019
  8. 'https://en.wikipedia.org/wiki/Richard_V._Southwell
  9. 'https://en.wikipedia.org/wiki/Relaxation_(iterative_method)
  10.  
  11. Const Xmax = 800, Ymax = 600
  12. ' normal x, y axis plane starts back and workes forward as 45 degrees with z height
  13. Const fSide = 15.73 / 8 ' half
  14. Const lSide = fSide * Sqr(2) / 3.9
  15. Const yOff = .67 * Ymax
  16. Dim Shared As _Unsigned Long fFace, sFace
  17. fFace = _RGB32(80, 0, 0)
  18. sFace = _RGB32(255, 40, 0)
  19.  
  20.  
  21. Screen _NewImage(Xmax, Ymax, 32)
  22. _ScreenMove 200, 60
  23.  
  24. 'blue green
  25. For i = 1 To 16
  26.     If i > 5 Then r = (i - 5) * 20 Else r = 0
  27.     pal(i) = _RGB32(r, i * 16, 5 * (14 - i))
  28. '============================================
  29. ' number of height/potential points in data
  30. '============================================
  31. stepper = 50 ' 40, 80, 160    20 is going to
  32. nP = (.5 * Xmax / stepper) * (.5 * Ymax / stepper)
  33. '============================================
  34. ReDim x(nP), y(nP), z(nP), E(nP), h(.5 * Xmax, .5 * Ymax)
  35. xoo = .25 * Xmax: yoo = 50 * Ymax
  36. restart:
  37. ReDim h(.5 * Xmax, .5 * Ymax)
  38.  
  39. i = 0
  40. For yy = stepper / 2 To .5 * Ymax Step stepper ' spread points evenly
  41.     For xx = stepper / 2 To .5 * Xmax Step stepper
  42.         i = i + 1
  43.         x(i) = xx: y(i) = yy: z(i) = Int(Rnd * 12 + 2)
  44.     Next
  45. 'Print nP, i     'check
  46. 'End
  47.  
  48. ztot = 0
  49. 'display height or potential data points on screen
  50. For i = 1 To nP
  51.     'fcirc x(i), .5*Ymax - y(i) , 3, pal(z(i))
  52.     ztot = ztot + z(i)
  53.  
  54. Color &HFFFFFF40
  55.  
  56. s$ = "Calculating Next Contour Map"
  57. x = (_Width - _PrintWidth(s$)) / 2
  58. _PrintString (x, _Height - 14), s$
  59. Line (x - 4, _Height - 18)-(x + _PrintWidth(s$) + 2, _Height - 2), , B
  60.  
  61.  
  62. 'initialize some variables
  63. zmean = ztot / nP
  64. wo = .5 * Xmax * .5 * Ymax / nP
  65. w = wo / 2
  66.  
  67. 'generate initial Error estimates
  68. For i = 1 To nP
  69.     E(i) = z(i) - zmean
  70. Legend
  71.  
  72. 'begin Relaxation (iterative method)
  73. For jj = 1 To 9 * nP 'find max error point  >> 9 to 16?
  74.     emax = 0
  75.     For i = 1 To nP
  76.         If Abs(E(i)) > emax Then
  77.             emax = Abs(E(i))
  78.             ii = i
  79.         End If
  80.         k = E(ii)
  81.     Next
  82.     'fixit
  83.     For i = 1 To nP
  84.         dx = x(i) - x(ii)
  85.         dy = y(i) - y(ii)
  86.         dsq = dx * dx + dy * dy
  87.         E(i) = E(i) - k * Exp(-(dsq / w))
  88.         If i = ii Then
  89.             'update map with revised height or potential estimates for each pixel
  90.             For fy = 1 To (.5 * Ymax - 1)
  91.                 For fx = 1 To (.5 * Xmax - 1)
  92.                     dx = fx - x(ii)
  93.                     dy = fy - y(ii)
  94.                     dsq2 = dx * dx + dy * dy
  95.                     dy = .5 * Ymax - fy
  96.                     h(fx, dy) = h(fx, dy) + k * Exp(-(dsq2 / w))
  97.                 Next
  98.             Next
  99.         End If
  100.     Next
  101. 'Draw calculated contour map
  102. Color pal(h(1, 1) + zmean)
  103. 'Line (200, 10)-(.5 * Xmax - 1, .5 * Ymax - 1), , B
  104. 'Line (0, 0)-(0, .5 * Ymax - 1)
  105.  
  106. For fy = 1 To .5 * Ymax - 1
  107.     For fx = 1 To .5 * Xmax - 1
  108.         If Int(h(fx, fy) + zmean) > 15 Then
  109.             c~& = pal(15)
  110.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  111.             c~& = pal(1)
  112.         Else
  113.             c~& = pal(Int(h(fx, fy) + zmean))
  114.         End If
  115.         PSet (fx + 200, fy + 10), c~&
  116.     Next
  117.  
  118. 'display height or potential data points on contour map
  119. For i = 1 To nP
  120.     'fcirc x(i), Ymax - y(i), 3, &HFF000000
  121.     fcirc x(i) + 200, .5 * Ymax - y(i) + 10, 2, pal(z(i))
  122. Legend
  123. For r = 0 To 239
  124.     For c = 319 To 0 Step -1
  125.         fx = c * 1.25 + .625: fy = r * 1.25 + .625
  126.         If Int(h(fx, fy) + zmean) > 15 Then
  127.             ch = 15: c~& = pal(15)
  128.         ElseIf Int(h(fx, fy) + zmean) < 1 Then
  129.             ch = 1: c~& = pal(1)
  130.         Else
  131.             ch = Int(h(fx, fy) + zmean): c~& = pal(Int(h(fx, fy) + zmean))
  132.         End If
  133.         cube c, r, 1.25 * ch, c~&
  134.     Next
  135. GoTo restart
  136.  
  137. Sub cube (col, row, height, tFace As _Unsigned Long)
  138.     x0 = (row + 1) * lSide + col * fSide ' bottom left corner of front face
  139.     y0 = yOff + (row + 1) * lSide
  140.     x1 = x0 '                               left corner above
  141.     y1 = y0 - height
  142.     x2 = x0 + fSide '                    right bottom corner
  143.     y2 = y0
  144.     x4 = x2 '                      right top corner of front face
  145.     y4 = y1
  146.     Line (x1, y1)-(x2, y2), fFace, BF
  147.  
  148.     ' side face
  149.     x1 = x0 - lSide
  150.     y1 = y0 - lSide
  151.     x2 = x1
  152.     y2 = y1 - height
  153.     x3 = x0
  154.     y3 = y0 - height
  155.     fquad x0, y0, x1, y1, x2, y2, x3, y3, sFace
  156.  
  157.     ' top face
  158.     x0 = x3
  159.     y0 = y3
  160.     x1 = x2
  161.     y1 = y2
  162.     x2 = x1 + fSide
  163.     y2 = y1
  164.     fquad x0, y0, x1, y1, x2, y2, x4, y4, tFace
  165.  
  166. Sub ftri (x1, y1, x2, y2, x3, y3, K As _Unsigned Long)
  167.     Dim D As Long
  168.     Static a&
  169.     D = _Dest
  170.     If a& = 0 Then a& = _NewImage(1, 1, 32)
  171.     _Dest a&
  172.     PSet (0, 0), K
  173.     _Dest D
  174.     _MapTriangle _Seamless(0, 0)-(0, 0)-(0, 0), a& To(x1, y1)-(x2, y2)-(x3, y3)
  175.  
  176. 'need 4 non linear points (not all on 1 line) list then clockwise so x2, y2 is opposite of x4, y4
  177. Sub fquad (x1, y1, x2, y2, x3, y3, x4, y4, K As _Unsigned Long)
  178.     ftri x1, y1, x2, y2, x4, y4, K
  179.     ftri x3, y3, x2, y2, x4, y4, K
  180.  
  181. 'show height or potential value for each color in map
  182. Sub Legend
  183.     posy = 20: xoff = _Width - 100: yoo = 50
  184.     ' just draw a balck box where legend is going
  185.     Line (xoff, yoo)-(xoff + 37, yoo + 200), &HFFFF0000, BF
  186.     Line (xoff + 1, yoo + 1)-(xoff + 36, yoo + 199), &HFFFFFF00, B
  187.     Color &HFFFFFF00
  188.     _Font 8
  189.     _PrintString (xoff + 7, yoo + 3), "Key"
  190.     For nn = 1 To 15
  191.         _PrintString (xoff + 4, yoo + posy), Right$(" " + _Trim$(Str$(nn)), 2)
  192.         Line (xoff + 21, yoo + posy + 1)-(xoff + 32, yoo + posy + 5), pal(nn), BF
  193.         posy = posy + 12
  194.     Next
  195.  
  196. Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
  197.     Dim Radius As Long, RadiusError As Long
  198.     Dim X As Long, Y As Long
  199.     Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
  200.     If Radius = 0 Then PSet (CX, CY), C: Exit Sub
  201.     Line (CX - X, CY)-(CX + X, CY), C, BF
  202.     While X > Y
  203.         RadiusError = RadiusError + Y * 2 + 1
  204.         If RadiusError >= 0 Then
  205.             If X <> Y + 1 Then
  206.                 Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
  207.                 Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
  208.             End If
  209.             X = X - 1
  210.             RadiusError = RadiusError - X * 2
  211.         End If
  212.         Y = Y + 1
  213.         Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
  214.         Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
  215.     Wend
  216.  
  217.  
3D Chart Height #2.PNG


Offline Richard Frost

  • Seasoned Forum Regular
  • Posts: 316
  • Needle nardle noo. - Peter Sellers
    • View Profile
Re: Contour Plot
« Reply #14 on: November 20, 2021, 03:54:19 pm »
I've wanted code to do that for over 20 years.  I did ONE hill in a sad, 2D fashion.

Lovely stuff, slightly marred by using FONT 8, so I modified your code to use Liberati 9.

My feeble attempt:

Code: QB64: [Select]
  1. Dim distance(360), elevation(360), active(20), angle(20)
  2. Screen _NewImage(800, 600, 12)
  3. xc = _Width \ 2: yc = _Height \ 2
  4. asp! = yc / xc
  5.  
  6.     Cls
  7.     mind = 9999: maxd = -mind
  8.     mine = 9999: maxe = -mine
  9.  
  10.     distance(0) = 0
  11.     e0 = 300 + Rnd * 20
  12.     active(0) = 360
  13.     elevation(0) = e0
  14.     elevation(360) = e0
  15.     maxe = e0
  16.  
  17.     n = 8
  18.     genang:
  19.     For i = 1 To n
  20.         angle(i) = Rnd * 360
  21.     Next i
  22.     sort:
  23.     sorted = 1
  24.     For i = 1 To n - 1
  25.         If angle(i) > angle(i + 1) Then sorted = 0: Swap angle(i), angle(i + 1)
  26.     Next i
  27.     If sorted = 0 Then GoTo sort
  28.     For i = 1 To n - 1
  29.         If (angle(i + 1) - angle(i)) < 20 Then GoTo genang
  30.     Next i
  31.  
  32.     For i = 1 To n
  33.         angle = angle(i)
  34.         angle = (angle + 360) Mod 360
  35.         active(i) = angle Mod 360
  36.         distance = 50 + Rnd * 250
  37.         elevation = 100 + Rnd * 150
  38.         If distance > maxd Then maxd = distance
  39.         If distance < mind Then mind = distance
  40.         If elevation > maxe Then maxe = elevation
  41.         If elevation < mine Then mine = elevation
  42.         distance(angle) = distance
  43.         elevation(angle) = elevation
  44.     Next i
  45.     xm = xc / (maxd - mind)
  46.     ym = yc / (maxd - mind)
  47.     active(n + 1) = active(1) '                            repeat 1st coord
  48.     distance(active(n + 1)) = distance(active(1))
  49.     elevation(active(n + 1)) = elevation(active(1))
  50.  
  51.     For i = 1 To n + 1 '                                   approx missing data each degree
  52.         angle1 = active(i - 1)
  53.         angle2 = active(i)
  54.         ddif = distance(angle2) - distance(angle1)
  55.         edif = elevation(angle2) - elevation(angle1)
  56.         If i = (n + 1) Then angle2 = angle2 + 360
  57.         a = 0: ai = 90 / (angle2 - angle1) '               angle, angle increment
  58.         For z = Int(angle1) To angle2
  59.             na = z Mod 360 '                               ensure array bounds not exceeded
  60.             q = Sin(_D2R(a))
  61.             q = q * q '                                    smoothing
  62.             a = a + ai
  63.             distance(na) = distance(angle1) + ddif * q
  64.             elevation(na) = elevation(angle1) + edif * q
  65.         Next z
  66.     Next i
  67.     EndinEl = Int(maxe / q) * q
  68.     For el = -300 To EndinEl Step 5
  69.         col = 1
  70.         If (el Mod 20) = 0 Then col = 9
  71.         If (el Mod 40) = 0 Then col = 4
  72.         pf = 0 '                                           plot flag (move/draw)
  73.         For mangle = 0 To 360
  74.             angle = mangle Mod 360
  75.             distance = distance(angle)
  76.             elevation = elevation(angle)
  77.             epf = distance / (e0 - elevation)
  78.             d = distance - ((el - elevation) * epf)
  79.             x = xc + d * xm * Cos(_D2R(angle)) * asp!
  80.             y = yc - d * ym * Sin(_D2R(angle))
  81.             If pf = 1 Then Line -(x, y), col Else pf = 1: PSet (x, y), col
  82.         Next mangle
  83.     Next el
  84.  
  85.     For i = 0 To n
  86.         d = 0: e = e0
  87.         If i Then d = distance(i): e = elevation(angle(i))
  88.         x = xc + d * Cos(_D2R(angle(i))) * asp!
  89.         y = yc - d * Sin(_D2R(angle(i)))
  90.         Circle (x, y), 1, 15
  91.         _PrintString (x, y + 5), Str$(Int(e))
  92.     Next i
  93.     Sleep
  94.  
« Last Edit: November 20, 2021, 04:08:38 pm by Richard Frost »
It works better if you plug it in.