Трассировщик лучей на визитке

Рассказывает Фабиен Санглард, автор блога fabiensanglard.net


Недавно в интернете я наткнулся на трассировщик лучей на визитке Пола Гекберта. Для тех, кто не в курсе: это очень известная задача, изначально предложенная Полом Гекбертом 4-ого мая 1984 на comp.graphics. Ее суть в том, чтобы написать демонстрацию метода бросания лучей, которая бы… умещалась на визитной карточке (больше об этом читайте в разделе «Трассировка лучей» из книги «Графические драгоценности IV»)!

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

Обратная сторона визитки

Вот так выглядит сам код:

 #include <stdlib.h>   // card > aek.ppm
    #include <stdio.h>
    #include <math.h>
    typedef int i;typedef float f;struct v{
    f x,y,z;v operator+(v r){return v(x+r.x
    ,y+r.y,z+r.z);}v operator*(f r){return
    v(x*r,y*r,z*r);}f operator%(v r){return
    x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r
    ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.
    y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v
    operator!(){return*this*(1/sqrt(*this%*
    this));}};i G[]={247570,280596,280600,
    249748,18578,18577,231184,16,16};f R(){
    return(f)rand()/RAND_MAX;}i T(v o,v d,f
    &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01
    <p)t=p,n=v(0,0,1),m=1;for(i k=19;k--;)
    for(i j=9;j--;)if(G[j]&1<<k){v p=o+v(-k
    ,0,-j-4);f b=p%d,c=p%p-1,q=b*b-c;if(q>0
    ){f s=-b-sqrt(q);if(s<t&&s>.01)t=s,n=!(
    p+d*t),m=2;}}return m;}v S(v o,v d){f t
    ;v n;i m=T(o,d,t,n);if(!m)return v(.7,
    .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R(
    ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l%
    n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b
    >0),99);if(m&1){h=h*.2;return((i)(ceil(
    h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b
    *.2+.1);}return v(p,p,p)+S(h,r)*.5;}i
    main(){printf("P6 512 512 255 ");v g=!v
    (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a
    )*.002,c=(a+b)*-256+g;for(i y=512;y--;)
    for(i x=512;x--;){v p(13,13,13);for(i r
    =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)*
    99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b
    *(y+R())+c)*16))*3.5+p;}printf("%c%c%c"
    ,(i)p.x,(i)p.y,(i)p.z);}}

Код выше выглядит… пугающе, но компилируется и запускается без проблем! Вы можете сохранить его на рабочем столе как card.cpp, открыть консоль и ввести:

    c++ -O3 -o card card.cpp
    ./card > card.ppm

Через 27 секунд на экране появится следующее изображение:

Возможности визитки-трассировщика лучей

Возможности просто поражают!

  • мир, состоящий из строго организованных сфер;
  • текстурированный пол;
  • небо с градиентом;
  • мягкие тени;
  • OMG, глубина резкости! Вы шутите?!

И все это на одной стороне визитной карточки! Посмотрим, как это работает.

Класс Vector

Рассмотрим первую часть кода:

    #include <stdlib.h>   // card > aek.ppm
    #include <stdio.h>
    #include <math.h>
    typedef int i;typedef float f;struct v{
    f x,y,z;v operator+(v r){return v(x+r.x
    ,y+r.y,z+r.z);}v operator*(f r){return
    v(x*r,y*r,z*r);}f operator%(v r){return
    x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r
    ){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.
    y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v
    operator!(){return*this*(1/sqrt(*this%*
    this));}};

Главная хитрость здесь — это сокращение ключевых слов типов int и float до i и f с помощью typedef. Другой ход, с помощью которого можно можно уменьшить количество кода — это класс v, используемый не только в качестве вектора, но и для обработки пикселей.

Вот код, приведенный выше, но отформатированный и с комментариями:

   #include <stdlib.h>   // card > aek.ppm
   #include <stdio.h>
   #include <math.h>

   typedef int i;       // Экономим место с помощью сокращения int до i
   typedef float f;     // Экономим еще больше места с f вместо float

   // Класс вектора с конструктором и операторами
   struct v{
      f x,y,z;  // Три координаты вектора
      v operator+(v r){return v(x+r.x,y+r.y,z+r.z);} // Сумма векторов
      v operator*(f r){return v(x*r,y*r,z*r);}       // Масштабирование векторов
      f operator%(v r){return x*r.x+y*r.y+z*r.z;}    // Скалярное произведение векторов
      v(){}                                  // Пустой конструктор
      v operator^(v r){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);} // Векторное произведение векторов
      v(f a,f b,f c){x=a;y=b;z=c;}            // Конструктор
      v operator!(){return *this*(1 /sqrt(*this%*this));} // Нормализация вектора
   };

Rand() и данные для генерации мира

i G[]={247570,280596,280600,
    249748,18578,18577,231184,16,16};f R(){
    return(f)rand()/RAND_MAX;}

Следующий код также экономит много места с помощью объявления функции R, которая возвращает случайное значение от 0 до 1 типа float. Это полезно при стохастическом сэмплировании, использующемся для blur-эффекта и мягких теней.

Массив G содержит в себе закодированное целыми числами положение сфер в мире. Совокупность всех чисел — это битовый вектор из 9 строк и 19 столбцов.

Вот код, приведенный выше, но отформатированный и с комментариями:

    // Набор позиций сфер, описывающий мир
    // Все эти числа, по сути, являются по сути битовым вектором
   i G[]={247570,280596,280600,249748,18578,18577,231184,16,16};
   
   // Генератор случайных чисел, возвращающий число с плавающей точкой в диапазоне 0-1
   f R(){return(f)rand()/RAND_MAX;}

Главный метод

i main(){printf("P6 512 512 255 ");v g=!v
    (-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a
    )*.002,c=(a+b)*-256+g;for(i y=512;y--;)
    for(i x=512;x--;){v p(13,13,13);for(i r
    =64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)*
    99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b
    *(y+R())+c)*16))*3.5+p;}printf("%c%c%c"
    ,(i)p.x,(i)p.y,(i)p.z);}}

Главный метод использует простой известный основанный на тексте формат изображений PPM. Изображение состоит из заголовка вида P6 [Ширина] [Высота] [Максимальное значение], за которым следует RGB-значение каждого пикселя.

Для каждого пикселя на изображении программа сэмплирует (S) цвет 64 лучей, аккумулирует результат и выводит его в stdout.

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

Вот код, приведенный выше, но отформатированный и с комментариями:

  // Главная функция. Выводит изображение.
  // Использовать программу просто: ./card > erk.ppm
  i main(){
    
     printf("P6 512 512 255 "); // Заголовок PPM
    
     // Оператор "!" осуществляет нормализацию вектора
     v g=!v(-6,-16,0),       // Направление камеры
       a=!(v(0,0,1)^g)*.002, // Вектор, отвечающий за высоту камеры...
       b=!(g^a)*.002,        // Правый вектор, получаемый с помощью векторного произведения
       c=(a+b)*-256+g;       // WTF? Вот здесь https:// news.ycombinator.com/item?id=6425965 написано про это подробнее.

     for(i y=512;y--;)    // Для каждого столбца
     for(i x=512;x--;){   // Для каждого пикселя в строке
       
        // Используем класс вектора, чтобы хранить цвет в RGB
        v p(13,13,13);     // Стандартный цвет пикселя — почти черный
       
        // Бросаем по 64 луча из каждого пикселя 
        for(i r=64;r--;){ 
          
            // Немного меняем влево/вправо и вверх/вниз координаты начала луча (для эффекта глубины резкости)
            v t=a*(R()-.5)*99+b*(R()-.5)*99;
                                           
            // Назначаем фокальной точкой камеры v(17,16,8) и бросаем луч 
            // Аккумулируем цвет, возвращенный в переменной t
            p=S(v(17,16,8)+t, // Начало луча
                !(t*-1+(a*(R()+x)+b*(y+R())+c)*16) // Направление луча с небольшим искажением
                                                   // ради эффекта стохастического сэмплирования
                )*3.5+p; // +p для аккумуляции цвета
        }
       
        printf("%c%c%c",(i)p.x,(i)p.y,(i)p.z);
       
     }
  }

Сэмплер

v S(v o,v d){f t
    ;v n;i m=T(o,d,t,n);if(!m)return v(.7,
    .6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R(
    ),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l%
    n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b
    >0),99);if(m&1){h=h*.2;return((i)(ceil(
    h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b
    *.2+.1);}return v(p,p,p)+S(h,r)*.5;}

Сэмплер S — это функция, возвращающая цвет пикселя по данным координатам точки начала луча о и его направлению d. Если она натыкается на сферу, то она вызывает себя рекурсивно, а в ином случае (если луч не имеет препятствий на своем пути) в зависимости от направления возвращает либо цвет неба, либо цвет пола (базируясь на его клетчатой текстуре).

Обратите внимание на вызов функции R при расчете направления света. Таким образом создается эффект «мягких теней».

Вот код, приведенный выше, но отформатированный и с комментариями:

   // (S)эмплируем мир и возвращаем цвет пикселя по
   // по лучу, начинающемуся в точке o (Origin) и имеющему направление d (Direction)
   v S(v o,v d){
      f t;
      v n;
    
      // Проверяем, натыкается ли луч на что-нибудь
      i m=T(o,d,t,n);

    
      if(!m) // m==0
      // Сфера не была найдена, и луч идет вверх: генерируем цвет неба  
      return v(.7,.6,1)*pow(1-d.z,4);

      // Возможно, луч задевает сферу
    
      v h=o+d*t,                    // h - координата пересечения
      l=!(v(9+R(),9+R(),16)+h*-1),  // 'l' = направление света (с небольшим искажеем для эффекта мягких теней)
      r=d+n*(n%d*-2);               // r = полувектор
 
      // Расчитываем коэффицент Ламберта
      f b=l%n;
    
      // Рассчитываем фактор освещения (коэффицент Ламберта > 0 или находимся в тени)?
      if(b<0||T(h,l,t,n))
         b=0;
   
      // Рассчитываем цвет p (с учетом диффузии и отражения света) 
      f p=pow(l%r*(b>0),99);
    
      if(m&1){   // m == 1
         h=h*.2; // Сфера не была задета, и луч уходит вниз, в пол: генерируем цвет пола
         return((i)(ceil(h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b*.2+.1);
      }
   
      // m == 2 Была задета сфера: генерируем луч, отскакивающий от поверхности сфера
      return v(p,p,p)+S(h,r)*.5; // Ослабляем цвет на 50%, так как он отскакивает от поверхности (* .5)
  }

Трэйсер

i T(v o,v d,f
    &t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01
    <p)t=p,n=v(0,0,1),m=1;for(i k=19;k--;)
    for(i j=9;j--;)if(G[j]&1<<k){v p=o+v(-k
    ,0,-j-4);f b=p%d,c=p%p-1,q=b*b-c;if(q>0
    ){f s=-b-sqrt(q);if(s<t&&s>.01)t=s,n=!(
    p+d*t),m=2;}}return m;}

Функция T (Tracer) отвечает за бросание луча из данной точки (o) в данном направлении (d). Она возвращает целое число, которое является кодом для результата бросания луча. 0 — луч ушел в небо, 1 — луч ушел в пол, 2 — луч наткнулся на сферу. Если была задета сфера, то функция обновляет переменные t (параметр, используемый для вычисления дистанции пересения) и n (полу-вектор при отскакивании от сферы).

Вот код, приведенный выше, но отформатированный и с комментариями:

   // Тест на пересечение для линии [o,v]
   // Возвращаем 2, если была задета сфера (а также дистанцию пересечения t и полу-вектор n).
   // Возвращаем 0, если луч ничего не задевает и идет вверх, в небо
   // Возвращаем 1, если луч ничего не задевает и идет вниз, в пол
   i T(v o,v d,f& t,v& n){ 
     t=1e9;
     i m=0;
     f p=-o.z/d.z;
      if(.01<p)›
       t=p,n=v(0,0,1),m=1;

     // Мир зашифрован в G, в 9 линий и 19 столбцов
     for(i k=19;k--;)  // Для каждого столбца
     for(i j=9;j--;)   // Для каждой строки
      
     if(G[j]&1<<k){ // Для этой линии j есть ли в столбце i cфера?
        
        // Сфера есть, но задевает ли ее луч?
         
        v p=o+v(-k,0,-j-4);
        f b=p%d,c=p%p-1,q=b*b-c;
     
        // Задевает ли луч сферу?
        if(q>0){
           // Да. Считаем расстояние от камеры до сферы
           f s=-b-sqrt(q);
             
           if(s<t && s>.01)
             // Это минимальное расстояние, сохраняем его. А также
             // вычиваем вектор отскакивающего луча и записываем его в 'n'  
             t=s,
             n=!(p+d*t),
             m=2;
         }
     }
    
     return m;
  }

Число Leet

Многие программисты пытались сократить код еще больше. Сам автор остановился на версии, предоставленной в этой статье. Знаете, почему?

      fabien$ wc card.cpp

      35      95    1337 card.cpp

Его размер — 1337 байт! Число Leet! Ха!

Материалы для ознакомления

Так как я провел очень много времени, читая про сам метод бросания лучей, вот несколько интересных ссылок и книг:

Перевод статьи «Decyphering the business card racetracer»