Shirone

1
#include <cmath>
2
#include <iostream>
3
#include <random>
4
#include <vector>
5
constexpr float EPSILON = 1e-5f;
6
constexpr float PI = 3.14159265358979323846f;
7
struct 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
};
25
Vec3 normalize(Vec3 v) { return v * (1.0f / sqrt(v.x * v.x + v.y * v.y + v.z * v.z)); }
26
float dot(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
27
float clamp(float x, float min, float max) {
28
  return x < min ? min : (x > max ? max : x);
29
}
30
float 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
}
36
struct 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
};
42
struct 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
};
49
struct HitInfo {
50
  bool is_hit;
51
  float t;
52
  Vec3 position;
53
  Vec3 normal;
54
  Sphere const* sphere = nullptr;
55
};
56
HitInfo 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
}
72
float 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
}
77
Vec3 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
}
90
Vec3 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
}
118
int 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
}