+added CUDA code

This commit is contained in:
Evghenii
2013-11-09 21:05:28 +01:00
parent 946530019a
commit dbd0581cb3

347
examples_cuda/rt/rt.cu Normal file
View File

@@ -0,0 +1,347 @@
#define programCount 32
#define programIndex (threadIdx.x & 31)
#define taskIndex (blockIdx.x*4 + (threadIdx.x >> 5))
#define taskCount (gridDim.x*4)
#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;
}
};
#if 0
#define DIRISNEG
#endif
struct Ray {
float3 origin, dir, invDir;
#ifdef DIRISNEG /* this fails to compile with nvvm */
unsigned int dirIsNeg[3];
#else
unsigned int dirIsNeg0, dirIsNeg1, dirIsNeg2;
#endif
float mint, maxt;
int hitId;
};
struct Triangle {
float p[3][4];
int id;
int pad[3];
};
#define int8 char
#define int16 short
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__
static inline 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;
#ifdef DIRISNEG
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__
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;
// 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__
static 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);
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__
static inline 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
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
#ifdef DIRISNEG
const int dirIsNeg = r.dirIsNeg[node.splitAxis];
#else
int dirIsNeg;
if (node.splitAxis == 0) dirIsNeg = r.dirIsNeg0;
if (node.splitAxis == 1) dirIsNeg = r.dirIsNeg1;
if (node.splitAxis == 2) dirIsNeg = r.dirIsNeg2;
#endif
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__
static inline 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);
// foreach_tiled (y = y0 ... y1, x = x0 ... x1)
for ( int y = y0; y < y1; y++)
for ( int xb = x0; xb < x1; xb += programCount)
{
const int x = xb + programIndex;
Ray ray;
generateRay(raster2camera, camera2world, x*widthScale,
y*heightScale, ray);
BVHIntersect(nodes, triangles, ray);
int offset = y * width + x;
if (x < x1)
{
image[offset] = ray.maxt;
id[offset] = ray.hitId;
}
}
}
extern "C"
__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[]) {
if (taskIndex >= taskCount) return;
int dx = 32, dy = 16; // 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);
}