Прошло много лет с тех пор, как я написал donut.c, и всё это время я не раз задумывался, можно ли как-то упростить этот проект. Например, может быть, нашёлся бы способ очертить пончик лучами, дописав для этого немного кода. В октябре 2023 года я написал твит о совершенно внезапном просветлении, позволившем мне найти новый подход к этой проблеме — без привлечения памяти, без каких-либо синусов, косинусов, без квадратных корней, деления, строго говоря, даже без умножения. Всё нужное можно отобразить с помощью одних только сдвигов и сложений. Вот обновлённая версия на C.

Затем мне потребовался почти целый год, чтобы воплотить эту идею на практике, воспользовавшись конкретной аппаратной реализацией. В начале сентября 2024 года я разработал 4-плиточную схему для Tiny Tapeout 8. Она занимает 0,8% 130-нм чипа, используемого и во многих других схемах. По состоянию на январь 2025 года этот чип производит компания SkyWater Technologies, и, если моя схема сработает, то на VGA-монитор будет выведена следующая картинка:

В данной схеме не предусмотрено отображение многоугольников! Вместо этого выполняется поступательное приближение пончика, контур которого вырисовывается при помощи лучевой трассировки. Причём возникает гонка с VGA-лучом, с периодичностью в 39 наносекунд монитор будет отображать следующий пиксель, независимо от того, целиком ли он отрисован. С практической точки зрения мы ограничены таковой частотой ~50 МГц. При этом у нас нет вообще никакого буфера памяти, а значит, нет и времени, чтобы выполнить работу качественно. Если фигура будет напоминать многоугольник, то это получится чисто случайно!

Предвосхищая ваш вопрос, отвечу на него — нет, сама кремниевая подложка не похожа по форме на пончик. Допустим, это задача на следующий раз. На проект требуется примерно 7000 стандартных ячеек (т.e., триггеров или логических вентилей). Кстати, здесь можно посмотреть объёмную модель этого кристалла благодаря замечательным инструментам, созданным командой Tiny Tapeout.

Визуализация компоновки того чипа, который даёт картинку пончика на Tiny Tapeout 8; это небольшой фрагмент гораздо более крупного кристалла

Располагая более обширным участком кристалла и/или повышенной тактовой частотой, эту работу можно было бы конвейеризовать и добиться идеального отображения каждого пикселя. Но для этого потребовалась бы более современная платформа, а не 130-нанометровый SkyWater. Впрочем, этот проект я пытался сделать настолько минималистическим, насколько это возможно.

Исходный код

Проект для TinyTapeout на языке Verilog выложен на github здесь вместе с тестовым стендом Verilator, с применением которого было подготовлено вышеприведённое видео. Там же выложена версия OrangeCrab FPGA, на которой я вёл разработку.

Поскольку я вёл разработку на OrangeCrab, в итоге я приспособился к частоте 48 МГц, нативной для этой платформы. Весь тайминг задан с такой настройкой. Обычно в VGA используется тактовая частота 25,175 МГц, так что я решил применить слегка странный VGA-режим 1220x480, чтобы в полной мере раскрыть возможности тайминга. Картинка отлично выглядит на аналоговом кинескопном проекторе (CRT), но на ЖК-мониторах могут возникать странные дефекты, связанные со смазыванием изображения. Будем знать.

Кроме этой я спроектировал на Tiny Tapeout 8 ещё две схемы, каждая из которых заслуживает отдельной статьи! Это мои первые проекты на Verilog. Код, честно говоря, ужасен, но зато мне удалось уложить проект в те жёсткие рамки, которые я сам себе задал.

Рендеринг пончиков только при помощи операций сдвига и сложения

Здесь нам нужно ухитриться представить тор как функцию расстояния со знаком, а потом трассировать её методом «обхода лучей». Это поразительно мощная техника, которую популяризовал Иньиго Килес, невероятно хороша на практике. Его сайт — настоящий кладезь, но здесь я покажу только небольшую выдержку со страницы, рассказывающей о функциях расчёта расстояния:

Расстояние со знаком до тора можно вычислить при помощи двух последовательных операций по расчёт длины на плоскости. Угадайте, как вычислить длину 2D-вектора? При помощи CORDIC, о другой вариации такого же фокуса я уже рассказывал здесь!

Вот что из этого следует: методом обхода лучей можно найти точку пересечения луча с пончиком. Начинаем с точки на экране, заданной в 3D-координатах, далее вычисляем вышеприведённую функцию, чтобы найти расстояние. Далее перемещаем отрезок, равный этому расстоянию, между началом координат и отображаемым пикселем. Повторяем этот процесс, пока не добьёмся сходимости расстояния около 0.

Так мы узнаем точку, но какой именно цвет будем отображать? Чтобы вычислить освещение, нам понадобится нормаль к поверхности, а для получения этой нормали потребуется вращать вектор направления света вдоль двух осей тора под таким углом, под которым пересекутся луч и тор. На первый взгляд кажется, что это сложно.

Вот в чём заключалось моё озарение: та величина, на которую требуется повернуть вектор освещения, чтобы найти нормаль поверхности, в точности равна ранее отменённым поворотам, при помощи которых мы нашли расстояние до тора! Более того, такой поворот практически ничего не будет нам стоить, поскольку его можно реализовать как побочный эффект при вычислении этих длин.

Давайте попробую объяснить.

Векторный режим CORDIC

Существует фокус, позволяющий вычислить при помощи только лишь сдвигов и сложений. Для этого нужно очень хитро воспользоваться CORDIC: поступательно поворачивать вектор по направлению к оси , один «бит» за раз (сначала , на , <, т.д). Так нам удастся восстановить не только длину, но и угол, фактически, выполнить преобразование из декартовых координат в полярные.

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

function cordicLength(x, y, n) {
    // x должна находиться в правой полуплоскости; поскольку нас интересует только расстояние,
    // можно проигнорировать это действие, но можно и установить здесь флаг, 
    // по которому впоследствии сможем восстановить правильный угол.
    if (x < 0) {
        x = -x;
    }
    for (var i=0; i < n; i++) {
        var t = x;
        if (y < 0) {
            x -= y>>i;
            y += t>>i;
        } else {
            x += y>>i;
            y -= t>>i;
        }
    }
    return x;
}

В данном случае при помощи сдвигов и сложений выполняется вращение и масштабирование. Фактически перемножаются следующие матрицы:

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

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

Векторизация дополнительного вектора

Немного откорректировав эту картину, можно повернуть (и масштабировать с тем же коэффициентом .607) второй вектор на угол исходного вектора, когда будем искать длину исходного вектора.

Теперь, если считать синий вектор «дополнительным», то именно длину чёрного вектора мы хотим найти. При этом зелёная линия находит относительный угол между двумя векторами, не зависящий от длины .

function rotateCordic(x, y, x2, y2, n) {
    xsign = 1;
    // если мы находимся в левой полуплоскости, то, задав отрицательную по модулю величину для x, y и получив на выходе
    // x2, получим значение y2, которое и будет верным результатом
    if (x < 0) {
        x = -x;
        y = -y;
        xsign = -1;
    }
    for (var i=0; i < n; i++) {
        var t = x;
        var t2 = x2;
        if (y < 0) {
            x -= y>>i;
            y += t>>i;
            x2 -= y2>>i;
            y2 += t2>>i;
        } else {
            x += y>>i;
            y -= t>>i;
            x2 += y2>>i;
            y2 -= t2>>i;
        }
    }
    // Теперь x — это длина x,y а [x2*xsign, y2*xsign] — это вектор после поворота
    return [x, x2*xsign, y2*xsign];
}

Откорректированная функция CORDIC для поворота вспомогательного вектора

Именно это нам и требовалось, чтобы вычислить освещение. Если сообщить программе единичный вектор направления движения и нормаль к плоскости, то координата зелёной линии будет равна косинусу угла между ними. Мы ухитрились получить скалярное произведение векторов, вообще не прибегая к умножению, а в качестве побочного продукта при вычислении расстояния.</p.

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

Вот небольшая анимация, иллюстрирующая, как мы находим ближайшую точку в торе. Сначала находим на внутреннем кольце тора точку, ближайшую к нашей точке. Всё равно, что спроецировать нашу точку на плоскость тора x-z и масштабировать расстояние до . Эта новая точка станет центром кольца с радиусом в плоскости, определяемой этой точкой, а также будет точкой начала координат и центром тора. Затем ещё раз спроецируем точку начала координат в этой плоскости, на этот раз на расстояние .

float distance_to_torus(float x, float y, float z, float r1, float r2) {
  // Сначала находим расстояние от точки p до круга с радиусом r1 в плоскости x-z 
  float dist_to_r1_circle = length(x, z) - r1;
  // Затем находим расстояние на сечении тора от центрального кольца до точки p
  return length(dist_to_r1_circle, y) - r2;
}

Эта последняя проекционная линия проводится для окончательного расчёта расстояния между тором и точкой, и эта линия в точности соответствует нормали к плоскости.

Такой приём возможен исключительно благодаря свойству торов, в соответствии с которым в двух измерениях у них две длины. Например, мне не удалось таким же образом отобразить сферу, хотя, казалось бы, это проще!

donut.c с обходом лучей

Вооружившись этой техникой, я написал ещё один donut.c (вот исходный код), опять же, только при помощи операций сдвига и сложения. Проект совсем не такой компактный, поскольку в нём гораздо больше кода на C, зато гораздо меньше логических вентилей. Я написал об этом в Твиттере, и вскоре мне ответил Бруно Леви, поделившийся интересной находкой: если ограничить количество итераций CORDIC, то пончик получается не гладким, а ребристым!

В окончательной версии ASIC я использовал 3 итерации и 2 для , поскольку именно столько комбинаторной логики мне удалось вложить в единственный цикл при частоте 48 МГц и при этом уложиться в требования по времени. Вот почему рендеринг выглядит именно так.

Составное вращение без умножения

Чтобы реализовать вращение целостного объекта, нужно отслеживать восемь переменных: sin/cos угла A, sin/cos угла B и четыре их произведения. В коде C это делается сверху вниз с применением макроса вращения R():

// HAKMEM 149 – Круговой алгоритм Минского
// Вращение происходит вокруг точки «вблизи» от начала координат, без потери магнитуды
// в течение длительных периодов времени, при условии, что хватает разрядов точности в 
// x и y. Здесь я использую 14. Экономичный способ вычислять приближения синусов и косинусов 
#define R(s,x,y) x-=(y>>s); y+=(x>>s)

    R(5, cA, sA);
    R(5, cAsB, sAsB);
    R(5, cAcB, sAcB);
    R(6, cB, sB);
    R(6, cAcB, cAsB);
    R(6, sAcB, sAsB);

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

К сожалению, со временем эти настройки могут сдвинуться¹, в результате чего рендер явственно скособочится. Решается так: можно заново инициализировать все эти значения в нулевые (то есть, в те, с которых начинается вращение), как только завершится полный цикл, захватывающий A и B.

Соображения о TinyTapeout

TinyTapeout рекомендует использовать значение тактовой частоты до 50 МГц, поэтому я остановился на 48 МГц (это значение действует по умолчанию на той плате FPGA, на которой я веду разработку). Таким образом, я получил примерно 1220 тактов на видимую часть растровой строки по горизонтали. Пожертвовав точностью, я смог уместить все итерации CORDIC в один цикл, но при обходе лучей итерации требуется продолжать до самого схождения, а при любом значении менее 8 итераций результат выглядел совсем плохо. Поэтому я решил остановиться на 8 итерациях, затем считывать вывод аппаратного обхода лучей и загружать новые направления лучей на следующие 8 тактов. Но, чтобы результат не выглядел слишком угловатым, я добавил дизеринг горизонтального смещения как в пространстве, так и во времени. Поэтому на видео пончик выглядит немного глючным/зашумленным.

Видеовывод обеспечивается через Tiny VGA Pmod, где предусмотрено всего два разряда на каждый канал RGB. Чтобы справиться с этим, я решил применить упорядоченный дизеринг, потратив 6 разрядов на внутреннюю глубину цвета — оказывается, это очень легко делается в Verilog. В результате генерируется матрица сглаживания, которая переключается на кадре и таким образом также обеспечивает дизеринг во времени.

// Дизеринг Байера
// это матрица Байера 8x4, которая переключается на каждом следующем кадре
wire [2:0] bayer_i = h_count[2:0] ^ {3{frame[0]}};
wire [1:0] bayer_j = v_count[1:0];
wire [2:0] bayer_x = {bayer_i[2], bayer_i[1]^bayer_j[1], bayer_i[0]^bayer_j[0]};
wire [4:0] bayer = {bayer_x[0], bayer_i[0], bayer_x[1], bayer_i[1], bayer_x[2]};

В результате генерируется матрица, имеющая следующий вид и отзеркаливаемая на каждом следующем кадре:

То же касается и 6-разрядных цветов, вот так:

// На выходе получаем сглаженный 2-разрядный цвет на основе 6-разрядного цвета и 5-разрядной матрицы Байера 
function [1:0] dither2;
    input [5:0] color6;
    input [4:0] bayer5;
    begin
        dither2 = ({1'b0, color6} + {2'b0, bayer5} + color6[0] + color6[5] + color6[5:1]) >> 5;
    end
endfunction

wire [1:0] rdither = dither2(r, bayer);
wire [1:0] gdither = dither2(g, bayer);
wire [1:0] bdither = dither2(b, bayer);

Как обойтись без умножения на последнем шаге

Также было непросто обеспечить отступ единичного 3D-вектора на указанное расстояние, чтобы далее выполнить обход лучей. В идеале для этого потребовались бы 3 параллельные 14-разрядные операции умножения, но в таком случае нам не хватило бы доступной площади кристалла, а я не хотел разматывать эту операцию ещё на множество циклов. Решил, что лучше применю приблизительный множитель, который находит ведущий разряд 1 и возвращает сдвинутую копию единичного вектора, тем самым задавая направление отступа.

Заключение

Этот проект во всём его разнообразии — от обфускации кода C, печатающего ASCII-картинки и до побитовых операций, применяемых для отрисовки 3D-графики и, наконец, до работы с крошечным срезом кремния, по которому бегает электронный луч — показывает, как можно переосмыслить и отточить одну и ту же простейшую идею. Я начал этот проект как забавную задачку по программированию, а в итоге открыл с её помощью элегантный способ отображения 3D-графики, требующий минимальных вычислительных ресурсов. Притом, что первый donut.c был для меня упражнением в творческом самоограничении — уместить механизм для 3D-рендеринга в минимальный объём кода — аппаратная реализация открыла передо мной совершенно иные горизонты. Каков минимум логических вентилей, минимум операций на пиксель, располагая которым всё равно можно запрограммировать узнаваемую анимацию — вращающийся пончик?

Рубленая форма, получившаяся у меня из-за ограничения количества итераций CORDIC — результат чистого везения. Она лишний раз напоминает, что иногда самые интересные находки попадаются именно при работе в условиях жесточайших ограничений.
Ожидается, что Tiny Tapeout 8 начнёт поставку чипов в мае 2025 года; а мне уже не терпится. На этом чипе также записаны условия конкурса Tiny Tapeout Demoscene Competition; возможно, это крутейший ASIC в истории.

---

1. Даже притом, что алгоритм HAKMEM 149 позволяет вращать точку в круге без потери магнитуды, произведения синусов и косинусов при этом вращаются не по кругу, а по фигуре Лиссажу, поэтому магнитуда всё-таки теряется. В результате получается очень испорченная матрица поворота без явных точек сброса. ↩