纯c++
效果图：
代码：
第一张图的代码
#include <time.h>
#include <iostream>
#include <math.h> // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h>
#include <stdio.h>
//#define M_PI 3.1415 //如使用VS编译去掉该行注释
//补充随机数实现（如使用LCG随机数可加快速度）
double erand48(unsigned short xsubi[3]) {
return (double)rand() / (double)RAND_MAX;
}
struct Vec { // Usage: time ./smallpt 5000 && xv image.ppm
double x, y, z; // position, also color (r,g,b)
Vec(double x_ = 0, double y_ = 0, double z_ = 0) { x = x_; y = y_; z = z_; }
Vec operator+(const Vec &b) const { return Vec(x + b.x, y + b.y, z + b.z); }
Vec operator-(const Vec &b) const { return Vec(x - b.x, y - b.y, z - b.z); }
Vec operator*(double b) const { return Vec(x*b, y*b, z*b); }
Vec mult(const Vec &b) const { return Vec(x*b.x, y*b.y, z*b.z); }
Vec& norm() { return *this = *this * (1 / sqrt(x*x + y * y + z * z)); }
double dot(const Vec &b) const { return x * b.x + y * b.y + z * b.z; } // cross:
Vec operator%(Vec&b) { return Vec(y*b.z - z * b.y, z*b.x - x * b.z, x*b.y - y * b.x); }
Vec operator-() const { return Vec(-x, -y, -z); }
};
struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR }; // material types, used in radiance()
struct Sphere {
double rad; // radius
Vec p, e, c; // position, emission, color
Refl_t refl; // reflection type (DIFFuse, SPECular, REFRactive)
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) :
rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
double intersect(const Ray &r) const { // returns distance, 0 if nohit
Vec op = p - r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
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);
}
};
struct Cube
{
Vec min, max;
Vec e, c; //emission, color
Refl_t refl;
Cube(Vec min_, Vec max_, Vec e_, Vec c_, Refl_t refl_) :
min(min_), max(max_), e(e_), c(c_), refl(refl_) {}
double intersect(const Ray &r, Vec &returnNormal) const //AABB检测
{
const Vec rayOrg = r.o;
const Vec rayDelta = r.d;
const double kNoIntersection = 1e30f;
bool inside = true;
double xt, xn;
if (rayOrg.x < min.x)
{
xt = min.x - rayOrg.x;
if (rayDelta.x < 1e-4 && rayDelta.x > -1e-4)
{
return kNoIntersection;
}
xt /= rayDelta.x;
inside = false;
xn = -1.0f;
}
else if (rayOrg.x > max.x)
{
xt = max.x - rayOrg.x;
if (rayDelta.x < 1e-4 && rayDelta.x > -1e-4)
{
return kNoIntersection;
}
xt /= rayDelta.x;
inside = false;
xn = 1.0f;
}
else
{
xt = -1.0f;
}
double yt, yn;
if (rayOrg.y < min.y)
{
yt = min.y - rayOrg.y;
if (rayDelta.y < 1e-4 && rayDelta.y > -1e-4)
{
return kNoIntersection;
}
yt /= rayDelta.y;
inside = false;
yn = -1.0f;
}
else if (rayOrg.y > max.y)
{
yt = max.y - rayOrg.y;
if (rayDelta.y < 1e-4 && rayDelta.y > -1e-4)
{
return kNoIntersection;
}
yt /= rayDelta.y;
inside = false;
yn = 1.0f;
}
else
{
yt = -1.0f;
}
double zt, zn;
if (rayOrg.z < min.z)
{
zt = min.z - rayOrg.z;
if (rayDelta.z < 1e-4 && rayDelta.z > -1e-4)
{
return kNoIntersection;
}
zt /= rayDelta.z;
inside = false;
zn = -1.0f;
}
else if (rayOrg.z > max.z)
{
zt = max.z - rayOrg.z;
if (rayDelta.z < 1e-4 && rayDelta.z > -1e-4)
{
return kNoIntersection;
}
zt /= rayDelta.z;
inside = false;
zn = 1.0f;
}
else
{
zt = -1.0f;
}
//inside?
if (inside)
{
returnNormal = -rayDelta;
returnNormal.norm();
return 0.0f;
}
//far
int which = 0;
double t = xt;
if (yt > t)
{
which = 1;
t = yt;
}
if (zt > t)
{
which = 2;
t = zt;
}
switch (which)
{
case 0:
{
double y = rayOrg.y + rayDelta.y*t;
if (y < min.y || y > max.y)
{
return kNoIntersection;
}
double z = rayOrg.z + rayDelta.z*t;
if (z < min.z || z > max.z)
{
return kNoIntersection;
}
returnNormal.x = xn;
returnNormal.y = 0.0f;
returnNormal.z = 0.0f;
}break;
case 1:
{
double x = rayOrg.x + rayDelta.x*t;
if (x < min.x || x > max.x)
{
return kNoIntersection;
}
double z = rayOrg.z + rayDelta.z*t;
if (z < min.z || z > max.z)
{
return kNoIntersection;
}
returnNormal.x = 0.0f;
returnNormal.y = yn;
returnNormal.z = 0.0f;
}break;
case 2:
{
double y = rayOrg.y + rayDelta.y*t;
if (y < min.y || y > max.y)
{
return kNoIntersection;
}
double x = rayOrg.x + rayDelta.x*t;
if (x < min.x || x > max.x)
{
return kNoIntersection;
}
returnNormal.x = 0.0f;
returnNormal.y = 0.0f;
returnNormal.z = zn;
}break;
}
return t;
}
};
//平行六面体
const int l = 3; //单位正方体的边长
Vec Ven(0, 0, 4 * l); //整体偏移
double z = 100; //z轴起始
double xoffset = 20.5; //x轴偏移
double yoffset = 0; //y轴偏移
//颜色
Vec Pink = Vec(1, .7529, .7961);
Vec Red = Vec(1, 0, 0);
Vec Lavender = Vec(.90196, .90196, .980392);
Vec CubeColor = Red;
Vec CubeEmission = Vec();
Cube cubes[] = {
/*
//I
Cube(Vec(-46 + xoffset,0 + yoffset,z),Vec(-46 + xoffset + l * 2, 0 + yoffset + l * 9, z + l) + Ven,CubeEmission,CubeColor,REFR),
*/
//L
Cube(Vec(-19 + xoffset,0 + yoffset,z),Vec(-19 + xoffset + l * 2, 0 + yoffset + l * 9, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(-19 + xoffset + 2 * l,0 + yoffset,z),Vec(-19 + xoffset + l * 2 + 5 * l, 0 + yoffset + l, z + l) + Ven,CubeEmission,CubeColor,REFR),
//O
Cube(Vec(7 + xoffset,l + yoffset,z),Vec(7 + xoffset + 2 * l, l + yoffset + 7 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(7 + xoffset + 5 * l,l + yoffset,z),Vec(7 + xoffset + 2 * l + 5 * l, l + yoffset + 7 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(7 + xoffset + l,0 + yoffset,z),Vec(7 + xoffset + l + 5 * l, l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(7 + xoffset + l,8 * l + yoffset,z),Vec(7 + xoffset + l + 5 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
//V
Cube(Vec(33 - l * .5 + xoffset,7 * l + yoffset,z),Vec(33 - l * .5 + xoffset + 2 * l, 7 * l + yoffset + 2 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(33 - l * .5 + xoffset + 6 * l,7 * l + yoffset,z),Vec(33 - l * .5 + xoffset + 6 * l + 2 * l, 7 * l + yoffset + 2 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(33 - l * .5 + xoffset + l,4 * l + yoffset,z),Vec(33 - l * .5 + xoffset + l + 2 * l, 4 * l + yoffset + 3 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(33 - l * .5 + xoffset + l + 4 * l,4 * l + yoffset,z),Vec(33 - l * .5 + xoffset + l + 2 * l + 4 * l, 4 * l + yoffset + 3 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(33 - l * .5 + xoffset + 2 * l,2 * l + yoffset,z),Vec(33 - l * .5 + xoffset + 4 * l + 2 * l, 2 * l + yoffset + 2 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(33 - l * .5 + xoffset + 3 * l,0 * l + yoffset,z),Vec(33 - l * .5 + xoffset + 2 * l + 3 * l, 0 * l + yoffset + 2 * l, z + l) + Ven,CubeEmission,CubeColor,REFR),
//E
Cube(Vec(59 + xoffset + 2 * l,0 + yoffset,z),Vec(59 + xoffset + 2 * l + 5 * l, l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(59 + xoffset + 2 * l,8 * l + yoffset,z),Vec(59 + xoffset + 2 * l + 5 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(59 + xoffset,0 + yoffset,z),Vec(59 + xoffset + 2 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(59 + xoffset + 2 * l,4 * l + yoffset,z),Vec(59 + xoffset + 2 * l + 4 * l, 5 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
/*
//Y
Cube(Vec(95 - l * .5 + xoffset,7 * l + yoffset,z),Vec(95 + xoffset + 2 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(95 - l * .5 + xoffset + 6 * l,7 * l + yoffset,z),Vec(95 + xoffset + 2 * l + 6 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(95 - l * .5 + xoffset + l,6 * l + yoffset,z),Vec(95 + xoffset + l + 2 * l, 7 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(95 - l * .5 + xoffset + l + 4 * l,6 * l + yoffset,z),Vec(95 + xoffset + l + 2 * l + 4 * l, 7 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(95 - l * .5 + xoffset + 2 * l,5 * l + yoffset,z),Vec(95 + xoffset + 2 * l + 4 * l, 6 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(95 - l * .5 + xoffset + 3 * l,0 * l + yoffset,z),Vec(95 + xoffset + 3 * l + 2 * l, 6 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
//O
Cube(Vec(121 + xoffset,l + yoffset,z),Vec(121 + xoffset + 2 * l, l + 7 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(121 + xoffset + 5 * l,l + yoffset,z),Vec(121 + xoffset + 2 * l + 5 * l, l + 7 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(121 + xoffset + l,0 + yoffset,z),Vec(121 + xoffset + l + 5 * l, l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(121 + xoffset + l,8 * l + yoffset,z),Vec(121 + xoffset + l + 5 * l, 8 * l + l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
//U
Cube(Vec(147 + xoffset,1 * l + yoffset,z),Vec(147 + xoffset + 2 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(147 + xoffset + 5 * l,1 * l + yoffset,z),Vec(147 + xoffset + 2 * l + 5 * l, 9 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
Cube(Vec(147 + xoffset + l,0 * l + yoffset,z),Vec(147 + xoffset + l + 5 * l, 1 * l + yoffset, z + l) + Ven,CubeEmission,CubeColor,REFR),
*/
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
Sphere(1e5, Vec(1e5 + 1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
Sphere(1e5, Vec(-1e5 + 99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
Sphere(1e5, Vec(50,40.8, 1e5), Vec(),Vec(.75,.75,.75),DIFF),//Back
Sphere(1e5, Vec(50,40.8,-1e5 + 170), Vec(),Vec(), DIFF),//Frnt
Sphere(1e5, Vec(50, 1e5, 81.6), Vec(),Vec(1,1,0),DIFF),//Botm
Sphere(1e5, Vec(50,-1e5 + 81.6,81.6),Vec(),Vec(.70,.95,.23),DIFF),//Top
Sphere(16.5,Vec(27,16.5,47), Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
// Sphere(16.5,Vec(73,16.5,47), Vec(),Vec(1,1,1)*.999, REFR),//Glas
Sphere(600, Vec(50,681.6 - .27,81.6),Vec(12,12,12), Vec(), DIFF) //Lite
};
const int n1 = sizeof(spheres) / sizeof(Sphere);
const int n2 = sizeof(cubes) / sizeof(Cube);
inline double clamp(double x) { return x < 0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x) { return int(pow(clamp(x), 1 / 2.2) * 255 + .5); }
inline bool intersect(const Ray &r, double &t, int &id, Vec &returnNormal) {
double d, inf = t = 1e20;
for (int i = n1 + n2; i--;)
{
if (i >= n2)
{
if ((d = spheres[i - n2].intersect(r)) && d < t) { t = d; id = i; }
}
else
{
if ((d = cubes[i].intersect(r, returnNormal)) && d < t)
{
t = d; id = i;
}
}
}
return t < inf;
}
Vec radiance(const Ray &r, int depth, unsigned short *Xi) {
double t; // distance to intersection
int id = 0; // id of intersected object
Vec returnNormal;
if (!intersect(r, t, id, returnNormal)) return Vec(); // if miss, return black
if (id >= n2)
{
const Sphere &obj = spheres[id - n2]; // the hit object
Vec x = r.o + r.d*t, n = (x - obj.p).norm(), nl = n.dot(r.d) < 0 ? n : n * -1, f = obj.c;
double p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z; // max refl
if (++depth > 5) if (erand48(Xi) < p) f = f * (1 / p); else return obj.e; //R.R.
if (obj.refl == DIFF) { // Ideal DIFFUSE reflection
double r1 = 2 * M_PI*erand48(Xi), r2 = erand48(Xi), r2s = sqrt(r2);
Vec w = nl, u = ((fabs(w.x) > .1 ? Vec(0, 1, 0) : Vec(1, 0, 0)) % w).norm(), v = w % u;
Vec d = (u*cos(r1)*r2s + v * sin(r1)*r2s + w * sqrt(1 - r2)).norm();
return obj.e + f.mult(radiance(Ray(x, d), depth, Xi));
}
else if (obj.refl == SPEC) // Ideal SPECULAR reflection
return obj.e + f.mult(radiance(Ray(x, r.d - n * 2 * n.dot(r.d)), depth, Xi));
Ray reflRay(x, r.d - n * 2 * n.dot(r.d)); // Ideal dielectric REFRACTION
bool into = n.dot(nl) > 0; // Ray from outside going in?
double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(nl), cos2t;
if ((cos2t = 1 - nnt * nnt*(1 - ddn * ddn)) < 0) // Total internal reflection
return obj.e + f.mult(radiance(reflRay, depth, Xi));
Vec 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 obj.e + f.mult(depth > 2 ? (erand48(Xi) < P ? // Russian roulette
radiance(reflRay, depth, Xi)*RP : radiance(Ray(x, tdir), depth, Xi)*TP) :
radiance(reflRay, depth, Xi)*Re + radiance(Ray(x, tdir), depth, Xi)*Tr);
}
else
{
const Cube &obj = cubes[id];
Vec x = r.o + r.d*t, n = returnNormal, nl = n.dot(r.d) < 0 ? n : n * -1, f = obj.c;
double p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z; // max refl
if (++depth > 5) if (erand48(Xi) < p) f = f * (1 / p); else return obj.e; //R.R.
Ray reflRay(x, r.d - n * 2 * n.dot(r.d)); // Ideal dielectric REFRACTION
bool into = n.dot(nl) > 0; // Ray from outside going in?
double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(nl), cos2t;
if ((cos2t = 1 - nnt * nnt*(1 - ddn * ddn)) < 0) // Total internal reflection
return obj.e + f.mult(radiance(reflRay, depth, Xi));
Vec 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 obj.e + f.mult(depth > 2 ? (erand48(Xi) < P ? // Russian roulette
radiance(reflRay, depth, Xi)*RP : radiance(Ray(x, tdir), depth, Xi)*TP) :
radiance(reflRay, depth, Xi)*Re + radiance(Ray(x, tdir), depth, Xi)*Tr);
}
}
int main(int argc, char *argv[]) {
int w = 1024 * .25, h = 768 * .25, samps = argc == 2 ? atoi(argv[1]) / 4 : 1; // # samples
Ray cam(Vec(50, 52, 295.6), Vec(0, -0.042612, -1).norm()); // cam pos, dir
Vec cx = Vec(w*.5135 / h, 0, 0), cy = (cx%cam.d).norm()*.5135, r, *c = new Vec[w*h];
clock_t start_time = clock();
#pragma omp parallel for schedule(dynamic, 1) private(r) // OpenMP
for (int y = 0; y < h; y++) { // Loop over image rows
fprintf(stderr, "\rRendering (%d spp) %5.2f%%", samps * 4, 100.*y / (h - 1));
for (unsigned short x = 0, Xi[3] = { 0,0,y*y*y }; x < w; x++) // Loop cols
for (int sy = 0, i = (h - y - 1)*w + x; sy < 2; sy++) // 2x2 subpixel rows
for (int sx = 0; sx < 2; sx++, r = Vec()) { // 2x2 subpixel cols
for (int s = 0; s < samps; s++)
{
double r1 = 2 * erand48(Xi), dx = r1 < 1 ? sqrt(r1) - 1 : 1 - sqrt(2 - r1);
double r2 = 2 * erand48(Xi), dy = r2 < 1 ? sqrt(r2) - 1 : 1 - sqrt(2 - r2);
Vec d = cx * (((sx + .5 + dx) / 2 + x) / w - .5) +
cy * (((sy + .5 + dy) / 2 + y) / h - .5) + cam.d;
r = r + radiance(Ray(cam.o + d * 140, d.norm()), 0, Xi)*(1. / samps);
} // Camera rays are pushed ^^^^^ forward to start in interior
c[i] = c[i] + Vec(clamp(r.x), clamp(r.y), clamp(r.z))*.25;
}
}
clock_t end_time = clock();
std::cout << "\nRunning time is: " << static_cast<double>(end_time - start_time) << "ms" << std::endl;
FILE *f = fopen("image.ppm", "w"); // Write image to PPM file.
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
for (int i = 0; i < w*h; i++)
fprintf(f, "%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
}
使用Dev编译：
1.修改编译选项：启用多线程，修改堆栈大小。
2.运行时参数为每像素采样点，默认为4，数值越大效果越好。
参考：smallpt，3D Math Primer for Graphics and Game Development.
扫一扫
1262
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
