Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Shaded / contour plot routine for 2D arrays

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Shaded / contour plot routine for 2D arrays

    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.


    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
Working...
X