This repository has been archived by the owner on Jun 16, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
blackhole-high-res.cpp
213 lines (175 loc) · 6.04 KB
/
blackhole-high-res.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
#include <chrono>
#include <cmath>
#include <cstring>
#include <ctime>
#include <iostream>
#include <thread>
#include "utils.h"
// FIXME: per sub-region should be doable
// use thread per fragment
#define USE_THREAD 1
constexpr float width = 130;
constexpr float height = 80;
// black hole configuration
constexpr float RS = .125; // singularity radius
constexpr float PS_RAD = 1.5 * RS; // photon sphere radius
constexpr float ACC_RAD = 3. * RS; // accretion disc radius
float gtime = .1;
// ray marching configuration
constexpr float MAX_DEPTH = 100.;
constexpr int MAX_STEPS = 300;
constexpr float DP = 0.05;
// core
inline float sdAccretionDisc(vec3 p) {
p = applyTransf(rotateX(-std::sin((PI / 3.) * 3.2)), p);
p = applyTransf(rotateZ(-std::sin((PI / 3.) * 3.2)), p);
const float p1 = sdRoundedCylinder(p, ACC_RAD, ACC_RAD / 6, .001);
const float p2 = sdSphere(p, ACC_RAD);
return sdSmoothSubtraction(p2, p1, .5);
}
// rough approximation of how light bends
inline float interpSpaceDistortion(float sz_rad, float dist_singularity) {
const float distortion_factor = 2.79;
const float dist_ratio = sz_rad / dist_singularity;
return pow(dist_ratio, distortion_factor);
}
vec3 bendLightDirection(vec3 bl_pos, vec3 ray_pos, vec3 ray_dir) {
// points at the current direction
const vec3 u = normalize(ray_dir);
// points at the singularity i.e. center of the black hole
// (u, v) angle is the maximum deviation angle starting from the current
// unchanged light direction
const vec3 v = normalize(sub(bl_pos, ray_pos));
// in this setup we use v directly
// let's bend u in such a way that it follows v
// also let's consider how close it is to bend it properly
float dist_how_close = length(sub(bl_pos, ray_pos));
float lerp_val = interpSpaceDistortion(RS, dist_how_close);
return normalize(lerp3(u, v, lerp_val));
}
inline float normalizedNoiseTexture(float x, float y) {
// noise-ish texture... a normalized chess board
const int sx = (int)(x * 10.f);
const int sy = (int)(y * 30.f);
if ((sx + sy) % 2 == 0)
return 0.3;
return 1.;
}
// coord is uv normalized !
inline float gridTexture(float coord_x, float coord_y) {
constexpr float scale = 3.;
const float dp = gtime / 5.;
const float gx = fixed_fmod((coord_x + dp) * scale, .5f);
const float gy = fixed_fmod(coord_y * scale, .5f);
// dot
// return min(1.0, step (gx + gy, .1));
// row col
return std::min(1.0f, step(gx, 0.1f) + step(gy, 0.1f));
}
void fragmentHandler(float x, float y, t_screen *screen) {
// we define a direction for each pixel
const vec3 camera{0., 0., 2.}; // the camera must be above
const vec3 cam_dir = normalize({x, y, -1.});
// ray marching
float dtravel = 0.05;
const vec3 bl_pos(0., 0., 0.);
vec3 ray_pos = camera;
vec3 ray_dir = cam_dir;
for (int i = 0; i < MAX_STEPS; i++) {
ray_dir = bendLightDirection(bl_pos, ray_pos, ray_dir);
ray_pos = add(ray_pos, scaleReal(ray_dir, DP));
float d = sdAccretionDisc(ray_pos);
dtravel += d;
if (d < EPSILON)
break;
if (d > MAX_DEPTH)
break;
}
// base color
float diffuse = 0.;
// just after the photon sphere
if (length(vec3(x, y, 0.)) < PS_RAD)
diffuse = 0.;
// construct the accretion disk and the photon sphere
if (dtravel < MAX_DEPTH) {
const vec3 delta = sub(ray_pos, bl_pos);
// make it glow more as it approaches the center
// create glow and diminish it with distance
float glow = .3 / std::pow(length(delta), 2.);
glow = clamp(glow, 0., 1.); // remove artifacts
// rotation effect
// the noise-ish texture is normalized and have deterministic values
// rotate the uv coordinate
constexpr float rot_vel = 3.;
const float s = std::sin(gtime * rot_vel), c = std::cos(gtime * rot_vel);
const float tuv_x = ray_pos.y;
const float tuv_y = ray_pos.z;
const float rot_x = c * tuv_x - s * tuv_y;
const float rot_y = s * tuv_x + c * tuv_y;
float tex_color = normalizedNoiseTexture(rot_x, rot_y);
// total color
diffuse += glow * tex_color;
}
// background space deformation
if (length(vec3(x, y, 0.)) >= PS_RAD) {
float back_col = gridTexture(ray_pos.x / 5., ray_pos.y / 5.);
diffuse += clamp(back_col, 0., .2);
}
// end ray marching
const char pixel =
screen->computeColorGivenDiffuseLight(diffuse, COLOR_STRONG);
// texture coords ---> screen coords
// -0.5 0.5 ---> -width/2 width/2
// -0.5 0.5 ---> -height/2 height/2
const int screen_x = (int)((x + 0.5) * width),
screen_y = (int)((y + 0.5) * height);
screen->put(screen_y, screen_x, pixel);
}
void computeScreenBufferConcurrent(t_screen *screen) {
// motivation : the bigger the screensize, the more the steps
// we can also assign an arbitrary step
const float dp = 1.f / std::max(height, width);
// normalized coordinates centered at (0., 0.)
constexpr float sy = -0.5, ey = 0.5;
constexpr float sx = -0.5, ex = 0.5;
const float ratio = width / height;
// iterate through the texture coordinate
std::vector<std::thread> frag_workers;
for (float y = sy; y < ey; y += dp) {
for (float x = sx; x < ex; x += dp) {
// FIXME: spread work accross max threads (thread spawn has significant
// overhead)
#if USE_THREAD == 1
std::thread worker(fragmentHandler, x, y * ratio, screen);
frag_workers.push_back(std::move(worker));
#else
fragmentHandler(x * ratio, y * ratio, screen);
#endif
}
}
#if USE_THREAD == 1
// for each subtask, sync with the main thread i.e. wait for all threads to
// terminate
for (auto &worker : frag_workers)
worker.join();
#endif
}
int main() {
t_screen screen(width, height);
screen.clear();
#if USE_THREAD == 1
std::cout << "/!\\ Preparing ~" << pow(std::max(height, width), 2.f)
<< " parallel tasks\n";
#endif
while (true) {
screen.saveCursor();
computeScreenBufferConcurrent(&screen);
screen.show();
screen.restoreCursor();
gtime += .2;
// attempt 60fps
using namespace std::literals::chrono_literals;
std::this_thread::sleep_for(16ms);
}
return 0;
}