Python Generating a Shadowmap from a Heightmap

Question

I've recently been experimenting with Python and PIL to generate heightmap images using simplex noise. The generated heightmaps are simply 2D numPy arrays with height values ranging from 0-1. I expanded upon this by implementing biomes based on the height values and thresholds, and converting them to RGB values. I want to experiment further by implementing rudimentary lighting, although I'm unsure how to implement the steps I believe need to be performed:

• Create a shadowmap of size, initialised to 1, assuming all elements are lit
• Create a sun at an arbitrary position, such as (0, 0, 2). This would be the upper left corner at double the height of the maximum value within the heightmap
• Loop through the heightmap, and for each position (P):
  • Draw a line (L) from the sun to P and obtain all elements within the heightmap that intersect L
  • Loop through each intersecting element and determine if the height of the element is greater than the height of L at the current position. If it is, then the sun is blocked, iteration can be stopped, and P within the shadowmap can be set to 0

I understand how to implement the shadows once I have a shadowmap based on the heightmap, but i'm struggling with how I would get a list of elements between the sun position and a target position within the heightmap. Alongside that, how would I determine the height of the line between the sun and target position, at an intersecting position?

I appreciate any help! Below is the script to generate heightmaps:

from PIL import Image
import random
import opensimplex as simplex
import numpy as np

#noise
#Produce a noise value depending on an amount of octaves and its corresponding amplitude via @amps, and an exponent @exp
def noise(nx, ny, amps, exp):
    #Rescale simplex noise output from -1.0:1.0 to 0.0:1.0
    def sn(nx, ny):
        return simplex.noise2(nx, ny) / 2.0 + 0.5
    
    e = ampTotal = 0

    #Iterate amplitude values for each octave
    for amp in amps:
        ampDiv = 1 / amp
        ampTotal += ampDiv

        #Sample different areas of noise to avoid correlation
        e += ampDiv * sn(nx * amp, ny * amp)

    return (e / ampTotal) ** exp

#normalise_heightmap
#Rescales noise values depending on @min and @max to the 0 to 1 range
def normalise_heightmap(min, max, size, heightmap):
    scale = max - min

    for y in range(0, size[1]):
        for x in range(0, size[0]):
            heightmap[y, x] = (heightmap[y, x] - min) / scale

    return heightmap

#generate_heightmap
#Creates a numPy array of noise values
def generate_heightmap(seed, size, freq = 4, amps = [1.0, 2.0, 4.0, 8.0], exp = 2):
    simplex.seed(seed)
    heightmap = np.zeros(size)

    min = 1
    max = 0
    freqX = size[0] / freq
    freqY = size[1] / freq

    for y in range(0, size[1]):
        for x in range(0, size[0]):
            nx = x / freqX
            ny = y / freqY
            nv = noise(nx, ny, amps, exp)
            heightmap[y, x] = nv

            if nv < min: min = nv
            if nv > max: max = nv
   
    return normalise_heightmap(min, max, size, heightmap)

#generate_image
#Creates an image from a heightmap
def generate_image(heightmap):
    img = Image.fromarray(heightmap.astype("uint8"), "L")
    img.show()

heightmap = generate_heightmap(int(random.random() * 1000), [256, 256])
generate_image(np.multiply(heightmap, 255))