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