diff --git a/examples_ptx/rt/.gitignore b/examples_ptx/rt/.gitignore new file mode 100644 index 00000000..5a95423b --- /dev/null +++ b/examples_ptx/rt/.gitignore @@ -0,0 +1,2 @@ +rt +*.ppm diff --git a/examples_ptx/rt/Makefile_cpu b/examples_ptx/rt/Makefile_cpu new file mode 100644 index 00000000..0c72f104 --- /dev/null +++ b/examples_ptx/rt/Makefile_cpu @@ -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 diff --git a/examples_ptx/rt/cornell.bvh b/examples_ptx/rt/cornell.bvh new file mode 100644 index 00000000..f7e0f3dd Binary files /dev/null and b/examples_ptx/rt/cornell.bvh differ diff --git a/examples_ptx/rt/cornell.camera b/examples_ptx/rt/cornell.camera new file mode 100644 index 00000000..0fec1642 Binary files /dev/null and b/examples_ptx/rt/cornell.camera differ diff --git a/examples_ptx/rt/rt.cpp b/examples_ptx/rt/rt.cpp new file mode 100644 index 00000000..ea1594c5 --- /dev/null +++ b/examples_ptx/rt/rt.cpp @@ -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 +#include +#include +#include +#include +#include +#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 [--scale=] [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; +} diff --git a/examples_ptx/rt/rt.cu b/examples_ptx/rt/rt.cu new file mode 100644 index 00000000..8decd03a --- /dev/null +++ b/examples_ptx/rt/rt.cu @@ -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(); +} + diff --git a/examples_ptx/rt/rt.ispc b/examples_ptx/rt/rt.ispc new file mode 100644 index 00000000..dd655085 --- /dev/null +++ b/examples_ptx/rt/rt.ispc @@ -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); +} + diff --git a/examples_ptx/rt/rt.vcxproj b/examples_ptx/rt/rt.vcxproj new file mode 100644 index 00000000..00b6dd3a --- /dev/null +++ b/examples_ptx/rt/rt.vcxproj @@ -0,0 +1,34 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + {E787BC3F-2D2E-425E-A64D-4721E2FF3DC9} + Win32Proj + rt + rt + sse2,sse4-x2,avx1-i32x8 + + + + + + + + diff --git a/examples_ptx/rt/rt1.ispc b/examples_ptx/rt/rt1.ispc new file mode 100644 index 00000000..4461bdb9 --- /dev/null +++ b/examples_ptx/rt/rt1.ispc @@ -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; +} + diff --git a/examples_ptx/rt/rt_serial.cpp b/examples_ptx/rt/rt_serial.cpp new file mode 100644 index 00000000..535f25e4 --- /dev/null +++ b/examples_ptx/rt/rt_serial.cpp @@ -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 +#include + +// 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; + } + } +} diff --git a/examples_ptx/rt/sponza.bvh b/examples_ptx/rt/sponza.bvh new file mode 100644 index 00000000..e59bde24 Binary files /dev/null and b/examples_ptx/rt/sponza.bvh differ diff --git a/examples_ptx/rt/sponza.camera b/examples_ptx/rt/sponza.camera new file mode 100644 index 00000000..7d44ec23 Binary files /dev/null and b/examples_ptx/rt/sponza.camera differ diff --git a/examples_ptx/rt/teapot.bvh b/examples_ptx/rt/teapot.bvh new file mode 100644 index 00000000..efcd7807 Binary files /dev/null and b/examples_ptx/rt/teapot.bvh differ diff --git a/examples_ptx/rt/teapot.camera b/examples_ptx/rt/teapot.camera new file mode 100644 index 00000000..9a98e3f6 Binary files /dev/null and b/examples_ptx/rt/teapot.camera differ