Skip to content

Instantly share code, notes, and snippets.

@horitaku1124

horitaku1124/Makefile

Last active June 12, 2018 12:48
Show Gist options
  • Save horitaku1124/5ca4853d12ad4b2c8ac64cc27fd2115c to your computer and use it in GitHub Desktop.
Save horitaku1124/5ca4853d12ad4b2c8ac64cc27fd2115c to your computer and use it in GitHub Desktop.
Download ZIP
smallpt for OpenMP original is here http://www.kevinbeason.com/smallpt/
Raw
libs.cc
struct Vec;
class Texture;
class SolidTexture;
const double EPS = 1e-4;
enum Refl_t { DIFF, SPEC, REFR }; // material types, used in radiance()
struct Ray {
Vec obj;
Vec dist;
Ray(Vec o_, Vec d_)
: obj(o_), dist(d_) {}
};
struct Col {
Vec emission;
Vec color;
Refl_t reflection;
Col();
Col(Vec _emission, Vec _color, Refl_t _reflection):
emission(_emission), color(_color), reflection(_reflection) {}
};
struct Point {
double x, y;
Point(double _x, double _y) {
x = _x;
y = _y;
}
};
/*
struct _Sphere {
double rad; // radius
Vec pos, emission, color; // 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_), pos(p_), emission(e_), color(c_), refl(refl_) {}
double intersect(const Ray &r) const { // returns distance, 0 if nohit
Vec op = pos - r.obj; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
double t,
b = op.dot(r.dist),
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 Surface {
public:
Vec pos;
Texture* texture;
Surface() {
}
Surface(Vec _pos, Texture* _texture):
pos(_pos), texture(_texture) {}
Surface(Vec _pos, Vec emission, Vec color, Refl_t reflection) {
pos = _pos;
texture = (Texture*)SolidTexture(emission, color, reflection);
}
// {
// // Surface(pos, SolidTexture(emission, color, reflection));
// pos(_pos);
// }
double intersect(Ray y);
void position(Vec p, Ray r, Vec* n, Col* c);
Point* makeXY(Vec p);
};
class Texture {
public:
Col getCol(Surface s, Vec x);
bool isHit(Surface s, Vec x);
};
class Sphere: public Surface{
public:
double rad; // radius
Sphere(double _rad, Vec p, Vec e, Vec c, Refl_t refl) {
// Surface(p, e, c, refl);
// rad = _rad;
}
Sphere(double _rad, Vec p, Texture texture) {
// Surface(p, texture);
// rad = _rad;
}
double intersect(const Ray &r) const { // returns distance, 0 if nohit
// Vec op = pos - r.obj; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
// double t, b=op.dot(r.dist), det=b*b-op.dot(op)+rad*rad;
// if (det<0) return 0;
// det = sqrt(det);
// return (t=b-det) > EPS ? t : ((t=b+det) > EPS ? t : 0);
}
void position(Vec x, Ray r, Vec n, Col* c) {
// n[0] = (x - pos).norm();
// c[0] = texture.getCol(*this, x);
}
Point* makeXY(Vec x) {
// Vec position = (x - pos) * (1 / rad);
// double phi = atan2(position.z, position.x);
// double theta = asin(position.y);
// double X = 1 - (phi + M_PI) / (2 * M_PI);
// double Y = (theta + M_PI / 2) / M_PI;
// return Point(X, Y);
return nullptr;
}
};
class Plane:public Surface {
double width, height;
public:
Plane(double x, double y, Vec pos, Vec emission, Vec color, Refl_t reflection) {
// Surface(pos, emission, color, reflection);
// width = x;
// height = y;
}
Plane(double x, double y, Vec pos, Texture tex) {
// Surface(pos, tex);
// width = x;
// height = y;
}
double intersect(Ray* ray) {
// if (ray.dist.z < EPS && ray.dist.z > -EPS) {
// return 0;
// }
// double d = (pos.z - ray.obj.z) / ray.dist.z;
// if (d < 0) {
// return 0;
// }
// Vec x = ray.obj + (ray.dist - d);
// Vec p = x - pos;
// if (p.x < 0 || p.x > width || p.y < 0 || p.y > height) {
// return 0;
// }
// if (!texture.isHit(*this, x)) {
// return 0;
// }
// return d;
}
void position(Vec p, Ray* r, Vec n, Col* c) {
// n[0] = Vec(0, 0, r.dist.z > 0 ? -1 : 1);
// c[0] = texture.getCol(*this, p);
}
Point makeXY(Vec p) {
// return Point((p.x - pos.x) / width, (p.y - pos.y) / height);
}
};
class SolidTexture: public Texture {
Col col;
public:
SolidTexture(Vec emission, Vec color, Refl_t ref) {
col = Col(emission, color, ref);
}
Col getCol(Surface s, Vec x) {
return col;
}
};
/*
class CheckTexture: public Texture {
Col col1, col2;
double freq;
public:
CheckTexture(Vec _col1, Vec _col2, double _freq) {
// col1 = Col(Vec(), _col1, DIFF);
// col2 = Col(Vec(), _col2, DIFF);
// freq = freq;
}
Col getCol(Surface s, Vec x) {
// Point p = s.makeXY(x);
// return (under(p.x / freq) - 0.5) * (under(p.y / freq) - 0.5) > 0 ? col1 : col2;
}
double under(double d) {
return d - floor(d);
}
};
class BitmapTexture:public Texture {
public:
BitmapTexture();
BitmapTexture(char* file);
// BufferedImage img;
// int width, height;
// Vec emission = Vec();
// public:
// BitmapTexture();
// BitmapTexture(std::string file) {
// try {
// img = ImageIO.read(SmallPT.class.getResourceAsStream(file));
// width = img.getWidth(null);
// height = img.getHeight(null);
// } catch (IOException ex) {
// throw new UncheckedIOException(ex);
// }
// }
// Col getCol(Surface s, Vec x) {
// int rgb = getRgb(s, x);
// return Col(emission, Vec((rgb >> 16 & 255) / 255., (rgb >> 8 & 255) / 255., (rgb & 255) / 255.), DIFF);
// }
// bool isHit(Surface s, Vec x) {
// int rgb = getRgb(s, x);
// return rgb >> 24 != 0;
// }
// int getRgb(Surface s, Vec x) {
// Point pos = s.makeXY(x);
// return img.getRGB((int)(pos.x * width), (int)((1 - pos.y) * height));
// }
};
class EmissionTexture:public BitmapTexture {
// Vec emission = Vec(12, 12, 12);
// Vec color = Vec();
public:
EmissionTexture(char* file) {
(file);
}
Col getCol(Surface s, Vec x) {
// return Col(emission, color, DIFF);
}
bool isHit(Surface s, Vec x) {
// int rgb = getRgb(s, x);
// return (rgb >> 24 != 0) && (rgb >> 16 & 255) < 80;
}
};
*/
Raw
Makefile
all: smallpt
smallpt:
g++-6 -O3 -fopenmp smallpt.cc -o smallpt
clean:
$(RM) smallpt
Raw
run.md 
g++-6 -O3 -fopenmp test2.cc -o test2 -mavx
g++-6 -O3 -fopenmp smallpt.cc -o smallpt
Raw
smallpt.cc
#include <math.h> // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h> // Remove "-fopenmp" for g++ version < 4.2
#include <time.h>
#include <pthread.h>
#include <sys/time.h>
#include "src/vec.cc"
#include "src/libs.cc"
static Vec UNIT_X = Vec(1, 0, 0);
static Vec UNIT_Y = Vec(0, 1, 0);
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(.75,.75,.75),DIFF),//Botm
Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),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,78), Vec(),Vec(1,1,1)*.999, REFR),//Glas
Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF) //Lite
};
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 double max(double x, double y, double z) {
double m = x;
m = m < y ? y : m;
m = m < z ? z : m;
return m;
}
inline bool intersect(const Ray &r, double &t, Surface* robj){
double n = sizeof(spheres) / sizeof(Sphere);
double d;
double inf = t = 1e20;
for(int i = int(n);i--;) {
Surface obj = spheres[i];
d = spheres[i].intersect(r);
if(d && d < t){
t = d;
// id = i;
robj = &obj;
}
}
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
Surface* obj;
Vec* rn;
Col* rc;
if (!intersect(r, t, obj))
return Vec(); // if miss, return black
Vec x = r.obj + r.dist * t;
obj->position(x, r, rn, rc);
Col tex = rc[0];
Vec n = rn[0];
Vec
nl = n.dot(r.dist) < 0 ? n : n * -1,
f = tex.color;
double p = max(f.x, f.y, f.z); // max refl
// 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 (depth < 50 && erand48(Xi) < p)
f = f * (1 / p);
else
return tex.emission; //R.R.
}
if (tex.reflection == DIFF){ // Ideal DIFFUSE reflection
double r1 = 2 * M_PI * erand48(Xi),
r2 = erand48(Xi),
r2s = sqrt(r2);
Vec vec_w = nl,
vec_u = ((fabs(vec_w.x) > .1 ? UNIT_Y : UNIT_X) % vec_w).norm(),
vec_v = vec_w % vec_u;
Vec d = (vec_u * cos(r1) * r2s + vec_v * sin(r1) * r2s + vec_w * sqrt(1 - r2)).norm();
return tex.emission + f.mult(radiance(Ray(x,d),depth,Xi));
} else if (tex.reflection == SPEC) {
// Ideal SPECULAR reflection
return tex.emission + f.mult(radiance(Ray(x,r.dist - n * 2 * n.dot(r.dist)),depth,Xi));
}
Ray reflRay(x, r.dist -n*2*n.dot(r.dist)); // 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.dist.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
return tex.emission + f.mult(radiance(reflRay,depth,Xi));
Vec tdir = (r.dist * 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 tex.emission + 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);
}
double getTimeStamp()
{
struct timeval tv;
gettimeofday(&tv, NULL);
return ((double)(tv.tv_sec) + (double)(tv.tv_usec) * 0.000001);
}
int main(int argc, char *argv[])
{
double start = getTimeStamp();
int w = 1024,
h = 768,
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),
cy = (cx % cam.dist).norm() * .5135,
r,
*c = new Vec[w * h];
printf("smallpt v:0.0.2 samps=%d\n", samps);
#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));
unsigned short yyy = y * y * y;
for (unsigned short x=0, Xi[3]={0,0,yyy}; 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.dist;
r = r + radiance(Ray(cam.obj + 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;
}
}
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));
double stop = getTimeStamp();
printf ( "time :%.3f sec\n", (stop - start));
}
Raw
test.cc
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <pthread.h>
#include <sys/time.h>
#include <immintrin.h>//AVX: -mavx
#include "vec.cc"
Vec test1(Vec vec_w, Vec vec_u)
{
Vec nl = Vec(0.1, 0.1, 0.1);
double r1 = 2 * M_PI * 0.1,
r2 = 0.7,
r2s = sqrt(r2);
Vec vec_v = vec_w % vec_u;
Vec d = (vec_u * cos(r1) * r2s + vec_v * sin(r1) * r2s + vec_w * sqrt(1 - r2)).norm();
return d;
}
Vec test2(Vec vec_w, Vec vec_u)
{
Vec nl = Vec(0.1, 0.1, 0.1);
double r1 = 2 * M_PI * 0.1,
r2 = 0.7,
r2s = sqrt(r2);
// vec_v = vec_w % vec_u;
double ux = vec_u.x,
uy = vec_u.y,
uz = vec_u.z;
// vec_u * cos(r1) * r2s
double cos_r1_r2s = cos(r1) * r2s;
ux *= cos_r1_r2s;
uy *= cos_r1_r2s;
uz *= cos_r1_r2s;
double vx = vec_w.y * vec_u.z - vec_w.z * vec_u.y,
vy = vec_w.z * vec_u.x - vec_w.x * vec_u.z,
vz = vec_w.x * vec_u.y - vec_w.y * vec_u.x;
// vec_v * sin(r1) * r2s
double sin_r2s = sin(r1) * r2s;
vx *= sin_r2s,
vy *= sin_r2s,
vz *= sin_r2s;
double sqrtMr2 = sqrt(1 - r2);
double wx = vec_w.x * sqrtMr2,
wy = vec_w.y * sqrtMr2,
wz = vec_w.z * sqrtMr2;
double dx = ux + vx + wx,
dy = uy + vy + wy,
dz = uz + vz + wz;
// d.norm();
double norm = (1 / sqrt(dx * dx + dy * dy + dz * dz));
dx *= norm;
dy *= norm;
dz *= norm;
return Vec(dx, dy, dz);
}
Vec test3(Vec vec_w, Vec vec_u)
{
double r1 = 2 * M_PI * 0.1,
r2 = 0.7,
r2s = sqrt(r2);
double __attribute__((aligned(32))) i_vec_u[4] = {0, 0, 0, 0};
double __attribute__((aligned(32))) i_vec_v[4] = {0, 0, 0, 0};
double __attribute__((aligned(32))) i_vec_w[4] = {0, 0, 0, 0};
double __attribute__((aligned(32))) multipulier[4] = {0, 0, 0, 0};
double __attribute__((aligned(32))) dist[4] = {0};
double vx = vec_w.y * vec_u.z - vec_w.z * vec_u.y,
vy = vec_w.z * vec_u.x - vec_w.x * vec_u.z,
vz = vec_w.x * vec_u.y - vec_w.y * vec_u.x;
// vec_u * cos(r1) * r2s
double cos_r1 = cos(r1) * r2s;
double sin_r2s = sin(r1) * r2s;
double sqrtMr2 = sqrt(1 - r2);
multipulier[0] = cos_r1;
multipulier[1] = cos_r1;
multipulier[2] = cos_r1;
i_vec_u[0] = vec_u.x;
i_vec_u[1] = vec_u.y;
i_vec_u[2] = vec_u.z;
i_vec_v[0] = vx;
i_vec_v[1] = vy;
i_vec_v[2] = vz;
i_vec_w[0] = vec_w.x;
i_vec_w[1] = vec_w.y;
i_vec_w[2] = vec_w.z;
__m256d v1, v2, v3;
v1 = _mm256_load_pd(i_vec_u);
v2 = _mm256_load_pd(multipulier);
v1 = _mm256_mul_pd(v1, v2);
multipulier[0] = sin_r2s;
multipulier[1] = sin_r2s;
multipulier[2] = sin_r2s;
v3 = _mm256_load_pd(i_vec_v);
v2 = _mm256_load_pd(multipulier);
v3 = _mm256_mul_pd(v3, v2);
// w + v + z
v1 = _mm256_add_pd(v1, v3);
multipulier[0] = sqrtMr2;
multipulier[1] = sqrtMr2;
multipulier[2] = sqrtMr2;
v3 = _mm256_load_pd(i_vec_w);
v2 = _mm256_load_pd(multipulier);
v3 = _mm256_mul_pd(v3, v2);
v1 = _mm256_add_pd(v1, v3);
_mm256_store_pd(i_vec_u, v1);
double dx = i_vec_u[0],
dy = i_vec_u[1],
dz = i_vec_u[2];
// d.norm();
double norm = (1 / sqrt(dx * dx + dy * dy + dz * dz));
dx *= norm;
dy *= norm;
dz *= norm;
return Vec(dx, dy, dz);
}
int main(int argc, char *argv[])
{
Vec nl = Vec(0.1, 0.1, 0.1);
Vec vec_w = nl,
vec_u = ((fabs(vec_w.x) > .1 ? Vec(0,1) : Vec(1)) % vec_w).norm();
Vec d1 = test1(vec_w, vec_u);
printf("d1.x=%.4f y=%.4f z=%.4f\n", d1.x, d1.y, d1.z);
Vec d2 = test2(vec_w, vec_u);
printf("d2.x=%.4f y=%.4f z=%.4f\n", d2.x, d2.y, d2.z);
Vec d3 = test3(vec_w, vec_u);
printf("d3.x=%.4f y=%.4f z=%.4f\n", d3.x, d3.y, d3.z);
return 0;
}
Raw
vec.cc
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
);
}
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment