c++实现光线跟踪算法

光线追踪算法是一种用于渲染3D图像的技术，它通过模拟光线在场景中的传播和反射来计算像素的颜色。以下是一个基于C++的光线追踪算法的简单实现：

#include <iostream>
#include <vector>
#include <cmath>
#include <limits>

// 定义向量类
class Vec3 {
public:
    double x, y, z;
    Vec3(double x_=0, double y_=0, double z_=0): x(x_), y(y_), z(z_) {}
    Vec3 operator+(const Vec3 &b) const { return Vec3(x+b.x, y+b.y, z+b.z); }
    Vec3 operator-(const Vec3 &b) const { return Vec3(x-b.x, y-b.y, z-b.z); }
    Vec3 operator*(double b) const { return Vec3(x*b, y*b, z*b); }
    Vec3 mult(const Vec3 &b) const { return Vec3(x*b.x, y*b.y, z*b.z); }
    Vec3& norm() { return *this = *this * (1 / sqrt(x*x+y*y+z*z)); }
    double dot(const Vec3 &b) const { return x*b.x+y*b.y+z*b.z; }
    Vec3 operator%(Vec3&b){return Vec3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};

// 定义光线类
class Ray {
public:
    Vec3 o, d;
    Ray(Vec3 o_, Vec3 d_): o(o_), d(d_) {}
};

// 定义材质类
class Material {
public:
    Vec3 color, emission;
    double refl, refr;
    Material(Vec3 color_, Vec3 emission_, double refl_=0, double refr_=0): color(color_), emission(emission_), refl(refl_), refr(refr_) {}
};

// 定义球体类
class Sphere {
public:
    double rad;
    Vec3 p, e, c;
    Material mat;
    Sphere(double rad_, Vec3 p_, Vec3 e_, Vec3 c_, Material mat_): rad(rad_), p(p_), e(e_), c(c_), mat(mat_) {}
    double intersect(const Ray &r) const {
        Vec3 op = p - r.o;
        double t, eps=1e-4, b = op.dot(r.d), det = b*b - op.dot(op) + rad*rad;
        if (det<0) return 0; else det=sqrt(det);
        return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
    }
};

// 定义三角形类
class Triangle {
public:
    Vec3 a, b, c;
    Material mat;
    Triangle(Vec3 a_, Vec3 b_, Vec3 c_, Material mat_): a(a_), b(b_), c(c_), mat(mat_) {}
    double intersect(const Ray &r) const {
        Vec3 ab = b - a, ac = c - a, pvec = r.d % ac;
        double det = ab.dot(pvec);
        if (fabs(det) < 1e-8) return 0;
        double invDet = 1 / det;
        Vec3 tvec = r.o - a;
        double u = (tvec.dot(pvec)) * invDet;
        if (u < 0 || u > 1) return 0;
        Vec3 qvec = tvec % ab;
        double v = (r.d.dot(qvec)) * invDet;
        if (v < 0 || u + v > 1) return 0;
        return ac.dot(qvec) * invDet;
    }
};

// 定义平面类
class Plane {
public:
    Vec3 n, p, e, c;
    Material mat;
    Plane(Vec3 n_, Vec3 p_, Vec3 e_, Vec3 c_, Material mat_): n(n_), p(p_), e(e_), c(c_), mat(mat_) {}
    double intersect(const Ray &r) const {
        double d = n.dot(p);
        double t = (d - n.dot(r.o)) / n.dot(r.d);
        if (t < 0) return 0;
        return t;
    }
};

// 定义场景类
class Scene {
public:
    std::vector<Sphere> spheres;
    std::vector<Triangle> triangles;
    std::vector<Plane> planes;
    Vec3 bg;
    Scene(Vec3 bg_): bg(bg_) {}
    void add(Sphere s) { spheres.push_back(s); }
    void add(Triangle t) { triangles.push_back(t); }
    void add(Plane p) { planes.push_back(p); }
    Vec3 intersect(const Ray &r, int depth) const {
        double t = std::numeric_limits<double>::infinity();
        const Sphere *s = NULL;
        const Triangle *tr = NULL;
        const Plane *pl = NULL;
        for (unsigned i = 0; i < spheres.size(); ++i) {
            double d = spheres[i].intersect(r);
            if (d && d < t) { t = d; s = &spheres[i]; }
        }
        for (unsigned i = 0; i < triangles.size(); ++i) {
            double d = triangles[i].intersect(r);
            if (d && d < t) { t = d; tr = &triangles[i]; }
        }
        for (unsigned i = 0; i < planes.size(); ++i) {
            double d = planes[i].intersect(r);
            if (d && d < t) { t = d; pl = &planes[i]; }
        }
        if (!s && !tr && !pl) return bg;
        Vec3 x = r.o + r.d * t, n;
        if (s) { n = (x - s->p).norm(); }
        else if (tr) { n = ((tr->b - tr->a) % (tr->c - tr->a)).norm(); }
        else { n = pl->n; }
        Vec3 f = s ? s->mat.color : (tr ? tr->mat.color : pl->mat.color);
        Vec3 e = s ? s->mat.emission : (tr ? tr->mat.emission : pl->mat.emission);
        double refl = s ? s->mat.refl : (tr ? tr->mat.refl : pl->mat.refl);
        double refr = s ? s->mat.refr : (tr ? tr->mat.refr : pl->mat.refr);
        if (depth > 5) return e;
        if (!refl && !refr) return e + f.mult(radiance(Ray(x, r.d - n * 2 * n.dot(r.d)), depth + 1));
        if (refr) {
            Ray reflRay(x, r.d - n * 2 * n.dot(r.d));
            bool into = n.dot(Vec3(0, 0, 1)) > 0;
            double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(n), cos2t;
            if ((cos2t = 1 - nnt * nnt * (1 - ddn * ddn)) < 0) return e + f.mult(radiance(reflRay, depth + 1));
            Vec3 tdir = (r.d * nnt - n * ((into ? 1 : -1) * (ddn * nnt + sqrt(cos2t)))).norm();
            double a = nt - nc, b = nt + nc, R0 = a * a / (b * b), c = 1 - (into ? -ddn : tdir.dot(n));
            double Re = R0 + (1 - R0) * c * c * c * c * c, Tr = 1 - Re, P = .25 + .5 * Re, RP = Re / P, TP = Tr / (1 - P);
            return e + f.mult(depth > 2 ? (erand48() < P ? radiance(reflRay, depth + 1) * RP : radiance(Ray(x, tdir), depth + 1) * TP) : radiance(reflRay, depth + 1) * Re + radiance(Ray(x, tdir), depth + 1) * Tr);
        }
        Ray reflRay(x, r.d - n * 2 * n.dot(r.d));
        bool into = n.dot(Vec3(0, 0, 1)) > 0;
        Vec3 nc = Vec3(1, 1, 1), nt = Vec3(1.5, 1.5, 1.5), nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(n), cos2t;
        double a = nt.dot(nnt), b = (nt * nnt).dot(n), R0 = (a - b) * (a - b) / ((a + b) * (a + b)), c = 1 - (into ? -ddn : r.d.dot(n));
        double Re = R0 + (1 - R0) * c * c * c * c * c, Tr = 1 - Re, P = .25 + .5 * Re, RP = Re / P, TP = Tr / (1 - P);
        return e + f.mult(depth > 2 ? (erand48() < P ? radiance(reflRay, depth + 1) * RP : radiance(Ray(x, r.d * nnt - n * ((into ? 1 : -1) * (ddn * nnt + sqrt(cos2t)))) * TP) : radiance(reflRay, depth + 1) * Re + radiance(Ray(x, r.d * nnt - n * ((into ? 1 : -1) * (ddn * nnt + sqrt(cos2t)))) * Tr));
    }
    Vec3 radiance(const Ray &r, int depth) const {
        double t = std::numeric_limits<double>::infinity();
        const Sphere *s = NULL;
        const Triangle *tr = NULL;
        const Plane *pl = NULL;
        for (unsigned i = 0; i < spheres.size(); ++i) {
            double d = spheres[i].intersect(r);
            if (d && d < t) { t = d; s = &spheres[i]; }
        }
        for (unsigned i = 0; i < triangles.size(); ++i) {
            double d = triangles[i].intersect(r);
            if (d && d < t) { t = d; tr = &triangles[i]; }
        }
        for (unsigned i = 0; i < planes.size(); ++i) {
            double d = planes[i].intersect(r);
            if (d && d < t) { t = d; pl = &planes[i]; }
        }
        if (!s && !tr && !pl) return bg;
        Vec3 x = r.o + r.d * t, n;
        if (s) { n = (x - s->p).norm(); }
        else if (tr) { n = ((tr->b - tr->a) % (tr->c - tr->a)).norm(); }
        else { n = pl->n; }
        Vec3 f = s ? s->mat.color : (tr ? tr->mat.color : pl->mat.color);
        Vec3 e = s ? s->mat.emission : (tr ? tr->mat.emission : pl->mat.emission);
        double refl = s ? s->mat.refl : (tr ? tr->mat.refl : pl->mat.refl);
        double refr = s ? s->mat.refr : (tr ? tr->mat.refr : pl->mat.refr);
        if (depth > 5) return e;
        if (!refl && !refr) return e + f.mult(intersect(Ray(x, r.d - n * 2 * n.dot(r.d)), depth + 1));
        if (refr) {
            Ray reflRay(x, r.d - n * 2 * n.dot(r.d));
            bool into = n.dot(Vec3(0, 0, 1)) > 0;
            double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(n), cos2t;
            if ((cos2t = 1 - nnt * nnt * (1 - ddn * ddn)) < 0) return e + f.mult(intersect(reflRay, depth + 1));
            Vec3 tdir = (r.d * nnt - n * ((into ? 1 : -1) * (ddn * nnt + sqrt(cos2t)))).norm();
            double a = nt - nc, b = nt + nc, R0 = a * a / (b * b), c = 1 - (into ? -ddn : tdir.dot(n));
            double Re = R0 + (1 - R0) * c * c * c * c * c, Tr = 1 - Re, P = .25 + .5 * Re, RP = Re / P, TP = Tr / (1 - P);
            return e + f.mult(depth > 2 ? (erand48() < P ? intersect(reflRay, depth + 1) * RP : intersect(Ray(x, tdir), depth + 1) * TP) : intersect(reflRay, depth + 1) * Re + intersect(Ray(x, tdir), depth + 1) * Tr);
        }
        Ray reflRay(x, r.d - n * 2 * n.dot(r.d));
        bool into = n.dot(Vec3(0, 0, 1)) > 0