【CG】150行C++路径追踪

  1#include <cmath>
  2#include <iostream>
  3#include <random>
  4#include <vector>
  5constexpr float EPSILON = 1e-5f;
  6constexpr float PI = 3.14159265358979323846f;
  7struct Vec3 {
  8  float x, y, z;
  9  Vec3(float x, float y, float z) : x(x), y(y), z(z) {}
 10  Vec3() : x(0), y(0), z(0) {}
 11  Vec3 operator+(const Vec3 &other) const {
 12    return Vec3(x + other.x, y + other.y, z + other.z);
 13  }
 14  Vec3 operator-(const Vec3 &other) const {
 15    return Vec3(x - other.x, y - other.y, z - other.z);
 16  }
 17  Vec3 operator*(float other) const {
 18    return Vec3(x * other, y * other, z * other);
 19  }
 20  Vec3 operator*(Vec3 other) {
 21    return Vec3(x * other.x, y * other.y, z * other.z);
 22  }
 23  Vec3 operator-() const { return Vec3(-x, -y, -z); }
 24};
 25Vec3 normalize(Vec3 v) { return v * (1.0f / sqrt(v.x * v.x + v.y * v.y + v.z * v.z)); }
 26float dot(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
 27float clamp(float x, float min, float max) {
 28  return x < min ? min : (x > max ? max : x);
 29}
 30float linear_to_srgb(float x) {
 31  x = clamp(x, 0.0f, 1.0f);
 32  if (x <= 0.0031308f)
 33    return 12.92f * x;
 34  return 1.055f * std::pow(x, 1.0f / 2.4f) - 0.055f;
 35}
 36struct Ray {
 37  Vec3 origin;
 38  Vec3 direction;
 39  Ray(Vec3 origin, Vec3 direction) : origin(origin), direction(direction) {}
 40  Vec3 at(float t) const { return origin + direction * t; }
 41};
 42struct Sphere {
 43  Sphere(float radius, Vec3 center, Vec3 color, bool light = false): radius(radius), center(center), color(color), is_light(light) {}
 44  float radius;
 45  Vec3 center;
 46  Vec3 color;
 47  bool is_light = false;
 48};
 49struct HitInfo {
 50  bool is_hit;
 51  float t;
 52  Vec3 position;
 53  Vec3 normal;
 54  Sphere const* sphere = nullptr;
 55};
 56HitInfo intersect(const Ray &ray, const Sphere &sphere) {
 57  Vec3 oc = ray.origin - sphere.center;
 58  float b = 2.0f * dot(oc, ray.direction);
 59  float c = dot(oc, oc) - sphere.radius * sphere.radius;
 60  float discriminant = b * b - 4.0f * c;
 61  if (discriminant < 0.0f)
 62    return {false, 0.0f, Vec3(), Vec3(), &sphere};
 63  float sqrt_d = sqrt(discriminant);
 64  float t0 = (-b - sqrt_d) * 0.5f,t1 = (-b + sqrt_d) * 0.5f;
 65  float t = t0;
 66  if (t < EPSILON)
 67    t = t1;
 68  if (t < EPSILON)
 69    return {false, 0.0f, Vec3(), Vec3(), &sphere};
 70  return {true, t, ray.at(t), normalize(ray.at(t) - sphere.center), &sphere};
 71}
 72float random_float() {
 73  thread_local std::mt19937 gen(std::minstd_rand{}());
 74  thread_local std::uniform_real_distribution<float> dis(0.0f, 1.0f);
 75  return dis(gen);
 76}
 77Vec3 sampling_hemisphere(const Vec3 &normal) {
 78  float z = random_float();
 79  float r = std::sqrt(1.0f - z * z);
 80  float phi = 2.0f * PI * random_float();
 81  Vec3 local(r * std::cos(phi), r * std::sin(phi), z);
 82  Vec3 tangent;
 83  if (std::abs(normal.x) > 0.9f)
 84    tangent = normalize(Vec3(0.0f, 1.0f, 0.0f) - normal * normal.y);
 85  else
 86    tangent = normalize(Vec3(1.0f, 0.0f, 0.0f) - normal * normal.x);
 87  Vec3 bitangent = Vec3(normal.y * tangent.z - normal.z * tangent.y,normal.z * tangent.x - normal.x * tangent.z,normal.x * tangent.y - normal.y * tangent.x);
 88  return normalize(tangent * local.x + bitangent * local.y + normal * local.z);
 89}
 90Vec3 Lo(Ray ray, const std::vector<Sphere> &spheres, int depth) {
 91  if (depth <= 0)
 92    return Vec3();
 93  float rr_prob = 0.9f;
 94  if (random_float() > rr_prob)
 95    return Vec3();
 96  HitInfo nearest_hit;
 97  nearest_hit.t = 1e10;
 98  nearest_hit.is_hit = false;
 99  for (auto &sphere : spheres) {
100    HitInfo hit = intersect(ray, sphere);
101    if (hit.is_hit && hit.t < nearest_hit.t)
102      nearest_hit = hit;
103  }
104  if (!nearest_hit.is_hit)
105    return Vec3();
106  if (nearest_hit.sphere->is_light) {
107    Vec3 Le = nearest_hit.sphere->color * 30.0f;
108    return Le;
109  }
110  Vec3 wo = -ray.direction,wi = sampling_hemisphere(nearest_hit.normal);
111  float pdf = 1.0f / (2.0f * PI);
112  Vec3 f_r = nearest_hit.sphere->color * (1.0f / PI);
113  float cos_theta = dot(wi, nearest_hit.normal);
114  Ray new_ray = Ray(nearest_hit.position + nearest_hit.normal * EPSILON, wi);
115  Vec3 Li = Lo(new_ray, spheres, depth - 1);
116  return f_r * Li * cos_theta * (1.0f / (rr_prob * pdf));
117}
118int main() {
119  int image_width = 256, image_height = 256,samples_per_pixel = 1024;
120  std::vector<Vec3> framebuffer(image_width * image_height);
121  Sphere left_wall = {1000.0f, Vec3(-1003.0f, 0.0f, -6.0f), Vec3(0.75f, 0.12f, 0.10f)};
122  Sphere right_wall = {1000.0f, Vec3(1003.0f, 0.0f, -6.0f),Vec3(0.12f, 0.45f, 0.12f)};
123  Sphere back_wall = {1000.0f, Vec3(0.0f, 0.0f, -1006.0f), Vec3(1.0f, 1.0f, 1.0f)};
124  Sphere floor = {1000.0f, Vec3(0.0f, -1002.6f, -6.0f), Vec3(1.0f, 1.0f, 1.0f)};
125  Sphere ceiling = {1000.0f, Vec3(0.0f, 1002.6f, -6.0f), Vec3(1.0f, 1.0f, 1.0f)};
126  Sphere light = {0.35f, Vec3(0.0f, 2.6f, -4.0f), Vec3(1.0f, 1.0f, 1.0f), true};
127  Sphere big_sphere = {1.3f, Vec3(-1.0f, -1.6f, -5.0f), Vec3(1.0f, 1.0f, 1.0f)};
128  Sphere small_sphere = {0.8f, Vec3(1.3f, -2.0f, -4.0f), Vec3(1.0f, 1.0f, 1.0f)};
129  std::vector<Sphere> spheres = {left_wall,  right_wall,  back_wall,floor,ceiling,light,big_sphere, small_sphere};
130  std::cout << "P3\n" << image_width << ' ' << image_height << "\n255\n";
131  for (int j = 0; j < image_height; j++) {
132    for (int i = 0; i < image_width; i++) {
133      Vec3 color;
134      for (int s = 0; s < samples_per_pixel; s++) {
135        float u = (((double(i) + random_float()) / image_width)-0.5f)*2;
136        float v = (((double(j) + random_float()) / image_height)-0.5f)*-2;
137        Vec3 ray_direction = normalize(Vec3(u, v, -1));
138        Ray ray = Ray(Vec3(), ray_direction);
139        color = color + Lo(ray, spheres, 10);
140      }
141      color = color * (1.0f / samples_per_pixel);
142      int ir = static_cast<int>(linear_to_srgb(color.x) * 255.0f);
143      int ig = static_cast<int>(linear_to_srgb(color.y) * 255.0f);
144      int ib = static_cast<int>(linear_to_srgb(color.z) * 255.0f);
145      std::cout << ir << ' ' << ig << ' ' << ib << '\n';
146    }
147  }
148}