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
