VB Snippet - Math - Shade a closed contour with parallel lines

[krtxmrtz]

Here are 2 functions (and a demo project):

InOrOut: a function that checks the position of a point relative to a closed contour (i.e. whether it lies inside, outside or on the border)
ContourShade: a subroutine that fills a contour (may be irregular) with parallel lines at an angle and separation given as parameters. It issues a call to InOrOut

I apologize because this is the translation into vb of some code I had written in Fortran a few years ago so, the style may look somewhat messy with all those go to's. But it works!

Public Function InOrOut(p, q, k, p0, q0)
'Function to determine the position of point (p0,q0) relative
'to the closed contour defined by the n points with
'coordinates p(i), q(i), i=1, 2, ..., k
'Point 1 must be the same as point k
'Function values():
'       1: point lies inside the contour
'       -1: point is outside the contour
'       0: point is on the contour border
'*********************************************************************
    Dim kross As Integer, i As Integer
    Dim pp As Single
    'Function initialization
    InOrOut = 0
    'Initialization of kross, a variable which keeps track
    'of how many times the horizontal semi-infinite straight
    'line starting at (p0,q0) and in the positive (right hand side)
    'x axis direction intercepts the contour
    '(The contour may have vertices with angles larger than 180 degrees)
    kross = 0
    'Loop over all contour sides
    For i = 1 To k - 1
        'If the side between the i and i+1 vertices lies entirely
        'above or below the point (p0,q0), then continue on
        'with the next side (there is no intercept)
        If (q(i) > q0 And q(i + 1) > q0) Or (q(i) < q0 And q(i + 1) < q0) Then GoTo NextItem
        'If the side is horizontal avoid the calculation of the intercept by interpolation
        'as there would be a division by zero
        If q(i) = q(i + 1) Then
            'It has to be determined if the point (p0,q0) lies on this horizontal segment
            If (p(i) > p0 And p(i + 1) > p0) Or (p(i) < p0 And p(i + 1) < p0) Then GoTo NextItem
            'If it doesn't, we're done!
            Exit Function
        End If
        'Calculation of pp, the coordinate of the point where the segment connecting the
        'sides i and i+1 and the horizontal straight line that goes through (p0,q0) intercept
        pp = p(i) + (q0 - q(i)) * ((p(i + 1) - p(i)) / (q(i + 1) - q(i)))
        'The sign of pp-p0 determines the position of the intercept relative to (p0,q0)
        If pp - p0 > 0 Then
            'Intercept to the right: increment counter
            kross = kross + 1
        ElseIf pp - p0 = 0 Then
            Exit Function
        End If
        'If the intercept lies to the left, take no action and continue on
NextItem:
    Next
    'End of loop, reinitialization of InOrOut
    InOrOut = 1
    'If the number of intercepts is even the point (p0,q0) lies out of the contour
    If kross Mod 2 = 0 Then InOrOut = -1
End Function
 
'*************************************************
'*************************************************
 
Public Sub ContourShade(X, Y, xc, n, sep, angle)
'It fills a closed contour defined by n points with coordinates
'x(i), y(i), i=1, 2, ..., n
'Points #1 and #n are the same, i.e. x(1)=x(n) & y(1)=y(n)
'xc is a workspace
'sep is the separation between the parallel shading lines
'angle is the angle in degrees between the lines and the x axis (horizontal)
'************************************************************************************
    Const Pi = 3.141592654
    Dim radians As Single
    Dim yLow As Single, yHigh As Single
    Dim dummyX As Single, dummyY As Single
    Dim x1 As Single, y1 As Single
    Dim Small As Single, Large As Single
    Dim v1 As Single, w1 As Single, v2 As Single, w2 As Single
    Dim i As Integer, ii As Integer, j As Integer, nc As Integer
    'Convert angle to radians
    radians = angle * Pi / 180
    'Rotate contour by angle -angle so that we will deal with horizontal parallel
    'lines, much easier to handle
    For i = 1 To n
        dummyX = X(i)
        dummyY = Y(i)
        X(i) = dummyX * Cos(radians) + dummyY * Sin(radians)
        Y(i) = -dummyX * Sin(radians) + dummyY * Cos(radians)
    Next
    'Now find the minimum and maximum y values
    yLow = Y(1)
    yHigh = Y(1)
    For i = 2 To n - 1
        If Y(i) < yLow Then yLow = Y(i)
        If Y(i) > yHigh Then yHigh = Y(i)
    Next
    'Initialize y1, the vertical coordinate of the line to be drawn
    y1 = yLow + sep
    'Loop over the lines
NextLine:
    'Initialize nc, the number of points where the line intercepts the contour
    nc = 0
    'Loop over the contour vertices. For each one of them, may or may not cross it.
    'If it does, the point where this occurs is stored in the workspace xc and nc is incremented
    For i = 1 To n - 1
        'Skip horizontal sides
        If Y(i) = Y(i + 1) Then GoTo NextItem
        'Calculate x1, the horizontal coordinate where the line intercepts the contour
        'side being considered
        x1 = X(i) + (y1 - Y(i)) * (X(i + 1) - X(i)) / (Y(i + 1) - Y(i))
        'Calculate xLarge and xSmall, the larger and smaller of the pair x(i) & x(i+1)
        If X(i) < X(i + 1) Then
            Small = X(i)
            Large = X(i + 1)
        Else
            Small = X(i + 1)
            Large = X(i)
        End If
        'Handle the vertical sides as a special case
        If X(i) = X(i + 1) Then
            'First calculate which of y(i) and y(i+1) is the larger and the smaller
            If Y(i) < Y(i + 1) Then
                Small = Y(i)
                Large = Y(i + 1)
            Else
                Small = Y(i + 1)
                Large = Y(i)
            End If
            'If the intercept point lies out of the side, go for the next side
            If Small > x1 Or Large < x1 Then GoTo NextItem
            'If it lies completely within the side, count it in as a real intercept
            'point and store it
            If (x1 <> X(i) Or y1 <> Y(i)) And (x1 <> X(i + 1) Or y1 <> Y(i + 1)) Then GoTo Increment
        Else
            'If the intercept lies out of the side, go to the next one
            If Small > x1 Or Large < x1 Then GoTo NextItem
            'If it lies completely within the side, count it really as an intercept and store it
            If (x1 <> X(i) Or y1 <> Y(i)) And (x1 <> X(i + 1) Or y1 <> Y(i + 1)) Then GoTo Increment
            'If it coincides with one of the vertices go to take special action
        End If
        'Now for the special action: we want to check whether both sides
        '(i, i+1) and (i+1, i+2) are above or below the horizontal line
        'at the same time
        ii = i + 2
        'Avoid having it out of range when i=n-1
        If i = n - 1 Then ii = 2
        'If on the same side then it is not an intercept
        If (Y(i) - y1) * (Y(ii) - y1) >= 0 Then GoTo NextItem
        'Increment the intercept counter
Increment:
        nc = nc + 1
        'Store the intercept x coordinate
        xc(nc) = x1
NextItem:
    Next
    'Arrange the xc from samller to larger
    For i = 1 To nc - 1
        ii = i
        For j = i + 1 To nc
            If xc(ii) > xc(j) Then ii = j
        Next
        dummyX = xc(ii)
        xc(ii) = xc(i)
        xc(i) = dummyX
    Next
    'Loop over the intercepts
    'Handle them in couples
    For i = 1 To nc - 1
        'dummyX is the coordinate of a point lying midway
        'between x(i) and x(i+1)
        dummyX = 0.5 * (xc(i) + xc(i + 1))
        'Call a function to determine whether this point (dummyX,y1)
        'lies inside the contour
        'If it does, plot the line. If it doesn't, check the next couple
        If InOrOut(X, Y, n, dummyX, y1) = 1 Then
            'Rotate back the 2 adjacent points to the original angle
            v1 = xc(i) * Cos(radians) - y1 * Sin(radians)
            w1 = xc(i) * Sin(radians) + y1 * Cos(radians)
            v2 = xc(i + 1) * Cos(radians) - y1 * Sin(radians)
            w2 = xc(i + 1) * Sin(radians) + y1 * Cos(radians)
            'And finally plot a line connecting them!
            Picture1.Line (v1, w1)-(v2, w2)
        End If
    Next
    'Now go for the next horizontal line
    y1 = y1 + sep
    'If we haven't yet reached the top go back to start with the next line
    If y1 < yHigh Then GoTo NextLine
    'Else, rotate back the contour
    For i = 1 To n
        dummyX = X(i)
        dummyY = Y(i)
        X(i) = dummyX * Cos(radians) - dummyY * Sin(radians)
        Y(i) = dummyX * Sin(radians) + dummyY * Cos(radians)
    Next
    'That's all folks!
End Sub

[krtxmrtz]

An improved version of the demo project (the basic code is the same).

[krtxmrtz]

Here's a minor change in the InOrOut function above, the purpose of which is to find out the position of a point relative to a closed contour.

There is a degenerate case that I had not dealt with: if the point has the same vertical (q) coordinate as one of the vertices of the contour, then the variable "kross" is counted twice. To avoid this I have added an if statement (see changes in bold red below).

Public Function InOrOut(p, q, k, p0, q0)
'Function to determine the position of point (p0,q0) relative
'to the closed contour defined by the n points with
'coordinates p(i), q(i), i=1, 2, ..., k
'Point 1 must be the same as point k
'Function values():
'       1: point lies inside the contour
'       -1: point is outside the contour
'       0: point is on the contour border
'*********************************************************************
    Dim kross As Integer, i As Integer
    Dim pp As Single
    'Function initialization
    InOrOut = 0
    'Initialization of kross, a variable which keeps track
    'of how many times the horizontal semi-infinite straight
    'line starting at (p0,q0) and in the positive (right hand side)
    'x axis direction intercepts the contour
    '(The contour may have vertices with angles larger than 180 degrees)
    kross = 0
    'Loop over all contour sides
    For i = 1 To k - 1
        'If the side between the i and i+1 vertices lies entirely
        'above or below the point (p0,q0), then continue on
        'with the next side (there is no intercept)
        If (q(i) > q0 And q(i + 1) > q0) Or (q(i) < q0 And q(i + 1) < q0) Then GoTo NextItem
        'If the side is horizontal avoid the calculation of the intercept by interpolation
        'as there would be a division by zero
        If q(i) = q(i + 1) Then
            'It has to be determined if the point (p0,q0) lies on this horizontal segment
            If (p(i) > p0 And p(i + 1) > p0) Or (p(i) < p0 And p(i + 1) < p0) Then GoTo NextItem
            'If it doesn't, we're done!
            Exit Function
        End If
        'Calculation of pp, the coordinate of the point where the segment connecting the
        'sides i and i+1 and the horizontal straight line that goes through (p0,q0) intercept
        pp = p(i) + (q0 - q(i)) * ((p(i + 1) - p(i)) / (q(i + 1) - q(i)))
        'The sign of pp-p0 determines the position of the intercept relative to (p0,q0)
        If pp - p0 > 0 Then

'Intercept to the right: increment counter (but only if the intercept does not lie
'on the first vertex, to avoid counting the same intercept twice!)
If q0 <> q(i) Then kross = kross + 1

ElseIf pp - p0 = 0 Then
            Exit Function
        End If
        'If the intercept lies to the left, take no action and continue on
NextItem:
    Next
    'If the number of intercepts is even, the point (p0,q0) lies out of the contour
    'If it's odd then it lies inside
    If kross Mod 2 = 0 Then
        InOrOut = -1
    Else
        InOrOut = 1
    End If
End Function

[krtxmrtz]

(Updated demo project)

[weicheng83]

i tried out your contourfill. but i hope that you can help me with this. cause the contour that i need to draw out is with 8 points and it will form a circle.

[krtxmrtz]

i tried out your contourfill. but i hope that you can help me with this. cause the contour that i need to draw out is with 8 points and it will form a circle.

Well, I hoped the demo project was sufficiently explicit, all you have to do is copy the code... In the demo itself you can change the points' coordinates in the Form_Load subroutine and place your eight points. Actually you must add an extra 9th point with the same coordinates as the first, for the purpose of closing the contour.
Otherwise let me know exactly what your problem is...