aaron/ispc
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

workingon rt example

Browse Source
This commit is contained in:
Evghenii
2014-01-05 11:45:19 +01:00
parent 10798cf2ae
commit 62deb82bf8
14 changed files with 1638 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
examples_ptx/rt/.gitignore vendored Normal file
View File

@@ -0,0 +1,2 @@
rt
*.ppm

8
examples_ptx/rt/Makefile_cpu Normal file
View File

@@ -0,0 +1,8 @@
EXAMPLE=rt
CPP_SRC=rt.cpp rt_serial.cpp
ISPC_SRC=rt.ispc
ISPC_IA_TARGETS=avx1-i32x8
ISPC_ARM_TARGETS=neon
include ../common.mk

BIN
examples_ptx/rt/cornell.bvh Normal file
View File

Binary file not shown.

BIN
examples_ptx/rt/cornell.camera Normal file
View File

Binary file not shown.

298
examples_ptx/rt/rt.cpp Normal file
View File

@@ -0,0 +1,298 @@
/*
Copyright (c) 2010-2011, Intel Corporation
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of Intel Corporation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#define NOMINMAX
#pragma warning (disable: 4244)
#pragma warning (disable: 4305)
#endif
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cassert>
#include <cstring>
#include <sys/types.h>
#include "../timing.h"
#include "../ispc_malloc.h"
#include "rt_ispc.h"
using namespace ispc;
typedef unsigned int uint;
extern void raytrace_serial(int width, int height, int baseWidth, int baseHeight,
const float raster2camera[4][4],
const float camera2world[4][4], float image[],
int id[], const LinearBVHNode nodes[],
const Triangle triangles[]);
static void writeImage(int *idImage, float *depthImage, int width, int height,
const char *filename) {
FILE *f = fopen(filename, "wb");
if (!f) {
perror(filename);
exit(1);
}
fprintf(f, "P6\n%d %d\n255\n", width, height);
for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
// use the bits from the object id of the hit object to make a
// random color
int id = idImage[y * width + x];
unsigned char r = 0, g = 0, b = 0;
for (int i = 0; i < 8; ++i) {
// extract bit 3*i for red, 3*i+1 for green, 3*i+2 for blue
int rbit = (id & (1 << (3*i))) >> (3*i);
int gbit = (id & (1 << (3*i+1))) >> (3*i+1);
int bbit = (id & (1 << (3*i+2))) >> (3*i+2);
// and then set the bits of the colors starting from the
// high bits...
r |= rbit << (7-i);
g |= gbit << (7-i);
b |= bbit << (7-i);
}
fputc(r, f);
fputc(g, f);
fputc(b, f);
}
}
fclose(f);
printf("Wrote image file %s\n", filename);
}
static void usage() {
fprintf(stderr, "rt <scene name base> [--scale=<factor>] [ispc iterations] [tasks iterations] [serial iterations]\n");
exit(1);
}
int main(int argc, char *argv[]) {
static unsigned int test_iterations[] = {3, 7, 1};
float scale = 1.f;
const char *filename = NULL;
if (argc < 2) usage();
filename = argv[1];
if (argc > 2) {
if (strncmp(argv[2], "--scale=", 8) == 0) {
scale = atof(argv[2] + 8);
}
}
if ((argc == 6) || (argc == 5)) {
for (int i = 0; i < 3; i++) {
test_iterations[i] = atoi(argv[argc - 3 + i]);
}
}
#define READ(var, n) \
if (fread(&(var), sizeof(var), n, f) != (unsigned int)n) { \
fprintf(stderr, "Unexpected EOF reading scene file\n"); \
return 1; \
} else /* eat ; */
//
// Read the camera specification information from the camera file
//
char fnbuf[1024];
sprintf(fnbuf, "%s.camera", filename);
FILE *f = fopen(fnbuf, "rb");
if (!f) {
perror(fnbuf);
return 1;
}
//
// Nothing fancy, and trouble if we run on a big-endian system, just
// fread in the bits
//
int baseWidth, baseHeight;
#if 0
float camera2world[4][4], raster2camera[4][4];
#else
float *camera2world_ispc, *raster2camera_ispc;
ispc_malloc((void**)&camera2world_ispc, 4*4*sizeof(float));
ispc_malloc((void**)&raster2camera_ispc, 4*4*sizeof(float));
float (*camera2world )[4] = (float (*)[4])camera2world_ispc;
float (*raster2camera)[4] = (float (*)[4])raster2camera_ispc;
#endif
READ(baseWidth, 1);
READ(baseHeight, 1);
READ(camera2world[0][0], 16);
READ(raster2camera[0][0], 16);
//
// Read in the serialized BVH
//
sprintf(fnbuf, "%s.bvh", filename);
f = fopen(fnbuf, "rb");
if (!f) {
perror(fnbuf);
return 1;
}
// The BVH file starts with an int that gives the total number of BVH
// nodes
uint nNodes;
READ(nNodes, 1);
#if 0
LinearBVHNode *nodes = new LinearBVHNode[nNodes];
#else
LinearBVHNode *nodes;
ispc_malloc((void**)&nodes, nNodes*sizeof(LinearBVHNode));
#endif
for (unsigned int i = 0; i < nNodes; ++i) {
// Each node is 6x floats for a boox, then an integer for an offset
// to the second child node, then an integer that encodes the type
// of node, the total number of int it if a leaf node, etc.
float b[6];
READ(b[0], 6);
nodes[i].bounds[0][0] = b[0];
nodes[i].bounds[0][1] = b[1];
nodes[i].bounds[0][2] = b[2];
nodes[i].bounds[1][0] = b[3];
nodes[i].bounds[1][1] = b[4];
nodes[i].bounds[1][2] = b[5];
READ(nodes[i].offset, 1);
READ(nodes[i].nPrimitives, 1);
READ(nodes[i].splitAxis, 1);
READ(nodes[i].pad, 1);
}
// And then read the triangles
uint nTris;
READ(nTris, 1);
Triangle *triangles = new Triangle[nTris];
for (uint i = 0; i < nTris; ++i) {
// 9x floats for the 3 vertices
float v[9];
READ(v[0], 9);
float *vp = v;
for (int j = 0; j < 3; ++j) {
triangles[i].p[j][0] = *vp++;
triangles[i].p[j][1] = *vp++;
triangles[i].p[j][2] = *vp++;
}
// And create an object id
triangles[i].id = i+1;
}
fclose(f);
int height = int(baseHeight * scale);
int width = int(baseWidth * scale);
// allocate images; one to hold hit object ids, one to hold depth to
// the first interseciton
#if 0
int *id = new int[width*height];
float *image = new float[width*height];
#else
int *id;
float *image;
ispc_malloc((void**)&id, sizeof( int)*width*height);
ispc_malloc((void**)&image, sizeof(float)*width*height);
#endif
//
// Run 3 iterations with ispc + 1 core, record the minimum time
//
#ifndef _CUDA_
double minTimeISPC = 1e30;
for (int i = 0; i < test_iterations[0]; ++i) {
reset_and_start_timer();
raytrace_ispc(width, height, baseWidth, baseHeight, raster2camera,
camera2world, image, id, nodes, triangles);
double dt = get_elapsed_msec();
printf("@time of ISPC run:\t\t\t[%.3f] msec\n", dt);
minTimeISPC = std::min(dt, minTimeISPC);
}
printf("[rt ispc, 1 core]:\t\t[%.3f] msec for %d x %d image\n",
minTimeISPC, width, height);
writeImage(id, image, width, height, "rt-ispc-1core.ppm");
#endif
memset(id, 0, width*height*sizeof(int));
memset(image, 0, width*height*sizeof(float));
//
// Run 3 iterations with ispc + 1 core, record the minimum time
//
double minTimeISPCtasks = 1e30;
for (int i = 0; i < test_iterations[1]; ++i) {
reset_and_start_timer();
raytrace_ispc_tasks(width, height, baseWidth, baseHeight, raster2camera,
camera2world, image, id, nodes, triangles);
double dt = get_elapsed_msec();
printf("@time of ISPC + TASKS run:\t\t\t[%.3f] msec\n", dt);
minTimeISPCtasks = std::min(dt, minTimeISPCtasks);
}
printf("[rt ispc + tasks]:\t\t[%.3f] msec for %d x %d image\n",
minTimeISPCtasks, width, height);
writeImage(id, image, width, height, "rt-ispc-tasks.ppm");
memset(id, 0, width*height*sizeof(int));
memset(image, 0, width*height*sizeof(float));
//
// And 3 iterations with the serial implementation, reporting the
// minimum time.
//
double minTimeSerial = 1e30;
for (int i = 0; i < test_iterations[2]; ++i) {
reset_and_start_timer();
raytrace_serial(width, height, baseWidth, baseHeight, raster2camera,
camera2world, image, id, nodes, triangles);
double dt = get_elapsed_msec();
printf("@time of serial run:\t\t\t[%.3f] msec\n", dt);
minTimeSerial = std::min(dt, minTimeSerial);
}
printf("[rt serial]:\t\t\t[%.3f] msec for %d x %d image\n",
minTimeSerial, width, height);
#ifndef _CUDA_
printf("\t\t\t\t(%.2fx speedup from ISPC, %.2fx speedup from ISPC + tasks)\n",
minTimeSerial / minTimeISPC, minTimeSerial / minTimeISPCtasks);
#else
printf("\t\t\t\t(%.2fx speedup from ISPC + tasks)\n",
minTimeSerial / minTimeISPCtasks);
#endif
writeImage(id, image, width, height, "rt-serial.ppm");
return 0;
}

359
examples_ptx/rt/rt.cu Normal file
View File

@@ -0,0 +1,359 @@
#define programCount 32
#define programIndex (threadIdx.x & 31)
#define taskIndex (blockIdx.x*4 + (threadIdx.x >> 5))
#define taskCount (gridDim.x*4)
#define warpIdx (threadIdx.x >> 5)
#define float3 Float3
struct Float3
{
float x,y,z;
__device__ friend Float3 operator+(const Float3 a, const Float3 b)
{
Float3 c;
c.x = a.x+b.x;
c.y = a.y+b.y;
c.z = a.z+b.z;
return c;
}
__device__ friend Float3 operator-(const Float3 a, const Float3 b)
{
Float3 c;
c.x = a.x-b.x;
c.y = a.y-b.y;
c.z = a.z-b.z;
return c;
}
__device__ friend Float3 operator/(const Float3 a, const Float3 b)
{
Float3 c;
c.x = a.x/b.x;
c.y = a.y/b.y;
c.z = a.z/b.z;
return c;
}
__device__ friend Float3 operator/(const float a, const Float3 b)
{
Float3 c;
c.x = a/b.x;
c.y = a/b.y;
c.z = a/b.z;
return c;
}
__device__ friend Float3 operator*(const Float3 a, const Float3 b)
{
Float3 c;
c.x = a.x*b.x;
c.y = a.y*b.y;
c.z = a.z*b.z;
return c;
}
__device__ friend Float3 operator*(const Float3 a, const float b)
{
Float3 c;
c.x = a.x*b;
c.y = a.y*b;
c.z = a.z*b;
return c;
}
};
#define int8 char
#define int16 short
struct Ray {
float3 origin, dir, invDir;
unsigned int dirIsNeg0, dirIsNeg1, dirIsNeg2;
float mint, maxt;
int hitId;
};
struct Triangle {
float p[3][4];
int id;
int pad[3];
};
struct LinearBVHNode {
float bounds[2][3];
unsigned int offset; // num primitives for leaf, second child for interior
unsigned int8 nPrimitives;
unsigned int8 splitAxis;
unsigned int16 pad;
};
__device__
static inline float3 Cross(const float3 v1, const float3 v2) {
float v1x = v1.x, v1y = v1.y, v1z = v1.z;
float v2x = v2.x, v2y = v2.y, v2z = v2.z;
float3 ret;
ret.x = (v1y * v2z) - (v1z * v2y);
ret.y = (v1z * v2x) - (v1x * v2z);
ret.z = (v1x * v2y) - (v1y * v2x);
return ret;
}
__device__
static inline float Dot(const float3 a, const float3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
__device__
inline
static void generateRay( const float raster2camera[4][4],
const float camera2world[4][4],
float x, float y, Ray &ray) {
ray.mint = 0.f;
ray.maxt = 1e30f;
ray.hitId = 0;
// transform raster coordinate (x, y, 0) to camera space
float camx = raster2camera[0][0] * x + raster2camera[0][1] * y + raster2camera[0][3];
float camy = raster2camera[1][0] * x + raster2camera[1][1] * y + raster2camera[1][3];
float camz = raster2camera[2][3];
float camw = raster2camera[3][3];
camx /= camw;
camy /= camw;
camz /= camw;
ray.dir.x = camera2world[0][0] * camx + camera2world[0][1] * camy +
camera2world[0][2] * camz;
ray.dir.y = camera2world[1][0] * camx + camera2world[1][1] * camy +
camera2world[1][2] * camz;
ray.dir.z = camera2world[2][0] * camx + camera2world[2][1] * camy +
camera2world[2][2] * camz;
ray.origin.x = camera2world[0][3] / camera2world[3][3];
ray.origin.y = camera2world[1][3] / camera2world[3][3];
ray.origin.z = camera2world[2][3] / camera2world[3][3];
ray.invDir = 1.f / ray.dir;
#if 0
ray.dirIsNeg[0] = any(ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg[1] = any(ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg[2] = any(ray.invDir.z < 0) ? 1 : 0;
#else
ray.dirIsNeg0 = any(ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg1 = any(ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg2 = any(ray.invDir.z < 0) ? 1 : 0;
#endif
}
__device__
inline
static bool BBoxIntersect(const float bounds[2][3],
const Ray &ray) {
float3 bounds0 = { bounds[0][0], bounds[0][1], bounds[0][2] };
float3 bounds1 = { bounds[1][0], bounds[1][1], bounds[1][2] };
float t0 = ray.mint, t1 = ray.maxt;
// Check all three axis-aligned slabs. Don't try to early out; it's
// not worth the trouble
float3 tNear = (bounds0 - ray.origin) * ray.invDir;
float3 tFar = (bounds1 - ray.origin) * ray.invDir;
if (tNear.x > tFar.x) {
float tmp = tNear.x;
tNear.x = tFar.x;
tFar.x = tmp;
}
t0 = max(tNear.x, t0);
t1 = min(tFar.x, t1);
if (tNear.y > tFar.y) {
float tmp = tNear.y;
tNear.y = tFar.y;
tFar.y = tmp;
}
t0 = max(tNear.y, t0);
t1 = min(tFar.y, t1);
if (tNear.z > tFar.z) {
float tmp = tNear.z;
tNear.z = tFar.z;
tFar.z = tmp;
}
t0 = max(tNear.z, t0);
t1 = min(tFar.z, t1);
return (t0 <= t1);
}
__device__
inline
static bool TriIntersect(const Triangle &tri, Ray &ray) {
float3 p0 = { tri.p[0][0], tri.p[0][1], tri.p[0][2] };
float3 p1 = { tri.p[1][0], tri.p[1][1], tri.p[1][2] };
float3 p2 = { tri.p[2][0], tri.p[2][1], tri.p[2][2] };
float3 e1 = p1 - p0;
float3 e2 = p2 - p0;
float3 s1 = Cross(ray.dir, e2);
float divisor = Dot(s1, e1);
bool hit = true;
if (divisor == 0.)
hit = false;
float invDivisor = 1.f / divisor;
// Compute first barycentric coordinate
float3 d = ray.origin - p0;
float b1 = Dot(d, s1) * invDivisor;
if (b1 < 0. || b1 > 1.)
hit = false;
// Compute second barycentric coordinate
float3 s2 = Cross(d, e1);
float b2 = Dot(ray.dir, s2) * invDivisor;
if (b2 < 0. || b1 + b2 > 1.)
hit = false;
// Compute _t_ to intersection point
float t = Dot(e2, s2) * invDivisor;
if (t < ray.mint || t > ray.maxt)
hit = false;
if (hit) {
ray.maxt = t;
ray.hitId = tri.id;
}
return hit;
}
__device__
inline
bool BVHIntersect(const LinearBVHNode nodes[],
const Triangle tris[], Ray &r,
int todo[]) {
Ray ray = r;
bool hit = false;
// Follow ray through BVH nodes to find primitive intersections
int todoOffset = 0, nodeNum = 0;
while (true) {
// Check ray against BVH node
LinearBVHNode node = nodes[nodeNum];
if (any(BBoxIntersect(node.bounds, ray))) {
unsigned int nPrimitives = node.nPrimitives;
if (nPrimitives > 0) {
// Intersect ray with primitives in leaf BVH node
unsigned int primitivesOffset = node.offset;
for ( unsigned int i = 0; i < nPrimitives; ++i) {
if (TriIntersect(tris[primitivesOffset+i], ray))
hit = true;
}
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
else {
// Put far BVH node on _todo_ stack, advance to near node
int dirIsNeg;
if (node.splitAxis == 0) dirIsNeg = r.dirIsNeg0;
if (node.splitAxis == 1) dirIsNeg = r.dirIsNeg1;
if (node.splitAxis == 2) dirIsNeg = r.dirIsNeg2;
if (dirIsNeg) {
todo[todoOffset++] = nodeNum + 1;
nodeNum = node.offset;
}
else {
todo[todoOffset++] = node.offset;
nodeNum = nodeNum + 1;
}
}
}
else {
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
}
r.maxt = ray.maxt;
r.hitId = ray.hitId;
return hit;
}
__device__
inline
static void raytrace_tile( int x0, int x1,
int y0, int y1,
int width, int height,
int baseWidth, int baseHeight,
const float raster2camera[4][4],
const float camera2world[4][4],
float image[], int id[],
const LinearBVHNode nodes[],
const Triangle triangles[]) {
float widthScale = (float)(baseWidth) / (float)(width);
float heightScale = (float)(baseHeight) / (float)(height);
#if 0
int * todo = new int[64];
#define ALLOC
#else
int todo[64];
#endif
for (int y = y0 ;y < y1; y++)
for (int x = x0 + programIndex; x < x1; x += programCount)
if (x < x1)
{
Ray ray;
generateRay(raster2camera, camera2world, x*widthScale,
y*heightScale, ray);
BVHIntersect(nodes, triangles, ray, todo);
int offset = y * width + x;
image[offset] = ray.maxt;
id[offset] = ray.hitId;
}
#ifdef ALLOC
delete todo;
#endif
}
__global__
void raytrace_tile_task( int width, int height,
int baseWidth, int baseHeight,
const float raster2camera[4][4],
const float camera2world[4][4],
float image[], int id[],
const LinearBVHNode nodes[],
const Triangle triangles[]) {
int dx = 64, dy = 8; // must match dx, dy below
int xBuckets = (width + (dx-1)) / dx;
int x0 = (taskIndex % xBuckets) * dx;
int x1 = min(x0 + dx, width);
int y0 = (taskIndex / xBuckets) * dy;
int y1 = min(y0 + dy, height);
raytrace_tile(x0, x1, y0, y1, width, height, baseWidth, baseHeight,
raster2camera, camera2world, image,
id, nodes, triangles);
}
extern "C" __global__ void raytrace_ispc_tasks( int width, int height,
int baseWidth, int baseHeight,
const float raster2camera[4][4],
const float camera2world[4][4],
float image[], int id[],
const LinearBVHNode nodes[],
const Triangle triangles[]) {
int dx = 64, dy = 8;
int xBuckets = (width + (dx-1)) / dx;
int yBuckets = (height + (dy-1)) / dy;
int nTasks = xBuckets * yBuckets;
if (programIndex == 0)
raytrace_tile_task<<<(nTasks+4-1)/4,128>>>(width, height, baseWidth, baseHeight,
raster2camera, camera2world,
image, id, nodes, triangles);
cudaDeviceSynchronize();
}

311
examples_ptx/rt/rt.ispc Normal file
View File

@@ -0,0 +1,311 @@
/*
Copyright (c) 2010-2011, Intel Corporation
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of Intel Corporation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#define bool int
typedef float<3> float3;
struct Ray {
float3 origin, dir, invDir;
uniform unsigned int dirIsNeg[3];
float mint, maxt;
int hitId;
};
struct Triangle {
float p[3][4];
int id;
int pad[3];
};
struct LinearBVHNode {
float bounds[2][3];
unsigned int offset; // num primitives for leaf, second child for interior
unsigned int8 nPrimitives;
unsigned int8 splitAxis;
unsigned int16 pad;
};
static inline float3 Cross(const float3 v1, const float3 v2) {
float v1x = v1.x, v1y = v1.y, v1z = v1.z;
float v2x = v2.x, v2y = v2.y, v2z = v2.z;
float3 ret;
ret.x = (v1y * v2z) - (v1z * v2y);
ret.y = (v1z * v2x) - (v1x * v2z);
ret.z = (v1x * v2y) - (v1y * v2x);
return ret;
}
static inline float Dot(const float3 a, const float3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
static void generateRay(uniform const float raster2camera[4][4],
uniform const float camera2world[4][4],
float x, float y, Ray &ray) {
ray.mint = 0.f;
ray.maxt = 1e30f;
ray.hitId = 0;
// transform raster coordinate (x, y, 0) to camera space
float camx = raster2camera[0][0] * x + raster2camera[0][1] * y + raster2camera[0][3];
float camy = raster2camera[1][0] * x + raster2camera[1][1] * y + raster2camera[1][3];
float camz = raster2camera[2][3];
float camw = raster2camera[3][3];
camx /= camw;
camy /= camw;
camz /= camw;
ray.dir.x = camera2world[0][0] * camx + camera2world[0][1] * camy +
camera2world[0][2] * camz;
ray.dir.y = camera2world[1][0] * camx + camera2world[1][1] * camy +
camera2world[1][2] * camz;
ray.dir.z = camera2world[2][0] * camx + camera2world[2][1] * camy +
camera2world[2][2] * camz;
ray.origin.x = camera2world[0][3] / camera2world[3][3];
ray.origin.y = camera2world[1][3] / camera2world[3][3];
ray.origin.z = camera2world[2][3] / camera2world[3][3];
ray.invDir = 1.f / ray.dir;
ray.dirIsNeg[0] = any(ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg[1] = any(ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg[2] = any(ray.invDir.z < 0) ? 1 : 0;
}
static bool BBoxIntersect(const uniform float bounds[2][3],
const Ray &ray) {
uniform float3 bounds0 = { bounds[0][0], bounds[0][1], bounds[0][2] };
uniform float3 bounds1 = { bounds[1][0], bounds[1][1], bounds[1][2] };
float t0 = ray.mint, t1 = ray.maxt;
// Check all three axis-aligned slabs. Don't try to early out; it's
// not worth the trouble
float3 tNear = (bounds0 - ray.origin) * ray.invDir;
float3 tFar = (bounds1 - ray.origin) * ray.invDir;
if (tNear.x > tFar.x) {
float tmp = tNear.x;
tNear.x = tFar.x;
tFar.x = tmp;
}
t0 = max(tNear.x, t0);
t1 = min(tFar.x, t1);
if (tNear.y > tFar.y) {
float tmp = tNear.y;
tNear.y = tFar.y;
tFar.y = tmp;
}
t0 = max(tNear.y, t0);
t1 = min(tFar.y, t1);
if (tNear.z > tFar.z) {
float tmp = tNear.z;
tNear.z = tFar.z;
tFar.z = tmp;
}
t0 = max(tNear.z, t0);
t1 = min(tFar.z, t1);
return (t0 <= t1);
}
static bool TriIntersect(const uniform Triangle &tri, Ray &ray) {
uniform float3 p0 = { tri.p[0][0], tri.p[0][1], tri.p[0][2] };
uniform float3 p1 = { tri.p[1][0], tri.p[1][1], tri.p[1][2] };
uniform float3 p2 = { tri.p[2][0], tri.p[2][1], tri.p[2][2] };
uniform float3 e1 = p1 - p0;
uniform float3 e2 = p2 - p0;
float3 s1 = Cross(ray.dir, e2);
float divisor = Dot(s1, e1);
bool hit = true;
if (divisor == 0.)
hit = false;
float invDivisor = 1.f / divisor;
// Compute first barycentric coordinate
float3 d = ray.origin - p0;
float b1 = Dot(d, s1) * invDivisor;
if (b1 < 0. || b1 > 1.)
hit = false;
// Compute second barycentric coordinate
float3 s2 = Cross(d, e1);
float b2 = Dot(ray.dir, s2) * invDivisor;
if (b2 < 0. || b1 + b2 > 1.)
hit = false;
// Compute _t_ to intersection point
float t = Dot(e2, s2) * invDivisor;
if (t < ray.mint || t > ray.maxt)
hit = false;
if (hit) {
ray.maxt = t;
ray.hitId = tri.id;
}
return hit;
}
bool BVHIntersect(const uniform LinearBVHNode nodes[],
const uniform Triangle tris[], Ray &r) {
Ray ray = r;
bool hit = false;
// Follow ray through BVH nodes to find primitive intersections
uniform int todoOffset = 0, nodeNum = 0;
uniform int todo[64];
while (true) {
// Check ray against BVH node
uniform LinearBVHNode node = nodes[nodeNum];
if (any(BBoxIntersect(node.bounds, ray))) {
uniform unsigned int nPrimitives = node.nPrimitives;
if (nPrimitives > 0) {
// Intersect ray with primitives in leaf BVH node
uniform unsigned int primitivesOffset = node.offset;
for (uniform unsigned int i = 0; i < nPrimitives; ++i) {
if (TriIntersect(tris[primitivesOffset+i], ray))
hit = true;
}
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
else {
// Put far BVH node on _todo_ stack, advance to near node
if (r.dirIsNeg[node.splitAxis]) {
todo[todoOffset++] = nodeNum + 1;
nodeNum = node.offset;
}
else {
todo[todoOffset++] = node.offset;
nodeNum = nodeNum + 1;
}
}
}
else {
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
}
r.maxt = ray.maxt;
r.hitId = ray.hitId;
return hit;
}
static void raytrace_tile(uniform int x0, uniform int x1,
uniform int y0, uniform int y1,
uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform float widthScale = (float)(baseWidth) / (float)(width);
uniform float heightScale = (float)(baseHeight) / (float)(height);
foreach_tiled (y = y0 ... y1, x = x0 ... x1) {
Ray ray;
generateRay(raster2camera, camera2world, x*widthScale,
y*heightScale, ray);
BVHIntersect(nodes, triangles, ray);
int offset = y * width + x;
image[offset] = ray.maxt;
id[offset] = ray.hitId;
}
}
export void raytrace_ispc(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
raytrace_tile(0, width, 0, height, width, height, baseWidth, baseHeight,
raster2camera, camera2world, image,
id, nodes, triangles);
}
task void raytrace_tile_task(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform int dx = 64, dy = 8; // must match dx, dy below
uniform int xBuckets = (width + (dx-1)) / dx;
uniform int x0 = (taskIndex % xBuckets) * dx;
uniform int x1 = min(x0 + dx, width);
uniform int y0 = (taskIndex / xBuckets) * dy;
uniform int y1 = min(y0 + dy, height);
raytrace_tile(x0, x1, y0, y1, width, height, baseWidth, baseHeight,
raster2camera, camera2world, image,
id, nodes, triangles);
}
export void raytrace_ispc_tasks(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform int dx = 64, dy = 8;
uniform int xBuckets = (width + (dx-1)) / dx;
uniform int yBuckets = (height + (dy-1)) / dy;
uniform int nTasks = xBuckets * yBuckets;
launch[nTasks] raytrace_tile_task(width, height, baseWidth, baseHeight,
raster2camera, camera2world,
image, id, nodes, triangles);
}

34
examples_ptx/rt/rt.vcxproj Normal file
View File

@@ -0,0 +1,34 @@
<?xml version="1.0" encoding="utf-8"?>
<Project DefaultTargets="Build" ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<ItemGroup Label="ProjectConfigurations">
<ProjectConfiguration Include="Debug|Win32">
<Configuration>Debug</Configuration>
<Platform>Win32</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Debug|x64">
<Configuration>Debug</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Release|Win32">
<Configuration>Release</Configuration>
<Platform>Win32</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Release|x64">
<Configuration>Release</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
</ItemGroup>
<PropertyGroup Label="Globals">
<ProjectGuid>{E787BC3F-2D2E-425E-A64D-4721E2FF3DC9}</ProjectGuid>
<Keyword>Win32Proj</Keyword>
<RootNamespace>rt</RootNamespace>
<ISPC_file>rt</ISPC_file>
<default_targets>sse2,sse4-x2,avx1-i32x8</default_targets>
</PropertyGroup>
<Import Project="..\common.props" />
<ItemGroup>
<ClCompile Include="rt.cpp" />
<ClCompile Include="rt_serial.cpp" />
<ClCompile Include="../tasksys.cpp" />
</ItemGroup>
</Project>

334
examples_ptx/rt/rt1.ispc Normal file
View File

@@ -0,0 +1,334 @@
/*
Copyright (c) 2010-2011, Intel Corporation
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of Intel Corporation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#define bool int
typedef float<3> float3;
struct Ray {
float3 origin, dir, invDir;
//uniform unsigned int dirIsNeg[3];
unsigned int dirIsNeg0, dirIsNeg1, dirIsNeg2;
float mint, maxt;
int hitId;
};
struct Triangle {
float p[3][4];
int id;
int pad[3];
};
struct LinearBVHNode {
float bounds[2][3];
unsigned int offset; // num primitives for leaf, second child for interior
unsigned int8 nPrimitives;
unsigned int8 splitAxis;
unsigned int16 pad;
};
static inline float3 Cross(const float3 v1, const float3 v2) {
float v1x = v1.x, v1y = v1.y, v1z = v1.z;
float v2x = v2.x, v2y = v2.y, v2z = v2.z;
float3 ret;
ret.x = (v1y * v2z) - (v1z * v2y);
ret.y = (v1z * v2x) - (v1x * v2z);
ret.z = (v1x * v2y) - (v1y * v2x);
return ret;
}
static inline float Dot(const float3 a, const float3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
inline
static void generateRay(uniform const float raster2camera[4][4],
uniform const float camera2world[4][4],
float x, float y, Ray &ray) {
ray.mint = 0.f;
ray.maxt = 1e30f;
ray.hitId = 0;
// transform raster coordinate (x, y, 0) to camera space
float camx = raster2camera[0][0] * x + raster2camera[0][1] * y + raster2camera[0][3];
float camy = raster2camera[1][0] * x + raster2camera[1][1] * y + raster2camera[1][3];
float camz = raster2camera[2][3];
float camw = raster2camera[3][3];
camx /= camw;
camy /= camw;
camz /= camw;
ray.dir.x = camera2world[0][0] * camx + camera2world[0][1] * camy +
camera2world[0][2] * camz;
ray.dir.y = camera2world[1][0] * camx + camera2world[1][1] * camy +
camera2world[1][2] * camz;
ray.dir.z = camera2world[2][0] * camx + camera2world[2][1] * camy +
camera2world[2][2] * camz;
ray.origin.x = camera2world[0][3] / camera2world[3][3];
ray.origin.y = camera2world[1][3] / camera2world[3][3];
ray.origin.z = camera2world[2][3] / camera2world[3][3];
ray.invDir = 1.f / ray.dir;
#if 0
ray.dirIsNeg[0] = any(ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg[1] = any(ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg[2] = any(ray.invDir.z < 0) ? 1 : 0;
#else
ray.dirIsNeg0 = any(ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg1 = any(ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg2 = any(ray.invDir.z < 0) ? 1 : 0;
#endif
}
inline
static bool BBoxIntersect(const uniform float bounds[2][3],
const Ray &ray) {
uniform float3 bounds0 = { bounds[0][0], bounds[0][1], bounds[0][2] };
uniform float3 bounds1 = { bounds[1][0], bounds[1][1], bounds[1][2] };
float t0 = ray.mint, t1 = ray.maxt;
// Check all three axis-aligned slabs. Don't try to early out; it's
// not worth the trouble
float3 tNear = (bounds0 - ray.origin) * ray.invDir;
float3 tFar = (bounds1 - ray.origin) * ray.invDir;
if (tNear.x > tFar.x) {
float tmp = tNear.x;
tNear.x = tFar.x;
tFar.x = tmp;
}
t0 = max(tNear.x, t0);
t1 = min(tFar.x, t1);
if (tNear.y > tFar.y) {
float tmp = tNear.y;
tNear.y = tFar.y;
tFar.y = tmp;
}
t0 = max(tNear.y, t0);
t1 = min(tFar.y, t1);
if (tNear.z > tFar.z) {
float tmp = tNear.z;
tNear.z = tFar.z;
tFar.z = tmp;
}
t0 = max(tNear.z, t0);
t1 = min(tFar.z, t1);
return (t0 <= t1);
}
inline
static bool TriIntersect(const uniform Triangle &tri, Ray &ray) {
uniform float3 p0 = { tri.p[0][0], tri.p[0][1], tri.p[0][2] };
uniform float3 p1 = { tri.p[1][0], tri.p[1][1], tri.p[1][2] };
uniform float3 p2 = { tri.p[2][0], tri.p[2][1], tri.p[2][2] };
uniform float3 e1 = p1 - p0;
uniform float3 e2 = p2 - p0;
float3 s1 = Cross(ray.dir, e2);
float divisor = Dot(s1, e1);
bool hit = true;
if (divisor == 0.)
hit = false;
float invDivisor = 1.f / divisor;
// Compute first barycentric coordinate
float3 d = ray.origin - p0;
float b1 = Dot(d, s1) * invDivisor;
if (b1 < 0. || b1 > 1.)
hit = false;
// Compute second barycentric coordinate
float3 s2 = Cross(d, e1);
float b2 = Dot(ray.dir, s2) * invDivisor;
if (b2 < 0. || b1 + b2 > 1.)
hit = false;
// Compute _t_ to intersection point
float t = Dot(e2, s2) * invDivisor;
if (t < ray.mint || t > ray.maxt)
hit = false;
if (hit) {
ray.maxt = t;
ray.hitId = tri.id;
}
return hit;
}
inline
bool BVHIntersect(const uniform LinearBVHNode nodes[],
const uniform Triangle tris[], Ray &r,
uniform int todo[]) {
Ray ray = r;
bool hit = false;
// Follow ray through BVH nodes to find primitive intersections
uniform int todoOffset = 0, nodeNum = 0;
while (true) {
// Check ray against BVH node
uniform LinearBVHNode node = nodes[nodeNum];
if (any(BBoxIntersect(node.bounds, ray))) {
uniform unsigned int nPrimitives = node.nPrimitives;
if (nPrimitives > 0) {
// Intersect ray with primitives in leaf BVH node
uniform unsigned int primitivesOffset = node.offset;
for (uniform unsigned int i = 0; i < nPrimitives; ++i) {
if (TriIntersect(tris[primitivesOffset+i], ray))
hit = true;
}
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
else {
// Put far BVH node on _todo_ stack, advance to near node
int dirIsNeg;
if (node.splitAxis == 0) dirIsNeg = r.dirIsNeg0;
if (node.splitAxis == 1) dirIsNeg = r.dirIsNeg1;
if (node.splitAxis == 2) dirIsNeg = r.dirIsNeg2;
if (dirIsNeg) {
todo[todoOffset++] = nodeNum + 1;
nodeNum = node.offset;
}
else {
todo[todoOffset++] = node.offset;
nodeNum = nodeNum + 1;
}
}
}
else {
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
}
r.maxt = ray.maxt;
r.hitId = ray.hitId;
return hit;
}
inline
static void raytrace_tile(uniform int x0, uniform int x1,
uniform int y0, uniform int y1,
uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform float widthScale = (float)(baseWidth) / (float)(width);
uniform float heightScale = (float)(baseHeight) / (float)(height);
#if 0
uniform int * uniform todo = uniform new uniform int[64];
#define ALLOC
#else
uniform int todo[64];
#endif
foreach_tiled (y = y0 ... y1, x = x0 ... x1) {
Ray ray;
generateRay(raster2camera, camera2world, x*widthScale,
y*heightScale, ray);
BVHIntersect(nodes, triangles, ray, todo);
int offset = y * width + x;
image[offset] = ray.maxt;
id[offset] = ray.hitId;
}
#ifdef ALLOC
delete todo;
#endif
}
export void raytrace_ispc(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
raytrace_tile(0, width, 0, height, width, height, baseWidth, baseHeight,
raster2camera, camera2world, image,
id, nodes, triangles);
}
task void raytrace_tile_task(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform int dx = 64, dy = 8; // must match dx, dy below
uniform int xBuckets = (width + (dx-1)) / dx;
uniform int x0 = (taskIndex % xBuckets) * dx;
uniform int x1 = min(x0 + dx, width);
uniform int y0 = (taskIndex / xBuckets) * dy;
uniform int y1 = min(y0 + dy, height);
raytrace_tile(x0, x1, y0, y1, width, height, baseWidth, baseHeight,
raster2camera, camera2world, image,
id, nodes, triangles);
}
export void raytrace_ispc_tasks(uniform int width, uniform int height,
uniform int baseWidth, uniform int baseHeight,
const uniform float raster2camera[4][4],
const uniform float camera2world[4][4],
uniform float image[], uniform int id[],
const uniform LinearBVHNode nodes[],
const uniform Triangle triangles[]) {
uniform int dx = 64, dy = 8;
uniform int xBuckets = (width + (dx-1)) / dx;
uniform int yBuckets = (height + (dy-1)) / dy;
uniform int nTasks = xBuckets * yBuckets;
launch[nTasks] raytrace_tile_task(width, height, baseWidth, baseHeight,
raster2camera, camera2world,
image, id, nodes, triangles);
sync;
}

292
examples_ptx/rt/rt_serial.cpp Normal file
View File

@@ -0,0 +1,292 @@
/*
Copyright (c) 2010-2011, Intel Corporation
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of Intel Corporation nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#define NOMINMAX
#pragma warning (disable: 4244)
#pragma warning (disable: 4305)
#endif
#include <algorithm>
#include <stdint.h>
// Just enough of a float3 class to do what we need in this file.
#ifdef _MSC_VER
__declspec(align(16))
#endif
struct float3 {
float3() { }
float3(float xx, float yy, float zz) { x = xx; y = yy; z = zz; }
float3 operator*(float f) const { return float3(x*f, y*f, z*f); }
float3 operator-(const float3 &f2) const {
return float3(x-f2.x, y-f2.y, z-f2.z);
}
float3 operator*(const float3 &f2) const {
return float3(x*f2.x, y*f2.y, z*f2.z);
}
float x, y, z;
float pad; // match padding/alignment of ispc version
}
#ifndef _MSC_VER
__attribute__ ((aligned(16)))
#endif
;
struct Ray {
float3 origin, dir, invDir;
unsigned int dirIsNeg[3];
float mint, maxt;
int hitId;
};
// Declare these in a namespace so the mangling matches
namespace ispc {
struct Triangle {
float p[3][4]; // extra float pad after each vertex
int32_t id;
int32_t pad[3]; // make 16 x 32-bits
};
struct LinearBVHNode {
float bounds[2][3];
int32_t offset; // primitives for leaf, second child for interior
uint8_t nPrimitives;
uint8_t splitAxis;
uint16_t pad;
};
}
using namespace ispc;
inline float3 Cross(const float3 &v1, const float3 &v2) {
float v1x = v1.x, v1y = v1.y, v1z = v1.z;
float v2x = v2.x, v2y = v2.y, v2z = v2.z;
float3 ret;
ret.x = (v1y * v2z) - (v1z * v2y);
ret.y = (v1z * v2x) - (v1x * v2z);
ret.z = (v1x * v2y) - (v1y * v2x);
return ret;
}
inline float Dot(const float3 &a, const float3 &b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
static void generateRay(const float raster2camera[4][4],
const float camera2world[4][4],
float x, float y, Ray &ray) {
ray.mint = 0.f;
ray.maxt = 1e30f;
ray.hitId = 0;
// transform raster coordinate (x, y, 0) to camera space
float camx = raster2camera[0][0] * x + raster2camera[0][1] * y + raster2camera[0][3];
float camy = raster2camera[1][0] * x + raster2camera[1][1] * y + raster2camera[1][3];
float camz = raster2camera[2][3];
float camw = raster2camera[3][3];
camx /= camw;
camy /= camw;
camz /= camw;
ray.dir.x = camera2world[0][0] * camx + camera2world[0][1] * camy +
camera2world[0][2] * camz;
ray.dir.y = camera2world[1][0] * camx + camera2world[1][1] * camy +
camera2world[1][2] * camz;
ray.dir.z = camera2world[2][0] * camx + camera2world[2][1] * camy +
camera2world[2][2] * camz;
ray.origin.x = camera2world[0][3] / camera2world[3][3];
ray.origin.y = camera2world[1][3] / camera2world[3][3];
ray.origin.z = camera2world[2][3] / camera2world[3][3];
ray.invDir.x = 1.f / ray.dir.x;
ray.invDir.y = 1.f / ray.dir.y;
ray.invDir.z = 1.f / ray.dir.z;
ray.dirIsNeg[0] = (ray.invDir.x < 0) ? 1 : 0;
ray.dirIsNeg[1] = (ray.invDir.y < 0) ? 1 : 0;
ray.dirIsNeg[2] = (ray.invDir.z < 0) ? 1 : 0;
}
static inline bool BBoxIntersect(const float bounds[2][3],
const Ray &ray) {
float3 bounds0(bounds[0][0], bounds[0][1], bounds[0][2]);
float3 bounds1(bounds[1][0], bounds[1][1], bounds[1][2]);
float t0 = ray.mint, t1 = ray.maxt;
float3 tNear = (bounds0 - ray.origin) * ray.invDir;
float3 tFar = (bounds1 - ray.origin) * ray.invDir;
if (tNear.x > tFar.x) {
float tmp = tNear.x;
tNear.x = tFar.x;
tFar.x = tmp;
}
t0 = std::max(tNear.x, t0);
t1 = std::min(tFar.x, t1);
if (tNear.y > tFar.y) {
float tmp = tNear.y;
tNear.y = tFar.y;
tFar.y = tmp;
}
t0 = std::max(tNear.y, t0);
t1 = std::min(tFar.y, t1);
if (tNear.z > tFar.z) {
float tmp = tNear.z;
tNear.z = tFar.z;
tFar.z = tmp;
}
t0 = std::max(tNear.z, t0);
t1 = std::min(tFar.z, t1);
return (t0 <= t1);
}
inline bool TriIntersect(const Triangle &tri, Ray &ray) {
float3 p0(tri.p[0][0], tri.p[0][1], tri.p[0][2]);
float3 p1(tri.p[1][0], tri.p[1][1], tri.p[1][2]);
float3 p2(tri.p[2][0], tri.p[2][1], tri.p[2][2]);
float3 e1 = p1 - p0;
float3 e2 = p2 - p0;
float3 s1 = Cross(ray.dir, e2);
float divisor = Dot(s1, e1);
if (divisor == 0.)
return false;
float invDivisor = 1.f / divisor;
// Compute first barycentric coordinate
float3 d = ray.origin - p0;
float b1 = Dot(d, s1) * invDivisor;
if (b1 < 0. || b1 > 1.)
return false;
// Compute second barycentric coordinate
float3 s2 = Cross(d, e1);
float b2 = Dot(ray.dir, s2) * invDivisor;
if (b2 < 0. || b1 + b2 > 1.)
return false;
// Compute _t_ to intersection point
float t = Dot(e2, s2) * invDivisor;
if (t < ray.mint || t > ray.maxt)
return false;
ray.maxt = t;
ray.hitId = tri.id;
return true;
}
bool BVHIntersect(const LinearBVHNode nodes[], const Triangle tris[],
Ray &r) {
Ray ray = r;
bool hit = false;
// Follow ray through BVH nodes to find primitive intersections
int todoOffset = 0, nodeNum = 0;
int todo[64];
while (true) {
// Check ray against BVH node
const LinearBVHNode &node = nodes[nodeNum];
if (BBoxIntersect(node.bounds, ray)) {
unsigned int nPrimitives = node.nPrimitives;
if (nPrimitives > 0) {
// Intersect ray with primitives in leaf BVH node
unsigned int primitivesOffset = node.offset;
for (unsigned int i = 0; i < nPrimitives; ++i) {
if (TriIntersect(tris[primitivesOffset+i], ray))
hit = true;
}
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
else {
// Put far BVH node on _todo_ stack, advance to near node
if (r.dirIsNeg[node.splitAxis]) {
todo[todoOffset++] = nodeNum + 1;
nodeNum = node.offset;
}
else {
todo[todoOffset++] = node.offset;
nodeNum = nodeNum + 1;
}
}
}
else {
if (todoOffset == 0)
break;
nodeNum = todo[--todoOffset];
}
}
r.maxt = ray.maxt;
r.hitId = ray.hitId;
return hit;
}
void raytrace_serial(int width, int height, int baseWidth, int baseHeight,
const float raster2camera[4][4],
const float camera2world[4][4],
float image[],
int id[],
const LinearBVHNode nodes[],
const Triangle triangles[]) {
float widthScale = float(baseWidth) / float(width);
float heightScale = float(baseHeight) / float(height);
for (int y = 0; y < height; ++y) {
for (int x = 0; x < width; ++x) {
Ray ray;
generateRay(raster2camera, camera2world, x * widthScale,
y * heightScale, ray);
BVHIntersect(nodes, triangles, ray);
int offset = y * width + x;
image[offset] = ray.maxt;
id[offset] = ray.hitId;
}
}
}

BIN
examples_ptx/rt/sponza.bvh Normal file
View File

Binary file not shown.

BIN
examples_ptx/rt/sponza.camera Normal file
View File

Binary file not shown.

BIN
examples_ptx/rt/teapot.bvh Normal file
View File

Binary file not shown.

BIN
examples_ptx/rt/teapot.camera Normal file
View File

Binary file not shown.