[Перевод] Трассировщик лучей с нуля за 100 строчек Python10.03.2021 14:17

Рисунок 1

В этом посте мы заглянем под капот алгоритмов компьютерной графики, пошагово разберем основные принципы трассировки лучей и напишем ее простую реализацию на Python. Никаких сторонних графических библиотек — только NumPy и голый код в компиляторе.

Примечание: Эта статья ни в коем случае не является полным руководством/объяснением трассировки лучей, поскольку эта тема слишком обширна, а скорее просто введением для любопытствующих.

Предпосылки

Нам понадобится только самая базовая векторная геометрия:

Пусть у нас есть две точки A и B — независимо от размерности: 1, 2, 3, …, n, — тогда вектор, идущий от A к B, может быть найден путем вычисления B — A (поэлементно);
Длину вектора — независимо от его размерности — можно найти, вычислив квадратный корень из суммы возведенных в квадрат компонентов. Длина вектора v обозначается ||v||;
Единичный вектор — это вектор длины 1: ||v|| = 1;
Для данного вектора другой вектор, идущий в том же направлении, но имеющий длину 1, может быть найден путем деления каждого компонента первого вектора на его длину — это называется нормализацией: u = v/||v||;
Точечное произведение для векторов вычисляется как: = ||v||².

Алгоритм трассировки лучей

Трассировка лучей — это метод рендеринга, который имитирует путь света и пересечения с объектами и позволяет создавать изображения с высокой степенью реализма. Более оптимизированные варианты этого алгоритма используются в видеоиграх.

Чтобы объяснить работу этого алгоритма, сначала нужно настроить сцену:

Трехмерное пространство (мы собираемся использовать три координаты для позиционирования объектов в пространстве);
Объекты в этом пространстве (поскольку мы собираемся воспроизвести рис. 1, то возьмем в качестве объектов сферы);
Источник света (в нашем случае это одна точка, излучающая свет во всех направлениях);
Камера для наблюдения за сценой;
Экран, через который камера смотрит на объекты (четыре точки для четырех углов прямоугольного экрана).

Рисунок 2

Экран — это некая определенная вами геометрическая фигура (например, прямоугольник 3×2). Но сами по себе цифры 3 и 2 ни о чем нам не говорят и действительно начинают приобретать какое-то значение только при сравнении их с размерами других объектов. Здесь важно то, как вы разделите ваш прямоугольник на более мелкие квадраты (пиксели), как на рисунке выше. Это определит размер конечного изображения. Другими словами, можно взять прямоугольник 3×2 и разделить его на 300×200 пикселей.

Напишем алгоритм трассировки лучей с учетом заданной сцены:

для каждого пикселя p(x, y, z) экрана:
   ассоциировать черный цвет с p
   если луч (линия), начинающийся от камеры и идущий к точке p, пересекает любой объект сцены, то:
       вычислить точку пересечения с ближайшим объектом
       если между точкой пересечения и источником света нет объекта, тогда:
           рассчитать цвет точки пересечения
           сопоставить цвет точки пересечения с p

Рисунок 3

Обратите внимание, что этот процесс на деле оказывается обратным реальному освещению. Ведь реальный свет выходит из источника во всех направлениях, отражается от объектов и попадает в камеру. Однако, поскольку не все лучи, выходящие из источника света, попадут в камеру, трассировка лучей выполняет обратный процесс для экономии времени вычислений (отслеживает лучи от камеры обратно к источнику света).

Настройка сцены

Перед тем, как начать писать код, нам нужно настроить сцену. В первую очередь определимся, где расположены камера и экран. Для этой цели примем некоторые упрощения, совместив объекты с координатными осями.

Рисунок 4

Допустим, камера расположена в точке (x = 0, y = 0, z = 1), а экран является частью плоскости, образованной осями x и y. Теперь мы можем написать скелет нашего кода.

Посмотреть код

import numpy as np
import matplotlib.pyplot as plt
width = 300
height = 200
camera = np.array([0, 0, 1])
ratio = float(width) / height
screen = (-1, 1 / ratio, 1, -1 / ratio) # слева, сверху, справа, снизу
image = np.zeros((height, width, 3))
for i, y in enumerate(np.linspace(screen[1], screen[3], height)):
    for j, x in enumerate(np.linspace(screen[0], screen[2], width)):
        # image[i, j] = ...
        print("progress: %d/%d" % (i + 1, height))
plt.imsave('image.png', image)

Камера — это просто точка, имеющая три координаты;
С другой стороны, экран определяется четырьмя числами (или двумя точками): слева, сверху, справа, снизу. Он находится в диапазоне от -1 до 1 в направлении x и от -1/ratio до 1/ratio в направлении y, где ratio — ширина изображения, деленная на его высоту. Это вытекает из того, что мы хотим, чтобы экран имел такое же соотношение сторон, что и фактическое изображение. При такой настройке экрана будет получено соотношение сторон (ширина к высоте): 2 /(2/ratio) = ratio, которое и является соотношением сторон желаемого изображения 300×200;
Цикл состоит из разделения экрана на точки в направлениях x и y, затем вычисляется цвет текущего пикселя;
Полученный код создаст — как и ожидалось на данном этапе — черное изображение.

Пересечение лучей

Следующий шаг алгоритма: если луч (линия), начинающийся от камеры и проходящий к точке p, пересекает объект сцены, тогда…

Разобьем его на две части. И начнем с определения того, какой луч (линия) начинается от камеры и идет к точке p?

Определение луча

Когда мы говорим »луч», на самом деле мы имеем в виду скорее »линию». Всякий раз при работе с геометрией лучше отдать предпочтение векторам, чем реальным линейным уравнениям: с ними легче работать, и они гораздо менее подвержены ошибкам, таким как деление на ноль.

Итак, поскольку луч начинается от камеры и идет в направлении текущего целевого пикселя, мы можем определить единичный вектор, указывающий в том же направлении. Поэтому мы определяем наш луч следующим уравнением:

Помните, что камера и пиксель — это 3D-точки. При t = 0 вы окажетесь в положении камеры, но с увеличением t будете удаляться от нее в направлении пикселя. Это параметрическое уравнение, которое даст точку вдоль линии для любого t.

Конечно, точно так же мы можем переписать уравнение и определить луч, который начинается в исходной точке (O) и идет к месту назначения (D) как:

Для удобства определим d как вектор направления.

Теперь мы можем добавить к нашему коду вычисление луча:

Посмотреть код

import numpy as np
import matplotlib.pyplot as plt
def normalize(vector):
    return vector / np.linalg.norm(vector)
width = 300
height = 200
camera = np.array([0, 0, 1])
ratio = float(width) / height
screen = (-1, 1 / ratio, 1, -1 / ratio) # слева, сверху, справа, снизу
image = np.zeros((height, width, 3))
for i, y in enumerate(np.linspace(screen[1], screen[3], height)):
    for j, x in enumerate(np.linspace(screen[0], screen[2], width)):
        pixel = np.array([x, y, 0])
        origin = camera
        direction = normalize(pixel - origin)
        # image[i, j] = ...
    print("progress: %d/%d" % (i + 1, height))
plt.imsave('image.png', image)

Мы добавили в код функцию normalize (vector), которая возвращает… собственно, нормализованный вектор;
Также мы добавили вычисление исходной точки и направления, которые вместе определяют луч. Обратите внимание, что пиксель имеет координату z = 0, поскольку он лежит на экране, который находится в плоскости, образованной осями x и y;

Теперь перейдем ко второй части, где луч пересекает объект сцены. Для простоты будем использовать только сферы.

Определение сферы

Сфера — довольно простой математический объект. Она определяется как набор точек, находящихся на одинаковом расстоянии r (радиус) от заданной точки (центра).

Следовательно, с учетом центра C сферы и ее радиуса r произвольная точка X лежит на сфере тогда, когда:

Для удобства возведем обе стороны в квадрат, чтобы избавиться от квадратного корня, обусловленного величиной X — C:

После этого мы сможем определить некоторые сферы сразу после объявления экрана:

objects = [
   { 'center': np.array([-0.2, 0, -1]), 'radius': 0.7 },
   { 'center': np.array([0.1, -0.3, 0]), 'radius': 0.1 },
   { 'center': np.array([-0.3, 0, 0]), 'radius': 0.15 }
]

Теперь вычислим пересечение луча и сферы.

Пересечение со сферой

Мы знаем уравнение лучей и знаем, какому условию должна удовлетворять точка, чтобы она лежала на сфере. Все, что нам нужно сделать, это подставить одно уравнение в другое и решить его относительно t. То есть, найти ответ на вопрос: для какого t точка луча ray (t) окажется на сфере?

Это обычное квадратное уравнение, которое просто решается относительно t. Мы будем вызывать коэффициенты, связанные с t², t¹, t⁰, a, b и c, соответственно. Вычислим дискриминант этого уравнения:

Поскольку направление d является единичным вектором, получим a = 1. После вычисления дискриминанта у нас есть три варианта:

Рисунок 5

Для обнаружения пересечений будем использовать только третий случай. Запишем функцию, отвечающую за обнаружение пересечения луча и сферы. Она возвращает расстояние от начала луча до ближайшей точки пересечения, если луч действительно пересекает сферу, иначе возвращает None:

Посмотреть код

def sphere_intersect(center, radius, ray_origin, ray_direction):
   b = 2 * np.dot(ray_direction, ray_origin — center)
   c = np.linalg.norm(ray_origin — center) ** 2 — radius ** 2
   delta = b ** 2 — 4 * c
   if delta > 0:
       t1 = (-b + np.sqrt(delta)) / 2
       t2 = (-b — np.sqrt(delta)) / 2
       if t1 > 0 and t2 > 0:
           return min(t1, t2)
   return None

Обратите внимание, что мы возвращаем только ближайшее пересечение из двух тогда, когда оба t1 и t2 положительны. Это связано с тем, что ответ уравнения может быть отрицательным, и в таком случае луч, пересекающий сферу, будет иметь не d в качестве вектора направления, а -d (например, если сфера находится за камерой и экраном).

Ближайший пересекаемый объект

Пока мы все еще не выполнили инструкцию из псевдокода: если луч (линия), начинающийся от камеры и идущий к точке p, пересекает любой объект сцены, то […]. Теперь нам нужно вычислить точку пересечения с ближайшим объектом.

Напишем функцию, которая использует sphere_intersect () для поиска ближайшего объекта, с которым пересекается луч, если он существует. Просто перебираем все сферы, ищем пересечения и сохраняем ближайшую сферу:

Посмотреть код

def nearest_intersected_object(objects, ray_origin, ray_direction):
   distances = [sphere_intersect(obj['center'], obj['radius'], ray_origin, ray_direction) for obj in objects]
   nearest_object = None
   min_distance = np.inf
   for index, distance in enumerate(distances):
       if distance and distance < min_distance:
           min_distance = distance
           nearest_object = objects[index]
   return nearest_object, min_distance

При вызове функции, если nearest_object = None, луч не пересекает никакого объекта, иначе его значением является ближайший пересекаемый объект, и мы получаем min_distance — расстояние от начала луча до точки пересечения.

Точка пересечения

Чтобы вычислить точку пересечения, используем предыдущую функцию:

nearest_object, distance = nearest_intersected_object(objects, o, d)
if nearest_object:
   intersection_point = o + d * distance

В результате получаем следующий код:

Посмотреть код

import numpy as np
import matplotlib.pyplot as plt
def normalize(vector):
   return vector / np.linalg.norm(vector)
def sphere_intersect(center, radius, ray_origin, ray_direction):
   b = 2 * np.dot(ray_direction, ray_origin — center)
   c = np.linalg.norm(ray_origin — center) ** 2 — radius ** 2
   delta = b ** 2 — 4 * c
   if delta > 0:
       t1 = (-b + np.sqrt(delta)) / 2
       t2 = (-b — np.sqrt(delta)) / 2
       if t1 > 0 and t2 > 0:
           return min(t1, t2)
   return None
def nearest_intersected_object(objects, ray_origin, ray_direction):
   distances = [sphere_intersect(obj['center'], obj['radius'], ray_origin, ray_direction) for obj in objects]
   nearest_object = None
   min_distance = np.inf
   for index, distance in enumerate(distances):
       if distance and distance < min_distance:
           min_distance = distance
           nearest_object = objects[index]
   return nearest_object, min_distance
width = 300
height = 200
camera = np.array([0, 0, 1])
ratio = float(width) / height
screen = (-1, 1 / ratio, 1, -1 / ratio) # слева, сверху, справа, снизу
objects = [
   { 'center': np.array([-0.2, 0, -1]), 'radius': 0.7 },
   { 'center': np.array([0.1, -0.3, 0]), 'radius': 0.1 },
   { 'center': np.array([-0.3, 0, 0]), 'radius': 0.15 }
]
image = np.zeros((height, width, 3))
for i, y in enumerate(np.linspace(screen[1], screen[3], height)):
   for j, x in enumerate(np.linspace(screen[0], screen[2], width)):
       pixel = np.array([x, y, 0])
       origin = camera
       direction = normalize(pixel — origin)
       # проверка пересечений
       nearest_object, min_distance = nearest_intersected_object(objects, origin, direction)
       if nearest_object is None:
           continue
       # вычисления пересечений между лучом и ближайшим объектом
       intersection = origin + min_distance * direction
       # image[i, j] = ...
       print("%d/%d" % (i + 1, height))
plt.imsave('image.png', image)

Пересечения света

Итак, мы нашли прямую линию, идущую от камеры к объекту, и знаем, что это за объект, а также на какую часть объекта мы смотрим. Но мы пока не знаем, освещена ли вообще эта конкретная точка. Возможно, свет не попадает конкретно на нее, и нет необходимости производить расчеты дальше, потому что мы ее не видим. Следовательно, следующий шаг — проверить, нет ли никаких посторонних объектов между точкой пересечения и источником света.

У нас уже есть функция, которая нам может помочь: near_intersected_object (). И мы хотим знать, пересекает ли луч, который начинается в точке пересечения и идет к свету, объект сцены перед тем, как пересечь свет. Это практически та же задача, что мы решали раньше: нам просто нужно изменить начало и направление луча. Во-первых, нам нужно определить свет. Можно сделать это сразу вместе с объявлением объектов:

light = { 'position': np.array([5, 5, 5]) }

Чтобы проверить, затеняет ли объект точку пересечения, мы должны пропустить луч, который начинается в точке пересечения и идет к свету, и посмотреть, действительно ли ближайший возвращенный объект оказывается ближе, чем свет, к точке пересечения (другими словами, находится ли он между ними).

# ...
intersection = origin + min_distance * direction
intersection_to_light = normalize(light['position'] — intersection)
_, min_distance = nearest_intersected_object(objects, intersection, intersection_to_light)
intersection_to_light_distance = np.linalg.norm(light['position'] — intersection)
is_shadowed = min_distance < intersection_to_light_distance

Что ж, это так не сработает. Нужно сделать небольшую корректировку. Если мы используем точку пересечения в качестве источника нового луча, мы можем в конечном итоге обнаружить сферу, на которой мы сейчас находимся, как объект между точкой пересечения и источником света. Быстрое и широко используемое решение этой проблемы — сделать небольшой шаг, который уводит нас от поверхности сферы. Обычно для этого вычисляется вектор нормали к поверхности и производится отступ в направлении этой нормали.

Рисунок 6

Этот трюк используется не только для сфер, но и для любых объектов.

Следовательно, правильный код будет таким:

Посмотреть код

# ...
intersection = origin + min_distance * direction
normal_to_surface = normalize(intersection — nearest_object['center'])
shifted_point = intersection + 1e-5 * normal_to_surface
intersection_to_light = normalize(light['position'] — shifted_point)
_, min_distance = nearest_intersected_object(objects, shifted_point, intersection_to_light)
intersection_to_light_distance = np.linalg.norm(light['position'] — intersection)
is_shadowed = min_distance < intersection_to_light_distance
if is_shadowed:
   continue

Модель отражения Блинна-Фонга

Итак, мы знаем, что луч света попал на объект, а отражение луча — прямо в камеру. Вопрос: что увидит камера? На него и пытается ответить модель Блинна-Фонга.

Модель Блинна-Фонга — это приближение к модели Фонга, требующее меньших вычислительных затрат.

Согласно этой модели, любой материал имеет четыре свойства:

Фоновый цвет (Ambient color): цвет, который имеет объект в отсутствие света;
Рассеянный цвет (Diffuse color): цвет, наиболее близкий к тому, что мы себе представляем;
Зеркальный цвет (Specular color): цвет блестящей части объекта, когда свет попадает на нее. В большинстве случаев это белый цвет;
Блеск (Shininess): коэффициент, показывающий, насколько блестящим является объект.

Примечание: Все цвета представлены в RGB в диапазоне 0 — 1.

Рисунок 7

Добавим эти свойства к сферам:

objects = [
   { 'center': np.array([-0.2, 0, -1]), 'radius': 0.7, 'ambient': np.array([0.1, 0, 0]), 'diffuse': np.array([0.7, 0, 0]), 'specular': np.array([1, 1, 1]), 'shininess': 100 },
   { 'center': np.array([0.1, -0.3, 0]), 'radius': 0.1, 'ambient': np.array([0.1, 0, 0.1]), 'diffuse': np.array([0.7, 0, 0.7]), 'specular': np.array([1, 1, 1]), 'shininess': 100 },
   { 'center': np.array([-0.3, 0, 0]), 'radius': 0.15, 'ambient': np.array([0, 0.1, 0]), 'diffuse': np.array([0, 0.6, 0]), 'specular': np.array([1, 1, 1]), 'shininess': 100 }
]

В нашем примере сферы имеют цвета красный, пурпурный и зеленый, соответственно.

Модель Блинн-Фонга утверждает, что свет также имеет три цветовых свойства: фоновый цвет, рассеянный и зеркальный. Их тоже добавим к модели:

light = { 'position': np.array([5, 5, 5]), 'ambient': np.array([1, 1, 1]), 'diffuse': np.array([1, 1, 1]), 'specular': np.array([1, 1, 1]) }

Учитывая эти свойства, модель Блинна-Фонга рассчитывает освещенность точки следующим образом:

где

ka, kd, ks — фоновое, рассеянное, зеркальное свойства объекта;
ia, id, is — фоновое, рассеянное, зеркальное свойства света;
L — единичный вектор направления от точки пересечения к свету;
N — единичный вектор нормали к поверхности объекта в точке пересечения;
V — единичный вектор направления от точки пересечения к камере;
α — блеск объекта.

Посмотреть код

# ...
if is_shadowed:
   break

# RGB
illumination = np.zeros((3))

# ambiant
illumination += nearest_object['ambient'] * light['ambient']

# diffuse
illumination += nearest_object['diffuse'] * light['diffuse'] * np.dot(intersection_to_light, normal_to_surface)

# specular
intersection_to_camera = normalize(camera — intersection)
H = normalize(intersection_to_light + intersection_to_camera)
illumination += nearest_object['specular'] * light['specular'] * np.dot(normal_to_surface, H) ** (nearest_object['shininess'] / 4)
image[i, j] = np.clip(illumination, 0, 1)

Обратите внимание, что в конце мы определяем цвет между 0 и 1, чтобы убедиться, что он находится в правильном диапазоне.

Запускаем код

Увеличим ширину и высоту для получения более высокого разрешения (ценой увеличения времени вычислений).

Рисунок 8

Можно заметить две вещи, которые отличают результат от первого изображения, показанного в начале:

Серый пол отсутствует;
Отсутствие отражений.

Фейковая плоскость

В идеале мы должны создать другой тип объекта — плоскость, но поскольку мы достаточно ленивы, то можем просто добавить другую сферу. Ведь если вы стоите на сфере, имеющей бесконечно большой радиус (по сравнению с вашим размером), тогда вам будет казаться, что вы стоите на плоской поверхности.

Добавим эту сферу в список объектов и снова выполним рендеринг:

{ 'center': np.array([0, -9000, 0]), 'radius': 9000 — 0.7, 'ambient': np.array([0.1, 0.1, 0.1]), 'diffuse': np.array([0.6, 0.6, 0.6]), 'specular': np.array([1, 1, 1]), 'shininess': 100 }

Отражение

Сейчас мы рендерим лучи, которые выходят из источника света, ударяются о поверхность объекта и отражаются в камеру. Но что, если луч попадет в несколько объектов, прежде чем попасть в камеру? Появится отражение.

Каждый объект имеет коэффициент отражения в диапазоне от 0 до 1. Здесь 0 означает, что объект матовый, 1 — что объект зеркальный. Добавим свойство отражения ко всем сферам:

{ 'center': np.array([-0.2, 0, -1]), ..., 'reflection': 0.5 }
{ 'center': np.array([0.1, -0.3, 0]), ..., 'reflection': 0.5 }
{ 'center': np.array([-0.3, 0, 0]), ..., 'reflection': 0.5 }
{ 'center': np.array([0, -9000, 0]), ..., 'reflection': 0.5 }

Алгоритм

Чтобы включить вычисление отражений, нужно отследить отраженный луч после того, как произошло пересечение, и учесть цветовой вклад каждой точки этого пересечения. Повторяем этот процесс несколько раз.

Рисунок 9

Расчет цвета

Чтобы получить цвет пикселя, нужно просуммировать вклад каждой точки пересечения луча:

где

c — конечный цвет пикселя;
i — освещенность, рассчитанная по модели Блинна-Фонга для точки пересечения;
r — отражение от пересекаемого объекта.

Когда прекратить вычисление этой суммы (и отслеживание отраженных лучей, соответственно), решаете вы сами.

Отраженный луч

Прежде чем мы сможем все это записать в виде кода, нам нужно определить направление отраженного луча. Можно вычислить его следующим образом:

где

R — нормализованный отраженный луч;
L — единичный вектор направления отражаемого луча;
N — единичный вектор направления нормали к поверхности хода луча.

Рисунок 10

Добавим этот метод в начало кода вместе с функцией normalize ():

def reflected(vector, axis):
   return vector — 2 * np.dot(vector, axis) * axis

Код

Посмотреть код

# глобальные переменные
max_depth = 3

# нужные данные для цикла

color = np.zeros((3))
reflection = 1

for k in range(max_depth):
   nearest_object, min_distance = # ...

   # ...
   illumination += # ...

   # отражение
   color += reflection * illumination
   reflection *= nearest_object['reflection']

   # начальная точка и направление нового луча
   origin = shifted_point
   direction = reflected(direction, normal_to_surface)
image[i, j] = np.clip(color, 0, 1)

Важно: теперь, когда мы поместили пересечения в другой цикл для отражений, нужно добавить операторы break во избежание бесполезных вычислений.

Вот и все. Запускаем код и наблюдаем результат:

Рисунок 11

Окончательный код

Итоговый код состоит из всего порядка сотни строк:

Посмотреть код

import numpy as np
import matplotlib.pyplot as plt

def normalize(vector):
    return vector / np.linalg.norm(vector)

def reflected(vector, axis):
    return vector - 2 * np.dot(vector, axis) * axis


def sphere_intersect(center, radius, ray_origin, ray_direction):
    b = 2 * np.dot(ray_direction, ray_origin - center)
    c = np.linalg.norm(ray_origin - center) ** 2 - radius ** 2
    delta = b ** 2 - 4 * c
    if delta > 0:
        t1 = (-b + np.sqrt(delta)) / 2
        t2 = (-b - np.sqrt(delta)) / 2
        if t1 > 0 and t2 > 0:
            return min(t1, t2)
    return None

def nearest_intersected_object(objects, ray_origin, ray_direction):
    distances = [sphere_intersect(obj['center'], obj['radius'], ray_origin, ray_direction) for obj in objects]
    nearest_object = None
    min_distance = np.inf
    for index, distance in enumerate(distances):
        if distance and distance < min_distance:
            min_distance = distance
            nearest_object = objects[index]
    return nearest_object, min_distance

width = 300
height = 200
max_depth = 3

camera = np.array([0, 0, 1])
ratio = float(width) / height
screen = (-1, 1 / ratio, 1, -1 / ratio) # left, top, right, bottom

light = { 'position': np.array([5, 5, 5]), 'ambient': np.array([1, 1, 1]), 'diffuse': np.array([1, 1, 1]), 'specular': np.array([1, 1, 1]) }

objects = [
    { 'center': np.array([-0.2, 0, -1]), 'radius': 0.7, 'ambient': np.array([0.1, 0, 0]), 'diffuse': np.array([0.7, 0, 0]), 'specular': np.array([1, 1, 1]), 'shininess': 100, 'reflection': 0.5 },
    { 'center': np.array([0.1, -0.3, 0]), 'radius': 0.1, 'ambient': np.array([0.1, 0, 0.1]), 'diffuse': np.array([0.7, 0, 0.7]), 'specular': np.array([1, 1, 1]), 'shininess': 100, 'reflection': 0.5 },
    { 'center': np.array([-0.3, 0, 0]), 'radius': 0.15, 'ambient': np.array([0, 0.1, 0]), 'diffuse': np.array([0, 0.6, 0]), 'specular': np.array([1, 1, 1]), 'shininess': 100, 'reflection': 0.5 },
    { 'center': np.array([0, -9000, 0]), 'radius': 9000 - 0.7, 'ambient': np.array([0.1, 0.1, 0.1]), 'diffuse': np.array([0.6, 0.6, 0.6]), 'specular': np.array([1, 1, 1]), 'shininess': 100, 'reflection': 0.5 }
]

image = np.zeros((height, width, 3))
for i, y in enumerate(np.linspace(screen[1], screen[3], height)):
    for j, x in enumerate(np.linspace(screen[0], screen[2], width)):

        # экран в начальной точке

        pixel = np.array([x, y, 0])
        origin = camera
        direction = normalize(pixel - origin)

        color = np.zeros((3))
        reflection = 1

        for k in range(max_depth):

            # проверка пересечений

            nearest_object, min_distance = nearest_intersected_object(objects, origin, direction)
            if nearest_object is None:

                break

            intersection = origin + min_distance * direction
            normal_to_surface = normalize(intersection - nearest_object['center'])
            shifted_point = intersection + 1e-5 * normal_to_surface
            intersection_to_light = normalize(light['position'] - shifted_point)

            _, min_distance = nearest_intersected_object(objects, shifted_point, intersection_to_light)

            intersection_to_light_distance = np.linalg.norm(light['position'] - intersection)
            is_shadowed = min_distance < intersection_to_light_distance

            if is_shadowed:
                break

            illumination = np.zeros((3))

            # ambiant

            illumination += nearest_object['ambient'] * light['ambient']

            # diffuse

            illumination += nearest_object['diffuse'] * light['diffuse'] * np.dot(intersection_to_light, normal_to_surface)

            # specular

            intersection_to_camera = normalize(camera - intersection)
            H = normalize(intersection_to_light + intersection_to_camera)
            illumination += nearest_object['specular'] * light['specular'] * np.dot(normal_to_surface, H) ** (nearest_object['shininess'] / 4)

            # reflection

            color += reflection * illumination
            reflection *= nearest_object['reflection']
            origin = shifted_point
            direction = reflected(direction, normal_to_surface)
        image[i, j] = np.clip(color, 0, 1)
    print("%d/%d" % (i + 1, height))
plt.imsave('image.png', image)

Что дальше?

Это очень упрощенная программа, предназначенная для ознакомления с основами трассировки лучей. Есть много способов улучшить ее и реализовать другие функции. Например:

Можно создать классы и выяснить, что является специфическим для сфер, а что нет, определить базовый класс и реализовать другие объекты, такие как плоскости или треугольники;
То же самое и со светом. Добавить сюда POO и сделать так, чтобы можно было добавить несколько источников света в сцену;
Отделить свойства материала от геометрических свойств, чтобы иметь возможность применять один материал (например, синий матовый) к любому типу объектов;
Найти способ правильно расположить экран при любом положении и направлении камеры;
Смоделировать свет по-другому. В настоящее время это одна точка, поэтому тени от объектов «жесткие» или четко очерченные. Чтобы получить «мягкие» тени , нужно смоделировать источник света как 2D- или 3D-объект: диск или сферу.

Бонус

Ниже приведена анимация трассировки лучей. По сути это просто несколько раз отрендеренная сцена с камерой в разных положениях:

Код написан на Kotlin (можно оценить, насколько медленный по сравнению с ним Python) и доступен на GitHub.