919 lines
28 KiB
C++
919 lines
28 KiB
C++
// main.cpp
|
||
#include <algorithm>
|
||
#include <array>
|
||
#include <cassert>
|
||
#include <cmath>
|
||
#include <cstdint>
|
||
#include <fstream>
|
||
#include <iostream>
|
||
#include <limits>
|
||
#include <optional>
|
||
#include <string>
|
||
#include <utility>
|
||
#include <vector>
|
||
|
||
static constexpr double kEps = 1e-9;
|
||
|
||
struct Vec2 {
|
||
double x = 0, y = 0;
|
||
Vec2() = default;
|
||
Vec2(double x_, double y_) : x(x_), y(y_) {}
|
||
Vec2 operator+(const Vec2& o) const { return {x + o.x, y + o.y}; }
|
||
Vec2 operator-(const Vec2& o) const { return {x - o.x, y - o.y}; }
|
||
Vec2 operator*(double s) const { return {x * s, y * s}; }
|
||
};
|
||
|
||
static inline double cross2(const Vec2& a, const Vec2& b) { return a.x * b.y - a.y * b.x; }
|
||
|
||
struct Vec3 {
|
||
double x = 0, y = 0, z = 0;
|
||
Vec3() = default;
|
||
Vec3(double x_, double y_, double z_) : x(x_), y(y_), z(z_) {}
|
||
Vec3 operator+(const Vec3& o) const { return {x + o.x, y + o.y, z + o.z}; }
|
||
Vec3 operator-(const Vec3& o) const { return {x - o.x, y - o.y, z - o.z}; }
|
||
Vec3 operator*(double s) const { return {x * s, y * s, z * s}; }
|
||
Vec3 operator/(double s) const { return {x / s, y / s, z / s}; }
|
||
};
|
||
|
||
static inline double dot(const Vec3& a, const Vec3& b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
|
||
|
||
static inline Vec3 cross(const Vec3& a, const Vec3& b) {
|
||
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
|
||
}
|
||
|
||
static inline double norm2(const Vec3& v) { return dot(v, v); }
|
||
static inline double norm(const Vec3& v) { return std::sqrt(norm2(v)); }
|
||
|
||
static inline Vec3 normalize(const Vec3& v) {
|
||
const double n = norm(v);
|
||
if (n < kEps) return v;
|
||
return v / n;
|
||
}
|
||
|
||
struct Color {
|
||
double r = 0, g = 0, b = 0;
|
||
Color() = default;
|
||
Color(double r_, double g_, double b_) : r(r_), g(g_), b(b_) {}
|
||
Color operator+(const Color& o) const { return {r + o.r, g + o.g, b + o.b}; }
|
||
Color operator*(double s) const { return {r * s, g * s, b * s}; }
|
||
Color operator*(const Color& o) const { return {r * o.r, g * o.g, b * o.b}; }
|
||
};
|
||
|
||
static inline Color clamp01(const Color& c) {
|
||
auto cl = [](double v) { return std::max(0.0, std::min(1.0, v)); };
|
||
return {cl(c.r), cl(c.g), cl(c.b)};
|
||
}
|
||
|
||
struct Plane {
|
||
Vec3 n; // normalized
|
||
double d = 0; // n·X + d = 0
|
||
|
||
double sd(const Vec3& p) const { return dot(n, p) + d; }
|
||
};
|
||
|
||
static inline Plane planeFromPointNormal(const Vec3& p, const Vec3& n_unit) {
|
||
Plane pl;
|
||
pl.n = n_unit;
|
||
pl.d = -dot(pl.n, p);
|
||
return pl;
|
||
}
|
||
|
||
static inline Plane planeThrough3(const Vec3& a, const Vec3& b, const Vec3& c) {
|
||
Vec3 n = normalize(cross(b - a, c - a));
|
||
return planeFromPointNormal(a, n);
|
||
}
|
||
|
||
struct Mat3 {
|
||
double m[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
|
||
|
||
static Mat3 identity() {
|
||
Mat3 I;
|
||
I.m[0][0] = I.m[1][1] = I.m[2][2] = 1;
|
||
return I;
|
||
}
|
||
};
|
||
|
||
static inline Mat3 mul(const Mat3& A, const Mat3& B) {
|
||
Mat3 C;
|
||
for (int i = 0; i < 3; ++i) {
|
||
for (int j = 0; j < 3; ++j) {
|
||
double s = 0;
|
||
for (int k = 0; k < 3; ++k) s += A.m[i][k] * B.m[k][j];
|
||
C.m[i][j] = s;
|
||
}
|
||
}
|
||
return C;
|
||
}
|
||
|
||
static inline Vec3 mul_h(const Mat3& H, const Vec3& p) {
|
||
Vec3 r;
|
||
r.x = H.m[0][0] * p.x + H.m[0][1] * p.y + H.m[0][2] * p.z;
|
||
r.y = H.m[1][0] * p.x + H.m[1][1] * p.y + H.m[1][2] * p.z;
|
||
r.z = H.m[2][0] * p.x + H.m[2][1] * p.y + H.m[2][2] * p.z;
|
||
return r;
|
||
}
|
||
|
||
static inline bool invertMat3(const Mat3& A, Mat3& invOut) {
|
||
const double a00 = A.m[0][0], a01 = A.m[0][1], a02 = A.m[0][2];
|
||
const double a10 = A.m[1][0], a11 = A.m[1][1], a12 = A.m[1][2];
|
||
const double a20 = A.m[2][0], a21 = A.m[2][1], a22 = A.m[2][2];
|
||
|
||
const double c00 = a11 * a22 - a12 * a21;
|
||
const double c01 = -(a10 * a22 - a12 * a20);
|
||
const double c02 = a10 * a21 - a11 * a20;
|
||
|
||
const double c10 = -(a01 * a22 - a02 * a21);
|
||
const double c11 = a00 * a22 - a02 * a20;
|
||
const double c12 = -(a00 * a21 - a01 * a20);
|
||
|
||
const double c20 = a01 * a12 - a02 * a11;
|
||
const double c21 = -(a00 * a12 - a02 * a10);
|
||
const double c22 = a00 * a11 - a01 * a10;
|
||
|
||
const double det = a00 * c00 + a01 * c01 + a02 * c02;
|
||
if (std::abs(det) < 1e-18) return false;
|
||
const double invDet = 1.0 / det;
|
||
|
||
invOut.m[0][0] = c00 * invDet;
|
||
invOut.m[0][1] = c10 * invDet;
|
||
invOut.m[0][2] = c20 * invDet;
|
||
|
||
invOut.m[1][0] = c01 * invDet;
|
||
invOut.m[1][1] = c11 * invDet;
|
||
invOut.m[1][2] = c21 * invDet;
|
||
|
||
invOut.m[2][0] = c02 * invDet;
|
||
invOut.m[2][1] = c12 * invDet;
|
||
invOut.m[2][2] = c22 * invDet;
|
||
return true;
|
||
}
|
||
|
||
// Solve 8x8 linear system via Gaussian elimination (small, fixed size)
|
||
static inline bool solve8(double A[8][8], double b[8], double x[8]) {
|
||
for (int col = 0; col < 8; ++col) {
|
||
int piv = col;
|
||
double best = std::abs(A[col][col]);
|
||
for (int r = col + 1; r < 8; ++r) {
|
||
double v = std::abs(A[r][col]);
|
||
if (v > best) { best = v; piv = r; }
|
||
}
|
||
if (best < 1e-18) return false;
|
||
if (piv != col) {
|
||
for (int c = col; c < 8; ++c) std::swap(A[piv][c], A[col][c]);
|
||
std::swap(b[piv], b[col]);
|
||
}
|
||
const double diag = A[col][col];
|
||
for (int c = col; c < 8; ++c) A[col][c] /= diag;
|
||
b[col] /= diag;
|
||
|
||
for (int r = 0; r < 8; ++r) {
|
||
if (r == col) continue;
|
||
const double f = A[r][col];
|
||
if (std::abs(f) < 1e-18) continue;
|
||
for (int c = col; c < 8; ++c) A[r][c] -= f * A[col][c];
|
||
b[r] -= f * b[col];
|
||
}
|
||
}
|
||
for (int i = 0; i < 8; ++i) x[i] = b[i];
|
||
return true;
|
||
}
|
||
|
||
// Homography from 4 point correspondences: (x,y)->(u,v), with h33=1
|
||
static inline bool homography4pt(const std::array<Vec2, 4>& src, const std::array<Vec2, 4>& dst, Mat3& H, Mat3& Hinv) {
|
||
double A[8][8];
|
||
double b[8];
|
||
for (int i = 0; i < 8; ++i) {
|
||
b[i] = 0;
|
||
for (int j = 0; j < 8; ++j) A[i][j] = 0;
|
||
}
|
||
|
||
for (int i = 0; i < 4; ++i) {
|
||
const double x = src[i].x, y = src[i].y;
|
||
const double u = dst[i].x, v = dst[i].y;
|
||
|
||
// row for u
|
||
A[2 * i + 0][0] = x;
|
||
A[2 * i + 0][1] = y;
|
||
A[2 * i + 0][2] = 1;
|
||
A[2 * i + 0][6] = -u * x;
|
||
A[2 * i + 0][7] = -u * y;
|
||
b[2 * i + 0] = u;
|
||
|
||
// row for v
|
||
A[2 * i + 1][3] = x;
|
||
A[2 * i + 1][4] = y;
|
||
A[2 * i + 1][5] = 1;
|
||
A[2 * i + 1][6] = -v * x;
|
||
A[2 * i + 1][7] = -v * y;
|
||
b[2 * i + 1] = v;
|
||
}
|
||
|
||
double xsol[8];
|
||
if (!solve8(A, b, xsol)) return false;
|
||
|
||
// fill H with h33=1
|
||
H = Mat3::identity();
|
||
H.m[0][0] = xsol[0]; H.m[0][1] = xsol[1]; H.m[0][2] = xsol[2];
|
||
H.m[1][0] = xsol[3]; H.m[1][1] = xsol[4]; H.m[1][2] = xsol[5];
|
||
H.m[2][0] = xsol[6]; H.m[2][1] = xsol[7]; H.m[2][2] = 1;
|
||
|
||
if (!invertMat3(H, Hinv)) return false;
|
||
return true;
|
||
}
|
||
|
||
enum class MaterialType { Curtain, Diffuse, Metal };
|
||
|
||
struct Material {
|
||
MaterialType type = MaterialType::Diffuse;
|
||
Color albedo = {1, 1, 1}; // Diffuse base or Metal reflectance
|
||
double roughness = 0.0; // unused in ideal mirror
|
||
};
|
||
|
||
struct Triangle {
|
||
Vec3 p0, p1, p2;
|
||
Material mat;
|
||
|
||
Plane plane() const { return planeThrough3(p0, p1, p2); }
|
||
|
||
Vec3 normalUnit() const { return normalize(cross(p1 - p0, p2 - p0)); }
|
||
|
||
Vec3 pointFromUV(const Vec2& uv) const { return p0 + (p1 - p0) * uv.x + (p2 - p0) * uv.y; }
|
||
|
||
// Solve X = p0 + u*(p1-p0) + v*(p2-p0) in least-square on plane basis via Gram 2x2.
|
||
Vec2 uvFromPoint(const Vec3& X) const {
|
||
const Vec3 e1 = p1 - p0;
|
||
const Vec3 e2 = p2 - p0;
|
||
const Vec3 w = X - p0;
|
||
|
||
const double a = dot(e1, e1);
|
||
const double b = dot(e1, e2);
|
||
const double c = dot(e2, e2);
|
||
const double d = dot(e1, w);
|
||
const double e = dot(e2, w);
|
||
|
||
const double det = a * c - b * b;
|
||
if (std::abs(det) < 1e-18) return {0, 0};
|
||
const double invDet = 1.0 / det;
|
||
const double u = (d * c - b * e) * invDet;
|
||
const double v = (a * e - b * d) * invDet;
|
||
return {u, v};
|
||
}
|
||
|
||
bool containsUV(const Vec2& uv) const {
|
||
return uv.x >= -1e-9 && uv.y >= -1e-9 && (uv.x + uv.y) <= 1.0 + 1e-9;
|
||
}
|
||
};
|
||
|
||
struct Image {
|
||
int w = 0, h = 0;
|
||
std::vector<Color> color;
|
||
std::vector<double> z;
|
||
|
||
Image(int w_, int h_) : w(w_), h(h_), color(w_ * h_), z(w_ * h_, std::numeric_limits<double>::infinity()) {}
|
||
|
||
void clear(const Color& c) {
|
||
std::fill(color.begin(), color.end(), c);
|
||
std::fill(z.begin(), z.end(), std::numeric_limits<double>::infinity());
|
||
}
|
||
|
||
void setIfCloser(int x, int y, double depth, const Color& c) {
|
||
if (x < 0 || x >= w || y < 0 || y >= h) return;
|
||
const int idx = y * w + x;
|
||
if (depth < z[idx]) {
|
||
z[idx] = depth;
|
||
color[idx] = c;
|
||
}
|
||
}
|
||
|
||
void writePPM(const std::string& path) const {
|
||
std::ofstream f(path, std::ios::binary);
|
||
f << "P6\n" << w << " " << h << "\n255\n";
|
||
for (int i = 0; i < w * h; ++i) {
|
||
Color c = clamp01(color[i]);
|
||
uint8_t r = (uint8_t)std::lround(c.r * 255.0);
|
||
uint8_t g = (uint8_t)std::lround(c.g * 255.0);
|
||
uint8_t b = (uint8_t)std::lround(c.b * 255.0);
|
||
f.write((char*)&r, 1);
|
||
f.write((char*)&g, 1);
|
||
f.write((char*)&b, 1);
|
||
}
|
||
}
|
||
};
|
||
|
||
struct Frustum {
|
||
std::array<Plane, 5> planes; // 3 side + near + far(optional as large)
|
||
int count = 0;
|
||
};
|
||
|
||
static inline std::optional<Vec3> intersectRayPlane(const Vec3& E, const Vec3& dir_unit, const Plane& pl, double& tOut) {
|
||
const double denom = dot(pl.n, dir_unit);
|
||
if (std::abs(denom) < 1e-12) return std::nullopt;
|
||
const double t = -(dot(pl.n, E) + pl.d) / denom;
|
||
tOut = t;
|
||
return E + dir_unit * t;
|
||
}
|
||
|
||
static inline std::vector<Vec3> clipPolyAgainstPlane(const std::vector<Vec3>& poly, const Plane& pl) {
|
||
std::vector<Vec3> out;
|
||
if (poly.empty()) return out;
|
||
auto inside = [&](const Vec3& p) { return pl.sd(p) >= -1e-9; };
|
||
|
||
for (size_t i = 0; i < poly.size(); ++i) {
|
||
const Vec3 S = poly[i];
|
||
const Vec3 E = poly[(i + 1) % poly.size()];
|
||
const double ds = pl.sd(S);
|
||
const double de = pl.sd(E);
|
||
const bool inS = ds >= -1e-9;
|
||
const bool inE = de >= -1e-9;
|
||
|
||
if (inS && inE) {
|
||
out.push_back(E);
|
||
} else if (inS && !inE) {
|
||
const double t = ds / (ds - de);
|
||
out.push_back(S + (E - S) * t);
|
||
} else if (!inS && inE) {
|
||
const double t = ds / (ds - de);
|
||
out.push_back(S + (E - S) * t);
|
||
out.push_back(E);
|
||
}
|
||
}
|
||
return out;
|
||
}
|
||
|
||
static inline std::vector<Triangle> triangulateFan(const std::vector<Vec3>& poly, const Material& mat) {
|
||
std::vector<Triangle> tris;
|
||
if (poly.size() < 3) return tris;
|
||
for (size_t i = 1; i + 1 < poly.size(); ++i) {
|
||
Triangle t;
|
||
t.p0 = poly[0];
|
||
t.p1 = poly[i];
|
||
t.p2 = poly[i + 1];
|
||
t.mat = mat;
|
||
tris.push_back(t);
|
||
}
|
||
return tris;
|
||
}
|
||
|
||
static inline Frustum buildFrustum(const Vec3& E, const Triangle& Tv, double near_eps, std::optional<double> far_opt) {
|
||
Frustum fr;
|
||
fr.count = 0;
|
||
|
||
// Ensure Tv normal points away from eye along view direction
|
||
Vec3 nTv = Tv.normalUnit();
|
||
Vec3 cTv = (Tv.p0 + Tv.p1 + Tv.p2) / 3.0;
|
||
Vec3 viewDir = normalize(cTv - E);
|
||
if (dot(nTv, viewDir) < 0) nTv = nTv * -1.0;
|
||
|
||
// Side planes: each passes through E and an edge of Tv; inward means Tv's opposite vertex is inside.
|
||
auto addSide = [&](const Vec3& A, const Vec3& B, const Vec3& C_opp) {
|
||
Vec3 v0 = A - E;
|
||
Vec3 v1 = B - E;
|
||
Vec3 n = normalize(cross(v0, v1)); // orientation depends
|
||
Plane pl = planeFromPointNormal(E, n);
|
||
if (pl.sd(C_opp) < 0) {
|
||
pl.n = pl.n * -1.0;
|
||
pl.d = -pl.d;
|
||
}
|
||
fr.planes[fr.count++] = pl;
|
||
};
|
||
addSide(Tv.p0, Tv.p1, Tv.p2);
|
||
addSide(Tv.p1, Tv.p2, Tv.p0);
|
||
addSide(Tv.p2, Tv.p0, Tv.p1);
|
||
|
||
// Near plane: perpendicular to viewDir (use nTv) at E + nTv * near_eps
|
||
fr.planes[fr.count++] = planeFromPointNormal(E + nTv * near_eps, nTv);
|
||
|
||
// Far plane if provided: at E + nTv * far, normal facing back (-nTv) to keep inside between near/far
|
||
if (far_opt.has_value()) {
|
||
fr.planes[fr.count++] = planeFromPointNormal(E + nTv * (*far_opt), nTv * -1.0);
|
||
}
|
||
return fr;
|
||
}
|
||
|
||
static inline std::vector<Triangle> clipTriangleToFrustum(const Triangle& tri, const Frustum& fr) {
|
||
std::vector<Vec3> poly = {tri.p0, tri.p1, tri.p2};
|
||
for (int i = 0; i < fr.count; ++i) {
|
||
poly = clipPolyAgainstPlane(poly, fr.planes[i]);
|
||
if (poly.size() < 3) return {};
|
||
}
|
||
return triangulateFan(poly, tri.mat);
|
||
}
|
||
|
||
// 2D Sutherland–Hodgman: clip polygon P by triangle T (given in 2D)
|
||
static inline std::vector<Vec2> clipPolyByTri2D(const std::vector<Vec2>& poly, const std::array<Vec2, 3>& tri) {
|
||
std::vector<Vec2> out = poly;
|
||
auto clipByEdge = [&](const Vec2& A, const Vec2& B, const Vec2& insideRef) {
|
||
std::vector<Vec2> res;
|
||
if (out.empty()) return res;
|
||
|
||
Vec2 e = B - A;
|
||
double sign = cross2(e, insideRef - A);
|
||
if (std::abs(sign) < 1e-18) sign = 1.0;
|
||
|
||
auto inside = [&](const Vec2& p) { return cross2(e, p - A) * sign >= -1e-9; };
|
||
|
||
for (size_t i = 0; i < out.size(); ++i) {
|
||
Vec2 S = out[i];
|
||
Vec2 E = out[(i + 1) % out.size()];
|
||
bool inS = inside(S);
|
||
bool inE = inside(E);
|
||
|
||
auto intersect = [&]() -> Vec2 {
|
||
Vec2 d = E - S;
|
||
double num = cross2(e, A - S);
|
||
double den = cross2(e, d);
|
||
if (std::abs(den) < 1e-18) return S;
|
||
double t = num / den;
|
||
return S + d * t;
|
||
};
|
||
|
||
if (inS && inE) {
|
||
res.push_back(E);
|
||
} else if (inS && !inE) {
|
||
res.push_back(intersect());
|
||
} else if (!inS && inE) {
|
||
res.push_back(intersect());
|
||
res.push_back(E);
|
||
}
|
||
}
|
||
return res;
|
||
};
|
||
|
||
// clip in order by each edge, using opposite vertex as inside reference
|
||
out = clipByEdge(tri[0], tri[1], tri[2]);
|
||
out = clipByEdge(tri[1], tri[2], tri[0]);
|
||
out = clipByEdge(tri[2], tri[0], tri[1]);
|
||
return out;
|
||
}
|
||
|
||
static inline double triArea2(const Vec2& a, const Vec2& b, const Vec2& c) {
|
||
return std::abs(cross2(b - a, c - a)) * 0.5;
|
||
}
|
||
|
||
static inline bool barycentric2D(const Vec2& p, const Vec2& a, const Vec2& b, const Vec2& c, double& w0, double& w1, double& w2) {
|
||
const Vec2 v0 = b - a;
|
||
const Vec2 v1 = c - a;
|
||
const Vec2 v2 = p - a;
|
||
const double den = cross2(v0, v1);
|
||
if (std::abs(den) < 1e-18) return false;
|
||
w1 = cross2(v2, v1) / den;
|
||
w2 = cross2(v0, v2) / den;
|
||
w0 = 1.0 - w1 - w2;
|
||
return true;
|
||
}
|
||
|
||
struct Screen {
|
||
// Root "curtain" is a rectangle in 3D: O + u*U + v*V, u,v in [0,1]
|
||
Vec3 O, U, V;
|
||
int W = 0, H = 0;
|
||
|
||
Vec3 pointFromUV01(const Vec2& uv01) const { return O + U * uv01.x + V * uv01.y; }
|
||
|
||
Vec2 uv01FromPixelCenter(int x, int y) const {
|
||
// y down
|
||
double u = (x + 0.5) / (double)W;
|
||
double v = (y + 0.5) / (double)H;
|
||
return {u, v};
|
||
}
|
||
};
|
||
|
||
struct DepthOverride {
|
||
Plane writePlane; // first mirror plane in root view
|
||
bool valid = false;
|
||
};
|
||
|
||
struct RenderConfig {
|
||
double near_eps = 1e-4;
|
||
std::optional<double> far = std::nullopt;
|
||
int L_max = 4;
|
||
double A_min_pixels = 2.0;
|
||
};
|
||
|
||
struct RenderContext {
|
||
const std::vector<Triangle>* scene = nullptr;
|
||
Image* img = nullptr;
|
||
Screen screen;
|
||
Vec3 rootEye;
|
||
RenderConfig cfg;
|
||
};
|
||
|
||
static inline Vec2 applyH(const Mat3& H, const Vec2& p) {
|
||
Vec3 hp = mul_h(H, Vec3{p.x, p.y, 1.0});
|
||
if (std::abs(hp.z) < 1e-18) return {hp.x, hp.y};
|
||
return {hp.x / hp.z, hp.y / hp.z};
|
||
}
|
||
|
||
static inline std::optional<std::pair<double, Vec3>> intersectRayTrianglePlane(const Vec3& E, const Vec3& dir_unit, const Triangle& tri) {
|
||
Plane pl = tri.plane();
|
||
double t = 0;
|
||
auto Xopt = intersectRayPlane(E, dir_unit, pl, t);
|
||
if (!Xopt.has_value()) return std::nullopt;
|
||
return std::make_pair(t, *Xopt);
|
||
}
|
||
|
||
static inline Vec3 reflectPointAboutPlane(const Vec3& P, const Plane& pl) {
|
||
// pl is normalized
|
||
const double s = pl.sd(P);
|
||
return P - pl.n * (2.0 * s);
|
||
}
|
||
|
||
static inline double estimateTrianglePixelArea(const Mat3& H_toRoot, const Triangle& Tv_local, const Screen& screen) {
|
||
// Tv_local canonical coords: (0,0),(1,0),(0,1)
|
||
Vec2 r0 = applyH(H_toRoot, Vec2{0, 0});
|
||
Vec2 r1 = applyH(H_toRoot, Vec2{1, 0});
|
||
Vec2 r2 = applyH(H_toRoot, Vec2{0, 1});
|
||
Vec2 p0 = {r0.x * screen.W, r0.y * screen.H};
|
||
Vec2 p1 = {r1.x * screen.W, r1.y * screen.H};
|
||
Vec2 p2 = {r2.x * screen.W, r2.y * screen.H};
|
||
return triArea2(p0, p1, p2);
|
||
}
|
||
|
||
static inline Color shadeDiffuseSimple(const Triangle& tri, const Vec3& X, const Vec3& dir_unit) {
|
||
Vec3 n = tri.normalUnit();
|
||
// Simple directional light + ambient
|
||
Vec3 L = normalize(Vec3{-0.3, -1.0, -0.2}); // light coming from above
|
||
double ndotl = std::max(0.0, dot(n, L * -1.0));
|
||
double ambient = 0.08;
|
||
double intensity = ambient + 0.92 * ndotl;
|
||
|
||
// optional: face-forward
|
||
if (dot(n, dir_unit) > 0) n = n * -1.0;
|
||
|
||
return tri.mat.albedo * intensity;
|
||
}
|
||
|
||
static void renderTriangleWithTriangle(
|
||
RenderContext& ctx,
|
||
const Vec3& E_cur,
|
||
const Triangle& Tv_local,
|
||
const Mat3& H_accum, const Mat3& H_accum_inv,
|
||
const Color& throughput,
|
||
const DepthOverride& depthOverride,
|
||
int depthLevel);
|
||
|
||
// Diffuse branch: clip piece already in frustum, project to Tv plane, rasterize in root pixels
|
||
static void renderDiffusePiece(
|
||
RenderContext& ctx,
|
||
const Vec3& E_cur,
|
||
const Triangle& Tv_local,
|
||
const Triangle& triPiece,
|
||
const Mat3& H_accum, const Mat3& H_accum_inv,
|
||
const Color& throughput,
|
||
const DepthOverride& depthOverride) {
|
||
|
||
// Project triPiece vertices to Tv plane via rays from E_cur
|
||
Plane tvPl = Tv_local.plane();
|
||
std::array<Vec2, 3> projUV;
|
||
std::array<Vec3, 3> projP;
|
||
|
||
std::array<Vec3, 3> V = {triPiece.p0, triPiece.p1, triPiece.p2};
|
||
for (int i = 0; i < 3; ++i) {
|
||
Vec3 dir = normalize(V[i] - E_cur);
|
||
double t = 0;
|
||
auto Xopt = intersectRayPlane(E_cur, dir, tvPl, t);
|
||
if (!Xopt.has_value() || t <= ctx.cfg.near_eps) return;
|
||
projP[i] = *Xopt;
|
||
projUV[i] = Tv_local.uvFromPoint(projP[i]); // in Tv_local uv-space (not necessarily inside Tv triangle, but should)
|
||
}
|
||
|
||
// Map projected triangle into root uv01 via H_accum
|
||
std::array<Vec2, 3> rootUV01 = {
|
||
applyH(H_accum, projUV[0]),
|
||
applyH(H_accum, projUV[1]),
|
||
applyH(H_accum, projUV[2]),
|
||
};
|
||
|
||
// Pixel bbox in root image
|
||
auto toPix = [&](const Vec2& uv01) -> Vec2 {
|
||
return Vec2{uv01.x * ctx.screen.W, uv01.y * ctx.screen.H};
|
||
};
|
||
Vec2 p0 = toPix(rootUV01[0]);
|
||
Vec2 p1 = toPix(rootUV01[1]);
|
||
Vec2 p2 = toPix(rootUV01[2]);
|
||
|
||
int xmin = (int)std::floor(std::min({p0.x, p1.x, p2.x}));
|
||
int xmax = (int)std::ceil (std::max({p0.x, p1.x, p2.x}));
|
||
int ymin = (int)std::floor(std::min({p0.y, p1.y, p2.y}));
|
||
int ymax = (int)std::ceil (std::max({p0.y, p1.y, p2.y}));
|
||
xmin = std::max(xmin, 0); ymin = std::max(ymin, 0);
|
||
xmax = std::min(xmax, ctx.img->w - 1); ymax = std::min(ymax, ctx.img->h - 1);
|
||
|
||
// For each pixel center: root uv01 -> tv uv via inverse H_accum -> inside test -> ray intersection
|
||
for (int y = ymin; y <= ymax; ++y) {
|
||
for (int x = xmin; x <= xmax; ++x) {
|
||
Vec2 uv01 = ctx.screen.uv01FromPixelCenter(x, y);
|
||
Vec2 uv_tv = applyH(H_accum_inv, uv01);
|
||
|
||
double w0, w1, w2;
|
||
if (!barycentric2D(uv_tv, projUV[0], projUV[1], projUV[2], w0, w1, w2)) continue;
|
||
if (w0 < -1e-8 || w1 < -1e-8 || w2 < -1e-8) continue;
|
||
|
||
// Build point on Tv plane for this uv_tv
|
||
Vec3 P_on_tv = Tv_local.pointFromUV(uv_tv);
|
||
Vec3 dir = normalize(P_on_tv - E_cur);
|
||
|
||
auto hit = intersectRayTrianglePlane(E_cur, dir, triPiece);
|
||
if (!hit.has_value()) continue;
|
||
double t_virtual = hit->first;
|
||
Vec3 X = hit->second;
|
||
if (t_virtual <= ctx.cfg.near_eps) continue;
|
||
|
||
// Ensure inside triPiece in its own uv
|
||
Vec2 uv_piece = triPiece.uvFromPoint(X);
|
||
if (!triPiece.containsUV(uv_piece)) continue;
|
||
|
||
// Depth for z-buffer: either real hit depth, or overridden mirror depth in root view
|
||
double depth = t_virtual;
|
||
if (depthOverride.valid) {
|
||
Vec3 P_root = ctx.screen.pointFromUV01(uv01);
|
||
Vec3 dir_root = normalize(P_root - ctx.rootEye);
|
||
double t_override = 0;
|
||
auto Xo = intersectRayPlane(ctx.rootEye, dir_root, depthOverride.writePlane, t_override);
|
||
if (!Xo.has_value() || t_override <= ctx.cfg.near_eps) continue;
|
||
depth = t_override;
|
||
}
|
||
|
||
Color base = shadeDiffuseSimple(triPiece, X, dir);
|
||
Color out = base * throughput;
|
||
ctx.img->setIfCloser(x, y, depth, out);
|
||
}
|
||
}
|
||
}
|
||
|
||
static void renderMetalTriangle(
|
||
RenderContext& ctx,
|
||
const Vec3& E_cur,
|
||
const Triangle& Tv_local, // parent curtain viewport
|
||
const Triangle& triMetal, // unclipped original metal tri OR clipped piece; both ok if plane same
|
||
const Mat3& H_accum, const Mat3& H_accum_inv,
|
||
const Color& throughput,
|
||
const DepthOverride& depthOverride,
|
||
int depthLevel) {
|
||
|
||
if (depthLevel >= ctx.cfg.L_max) return;
|
||
|
||
Plane metalPl = triMetal.plane(); // normalized
|
||
// Footprint of parent frustum on metal plane: intersect 3 boundary rays (E_cur -> Tv_local vertices)
|
||
std::array<Vec3, 3> Q;
|
||
std::array<Vec3, 3> Sv = {Tv_local.p0, Tv_local.p1, Tv_local.p2};
|
||
for (int i = 0; i < 3; ++i) {
|
||
Vec3 dir = normalize(Sv[i] - E_cur);
|
||
double t = 0;
|
||
auto Xopt = intersectRayPlane(E_cur, dir, metalPl, t);
|
||
if (!Xopt.has_value() || t <= ctx.cfg.near_eps) return;
|
||
Q[i] = *Xopt;
|
||
}
|
||
|
||
// Clip footprint triangle Q by the metal triangle in-plane (2D)
|
||
// Use metal triangle local 2D coords in its own basis
|
||
auto to2D_inMetal = [&](const Vec3& X) { return triMetal.uvFromPoint(X); };
|
||
std::array<Vec2, 3> metalTri2D = {Vec2{0, 0}, Vec2{1, 0}, Vec2{0, 1}};
|
||
std::vector<Vec2> poly2D = {to2D_inMetal(Q[0]), to2D_inMetal(Q[1]), to2D_inMetal(Q[2])};
|
||
poly2D = clipPolyByTri2D(poly2D, metalTri2D);
|
||
if (poly2D.size() < 3) return;
|
||
|
||
// Triangulate child viewport polygon on metal plane
|
||
std::vector<Triangle> childViewports;
|
||
for (size_t i = 1; i + 1 < poly2D.size(); ++i) {
|
||
Triangle tv;
|
||
tv.p0 = triMetal.pointFromUV(poly2D[0]);
|
||
tv.p1 = triMetal.pointFromUV(poly2D[i]);
|
||
tv.p2 = triMetal.pointFromUV(poly2D[i + 1]);
|
||
tv.mat.type = MaterialType::Curtain; // IMPORTANT: child is a viewport, not a drawable
|
||
tv.mat.albedo = {0, 0, 0};
|
||
childViewports.push_back(tv);
|
||
}
|
||
|
||
// Mirror camera
|
||
Vec3 E_mirror = reflectPointAboutPlane(E_cur, metalPl);
|
||
|
||
// Decide depth override: lock to the first mirror plane in root view
|
||
DepthOverride nextOverride = depthOverride;
|
||
if (!nextOverride.valid) {
|
||
nextOverride.valid = true;
|
||
nextOverride.writePlane = metalPl;
|
||
}
|
||
|
||
// Recurse for each child viewport
|
||
for (const Triangle& Tv_child : childViewports) {
|
||
// Stop condition by projected pixel area in root
|
||
// Need H_step: child uv -> parent uv, induced by central projection from E_cur onto parent plane
|
||
Plane parentPl = Tv_local.plane();
|
||
|
||
auto projToParent = [&](const Vec3& X_on_child) -> Vec3 {
|
||
Vec3 dir = normalize(X_on_child - E_cur);
|
||
double t = 0;
|
||
auto Xp = intersectRayPlane(E_cur, dir, parentPl, t);
|
||
if (!Xp.has_value()) return Tv_local.p0;
|
||
return *Xp;
|
||
};
|
||
|
||
// 4 correspondences: 3 vertices + centroid
|
||
std::array<Vec2, 4> src_uv_child = {
|
||
Vec2{0, 0},
|
||
Vec2{1, 0},
|
||
Vec2{0, 1},
|
||
Vec2{1.0 / 3.0, 1.0 / 3.0},
|
||
};
|
||
Vec3 C3 = (Tv_child.p0 + Tv_child.p1 + Tv_child.p2) / 3.0;
|
||
std::array<Vec3, 4> child3 = {Tv_child.p0, Tv_child.p1, Tv_child.p2, C3};
|
||
std::array<Vec2, 4> dst_uv_parent;
|
||
for (int i = 0; i < 4; ++i) {
|
||
Vec3 Xp = projToParent(child3[i]);
|
||
dst_uv_parent[i] = Tv_local.uvFromPoint(Xp);
|
||
}
|
||
|
||
Mat3 H_step, H_step_inv;
|
||
if (!homography4pt(src_uv_child, dst_uv_parent, H_step, H_step_inv)) continue;
|
||
|
||
Mat3 H_child_toRoot = mul(H_accum, H_step);
|
||
Mat3 H_child_toRoot_inv;
|
||
if (!invertMat3(H_child_toRoot, H_child_toRoot_inv)) continue;
|
||
|
||
double areaPx = estimateTrianglePixelArea(H_child_toRoot, Tv_child, ctx.screen);
|
||
if (areaPx < ctx.cfg.A_min_pixels) continue;
|
||
|
||
Color nextThroughput = throughput * triMetal.mat.albedo; // metal reflectance tint
|
||
renderTriangleWithTriangle(ctx, E_mirror, Tv_child, H_child_toRoot, H_child_toRoot_inv, nextThroughput, nextOverride, depthLevel + 1);
|
||
}
|
||
}
|
||
|
||
static void renderTriangleWithTriangle(
|
||
RenderContext& ctx,
|
||
const Vec3& E_cur,
|
||
const Triangle& Tv_local,
|
||
const Mat3& H_accum, const Mat3& H_accum_inv,
|
||
const Color& throughput,
|
||
const DepthOverride& depthOverride,
|
||
int depthLevel) {
|
||
|
||
if (Tv_local.mat.type != MaterialType::Curtain) {
|
||
// In this design, only Curtain triangles are viewports.
|
||
// (Metal recursion creates Tv_child as Curtain.)
|
||
return;
|
||
}
|
||
|
||
// Build frustum from E_cur and current viewport triangle
|
||
Frustum fr = buildFrustum(E_cur, Tv_local, ctx.cfg.near_eps, ctx.cfg.far);
|
||
|
||
// Collect candidates: clip each scene triangle against frustum
|
||
for (const Triangle& tri : *ctx.scene) {
|
||
if (tri.mat.type == MaterialType::Curtain) continue; // do not render curtain geometry
|
||
std::vector<Triangle> pieces = clipTriangleToFrustum(tri, fr);
|
||
for (const Triangle& piece : pieces) {
|
||
if (piece.mat.type == MaterialType::Diffuse) {
|
||
renderDiffusePiece(ctx, E_cur, Tv_local, piece, H_accum, H_accum_inv, throughput, depthOverride);
|
||
} else if (piece.mat.type == MaterialType::Metal) {
|
||
// Metal is recursive, no direct color; use parent frustum footprint method
|
||
renderMetalTriangle(ctx, E_cur, Tv_local, piece, H_accum, H_accum_inv, throughput, depthOverride, depthLevel);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// Build Cornell Box (5 faces) + 2 metal boxes
|
||
static std::vector<Triangle> buildScene() {
|
||
std::vector<Triangle> s;
|
||
|
||
auto addQuad = [&](Vec3 a, Vec3 b, Vec3 c, Vec3 d, const Material& m) {
|
||
// (a,b,c) and (a,c,d)
|
||
Triangle t1{a, b, c, m};
|
||
Triangle t2{a, c, d, m};
|
||
s.push_back(t1);
|
||
s.push_back(t2);
|
||
};
|
||
|
||
// Cornell dimensions
|
||
// Coordinate system: x right, y up, z forward (into the box)
|
||
const double W = 2.0;
|
||
const double H = 2.0;
|
||
const double D = 2.5;
|
||
|
||
Vec3 A{-W/2, 0, 0};
|
||
Vec3 B{ W/2, 0, 0};
|
||
Vec3 C{ W/2, H, 0};
|
||
Vec3 Dd{-W/2, H, 0};
|
||
|
||
Vec3 A2{-W/2, 0, D};
|
||
Vec3 B2{ W/2, 0, D};
|
||
Vec3 C2{ W/2, H, D};
|
||
Vec3 D2{-W/2, H, D};
|
||
|
||
Material white{MaterialType::Diffuse, {0.75, 0.75, 0.75}};
|
||
Material red{MaterialType::Diffuse, {0.75, 0.15, 0.15}};
|
||
Material green{MaterialType::Diffuse, {0.15, 0.75, 0.15}};
|
||
|
||
// Floor: A..B..B2..A2
|
||
addQuad(A, B, B2, A2, white);
|
||
// Ceiling: Dd..C..C2..D2
|
||
addQuad(Dd, C, C2, D2, white);
|
||
// Back wall: A2..B2..C2..D2
|
||
addQuad(A2, B2, C2, D2, white);
|
||
// Left wall: A..A2..D2..Dd
|
||
addQuad(A, A2, D2, Dd, red);
|
||
// Right wall: B..C..C2..B2
|
||
addQuad(B, C, C2, B2, green);
|
||
|
||
auto addBox = [&](Vec3 minP, Vec3 maxP, const Material& m) {
|
||
Vec3 p000{minP.x, minP.y, minP.z};
|
||
Vec3 p100{maxP.x, minP.y, minP.z};
|
||
Vec3 p110{maxP.x, maxP.y, minP.z};
|
||
Vec3 p010{minP.x, maxP.y, minP.z};
|
||
|
||
Vec3 p001{minP.x, minP.y, maxP.z};
|
||
Vec3 p101{maxP.x, minP.y, maxP.z};
|
||
Vec3 p111{maxP.x, maxP.y, maxP.z};
|
||
Vec3 p011{minP.x, maxP.y, maxP.z};
|
||
|
||
// 6 faces
|
||
addQuad(p000, p100, p110, p010, m); // front (z=min)
|
||
addQuad(p001, p011, p111, p101, m); // back (z=max)
|
||
addQuad(p000, p010, p011, p001, m); // left
|
||
addQuad(p100, p101, p111, p110, m); // right
|
||
addQuad(p010, p110, p111, p011, m); // top
|
||
addQuad(p000, p001, p101, p100, m); // bottom
|
||
};
|
||
|
||
Material metalBlue{MaterialType::Metal, {0.35, 0.45, 0.95}};
|
||
Material metalYellow{MaterialType::Metal, {0.95, 0.85, 0.25}};
|
||
|
||
// Two boxes inside (raise above floor a bit)
|
||
addBox(Vec3{-0.6, 0.0, 0.7}, Vec3{-0.1, 0.9, 1.2}, metalBlue);
|
||
addBox(Vec3{ 0.1, 0.0, 1.2}, Vec3{ 0.7, 1.4, 1.9}, metalYellow);
|
||
|
||
return s;
|
||
}
|
||
|
||
struct Camera {
|
||
Vec3 E;
|
||
Triangle Tv0, Tv1; // Curtain triangles covering screen
|
||
Screen screen;
|
||
|
||
static Camera makeDefault(int W, int H) {
|
||
Camera cam;
|
||
cam.E = Vec3{0, 1.0, -3.2};
|
||
|
||
// Screen rectangle in world: put in front of eye
|
||
// Define 4 corners of screen quad in world (z=-1.0 plane, spanning [-1,1]x[0,2])
|
||
Vec3 a{-1.0, 0.0, -1.0};
|
||
Vec3 b{ 1.0, 0.0, -1.0};
|
||
Vec3 c{ 1.0, 2.0, -1.0};
|
||
Vec3 d{-1.0, 2.0, -1.0};
|
||
|
||
Material curtain{MaterialType::Curtain, {0, 0, 0}};
|
||
cam.Tv0 = Triangle{a, b, c, curtain};
|
||
cam.Tv1 = Triangle{a, c, d, curtain};
|
||
|
||
cam.screen.O = a;
|
||
cam.screen.U = (b - a);
|
||
cam.screen.V = (d - a);
|
||
cam.screen.W = W;
|
||
cam.screen.H = H;
|
||
|
||
return cam;
|
||
}
|
||
|
||
void render(const std::vector<Triangle>& scene, Image& img, const RenderConfig& cfg) const {
|
||
RenderContext ctx;
|
||
ctx.scene = &scene;
|
||
ctx.img = &img;
|
||
ctx.screen = screen;
|
||
ctx.rootEye = E;
|
||
ctx.cfg = cfg;
|
||
|
||
Mat3 H0 = Mat3::identity();
|
||
Mat3 H0inv = Mat3::identity();
|
||
DepthOverride none;
|
||
none.valid = false;
|
||
|
||
Color throughput{1, 1, 1};
|
||
|
||
renderTriangleWithTriangle(ctx, E, Tv0, H0, H0inv, throughput, none, 0);
|
||
renderTriangleWithTriangle(ctx, E, Tv1, H0, H0inv, throughput, none, 0);
|
||
}
|
||
};
|
||
|
||
int main(int argc, char** argv) {
|
||
std::string out = "out.ppm";
|
||
if (argc >= 2) out = argv[1];
|
||
|
||
int W = 640, H = 480;
|
||
Image img(W, H);
|
||
img.clear(Color{0.0, 0.0, 0.0});
|
||
|
||
std::vector<Triangle> scene = buildScene();
|
||
Camera cam = Camera::makeDefault(W, H);
|
||
|
||
RenderConfig cfg;
|
||
cfg.near_eps = 1e-4;
|
||
cfg.far = std::nullopt; // or 10.0
|
||
cfg.L_max = 4;
|
||
cfg.A_min_pixels = 2.0;
|
||
|
||
cam.render(scene, img, cfg);
|
||
img.writePPM(out);
|
||
|
||
std::cerr << "Wrote " << out << "\n";
|
||
return 0;
|
||
}
|