function render_bezier(beg_point, ctrl_point, end_point, output_buffer, width, height):
    // set up curve structure
    curve = Curve(beg_point, end_point, ctrl_point)

    // set up raster buffer
    buffer = allocate_buffer(width, height)
    clear_buffer(buffer)

    // tessellate and draw
    tessellate_and_draw_curve(curve, buffer, width, height)

    // convert to final image
    post_process(buffer, output_buffer, width, height)

function is_flat(beg, ctrl, end):
    // calculate twice the area of the triangle formed by the three points
    // this gives a measure of how "curvy" the curve is
    vector1 = (ctrl.x - beg.x, ctrl.y - beg.y)
    vector2 = (end.x - beg.x, end.y - beg.y)
    area2 = abs(vector1.x * vector2.y - vector1.y * vector2.x)

    // if area is small enough, the curve can be approximated as a straight line
    return area2 <= 2.0  // threshold value

function tessellate_and_draw_curve(curve, buffer, width, height):
    stack = empty stack
    stack.push(curve)

    while stack is not empty:
        current_curve = stack.pop()

        if is_flat(points[current_curve.beg],
                  points[current_curve.ctrl],
                  points[current_curve.end]) or stack is full:
            // draw as line segment if flat enough
            draw_line(buffer, width, height,
                     points[current_curve.beg],
                     points[current_curve.end])
        else:
            // subdivide curve at midpoint using de Casteljau's algorithm
            left_ctrl = midpoint(points[current_curve.beg], points[current_curve.ctrl])
            right_ctrl = midpoint(points[current_curve.ctrl], points[current_curve.end])
            center = midpoint(left_ctrl, right_ctrl)

            // add new points to collection
            center_idx = add_point(center)
            left_ctrl_idx = add_point(left_ctrl)
            right_ctrl_idx = add_point(right_ctrl)

            // push subdivided curves onto stack (right then left for proper order)
            stack.push(Curve(center_idx, current_curve.end, right_ctrl_idx))
            stack.push(Curve(current_curve.beg, center_idx, left_ctrl_idx))

function draw_line(buffer, width, height, origin, goal):
    delta = (goal.x - origin.x, goal.y - origin.y)
    direction = (sign(delta.x), sign(delta.y))

    // skip horizontal lines (no area contribution)
    if direction.y == 0:
        return

    // init grid traversal
    init_pixel_coordinates(origin, direction, &pixel)
    init_crossing_distances(origin, delta, direction, pixel, &next_crossing, &crossing_incr)

    prev_distance = 0.0
    num_steps = calculate_steps(origin, goal, direction)
    half_delta_x = 0.5 * delta.x

    // walk through each pixel the line passes through
    for step = 0 to num_steps-1:
        next_distance = min(next_crossing.x, next_crossing.y)

        // calculate trapezoidal area contribution to current pixel
        x_average = origin.x + (prev_distance + next_distance) * half_delta_x
        y_difference = (next_distance - prev_distance) * delta.y

        // update coverage and area for this pixel
        cell = buffer[pixel.y * width + pixel.x]
        cell.cover += y_difference
        cell.area += (1.0 - (x_average - pixel.x)) * y_difference
        buffer[pixel.y * width + pixel.x] = cell

            // advance to next pixel
        if next_crossing.x < next_crossing.y:
            pixel.x += direction.x
            next_crossing.x += crossing_incr.x
        else:
            pixel.y += direction.y
            next_crossing.y += crossing_incr.y

        prev_distance = next_distance

    // handle final segment to line endpoint
    next_distance = 1.0
    x_average = origin.x + (prev_distance + next_distance) * half_delta_x
    y_difference = (next_distance - prev_distance) * delta.y

    cell = buffer[pixel.y * width + pixel.x]
    cell.cover += y_difference
    cell.area += (1.0 - (x_average - pixel.x)) * y_difference
    buffer[pixel.y * width + pixel.x] = cell

function post_process(buffer, output_image, width, height):
    accumulated = 0.0

    for each pixel in buffer (row-major order):
        // combine accumulated coverage with local area
        value = abs(accumulated + pixel.area)
        value = clamp(value, 0.0, 1.0)

        // convert to 8-bit grayscale
        output_image[pixel_index] = round(value * 255.0)

        // update accumulated coverage for next pixel
        accumulated += pixel.cover