Tweet
This routine produces topographic map style plots from a 2D array representing
measured or calculated values on a regular rectangular grid.
Contours are obtained by inserting contrasting colours in the color table at
a specified interval.
Although the resulting contours are sometimes not perfect lines, this simple
technique may be useful for many practical purposes.
The routine was written for PBCC.
measured or calculated values on a regular rectangular grid.
Contours are obtained by inserting contrasting colours in the color table at
a specified interval.
Although the resulting contours are sometimes not perfect lines, this simple
technique may be useful for many practical purposes.
The routine was written for PBCC.
Code:
' --------------------------------------------------------------------------- ' SHADED / CONTOUR PLOTS FROM 2D ARRAYS ' for PBCC ' Arie Verheul - october 2008 ' --------------------------------------------------------------------------- ' This routine produces topographic map style plots from a 2D array representing ' measured or calculated values on a regular rectangular grid. ' Contours are obtained by inserting contrasting colours in the color table at ' a specified interval. ' Although the resulting contours are sometimes not perfect lines, this simple ' technique may be useful for many practical purposes. ' The interpolated array ScreenArray() may be used to retrieve local values. ' For clarity this option, and some control options, were omitted from the code. ' The color range may be limited according to a specific interest. ' Therefore the color scale is extended at both sides with a dark green overflow color. ' The "zero" contour is dark gray. ' To use the plot routine in a different context requires the Subs DetermineExtremes, ' MapArrayToWindow and ColorPlot, plus the Type variable ContourParmset. ' --------------------------------------------------------------------------- ' For those interested in details about the example used here: ' Finite difference is a simple algorithm to (approximately) solve the Laplacian. ' It uses a 2D grid of points, usually representing some aspect of a physical structure. ' If grid spacing in both directions is equal, the algorithm involves repeatedly setting ' each gridpoint (except for boundary points) at the average of its four neighbours. ' The example uses a grid of 101 x 101 gridpoints. ' The plot shows an example of "field-mapping", a technique to use the field ' from point or line sources (which can be easily calcuated) to model something else. ' The idea behind this is that each equipotential contour may be interpreted ' as a conducting surface (or an electrode). ' Therefore, the calculated 2D field from five line sources perpendicular to the ' plane may as well represent the field from five rods. ' --------------------------------------------------------------------------- #Compiler PBCC 4.04 #Compile Exe #Dim All #Console Off ' --------------------------------------------------------------------------- Type ContourParmset NumColor As Long ' Number of shades in positive and negative half of range PlotRange As Single ' Maximum value to be rendered (positive value) ContInt As Single ' Contour interval, in same units as PlotRange (positive value) ArrayMax As Single ' Maximum and minimum values, determined internally ArrayMin As Single End Type ' --------------------------------------------------------------------------- Declare Sub FiniteDifference (PlotArray() As Single) Declare Sub DetermineExtremes (PlotArray() As Single, PlotParms As ContourParmset) Declare Sub MapArrayToWindow (PlotArray() As Single, ScreenArray() As Single) Declare Sub ColorPlot (ScreenArray() As Single, PlotParms As ContourParmset) Declare Function IsWindow Lib "USER32.DLL" Alias "IsWindow" (ByVal hWnd As Dword) As Long ' --------------------------------------------------------------------------- Function PBMain () As Long Local hWnd As Dword Dim ScreenArray () As Single Dim PlotArray() As Single Dim PlotParms As ContourParmset ' Setup window Graphic Window "Electric potential plot",200,100,600,600 To hWnd Graphic Attach hWnd, 0 Graphic Color %White, %Black Graphic Clear ' Plot settings (may be changed) PlotParms.NumColor = 250 PlotParms.PlotRange = 500 PlotParms.ContInt = 50 ' Main program FiniteDifference PlotArray() ' Produce some data to be plotted DetermineExtremes PlotArray(), PlotParms ' Determine Max/Min values MapArrayToWindow PlotArray(), ScreenArray() ' Calculate interpolated values for each pixel ColorPlot ScreenArray(), PlotParms ' Translate values to shades While IsWindow(hWnd): Sleep 1000: Wend ' Terminate program after window closure End Function ' --------------------------------------------------------------------------- Sub FiniteDifference (PlotArray() As Single) ' The purpose of this Sub is just to produce some data to be plotted ' It must not be considered as a good or general example to use this algorithm %NumInt = 100 Local N, M, I As Long Local XS1, XS2, XS3, YS1, YS2, YS3 As Long ReDim PlotArray(%NumInt,%NumInt) As Single ' --------------------------------------------------------------------------- ' Set internal boundaries XS1 = .2 * %NumInt ' Location of internal boundary points XS2 = %NumInt - XS1 YS1 = .2 * %NumInt YS2 = %NumInt - YS1 XS3 = %NumInt / 2 YS3 = %NumInt / 2 PlotArray(XS1,YS1) = 1250 ' Values assigned to internal boundary points PlotArray(XS2,YS2) = 1250 PlotArray(XS1,YS2) = - 1375 PlotArray(XS2,YS1) = - 1375 PlotArray(XS3,YS3) = 1125 ' --------------------------------------------------------------------------- For I = 1 To 2000 ' Iterations For M = 1 To %NumInt - 1 ' Vert. direction For N = 1 To %NumInt - 1 ' Hor. direction If M = YS1 Or M = YS2 Then If N = XS1 Or N = XS2 Then Iterate For ' Skip internal boundary points End If If M = YS3 And N = XS3 Then Iterate For PlotArray(N,M) = (PlotArray(N-1,M) + PlotArray(N+1,M)_ + PlotArray(N,M-1) + PlotArray(N,M+1)) / 4 Next Next Next End Sub ' --------------------------------------------------------------------------- Sub DetermineExtremes (PlotArray() As Single, PlotParms As ContourParmset) ' Determines Maximum and Minimum value in Plotarray() Local M, N As Long Local TestValue, ArrayMax, ArrayMin As Single ArrayMin = PlotArray(0,0) ArrayMax = ArrayMin For M = 0 To UBound(PlotArray, 2) For N = 0 To UBound(PlotArray, 1) TestValue = PlotArray(N,M) 'Skip poles, if present If (TestValue > 1E35) Or (TestValue < -1E35) Then Iterate For ArrayMax = Max(TestValue, ArrayMax) ArrayMin = Min(TestValue, ArrayMin) Next Next PlotParms.ArrayMax = ArrayMax PlotParms.ArrayMin = ArrayMin End Sub ' --------------------------------------------------------------------------- Sub MapArrayToWindow (PlotArray() As Single, ScreenArray() As Single) ' This sub calculates linearly interpolated values from PlotArray() ' for each pixel in ScreenArray() ' Note: in the following the rectangular area enclosed by ' four neighbouring gridpoints is designated as an element Local WindowWidth, WindowHeight As Long Local Xmax, Ymax As Long Local I, ElementRow, Pix, PixRow, CurrentRow As Long Local X, PixRowPos, VertIntFactor As Single Graphic Get Client To WindowWidth, WindowHeight Xmax = UBound(PlotArray(),1) Ymax = UBound(PlotArray(),2) Dim TopRow(Xmax) As Single ' Top nodes of current element row Dim BotRow(Xmax) As Single ' Bottom nodes of current element row Dim RowDiff(Xmax) As Single ' Differential (TopRow() - BotRow) Dim InterpRow(Xmax) As Single ' Interpolated values for position of pixel row Dim ElDiff(Xmax - 1) As Single ' Hor. differential in InterpRow() Dim HorIntNumber (WindowWidth) As Long ' Element to use for interpolation for specific pixel Dim HorIntFactor (WindowWidth) As Single ' Interpolation factor for specific pixel ReDim ScreenArray (1 To WindowWidth, 1 To WindowHeight) As Single ' --------------------------------------------------------------------------- ' Make tables for pixel locations and interpolation factors in element grid coordinates ' in horizontal direction; these tables are valid for each row of pixels. For Pix = 1 To WindowWidth ' For each pixel in horizontal row X = Pix * Xmax / WindowWidth ' > Pixel position in element grid coordinates HorIntNumber (Pix) = Int (X) ' > Element where pixel is located HorIntFactor (Pix) = Frac (X) ' > Interpolation factor within element Next ' --------------------------------------------------------------------------- ' Do the actual interpolation, working row by row CurrentRow = -1 ' Counter to keep track of current element row For PixRow = 1 To WindowHeight ' Work row by row PixRowPos = Ymax * PixRow / WindowHeight ' Vertical position of pixel row ' in element grid coordinates ElementRow = Int (PixRowPos) ' Element row to use for interpolation VertIntFactor = Frac (PixRowPos) ' Interpolation factor If ElementRow <> CurrentRow Then ' If interpolation has reached next ' element row, data must be refreshed ' Full row of data is extracted from PlotArray() by placing 1D array on top of it ' Using MAT these data are manipulated with a full row at the time ReDim TopRow(Xmax) As Single At VarPtr (PlotArray(0, ElementRow + 1)) ReDim BotRow(Xmax) As Single At VarPtr (PlotArray(0, ElementRow)) Mat RowDiff() = TopRow() - BotRow() ' Vertical differentials CurrentRow = ElementRow End If Mat InterpRow () = (VertIntFactor) * RowDiff() Mat InterpRow () = BotRow() + InterpRow() ' InterpRow() now contains an interpolated row of data from PlotArray() ' at PixRowPos, the position of the row of pixels that is treated currently. ' --------------------------------------------------------------------------- ' Next within this row the differentials are calculated, and stored in ElDiff() For I = 0 To Xmax - 1 ' Horizontal differential in InterpRow() ElDiff(I) = InterpRow(I + 1) - InterpRow(I) Next ' --------------------------------------------------------------------------- ' From InterpRow() and ElDiff() the values for each pixel are calculated For Pix = 1 To WindowWidth ' Work along a row of pixels I = HorIntNumber(Pix) ' Element where pixel is located ScreenArray(Pix, PixRow) = InterpRow(I) + ElDiff(I) * HorIntFactor(Pix) Next ' Next pixel Next ' Next row of pixels End Sub ' --------------------------------------------------------------------------- Sub ColorPlot (ScreenArray() As Single, PlotParms As ContourParmset) ' This sub translates values associated with pixels into shades Local R,G,B As Integer Local N, M, ColorIndex As Long Local ColorScaleUBound, ColorScaleLBound As Long Local WindowWidth, WindowHeight As Long Local X, NormalisedIndex, NormalisedIndexSqr As Single Local WindowMap As String Dim PixelPtr As Long Pointer ' --------------------------------------------------------------------------- ' Set default values if necessary If PlotParms.NumColor <= 0 Then PlotParms.NumColor = 300 PlotParms.PlotRange = Abs(PlotParms.PlotRange) If PlotParms.PlotRange = 0 Then PlotParms.PlotRange = Max(Abs(PlotParms.ArrayMax), Abs(PlotParms.ArrayMin)) End If PlotParms.ContInt = Abs(PlotParms.ContInt) If PlotParms.ContInt = 0 Then PlotParms.ContInt = PlotParms.PlotRange / 20 ' --------------------------------------------------------------------------- ColorScaleUBound = PlotParms.NumColor + 1 ColorScaleLBound = -PlotParms.NumColor - 1 Dim Shade(ColorScaleLBound To ColorScaleUBound) As Long ' --------------------------------------------------------------------------- ' make color table, involving the following steps: ' a. a fractional value <NormalisedIndex> is calculated from the color index. ' b. using the expression contained in ColorAlgorithm the "positive" shading. ' colors are calculated (BGR format). ' c. swapping R and B value gives "negative" shading colors. ' d. contour colors are obtained by adding .35 to <NormalisedIndex>. ' which shifts them to a brighter value. ' e. if <NormalisedIndex> > 1 then white is used for contours. ' f. "zero" contour is assigned dark gray. ' g. at both sides of the color scale an overflow color, dark green, is added. ' --------------------------------------------------------------------------- ' background colors For N = 1 To PlotParms.NumColor NormalisedIndex = N / PlotParms.NumColor GoSub ColorAlgorithm Next ' --------------------------------------------------------------------------- ' contour colors For X = 0 To PlotParms.PlotRange Step PlotParms.ContInt NormalisedIndex = X / PlotParms.PlotRange N = NormalisedIndex * PlotParms.NumColor NormalisedIndex = NormalisedIndex + .35 If NormalisedIndex <= 1 Then ' the above line may produce values > 1 GoSub ColorAlgorithm Else Shade (N) = %White ' use white instead Shade(-N) = %White End If Next ' --------------------------------------------------------------------------- Shade(0) = RGB(80,80,80) ' "zero" contour Shade( PlotParms.NumColor + 1) = RGB(0,80,0) ' overflow color Shade(-PlotParms.NumColor - 1) = RGB(0,80,0) ' --------------------------------------------------------------------------- ' GRAPHIC GET BITS / GRAPHIC SET BITS is used here because it is ' much faster than GRAPHIC SET PIXEL Graphic Get Client To WindowWidth, WindowHeight Graphic Get Bits To WindowMap ' Get copy of video memory PixelPtr = StrPtr(WindowMap) + 8 ' First two Long Int. contain size For M = WindowHeight To 1 Step -1 ' Origin at left bottom corner ' For origin at left top corner reverse scanning direction For N = 1 To WindowWidth ColorIndex = ScreenArray(N,M) * (PlotParms.NumColor/PlotParms.PlotRange) ' overflow check ColorIndex = Min(ColorIndex, ColorScaleUBound) ColorIndex = Max(ColorIndex, ColorScaleLBound) @PixelPtr = Shade(ColorIndex) ' Write color to Windowmap Incr PixelPtr Next Next Graphic Set Bits WindowMap ' Put modified WindowMap back Exit Sub ' --------------------------------------------------------------------------- ColorAlgorithm: NormalisedIndexSqr = NormalisedIndex^2 R = 600 * NormalisedIndex - 345 * NormalisedIndexSqr ' color expression G = 255 * NormalisedIndexSqr B = 0 Shade (N) = RGB(B,G,R) ' BGR color format Shade(-N) = RGB(R,G,B) ' for "negative" colors swap R and B Return End Sub