1 - raytracer 프로젝트 - 시작!

이 글은 raytracer 시리즈입니다.

목차

• 이미지 출력하기
• 수학적 배경
  • 벡터
  • 행렬
• 카메라
• 구

이미지 출력하기

이미지를 GUI로 보여주기보다는 그냥 BMP파일로 출력하는 식으로 진행할 것이다.

BMP 출력은 궁금하다면 이 글을 참고하자. bmp.c, bmp.h를 그대로 사용할 것이다.

#pragma once

#include <stddef.h>

#include "external/bmp.h"

typedef bmp_pixel_t (*renderer_t)(void *context, float x, float y);

void write_image(const char *output_path, renderer_t renderer, void *context,
                 size_t width, size_t height);

일단 간단한 이미지 출력 함수를 작성해보자

#include "write_image.h"

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

void write_image(const char *output_path, renderer_t renderer, void *context,
                 size_t width, size_t height) {
  const size_t pixel_count = width * height;
  bmp_t *bmp = malloc(sizeof(bmp_t) + sizeof(bmp_pixel_t) * pixel_count);
  if (!bmp) {
    fprintf(stderr, "memory allocation failure\n");
    exit(EXIT_FAILURE);
  }
  bmp->width = width;
  bmp->height = height;

  const size_t sqrt_pixel_count = (size_t)sqrtf((float)pixel_count);
  for (size_t y = 0; y < height; y++) {
    for (size_t x = 0; x < width; x++) {
      const size_t index = y * width + x;
      if (index % sqrt_pixel_count == 0) {
        printf("progress: %zu / %zu\n", index, pixel_count);
      }
      bmp->extra[index] =
          renderer(context, x / (width - 1.0f), y / (height - 1.0f));
    }
  }

  char *buffer;
  size_t length;
  if (serialize_bmp(bmp, &buffer, &length)) {
    fprintf(stderr, "memory allocation failure\n");
    exit(EXIT_FAILURE);
  }

  free(bmp);

  FILE *fp = fopen(output_path, "wb");
  if (!fp) {
    fprintf(stderr, "Failed to open %s\n", output_path);
    exit(EXIT_FAILURE);
  }
  if (fwrite(buffer, 1, length, fp) != length) {
    fprintf(stderr, "Failed to write to %s\n", output_path);
    exit(EXIT_FAILURE);
  }
  free(buffer);
  fclose(fp);

  printf("Successfully wrote image to %s\n", output_path);
}

이미지를 한 픽셀씩 그리는 함수 renderer_t를 받고, 이미지 파일로 출력하면서 어디까지 진행됐는지 출력하는 함수이다.

#include <stdint.h>
#include <stdio.h>

#include "../../write_image.h"

bmp_pixel_t renderer(void *context, float x, float y) {
  return (bmp_pixel_t){
    (uint8_t)(x * 255 + 0.5f),
    (uint8_t)(y * 255 + 0.5f),
    255,
  };
}

int main(int argc, char **argv) {
  if (argc != 2) {
    return 1;
  }

  write_image(argv[1], renderer, NULL, 256, 256);
}

일단 아무 이미지나 나오게 만들어서, 아래 파일들로 테스트하면 아래 이미지가 나온다

cmake_minimum_required(VERSION 3.10)
project(write_image)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_CURRENT_BINARY_DIR})
add_executable(main main.c ../../write_image.c ../../external/bmp.c)
if (NOT WIN32)
  target_link_libraries(main PUBLIC m)
endif()

#!/bin/bash

set -e

cd "$(dirname "$0")"

cmake -DCMAKE_BUILD_TYPE=Release -B builddir
cmake --build builddir --config Release
builddir/main ./result.bmp

수학적 배경

레이트레이싱은 픽셀마다 광선을 발사해서, 그 광선이 부딛힌 곳의 색을 구하는 방식으로 동작한다.

광선에는 두가지 속성이 있다.

1. 광선이 시작한 위치
2. 광선이 향하는 방향

여기서 위치와 방향이 모두 벡터이다.

벡터

벡터는 가변 갈이 배열이라는 의미도 있지만, 여기서의 벡터의 의미는 수학의 벡터이다.

3차원 공간상의 위치와 방향을 나타낼 것이기에 3차원 벡터를 사용할 것이다.

#pragma once

typedef struct f3 {
  float x;
  float y;
  float z;
} f3_t;

f3_t f3(float x, float y, float z);
f3_t f3_add(f3_t a, f3_t b);
f3_t f3_sub(f3_t a, f3_t b);
f3_t f3_mul1(f3_t a, float b);
f3_t f3_mul3(f3_t a, f3_t b);
f3_t f3_div1(f3_t a, float b);
f3_t f3_div3(f3_t a, f3_t b);
float f3_length(f3_t f3);
float f3_dot(f3_t a, f3_t b);
f3_t f3_cross(f3_t a, f3_t b);
f3_t f3_normalize(f3_t f3);
f3_t f3_rotate(f3_t point, f3_t rotation);

여기서 normalize(정규화)는 길이가 1이 아닌 벡터의 길이를 1로 만드는 것이다.

벡터를 아래의 세 가지 의미 중 하나로 쓸 것이다.

1. 위치
2. 이동 (위치의 차이, 방향 + 거리)
3. 방향

여기서 방향을 나타내는 벡터는 길이가 항상 1이어야 한다.

위치 $B$에서 위치 $A$를 빼면 $A$에서 $B$로의 이동을 얻을 수 있고, 이를 정규화하면 $A$에서 $B$로 향하는 방향을 얻을 수 있다.

dot(내적)과 cross(외적)는 벡터의 특수한 곱셈이다.

행렬

행렬은 벡터의 벡터라고 생각하면 될 것 같다. 벡터가 1차원이라면 행렬은 2차원이다.

여기서는 3*3 행렬을 f3x3으로 정의한다.

#pragma once

#include "f3.h"

typedef struct f3x3
{
  f3_t x;
  f3_t y;
  f3_t z;
} f3x3_t;

f3x3_t f3x3(f3_t x, f3_t y, f3_t z);
f3_t f3x3_col_x(f3x3_t mat);
f3_t f3x3_col_y(f3x3_t mat);
f3_t f3x3_col_z(f3x3_t mat);
f3x3_t f3x3_mul(f3x3_t a, f3x3_t b);
f3x3_t f3x3_transpose(f3x3_t mat);
f3x3_t f3x3_rotate(float yaw, float pitch, float roll);
f3x3_t f3x3_rotate_reverse(float yaw, float pitch, float roll);
f3_t f3x3_mul_f3(f3_t vec, f3x3_t mat);

extern f3x3_t f3x3_i;

행렬의 곱 $AB$의 각 원소는 $A$의 해당 행을 이루는 벡터와 $B$의 해당 행을 이루는 벡터의 내적이 된다.

즉, $A \times B$ 행렬은 $B \times C$ 행렬과 곱할 수 있으며, 결과는 $A \times C$ 행렬이 된다.

행렬의 곱셈은 교환법칙은 성립하지 않지만 결합법칙은 성립한다.

$AB != BA$일 수 있는데, $A(BC) == (AB)C$는 항상 성립한다.

transpose는 전치행렬을 의미한다. 행/열을 바꾼 것이다.
직교행렬(x(첫번째 줄), y(두번째 줄), z(세번째 줄)가 서로 직교하는 행렬)은 전치행렬이 역행렬이다.

단위행렬 $I$는 항등원으로, 어떤 행렬과 곱해도 곱한 행렬이 그대로 나온다.
어떤 행렬의 역행렬은 그 행렬의 역원, 즉 그 행렬과 곱해서 단위행렬이 되게 하는 행렬이다.

구현은 대충 이렇다.

#include "f3.h"

#include <math.h>

#include "f3x3.h"

f3_t f3(float x, float y, float z) { return (f3_t){x, y, z}; }

f3_t f3_add(f3_t a, f3_t b) { return (f3_t){a.x + b.x, a.y + b.y, a.z + b.z}; }

f3_t f3_sub(f3_t a, f3_t b) { return (f3_t){a.x - b.x, a.y - b.y, a.z - b.z}; }

f3_t f3_mul1(f3_t a, float b) { return (f3_t){a.x * b, a.y * b, a.z * b}; }

f3_t f3_mul3(f3_t a, f3_t b) { return (f3_t){a.x * b.x, a.y * b.y, a.z * b.z}; }

f3_t f3_div1(f3_t a, float b) { return (f3_t){a.x / b, a.y / b, a.z / b}; }

f3_t f3_div3(f3_t a, f3_t b) { return (f3_t){a.x / b.x, a.y / b.y, a.z / b.z}; }

float f3_length(f3_t f3) {
  return sqrtf(f3.x * f3.x + f3.y * f3.y + f3.z * f3.z);
}

float f3_dot(f3_t a, f3_t b) { return a.x * b.x + a.y * b.y + a.z * b.z; }

f3_t f3_cross(f3_t a, f3_t b) {
  return (f3_t){a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z,
                a.x * b.y - a.y * b.x};
}

f3_t f3_normalize(f3_t f3) { return f3_div1(f3, f3_length(f3)); }

f3_t f3_rotate(f3_t point, f3_t rotation) {
  const float yaw = rotation.x;
  const float pitch = rotation.y;
  const float roll = rotation.z;

  return (f3x3_mul_f3(point, f3x3_rotate(yaw, pitch, roll)));
}

#include "f3x3.h"

#include <math.h>

f3x3_t f3x3(f3_t x, f3_t y, f3_t z) { return (f3x3_t){x, y, z}; }

f3_t f3x3_col_x(f3x3_t f3x3) { return (f3_t){f3x3.x.x, f3x3.y.x, f3x3.z.x}; }

f3_t f3x3_col_y(f3x3_t f3x3) { return (f3_t){f3x3.x.y, f3x3.y.y, f3x3.z.y}; }

f3_t f3x3_col_z(f3x3_t f3x3) { return (f3_t){f3x3.x.z, f3x3.y.z, f3x3.z.z}; }

f3x3_t f3x3_mul(f3x3_t a, f3x3_t b) {
  const f3_t x = {
      f3_dot(a.x, f3x3_col_x(b)),
      f3_dot(a.x, f3x3_col_y(b)),
      f3_dot(a.x, f3x3_col_z(b)),
  };
  const f3_t y = {
      f3_dot(a.y, f3x3_col_x(b)),
      f3_dot(a.y, f3x3_col_y(b)),
      f3_dot(a.y, f3x3_col_z(b)),
  };
  const f3_t z = {
      f3_dot(a.z, f3x3_col_x(b)),
      f3_dot(a.z, f3x3_col_y(b)),
      f3_dot(a.z, f3x3_col_z(b)),
  };

  return f3x3(x, y, z);
}

f3x3_t f3x3_transpose(f3x3_t m) {
  return f3x3(f3x3_col_x(m), f3x3_col_y(m), f3x3_col_z(m));
}

static f3x3_t rotate_yaw(float rad) {
  const f3_t x = {
      cosf(rad),
      sinf(rad),
      0,
  };
  const f3_t y = {
      -sinf(rad),
      cosf(rad),
      0,
  };
  const f3_t z = {
      0,
      0,
      1,
  };

  return f3x3(x, y, z);
}

static f3x3_t rotate_pitch(float rad) {
  const f3_t x = {
      1,
      0,
      0,
  };
  const f3_t y = {
      0,
      cosf(rad),
      sinf(rad),
  };
  const f3_t z = {
      0,
      -sin(rad),
      cos(rad),
  };

  return f3x3(x, y, z);
}

static f3x3_t rotate_roll(float rad) {
  const f3_t x = {
      cosf(rad),
      0,
      sinf(rad),
  };
  const f3_t y = {
      0,
      1,
      0,
  };
  const f3_t z = {
      -sinf(rad),
      0,
      cosf(rad),
  };

  return f3x3(x, y, z);
}

f3x3_t f3x3_rotate(float yaw, float pitch, float roll) {
  const f3x3_t mat_yaw = rotate_yaw(yaw);
  const f3x3_t mat_pitch = rotate_pitch(pitch);
  const f3x3_t mat_roll = rotate_roll(roll);

  return (f3x3_mul(f3x3_mul(mat_yaw, mat_pitch), mat_roll));
}

f3x3_t f3x3_rotate_reverse(float yaw, float pitch, float roll) {
  return f3x3_transpose(f3x3_rotate(yaw, pitch, roll));
}

f3_t f3x3_mul_f3(f3_t vec, f3x3_t mat) {
  return f3(f3_dot(mat.x, vec), f3_dot(mat.y, vec), f3_dot(mat.z, vec));
}

f3x3_t f3x3_i = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};

벡터와 행렬에 대해서는 다음 글에서 더 자세히 알아보기로 하고 일단 넘어가자.

카메라

카메라에서 광선을 발사할 것이다. 그런데 화면상의 위치를 어떻게 광선의 위치, 방향으로 바꿀까?

카메라에는 일단 크게 두 종류가 있다.

1. 직교 카메라 (orthogonal)
2. 원근 카메라 (perspective)

상세한 구현은 나중에 하기로 하고, 우선 간이로 간단한 직교 카메라만 구현해 보자.

x가 오른쪽, y가 앞쪽, z가 위쪽인 좌표계를 생각해 보자.

#include <stdio.h>

#include "../../f3.h"
#include "../../write_image.h"

typedef struct ray {
  f3_t origin;
  f3_t direction;
} ray_t;

const float orthogonal_camera_scale = 2.0f;
const f3_t camera_position = {0, 0, 0};
const f3_t camera_direction = {0, 1, 0};
const bmp_pixel_t black = {0, 0, 0};
const bmp_pixel_t white = {255, 255, 255};

bmp_pixel_t renderer(void *context, float x_on_image, float y_on_image) {
  const float x = (2.0f * x_on_image - 1.0f) * orthogonal_camera_scale;
  // y_on_image는 높을수록 아래인데 z는 높을수록 위이므로, 반대로 처리
  const float z = -(2.0f * y_on_image - 1.0f) * orthogonal_camera_scale;
  const f3_t origin = f3_add(camera_position, f3(x, 0, z));
  const ray_t ray = {origin, camera_direction};

  if (ray.origin.z < 0.0f) {
    return black;
  }
  return white;
}

int main(int argc, char **argv) {
  if (argc != 2) {
    return 1;
  }

  write_image(argv[1], renderer, NULL, 480, 270);
}

ray(광선)를 만들고, origin(시작한 위치)가 0보다 작으면 검정색, 아니면 흰색의 이미지를 만들 것이다.

cmake_minimum_required(VERSION 3.10)
project(camera)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_CURRENT_BINARY_DIR})
add_compile_options("$<$<C_COMPILER_ID:MSVC>:/utf-8>")
add_executable(main main.c ../../write_image.c ../../f3.c ../../f3x3.c ../../external/bmp.c)
if (NOT WIN32)
  target_link_libraries(main PUBLIC m)
endif()

이하 build.sh 내용은 위의 build.sh와 동일.

구

구와 광선의 충돌 확인을 통해 구를 렌더링 해 보자.

#include <math.h>
#include <stdint.h>
#include <stdio.h>

#include "../../f3.h"
#include "../../write_image.h"

typedef struct ray {
  f3_t origin;
  f3_t direction;
} ray_t;

typedef struct sphere {
  f3_t center;
  float radius;
} sphere_t;

const float orthogonal_camera_scale = 2.0f;
const f3_t camera_position = {0, 0, 0};
const f3_t camera_direction = {0, 1, 0};
const bmp_pixel_t white = {255, 255, 255};
const sphere_t sphere = {{0, 5, 0}, 1};

static inline float sqr(float x) { return x * x; }

bool collide(ray_t ray, sphere_t sphere, f3_t *out_normal) {
  // 카메라가 구 안에 있는지
  if (f3_length(f3_sub(ray.origin, sphere.center)) < sphere.radius) {
    *out_normal = f3_mul1(ray.direction, -1);
    return true;
  }

  // 편의상 ray와 sphere의 위치에서 sphere의 위치를 빼서 sphere를 원점으로 만듦
  ray.origin = f3_sub(ray.origin, sphere.center);
  sphere.center = f3(0, 0, 0);

  /*
    이제 sphere를 x^2 + y^2 + z^2 - r^2 = 0으로 나타낼 수 있음
    x, y, z는 ray.origin + ray.direction * t의 x, y, z가 됨
    여기서 t는 광선이 이동한 거리, r은 ray.origin

    x^2 + y^2 + z^2 - r^2 = 0이라는 식에서 x^2는 이렇게 다시 쓸 수 있음
    (ray.origin.x + ray.direction.x * t)^2

    편의상 ray.origin을 (a, b, c), ray.direction을 (u, v, w)라고 하면
    (a + ut)^2 + (b + vt)^2 + (c + wt)^2 = r^2이 되고, 이를 풀어 쓰면
    a^2 + 2aut + (ut)^2 + b^2 + 2bvt + (vt)^2 + c^2 + 2cwt + (wt)^2 - r^2 = 0
    이를 만족하는 t를 찾아야 함

    근의 공식에 따라 A𝑥^2 + B𝑥 + C = 0이라고 할 때
    x는 (-B ± sqrt(B^2 - 4AC)) / 2A가 되고,
    여기서 𝑥가 t라고 할 때 a, b, c를 구해서 근의 공식대로 𝑥(=t)를 구하면 된다

    여기서의 A와 B, C를 각각 구하면 다음과 같다
    A = u^2 + v^2 + w^2
    B = 2au + 2bv + 2cw
    C = a^2 + b^2 + c^2
  */
  const float a =
      sqr(ray.direction.x) + sqr(ray.direction.y) + sqr(ray.direction.z);
  const float b =
      2 * (ray.origin.x * ray.direction.x + ray.origin.y * ray.direction.y +
           ray.origin.z * ray.direction.z);
  const float c = sqr(ray.origin.x) + sqr(ray.origin.y) + sqr(ray.origin.z) -
                  sqr(sphere.radius);
  const float discriminant = sqr(b) - 4 * a * c;

  if (discriminant < 0) {
    // 실수 x가 존재하지 않으면 교점이 없는 것, 즉 충돌하지 않음
    return false;
  }
  // 교점이 ± 때문에 두 개가 나오는데, ±가 +이면 뒷면이므로 여기서 ±는 아마도 -
  float t = (-b - sqrtf(discriminant)) / (2 * a);
  if (t < 0) { // 음수인 경우 카메라 뒤쪽이므로 +로 재시도
    t = (-b + sqrtf(discriminant)) / (2 * a);
    if (t < 0) { // 이래도 음수면 결국 충돌 안 한 거
      return false;
    }
  }

  // t를 통해 x^2 + y^2 + z^2 - r^2 = 0애서의 x, y, z를 구할 수 있음
  const float x = ray.origin.x + ray.direction.x * t;
  const float y = ray.origin.y + ray.direction.y * t;
  const float z = ray.origin.z + ray.direction.z * t;

  // x^2 + y^2 + z^2 - r^2를 각각 x, y, z에 대해 편미분하면 법선 벡터가 나옴
  out_normal->x = 2 * x; // f(x) = x^2 + y^2 + z^2, f'(x) = 2x
  out_normal->y = 2 * y; // f(y) = x^2 + y^2 + z^2, f'(y) = 2y
  out_normal->z = 2 * z; // f(z) = x^2 + y^2 + z^2, f'(z) = 2z
  // 하지만 길이가 1이 아닐 수 있으므로 정규화를 해야 함
  *out_normal = f3_normalize(*out_normal);

  return true;
}

bmp_pixel_t renderer(void *context, float x_on_image, float y_on_image) {
  const float x = (2.0f * x_on_image - 1.0f) * orthogonal_camera_scale;
  // y_on_image는 높을수록 아래인데 z는 높을수록 위이므로, 반대로 처리
  const float z = -(2.0f * y_on_image - 1.0f) * orthogonal_camera_scale;
  const f3_t origin = f3_add(camera_position, f3(x, 0, z));
  const ray_t ray = {origin, camera_direction};

  f3_t normal;
  if (!collide(ray, sphere, &normal)) {
    return white;
  }
  return (bmp_pixel_t){
      // normal의 각 요소는 -1에서 1 사이이므로, 0에서 255 사이로 만든다
      (uint8_t)((normal.x * 0.5f + 0.5f) * 255 + 0.5f),
      (uint8_t)((normal.y * 0.5f + 0.5f) * 255 + 0.5f),
      (uint8_t)((normal.z * 0.5f + 0.5f) * 255 + 0.5f),
  };
}

int main(int argc, char **argv) {
  if (argc != 2) {
    return 1;
  }

  write_image(argv[1], renderer, NULL, 1024, 1024);
}

자세한 계산 과정은 주석 참고.

이하 CMakeLists.txt 내용은 위의 CMakeLists.txt와 동일.

실행하면 이런 결과를 얻을 수 있다.

---

다음 글 보기: 클릭!