// пример рэйтрэйсинга (точнее, рэйкастинга)

#include <math.h>  // для sqrt()
#include <stdio.h> // для работы с файлами
#include "video.h"

// определяем тип "32-битное целое число" в зависимости от компилятора

#ifdef __WATCOMC__
typedef int int32;
#endif

#ifdef __BORLANDC__
typedef long int32;
#endif

// 3D-вектор
typedef struct {
  float x, y, z;
} vector;

// RGB-цвет
typedef struct {
  float r, g, b;
} color;

// сфера
typedef struct {
  vector center;     // центр
  color  col;        // цвет
  float  radius;     // радиус
  float  kd, ks, ns; // характеристики поверхности
} sphere;

// источник света
typedef struct {
  vector where; // местоположение
  color  col;   // цвет
} light;

#define NUMSPHERES 3        // число сфер
#define NUMLIGHTS  3        // число источников света
#define INFINITY   1000000L // достаточно большое число
#define EPSILON    0.001    // достаточно маленькое число

color ambient = {0.3, 0.2, 0.1}; // фоновая освещенность

// задаем сферы, составляющие сцену
sphere spheres[NUMSPHERES] = {
  { {  0,  0,    0}, {1, 1, 1}, 64, 0.7, 0.5, 4 },
  { { 10, 65,  -60}, {1, 1, 1}, 30, 0.8, 0.7, 2 },
  { {-40, 55, -100}, {1, 1, 1},  8, 0.6, 0.3, 3 }
};

// задаем источники света
light lights[NUMLIGHTS] = {
  { { 10,  100, -120}, {0, 1, 0} },
  { {200,    0,    0}, {1, 0, 0} },
  { {-50,   50, -120}, {0, 0, 1} }
};

// заголовок TGA-файла
char header[12] = {0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0};

// скалярное произведение векторов
float dot(vector a, vector b) {
  return (a.x * b.x + a.y * b.y + a.z * b.z);
}

// максимум из двух чисел
float fmax(float a, float b) {
  return (a > b) ? a : b;
}

// косинус угла между векторами или ноль, если угол больше 90 градусов
float vcos(vector a, vector b) {
  return fmax(dot(a, b), 0) / sqrt(dot(a, a) * dot(b, b));
}

// линейная комбинация векторов, linear = ka * a + kb * b
vector linear(vector a, vector b, float ka, float kb) {
  vector res;

  res.x = a.x * ka + b.x * kb;
  res.y = a.y * ka + b.y * kb;
  res.z = a.z * ka + b.z * kb;
  return res;
}

// приведение вектора к единичному
void normalize(vector *a) {
  float l = sqrt(dot(*a, *a));

  a->x /= l;
  a->y /= l;
  a->z /= l;
}

// поиск пересечения сферы и луча (o + k * t), k > 0, возвращает k, в случае
// отсутствия пересечения возвращает INFINITY
float isectSphere(vector o, vector t, sphere s) {
  float a, b, c, d, t1, t2;
  vector m;

  // точка пересечения лежит и на луче, и на сфере; таким образом,
  // имеем уравнение
  //
  //     length((o + k * t) - s.center) = radius,
  //
  // где length - длина вектора, известно, что length(a) = sqrt(dot(a,a)).
  // отсюда после преобразований получаем квадратное уравнение относительно
  // k, (a*k*k + b*k + c) = 0 с дискриминантом d.
  m = linear(o, s.center, 1, -1);
  a = dot(t, t);
  b = 2 * dot(m, t);
  c = dot(m, m) - s.radius * s.radius;
  d = b * b - 4 * a * c;

  // если корней нет, то и пересечений нет
  if (d < 0) return INFINITY;

  // ищем корни
  d = sqrt(d);
  t1 = (-b - d) / (2 * a);
  t2 = (-b + d) / (2 * a);

  // если больший корень отрицателен, то оба отрицательны и k > 0 нету
  if (t2 < 0) return INFINITY;

  // если меньший корень положителен, то это и есть наш k. если нет,
  // то остается только больший (он уже заведомо неотрицателен).
  if (t1 >= 0) return t1;
  return t2;
}

// поиск "ближайшего" пересечения луча со сценой, возвращает номер объекта, с
// которым пересекся луч, или -1, если пересечений нет. точка o может лежать
// на каком-либо объекте, это НЕ будет считаться за пересечение.
int isectAll(vector o, vector t, float *k) {
  int    i, closest = -1;
  float  tmpf;

  *k = INFINITY;
  for (i = 0; i < NUMSPHERES; i++) {
    tmpf = isectSphere(o, t, spheres[i]);
    // если пересечение с очередным объектом находится ближе, чем текущее,
    // то запоминаем его
    if (*k > tmpf && tmpf > EPSILON) {
      closest = i;
      *k = tmpf;
    }
  }
  return closest;
}

// собственно трассировка лучей. возвращает освещенность, приходящую в
// точку o по направлению t - то есть, трассирует луч (o + k * t), k > 0.
color trace(vector o, vector t) {
  color  res = {0, 0, 0}, tmpc;
  vector p, n, l, r;
  float  k, tmpf;
  int    i, closest;
  sphere sp;

  // ищем пересечение заданного со сценой
  closest = isectAll(o, t, &k);

  // нет пересечения, нет и освещенности
  if (closest == -1) return res;

  // sp - объект, с которым пересекся наш луч
  sp = spheres[closest];

  // считаем p и n, точку пересечения и нормаль к объекту в этой точке
  p = linear(o, t, 1, k);
  n = linear(p, sp.center, 1, -1);
  normalize(&n);

  // начинаем накапливать освещенность, первое слагаемое: intensity = Ka
  res = ambient;

  // для каждого из источников света применяем уравнение Фонга
  for (i = 0; i < NUMLIGHTS; i++) {
    // l - вектор, проведенный из точки в источник света
    l = linear(lights[i].where, p, 1, -1);

    // если луч из точки p в направлении l НЕ пересекается со сценой,
    // то точка p "видна" источнику света, считаем освещенность
    if (isectAll(p, l, &tmpf) == -1) {

      // диффузная (рассеянная компонента), intensity += kd * cos(n, l)
      tmpf = vcos(n, l) * sp.kd;
      res.r += lights[i].col.r * tmpf;
      res.g += lights[i].col.g * tmpf;
      res.b += lights[i].col.b * tmpf;

      // отраженная компонента, отблеск; intensity += ks * pow(cos(r, v), ns)
      // здесь r - симметричный l относительно n вектор, v - вектор из точки p
      // в направлении камеры; таким образом, в нашем случае v = t.
      r = linear(n, l, -2 * dot(n, l), 1);
      tmpf = pow(vcos(r, t), sp.ns) * sp.ks;
      res.r += lights[i].col.r * tmpf;
      res.g += lights[i].col.g * tmpf;
      res.b += lights[i].col.b * tmpf;
    }
  }

  // полученную освещенность осталось домножить на цвет объекта
  res.r *= sp.col.r;
  res.g *= sp.col.g;
  res.b *= sp.col.b;

  return res;
}

// функция упаковки полей структуры color в обычный 32-битный пиксел
int32 packColor(color c) {
  int32 r = (int)(((c.r > 1) ? 1 : c.r) * 255);
  int32 g = (int)(((c.g > 1) ? 1 : c.g) * 255);
  int32 b = (int)(((c.b > 1) ? 1 : c.b) * 255);

  return (b + (g << 8) + (r << 16));
}

void main() {
  int    i, j, k;
  int32  c;
  vector o, t;
  FILE   *f;
  char   *dest;

  // создаем traced.tga, пишем туда заголовок
  f = fopen("traced.tga", "wb+");
  fwrite(&header, 1, sizeof(header), f);
  i = XSIZE; // размер картинки по x
  j = YSIZE; // размер картинки по y
  fwrite(&i, 1, 2, f);
  fwrite(&j, 1, 2, f);
  i = 24;    // битность цвета
  fwrite(&i, 1, 2, f);

  // инициализируем видеосистему
  initVideo();

  // ставим палитру для предварительного просмотра (получаем как бы RGB332)
  outp(0x3C8, 0x00);
  for (i = 0; i < 8; i++) {
    for (j = 0; j < 8; j++) {
      for (k = 0; k < 4; k++) {
	outp(0x3C9, i << 3);
	outp(0x3C9, j << 3);
	outp(0x3C9, k << 4);
      }
    }
  }

  // трассируем каждый пиксел построчно снизу вверх (так как в TGA-файле
  // строки идут как раз снизу вверх)
  for (j = YSIZE - 1; j >= 0; j--) {
    dest = videoMem + j * XSIZE;
    for (i = 0; i < XSIZE; i++) {
      // точка o - камера, она же начало луча, (0,0,-DIST)
      o.x = o.y = 0;
      o.z = -DIST;

      // t - направление луча
      t.x = i - (XSIZE / 2);
      t.y = (YSIZE / 2) - j;
      t.z = DIST;

      // трассируем луч
      c = packColor(trace(o, t));

      // выводим на экран переведенный в RGB332 цвет
      *dest++ = ((c >> 6) & 0x03) + ((c >> 11) & 0x1C) + ((c >> 16) & 0xE0);

      // пишем очередной пиксел в файл
      fwrite(&c, 1, 3, f);
    }
  }

  fclose(f);   // закрываем файл
  getch();     // ждем нажатия клавиши
  doneVideo(); // возвращаемся в текстовый режим
}