mirror/OpenSpace
1
0
Fork 0
mirror of https://github.com/OpenSpace/OpenSpace.git synced 2026-02-26 14:58:51 -06:00
Code Issues Packages Projects Releases Wiki Activity

Remove old OpenCL kernels

Browse Source
This commit is contained in:
Alexander Bock
2016-05-01 12:38:16 +02:00
parent 8043963c20
commit f3a5b976ad
8 changed files with 0 additions and 1798 deletions
Download Patch File Download Diff File Expand all files Collapse all files

202
kernels/Raycaster.cl
View File

@@ -1,202 +0,0 @@
// Needs to mirror struct on host
struct KernelConstants {
float stepSize;
float intensity;
int aDim;
int bDim;
int cDim;
};
// Turn normalized [0..1] coordinates into array index
int CoordsToIndex(float3 _coordinates,
int3 _dimensions) {
// Put coords in [0 .. dim-1] range
int x = (float)(_dimensions.x-1) * _coordinates.x;
int y = (float)(_dimensions.y-1) * _coordinates.y;
int z = (float)(_dimensions.z-1) * _coordinates.z;
// Return index
return x + y*_dimensions.x + z*_dimensions.x*_dimensions.y;
}
// Linearly interpolate between two values. Distance
// is assumed to be normalized.
float Lerp(float _v0, float _v1, float _d) {
return _v0*(1.0 - _d) + _v1*_d;
}
// Sample a volume given spherical integer coordinates
float Sample(__global __read_only float *_data,
int3 _coords,
int3 _dims) {
int r = (_coords.x < _dims.x) ? _coords.x : _dims.x-1;
int t = (_coords.y < _dims.y) ? _coords.y : _dims.y-1;
int p = (_coords.z < _dims.z) ? _coords.z : _dims.z-1;
int idx = r + t*_dims.x + p*_dims.x*_dims.y;
return _data[idx];
}
// Trilinear interpolation and sampling
float Trilerp(__global __read_only float *_data,
float3 _spherical,
int3 _dims) {
// TODO fix wrap issue
// Get coordinates in [0..dim-1] range
float r = (float)(_dims.x-1) * _spherical.x;
float t = (float)(_dims.y-1) * _spherical.y;
float p = (float)(_dims.z-1) * _spherical.z;
// Lower values
int r0 = (int)floor(r);
int t0 = (int)floor(t);
int p0 = (int)floor(p);
// Upper values
int r1 = r0+1;
int t1 = t0+1;
int p1 = p0+1;
// Interpolation values
float dr = r - floor(r);
float dt = t - floor(t);
float dp = p - floor(p);
// Sample corners
float v000 = Sample(_data, (int3)(r0, t0, p0), _dims);
float v100 = Sample(_data, (int3)(r1, t0, p0), _dims);
float v010 = Sample(_data, (int3)(r0, t1, p0), _dims);
float v110 = Sample(_data, (int3)(r1, t1, p0), _dims);
float v001 = Sample(_data, (int3)(r0, t0, p1), _dims);
float v101 = Sample(_data, (int3)(r1, t0, p1), _dims);
float v011 = Sample(_data, (int3)(r0, t1, p1), _dims);
float v111 = Sample(_data, (int3)(r1, t1, p1), _dims);
// Interpolate
float v00 = Lerp(v000, v100, dr);
float v10 = Lerp(v010, v110, dr);
float v01 = Lerp(v001, v101, dr);
float v11 = Lerp(v011, v111, dr);
float v0 = Lerp(v00, v10, dt);
float v1 = Lerp(v01, v11, dt);
float v = Lerp(v0, v1, dp);
return min(v, 1.0);
}
// Turn normalized [0..1] cartesian coordinates
// to normalized spherical [0..1] coordinates
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acos(_cartesian.z/r) / M_PI;
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / sqrt(3.0);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
float4 TransferFunction(__global __read_only float *_tf, float _i) {
// TODO remove 512 hard-coded value and change to 1D texture
int i0 = (int)floor(511.0*_i);
int i1 = (i0 < 511) ? i0+1 : i0;
float di = _i - floor(_i);
float tfr0 = _tf[i0*4+0];
float tfr1 = _tf[i1*4+0];
float tfg0 = _tf[i0*4+1];
float tfg1 = _tf[i1*4+1];
float tfb0 = _tf[i0*4+2];
float tfb1 = _tf[i1*4+2];
float tfa0 = _tf[i0*4+3];
float tfa1 = _tf[i1*4+3];
float tfr = Lerp(tfr0, tfr1, di);
float tfg = Lerp(tfg0, tfg1, di);
float tfb = Lerp(tfb0, tfb1, di);
float tfa = Lerp(tfa0, tfa1, di);
return (float4)(tfr, tfg, tfb, tfa);
}
__kernel void
Raycaster(__global __read_only image2d_t _cubeFront,
__global __read_only image2d_t _cubeBack,
__global __write_only image2d_t _output,
__global __read_only float * _voxelData,
__constant struct KernelConstants *_constants,
__global __read_only float *_tf) {
int3 dimensions = (int3)(_constants->aDim,
_constants->bDim,
_constants->cDim);
// Kernel should be launched in 2D with one work item per pixel
int idx = get_global_id(0);
int idy = get_global_id(1);
int2 intCoords = (int2)(idx, idy);
// Sampler for texture reading
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE;
const sampler_t volumeSampler = CLK_FILTER_LINEAR |
CLK_NORMALIZED_COORDS_TRUE |
CLK_ADDRESS_CLAMP_TO_EDGE;
// Read from textures
float4 cubeFrontColor = read_imagef(_cubeFront, sampler, intCoords);
float4 cubeBackColor = read_imagef(_cubeBack, sampler, intCoords);
// Figure out the direction
float3 direction = (cubeBackColor-cubeFrontColor).xyz;
float maxDistance = length(direction);
direction = normalize(direction);
// Keep track of distance traversed
float traversed = 0.0;
// Sum colors
float stepSize = _constants->stepSize;
float3 samplePoint = cubeFrontColor.xyz;
float3 spherical;
float4 associatedColor;
float4 color = (float4)(0.0, 0.0, 0.0, 0.0);
while (traversed < maxDistance) {
spherical = CartesianToSpherical(samplePoint);
//int index = timestepOffset + CoordsToIndex(spherical, dimensions);
// Get intensity from data
//float i = _voxelData[index];
float i = Trilerp(_voxelData, spherical, dimensions);
// Texture stores intensity in R channel
//float i = read_imagef(_voxelData, volumeSampler, spherical).x;
// Front-to-back compositing
float4 tf = TransferFunction(_tf, i);
color += (1.0 - color.w)*tf;
// color += (float4)(i, i, i, 1.0);
samplePoint += direction * stepSize;
traversed += stepSize;
}
// Output
float intensity = _constants->intensity;
color *= intensity*stepSize;
// Write to image
write_imagef(_output, intCoords, color);
}

297
kernels/RaycasterBricks.cl
View File

@@ -1,297 +0,0 @@
// Needs to mirror struct on host
struct KernelConstants {
float stepSize;
float intensity;
int numBoxesPerAxis;
};
// Linearly interpolate between two values. Distance
// is assumed to be normalized.
float Lerp(float _v0, float _v1, float _d) {
return _v0*(1.0 - _d) + _v1*_d;
}
// Turn normalized [0..1] cartesian coordinates
// to normalized spherical [0..1] coordinates
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acospi(_cartesian.z/r);
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / native_sqrt(3.0);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
float4 TransferFunction(__global __read_only float *_tf, float _i) {
// TODO remove hard-coded value and change to 1D texture
int i0 = (int)floor(1023.0*_i);
int i1 = (i0 < 1023) ? i0+1 : i0;
float di = _i - floor(_i);
float tfr0 = _tf[i0*4+0];
float tfr1 = _tf[i1*4+0];
float tfg0 = _tf[i0*4+1];
float tfg1 = _tf[i1*4+1];
float tfb0 = _tf[i0*4+2];
float tfb1 = _tf[i1*4+2];
float tfa0 = _tf[i0*4+3];
float tfa1 = _tf[i1*4+3];
float tfr = Lerp(tfr0, tfr1, di);
float tfg = Lerp(tfg0, tfg1, di);
float tfb = Lerp(tfb0, tfb1, di);
float tfa = Lerp(tfa0, tfa1, di);
return (float4)(tfr, tfg, tfb, tfa);
}
// Translates a global volume coordinate [0..1] to a box coordinate,
// given a number of boxes that fit along each axis
int3 BoxCoords(float3 _globalCoords, int _boxesPerAxis) {
int3 boxCoords = convert_int3(floor(_globalCoords * (float)_boxesPerAxis));
return clamp(boxCoords, (int3)(0, 0, 0), (int3)(_boxesPerAxis-1));
}
// Calculate global volume coordinates for the corners of a box, given
// the box's coordinates [0..NumBoxesPerAxis] and the number of boxes
// that fit along each axis
void BoxCorners(int3 _boxCoords, float3 *_minCorner, float3 *_maxCorner,
int _boxesPerAxis) {
// TODO figure out a better way to handle offsets
*_minCorner = convert_float3(_boxCoords) / (float)_boxesPerAxis +
(float3)(0.00001);
*_maxCorner = convert_float3((_boxCoords+1)) / (float)_boxesPerAxis -
(float3)(0.00001);
}
// Intersect a ray specifiec by origin and direction
// with a box specified by opposing corners.
// Returns intersect/no intersect along with t values
// for intersections points.
bool IntersectBox(float3 _boundsMin, float3 _boundsMax,
float3 _rayO, float3 _rayD,
float *_tMinOut, float *_tMaxOut) {
float _tMin, _tMax, tYMin, tYMax, tZMin, tZMax;
float divx = 1.0/_rayD.x;
if (divx >= 0.0) {
_tMin = (_boundsMin.x - _rayO.x) * divx;
_tMax = (_boundsMax.x - _rayO.x) * divx;
} else {
_tMin = (_boundsMax.x - _rayO.x) * divx;
_tMax = (_boundsMin.x - _rayO.x) * divx;
}
float divy = 1.0/_rayD.y;
if (divy >= 0.0) {
tYMin = (_boundsMin.y - _rayO.y) * divy;
tYMax = (_boundsMax.y - _rayO.y) * divy;
} else {
tYMin = (_boundsMax.y - _rayO.y) * divy;
tYMax = (_boundsMin.y - _rayO.y) * divy;
}
if ( (_tMin > tYMax || tYMin > _tMax) ) return false;
if (tYMin > _tMin) _tMin = tYMin;
if (tYMax < _tMax) _tMax = tYMax;
float divz = 1.0/_rayD.z;
if (divz >= 0.0) {
tZMin = (_boundsMin.z - _rayO.z) * divz;
tZMax = (_boundsMax.z - _rayO.z) * divz;
} else {
tZMin = (_boundsMax.z - _rayO.z) * divz;
tZMax = (_boundsMin.z - _rayO.z) * divz;
}
if ( (_tMin > tZMax || tZMin > _tMax) ) return false;
if (tZMin > _tMin) _tMin = tZMin;
if (tZMax < _tMax) _tMax = tZMax;
*_tMinOut = _tMin;
*_tMaxOut = _tMax;
return ( (_tMin < 1e20 && _tMax > -1e20 ) );
}
int3 BrickAtlasCoords(int3 _boxCoords,
__global int *_boxList,
int _numBoxesPerAxis) {
int boxIndex = _boxCoords.x +
_boxCoords.y*_numBoxesPerAxis +
_boxCoords.z*_numBoxesPerAxis*_numBoxesPerAxis;
int x = _boxList[5*boxIndex+1];
int y = _boxList[5*boxIndex+2];
int z = _boxList[5*boxIndex+3];
return (int3)(x, y, z);
}
int BrickSize(int3 _boxCoords,
__global int * _brickList,
int _numBoxesPerAxis) {
int boxIndex = _boxCoords.x +
_boxCoords.y*_numBoxesPerAxis +
_boxCoords.z*_numBoxesPerAxis*_numBoxesPerAxis;
return _brickList[5*boxIndex+4];
}
// Traverse one brick, sampling individual voxels
float4 TraverseBrick(__global __read_only image3d_t _textureAtlas,
__global __read_only float *_tf,
int3 _brickAtlasCoords, int _numBoxesPerAxis,
int _brickSize, float _globalStepSize,
float3 _brickEntry, float3 _brickExit) {
// Number of bricks of this size for each axix
int numBricksPerAxis = _numBoxesPerAxis/_brickSize;
// Find brick coordinates [0 .. NumBricksPerAxis]
// Also taking brick size into account
int3 brickCoords = BoxCoords(_brickEntry, numBricksPerAxis);
float numBoxesPerAxisf = (float)_numBoxesPerAxis;
float numBricksPerAxisf = (float)numBricksPerAxis;
// Calculate local brick entry and exit coordinates [0..1]
float3 localEntryCoords = numBricksPerAxisf *
(_brickEntry - convert_float3(brickCoords)/numBricksPerAxisf);
float3 localExitCoords = numBricksPerAxisf *
(_brickExit - convert_float3(brickCoords)/numBricksPerAxisf);
// Calculate offset into texture atlas
float3 offset = convert_float3(_brickAtlasCoords)/numBoxesPerAxisf;
float3 atlasEntry = localEntryCoords/numBoxesPerAxisf + offset;
float3 atlasExit = localExitCoords/numBoxesPerAxisf + offset;
float3 direction = atlasExit - atlasEntry;
float maxDistance = length(atlasExit - atlasEntry);
direction = normalize(direction);
float localStepSize = _globalStepSize * (float)_brickSize;
float4 color = (float4)(0.0, 0.0, 0.0, 0.0);
float traversed = 0.0;
float3 samplePoint = atlasEntry;
float4 spherical;
const sampler_t atlasSampler = CLK_FILTER_LINEAR |
CLK_NORMALIZED_COORDS_TRUE |
CLK_ADDRESS_CLAMP_TO_EDGE;
while (traversed < maxDistance) {
spherical.xyz = CartesianToSpherical(samplePoint);
spherical.w = 1.0;
// Texture stores intensity in R channel
float i = read_imagef(_textureAtlas, atlasSampler, spherical).x;
// Front-to-back compositing
float4 tf = TransferFunction(_tf, i);
color += (1.0 - color.w)*tf;
traversed += localStepSize;
samplePoint += localStepSize * direction;
}
return color;
}
__kernel void
Raycaster(__global __read_only image2d_t _cubeFront,
__global __read_only image2d_t _cubeBack,
__global __write_only image2d_t _output,
__global __read_only image3d_t _textureAtlas,
__constant struct KernelConstants *_constants,
__global __read_only float *_tf,
__global int *_brickList) {
// Kernel should be launched in 2D with one work item per pixel
int idx = get_global_id(0);
int idy = get_global_id(1);
int2 intCoords = (int2)(idx, idy);
// Sampler for texture reading
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE;
// Read from textures
float4 cubeFrontColor = read_imagef(_cubeFront, sampler, intCoords);
float4 cubeBackColor = read_imagef(_cubeBack, sampler, intCoords);
// Figure out the direction
float3 direction = (cubeBackColor-cubeFrontColor).xyz;
float maxDistance = length(direction);
direction = normalize(direction);
float3 origin = cubeFrontColor.xyz - 0.001 * direction;
// Keep track of distance traversed
float traversed = 0.0;
// Sum colors
float stepSize = _constants->stepSize;
float3 samplePoint = cubeFrontColor.xyz;
float4 color = (float4)(0.0, 0.0, 0.0, 0.0);
// TODO temp
float a, b;
if (IntersectBox((float3)(0.0), (float3)(1.0), origin, direction, &a, &b)) {
color += (float4)(0.0, 0.0, 0.08, 1.0);
}
while (traversed < maxDistance) {
int numBoxesPerAxis = _constants->numBoxesPerAxis;
// Convert the sample point to coords
int3 boxCoords = BoxCoords(samplePoint, numBoxesPerAxis);
// Lookup brick coords and size
int3 brickAtlasCoords = BrickAtlasCoords(boxCoords,
_brickList,
numBoxesPerAxis);
int brickSize = BrickSize(boxCoords, _brickList, numBoxesPerAxis);
// Calculate the box's corners
float3 minCorner, maxCorner;
BoxCorners(boxCoords, &minCorner, &maxCorner, numBoxesPerAxis);
// Intersect ray with box
float tMin, tMax;
IntersectBox(minCorner, maxCorner, origin, direction, &tMin, &tMax);
float3 brickEntry = origin + tMin*direction;
// TODO calculate BRICK exit
float3 brickExit = origin + tMax*direction;
float brickDist = length(brickExit - brickEntry);
// Traverse brick
float4 brickColor = TraverseBrick(_textureAtlas, _tf,
brickAtlasCoords, numBoxesPerAxis,
brickSize, stepSize,
brickEntry, brickExit);
// Compositing
color += (1.0 - brickColor.w)*brickColor;
// Advance ray
samplePoint = brickExit;
traversed += brickDist;
// color += (float4)(i, i, i, 1.0);
samplePoint += direction * stepSize;
traversed += stepSize;
}
// Output
float intensity = _constants->intensity;
color *= intensity*stepSize;
// Write to image
write_imagef(_output, intCoords, color);
}

579
kernels/RaycasterTSP.cl
View File

@@ -1,579 +0,0 @@
struct KernelConstants {
int gridType_;
float stepsize_;
float intensity_;
int numTimesteps_;
int numValuesPerNode_;
int numOTNodes_;
int numBoxesPerAxis_;
float temporalTolerance_;
float spatialTolerance_;
int rootLevel_;
int paddedBrickDim_;
};
float3 CartesianToSpherical(float3 _cartesian);
float Lerp(float _v0, float _v1, float _d);
int LeftBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp);
int RightBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp);
int ChildNodeIndex(int _bstNodeIndex,
int *_timespanStart,
int *_timespanEnd,
int _timestep,
int _numValuesPerNode,
int _numOTNodes,
bool _bstRoot,
__global __read_only int *_tsp);
int BrickIndex(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp);
bool IsBSTLeaf(int _bstNodeIndex, int _numValuesPerNode,
bool _bstRoot, __global __read_only int *_tsp);
bool IsOctreeLeaf(int _otNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp);
int OTChildIndex(int _otNodeIndex, int _numValuesPerNode,
int _child,
__global __read_only int *_tsp) ;
float TemporalError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) ;
float SpatialError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp);
int3 BoxCoords(float3 _globalCoords, int _boxesPerAxis) ;
int3 AtlasBoxCoords(int _brickIndex,
__global __read_only int *_brickList);
float3 AtlasCoords(float3 _globalCoords, int _brickIndex, int _boxesPerAxis,
int _paddedBrickDim, int _level,
__global __read_only int *_brickList) ;
void SampleAtlas(float4 *_color, float3 _coords, int _brickIndex,
int _boxesPerAxis, int _paddedBrickDim, int _level,
const sampler_t _atlasSampler,
__read_only image3d_t _textureAtlas,
__read_only image2d_t _transferFunction,
const sampler_t _tfSampler,
__global __read_only int *_brickList);
bool TraverseBST(int _otNodeIndex, int *_brickIndex,
__constant __read_only struct KernelConstants *_constants,
__global __read_only int *_tsp, const int _timestep);
int EnclosingChild(float3 _P, float _boxMid, float3 _offset) ;
void UpdateOffset(float3 *_offset, float _boxDim, int _child) ;
float4 TraverseOctree(float3 _rayO, float3 _rayD, float _maxDist,
__read_only image3d_t _textureAtlas,
__constant struct KernelConstants *_constants,
__read_only image2d_t _transferFunction,
__global __read_only int *_tsp,
__global __read_only int *_brickList,
const int _timestep);
__kernel void RaycasterTSP(__read_only image2d_t _cubeFront,
__read_only image2d_t _cubeBack,
__write_only image2d_t _output,
__read_only image3d_t _textureAtlas,
__constant struct KernelConstants *_constants,
__read_only image2d_t _transferFunction,
//__global __read_only float *_transferFunction,
__global __read_only int *_tsp,
__global __read_only int *_brickList,
const int _timestep);
// Turn normalized [0..1] cartesian coordinates
// to normalized spherical [0..1] coordinates
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0f* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acospi(_cartesian.z/r);
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / native_sqrt(3.0f);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
// Linearly interpolate between two values. Distance
// is assumed to be normalized.
float Lerp(float _v0, float _v1, float _d) {
return _v0*(1.0 - _d) + _v1*_d;
}
/*
float4 TransferFunction(__global __read_only image2d_t _tf,
const sampler_t _tfSampler,
float _i) {
// TODO remove hard-coded value and change to 1D texture
int i0 = (int)floor(1023.0*_i);
int i1 = (i0 < 1023) ? i0+1 : i0;
float di = _i - floor(_i);
float tfr0 = _tf[i0*4+0];
float tfr1 = _tf[i1*4+0];
float tfg0 = _tf[i0*4+1];
float tfg1 = _tf[i1*4+1];
float tfb0 = _tf[i0*4+2];
float tfb1 = _tf[i1*4+2];
float tfa0 = _tf[i0*4+3];
float tfa1 = _tf[i1*4+3];
float tfr = Lerp(tfr0, tfr1, di);
float tfg = Lerp(tfg0, tfg1, di);
float tfb = Lerp(tfb0, tfb1, di);
float tfa = Lerp(tfa0, tfa1, di);
return clamp((float4)(tfr, tfg, tfb, tfa), (float4)(0.0), (float4)(1.0));
//return (float4)(tfr, tfg, tfb, tfa);
float2 sampleCoords = (float2)(1.0, _i);
return read_imagef(_tf, _tfSampler, sampleCoords);
}
*/
/*
// Return index to the octree root (same index as BST root at that OT node)
int OctreeRootNodeIndex() {
return 0;
}
*/
// Return index to left BST child (low timespan)
int LeftBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp) {
// If the BST node is a root, the child pointer is used for the OT.
// The child index is next to the root.
// If not root, look up in TSP structure.
if (_bstRoot) {
return _bstNodeIndex + _numOTNodes;
} else {
return _tsp[_bstNodeIndex*_numValuesPerNode + 1];
}
}
// Return index to right BST child (high timespan)
int RightBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp) {
if (_bstRoot) {
return _bstNodeIndex + _numOTNodes*2;
} else {
return _tsp[_bstNodeIndex*_numValuesPerNode + 1] + _numOTNodes;
}
}
// Return child node index given a BST node, a time span and a timestep
// Updates timespan
int ChildNodeIndex(int _bstNodeIndex,
int *_timespanStart,
int *_timespanEnd,
int _timestep,
int _numValuesPerNode,
int _numOTNodes,
bool _bstRoot,
__global __read_only int *_tsp) {
// Choose left or right child
int middle = *_timespanStart + (*_timespanEnd - *_timespanStart)/2;
if (_timestep <= middle) {
// Left subtree
*_timespanEnd = middle;
return LeftBST(_bstNodeIndex, _numValuesPerNode, _numOTNodes,
_bstRoot, _tsp);
} else {
// Right subtree
*_timespanStart = middle+1;
return RightBST(_bstNodeIndex, _numValuesPerNode, _numOTNodes,
_bstRoot, _tsp);
}
}
// Return the brick index that a BST node represents
int BrickIndex(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return _tsp[_bstNodeIndex*_numValuesPerNode + 0];
}
// Checks if a BST node is a leaf ot not
bool IsBSTLeaf(int _bstNodeIndex, int _numValuesPerNode,
bool _bstRoot, __global __read_only int *_tsp) {
if (_bstRoot) return false;
return (_tsp[_bstNodeIndex*_numValuesPerNode + 1] == -1);
}
// Checks if an OT node is a leaf or not
bool IsOctreeLeaf(int _otNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
// CHILD_INDEX is at offset 1, and -1 represents leaf
return (_tsp[_otNodeIndex*_numValuesPerNode + 1] == -1);
}
// Return OT child index given current node and child number (0-7)
int OTChildIndex(int _otNodeIndex, int _numValuesPerNode,
int _child,
__global __read_only int *_tsp) {
int firstChild = _tsp[_otNodeIndex*_numValuesPerNode + 1];
return firstChild + _child;
}
float TemporalError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return as_float(_tsp[_bstNodeIndex*_numValuesPerNode + 3]);
}
float SpatialError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return as_float(_tsp[_bstNodeIndex*_numValuesPerNode + 2]);
}
// Converts a global coordinate [0..1] to a box coordinate [0..boxesPerAxis]
int3 BoxCoords(float3 _globalCoords, int _boxesPerAxis) {
int3 boxCoords = convert_int3((_globalCoords * (float)_boxesPerAxis));
return clamp(boxCoords, (int3)(0, 0, 0), (int3)(_boxesPerAxis-1));
}
// Fetch atlas box coordinates from brick list
int3 AtlasBoxCoords(int _brickIndex,
__global __read_only int *_brickList) {
int x = _brickList[3*_brickIndex+0];
int y = _brickList[3*_brickIndex+1];
int z = _brickList[3*_brickIndex+2];
return (int3)(x, y, z);
}
// Convert a global coordinate to a local in-box coordinate, given
// the number of boxes (of this size) per axis and the box coordinates
float3 InBoxCoords(float3 _globalCoords, int3 _boxCoords,
int _boxesPerAxis, int _paddedBrickDim) {
// Calculate [0.0 1.0] box coordinates
float3 inbox = (_globalCoords - convert_float3(_boxCoords)/(float)_boxesPerAxis)
* (float)_boxesPerAxis;
// Map to padding range
float low = 1.0/(float)_paddedBrickDim;
float high = (float)(_paddedBrickDim-1)/(float)_paddedBrickDim;
return (float3)(low) + inbox * ((float3)(high)-(float3)(low));
}
float3 AtlasCoords(float3 _globalCoords, int _brickIndex, int _boxesPerAxis,
int _paddedBrickDim, int _level,
__global __read_only int *_brickList) {
// Use current octree level to calculate dividing factor for coordinates
int divisor = (int)pow(2.0, _level);
// Calculate box coordinates, taking current subdivision level into account
int3 boxCoords = BoxCoords(_globalCoords, _boxesPerAxis/divisor);
// Calculate local in-box coordinates for the point
float3 inBoxCoords = InBoxCoords(_globalCoords, boxCoords,
_boxesPerAxis/divisor,
_paddedBrickDim*divisor);
// Fetch atlas box coordinates
int3 atlasBoxCoords = AtlasBoxCoords(_brickIndex, _brickList);
// Transform coordinates to atlas coordinates
return inBoxCoords/(float)_boxesPerAxis +
convert_float3(atlasBoxCoords)/(float)_boxesPerAxis;
}
// Sample atlas
void SampleAtlas(float4 *_color, float3 _coords, int _brickIndex,
int _boxesPerAxis, int _paddedBrickDim, int _level,
const sampler_t _atlasSampler,
__read_only image3d_t _textureAtlas,
__read_only image2d_t _transferFunction,
const sampler_t _tfSampler,
__global __read_only int *_brickList) {
// Find the texture atlas coordinates for the point
float3 atlasCoords = AtlasCoords(_coords, _brickIndex,
_boxesPerAxis, _paddedBrickDim,
_level, _brickList);
int3 boxCoords = BoxCoords(_coords, _boxesPerAxis);
float4 a4 = (float4)(atlasCoords.x, atlasCoords.y, atlasCoords.z, 1.0);
// Sample the atlas
float sample = read_imagef(_textureAtlas, _atlasSampler, a4).x;
// Composition
float4 tf = read_imagef(_transferFunction, _tfSampler, (float2)(sample, 0.0));
*_color += (1.0f - _color->w)*tf;
}
bool TraverseBST(int _otNodeIndex, int *_brickIndex,
__constant __read_only struct KernelConstants *_constants,
__global __read_only int *_tsp, const int _timestep) {
// Start att the root of the current BST
int bstNodeIndex = _otNodeIndex;
bool bstRoot = true;
int timespanStart = 0;
int timespanEnd = _constants->numTimesteps_;
while (true) {
*_brickIndex = BrickIndex(bstNodeIndex,
_constants->numValuesPerNode_,
_tsp);
// Check temporal error
if (TemporalError(bstNodeIndex, _constants->numValuesPerNode_, _tsp) <=
_constants->temporalTolerance_) {
// If the OT node is a leaf, we cannot do any better spatially
if (IsOctreeLeaf(_otNodeIndex, _constants->numValuesPerNode_, _tsp)) {
return true;
} else if (SpatialError(bstNodeIndex, _constants->numValuesPerNode_,
_tsp) <= _constants->spatialTolerance_) {
return true;
} else if (IsBSTLeaf(bstNodeIndex, _constants->numValuesPerNode_,
bstRoot, _tsp)) {
return false;
} else {
// Keep traversing
bstNodeIndex = ChildNodeIndex(bstNodeIndex, &timespanStart,
&timespanEnd,
_timestep,
_constants->numValuesPerNode_,
_constants->numOTNodes_,
bstRoot, _tsp);
}
} else if (IsBSTLeaf(bstNodeIndex, _constants->numValuesPerNode_,
bstRoot, _tsp)) {
return false;
} else {
// Keep traversing
bstNodeIndex = ChildNodeIndex(bstNodeIndex, &timespanStart,
&timespanEnd,
_timestep,
_constants->numValuesPerNode_,
_constants->numOTNodes_,
bstRoot, _tsp);
}
bstRoot = false;
}
}
// Given a point, a box mid value and an offset, return enclosing child
int EnclosingChild(float3 _P, float _boxMid, float3 _offset) {
if (_P.x < _boxMid+_offset.x) {
if (_P.y < _boxMid+_offset.y) {
if (_P.z < _boxMid+_offset.z) {
return 0;
} else {
return 4;
}
} else {
if (_P.z < _boxMid+_offset.z) {
return 2;
} else {
return 6;
}
}
} else {
if (_P.y < _boxMid+_offset.y) {
if (_P.z < _boxMid+_offset.z) {
return 1;
} else {
return 5;
}
} else {
if (_P.z < _boxMid+_offset.z) {
return 3;
} else {
return 7;
}
}
}
}
void UpdateOffset(float3 *_offset, float _boxDim, int _child) {
if (_child == 0) {
// do nothing
} else if (_child == 1) {
_offset->x += _boxDim;
} else if (_child == 2) {
_offset->y += _boxDim;
} else if (_child == 3) {
_offset->x += _boxDim;
_offset->y += _boxDim;
} else if (_child == 4) {
_offset->z += _boxDim;
} else if (_child == 5) {
_offset->x += _boxDim;
_offset->z += _boxDim;
} else if (_child == 6) {
_offset->y += _boxDim;
_offset->z += _boxDim;
} else if (_child == 7) {
*_offset += (float3)(_boxDim);
}
}
float4 TraverseOctree(float3 _rayO, float3 _rayD, float _maxDist,
__read_only image3d_t _textureAtlas,
__constant struct KernelConstants *_constants,
__read_only image2d_t _transferFunction,
__global __read_only int *_tsp,
__global __read_only int *_brickList,
const int _timestep) {
float stepsize = _constants->stepsize_;
// Sample point
float3 cartesianP = _rayO;
// Keep track of distance traveled along ray
float traversed = 0;
// Cumulative color for ray to return
float4 color = (float4)(0.0);
// Sampler for texture atlas
const sampler_t atlasSampler = CLK_FILTER_LINEAR |
CLK_NORMALIZED_COORDS_TRUE |
CLK_ADDRESS_CLAMP_TO_EDGE;
// Sampler for transfer function texture
const sampler_t tfSampler = CLK_FILTER_LINEAR |
CLK_NORMALIZED_COORDS_TRUE |
CLK_ADDRESS_CLAMP_TO_EDGE;
// Traverse until sample point is outside of volume
while (traversed < _maxDist) {
// Reset octree traversal variables
float3 offset = (float3)(0.0);
float boxDim = 1.0;
bool foundBrick = false;
int child;
int level = _constants->rootLevel_;
int otNodeIndex = 0;
// Rely on finding a leaf for loop termination
while (true) {
// Traverse BST to get a brick index, and see if the found brick
// is good enough
int brickIndex;
bool bstSuccess = TraverseBST(otNodeIndex,
&brickIndex,
_constants,
_tsp,
_timestep);
// Convert to spherical if needed
float3 sampleP;
if (_constants->gridType_ == 0) { // cartesian
sampleP = cartesianP;
} else { // spherical ( == 1)
sampleP = CartesianToSpherical(cartesianP);
}
if (bstSuccess ||
IsOctreeLeaf(otNodeIndex, _constants->numValuesPerNode_, _tsp)) {
//float s = 0.008*SpatialError(brickIndex, 4, _tsp);
//color += (float4)(s);
// Sample the brick
SampleAtlas(&color, sampleP, brickIndex,
_constants->numBoxesPerAxis_,
_constants->paddedBrickDim_,
level,
atlasSampler, _textureAtlas,
_transferFunction,
tfSampler, _brickList);
break;
} else {
// Keep traversing the octree
// Next box dimension
boxDim /= 2.0f;
// Current mid point
float boxMid = boxDim;
// Check which child encloses the sample point
child = EnclosingChild(sampleP, boxMid, offset);
// Update offset for next level
UpdateOffset(&offset, boxDim, child);
// Update index to new node
otNodeIndex = OTChildIndex(otNodeIndex, _constants->numValuesPerNode_,
child, _tsp);
level--;
}
} // while traversing
// Update sample point
traversed += stepsize;
cartesianP += stepsize * _rayD;
} // while (traversed < maxDist)
return color;
}
__kernel void RaycasterTSP(__read_only image2d_t _cubeFront,
__read_only image2d_t _cubeBack,
__write_only image2d_t _output,
__read_only image3d_t _textureAtlas,
__constant struct KernelConstants *_constants,
__read_only image2d_t _transferFunction,
//__global __read_only float *_transferFunction,
__global __read_only int *_tsp,
__global __read_only int *_brickList,
const int _timestep) {
// Kernel should be launched in 2D with one work item per pixel
int2 intCoords = (int2)(get_global_id(0), get_global_id(1));
// Sampler for color cube reading
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE;
// Read from color cube textures
float4 cubeFrontColor = read_imagef(_cubeFront, sampler, intCoords);
float4 cubeBackColor = read_imagef(_cubeBack, sampler, intCoords);
// Figure out ray direction and distance to traverse
float3 direction = cubeBackColor.xyz - cubeFrontColor.xyz;
float maxDist = length(direction);
direction = normalize(direction);
float4 color = TraverseOctree(cubeFrontColor.xyz, // ray origin
direction, // direction
maxDist, // distance to traverse
_textureAtlas, // voxel data atlas
_constants, // kernel constants
_transferFunction, // transfer function
_tsp, // TSP tree struct
_brickList,
_timestep);
//color = 0.0001*color + cubeFrontColor;
write_imagef(_output, intCoords, _constants->intensity_*color);
return;
}

124
kernels/RaycasterTexture.cl
View File

@@ -1,124 +0,0 @@
// Needs to mirror struct on host
struct KernelConstants {
float stepSize;
float intensity;
int aDim;
int bDim;
int cDim;
};
// Linearly interpolate between two values. Distance
// is assumed to be normalized.
float Lerp(float _v0, float _v1, float _d) {
return _v0*(1.0 - _d) + _v1*_d;
}
// Turn normalized [0..1] cartesian coordinates
// to normalized spherical [0..1] coordinates
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acospi(_cartesian.z/r);
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / native_sqrt(3.0);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
float4 TransferFunction(__global __read_only float *_tf, float _i) {
// TODO remove hard-coded value and change to 1D texture
int i0 = (int)floor(1023.0*_i);
int i1 = (i0 < 1023) ? i0+1 : i0;
float di = _i - floor(_i);
float tfr0 = _tf[i0*4+0];
float tfr1 = _tf[i1*4+0];
float tfg0 = _tf[i0*4+1];
float tfg1 = _tf[i1*4+1];
float tfb0 = _tf[i0*4+2];
float tfb1 = _tf[i1*4+2];
float tfa0 = _tf[i0*4+3];
float tfa1 = _tf[i1*4+3];
float tfr = Lerp(tfr0, tfr1, di);
float tfg = Lerp(tfg0, tfg1, di);
float tfb = Lerp(tfb0, tfb1, di);
float tfa = Lerp(tfa0, tfa1, di);
return (float4)(tfr, tfg, tfb, tfa);
}
__kernel void
Raycaster(__global __read_only image2d_t _cubeFront,
__global __read_only image2d_t _cubeBack,
__global __write_only image2d_t _output,
__global __read_only image3d_t _voxelData,
__constant struct KernelConstants *_constants,
__global __read_only float *_tf) {
int3 dimensions = (int3)(_constants->aDim,
_constants->bDim,
_constants->cDim);
// Kernel should be launched in 2D with one work item per pixel
int idx = get_global_id(0);
int idy = get_global_id(1);
int2 intCoords = (int2)(idx, idy);
// Sampler for texture reading
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE;
const sampler_t volumeSampler = CLK_FILTER_LINEAR |
CLK_NORMALIZED_COORDS_TRUE |
CLK_ADDRESS_CLAMP_TO_EDGE;
// Read from textures
float4 cubeFrontColor = read_imagef(_cubeFront, sampler, intCoords);
float4 cubeBackColor = read_imagef(_cubeBack, sampler, intCoords);
// Figure out the direction
float3 direction = (cubeBackColor-cubeFrontColor).xyz;
float maxDistance = length(direction);
direction = normalize(direction);
// Keep track of distance traversed
float traversed = 0.0;
// Sum colors
float stepSize = _constants->stepSize;
float3 samplePoint = cubeFrontColor.xyz;
float4 spherical;
float4 color = (float4)(0.0, 0.0, 0.0, 0.0);
while (traversed < maxDistance) {
spherical.xyz = CartesianToSpherical(samplePoint);
spherical.w = 1.0;
// Texture stores intensity in R channel
float i = read_imagef(_voxelData, volumeSampler, spherical).x;
// Front-to-back compositing
float4 tf = TransferFunction(_tf, i);
color += (1.0 - color.w)*tf;
// color += (float4)(i, i, i, 1.0);
samplePoint += direction * stepSize;
traversed += stepSize;
}
// Output
float intensity = _constants->intensity;
color *= intensity*stepSize;
// Write to image
write_imagef(_output, intCoords, color);
}

445
kernels/TSPTraversal.cl
View File

@@ -1,445 +0,0 @@
// Mirrors struct on host side
struct TraversalConstants {
int gridType_;
float stepsize_;
int numTimesteps_;
int numValuesPerNode_;
int numOTNodes_;
float temporalTolerance_;
float spatialTolerance_;
};
float3 CartesianToSpherical(float3 _cartesian);
int OctreeRootNodeIndex();
int LeftBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp);
int RightBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp);
int ChildNodeIndex(int _bstNodeIndex,
int *_timespanStart,
int *_timespanEnd,
int _timestep,
int _numValuesPerNode,
int _numOTNodes,
bool _bstRoot,
__global __read_only int *_tsp);
int BrickIndex(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) ;
bool IsBSTLeaf(int _bstNodeIndex, int _numValuesPerNode,
bool _bstRoot, __global __read_only int *_tsp) ;
bool IsOctreeLeaf(int _otNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) ;
int OTChildIndex(int _otNodeIndex, int _numValuesPerNode,
int _child,
__global __read_only int *_tsp);
void AddToList(int _brickIndex,
__global volatile int *_reqList);
float TemporalError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) ;
float SpatialError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp);
bool TraverseBST(int _otNodeIndex,
int *_brickIndex,
int _timestep,
__constant struct TraversalConstants *_constants,
__global volatile int *_reqList,
__global __read_only int *_tsp);
int EnclosingChild(float3 _P, float _boxMid, float3 _offset);
void UpdateOffset(float3 *_offset, float _boxDim, int _child) ;
void TraverseOctree(float3 _rayO,
float3 _rayD,
float _maxDist,
__constant struct TraversalConstants *_constants,
__global volatile int *_reqList,
__global __read_only int *_tsp,
const int _timestep);
__kernel void TSPTraversal(__read_only image2d_t _cubeFront,
__read_only image2d_t _cubeBack,
__constant struct TraversalConstants *_constants,
__global __read_only int *_tsp,
__global int *_reqList,
const int _timestep) ;
// Turn normalized [0..1] cartesian coordinates
// to normalized spherical [0..1] coordinates
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0f* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acospi(_cartesian.z/r);
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / native_sqrt(3.0f);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
// Return index to the octree root (same index as BST root at that OT node)
int OctreeRootNodeIndex() {
return 0;
}
// Return index to left BST child (low timespan)
int LeftBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp) {
// If the BST node is a root, the child pointer is used for the OT.
// The child index is next to the root.
// If not root, look up in TSP structure.
if (_bstRoot) {
return _bstNodeIndex + _numOTNodes;
//return _bstNodeIndex + 1;
} else {
return _tsp[_bstNodeIndex*_numValuesPerNode + 1];
}
}
// Return index to right BST child (high timespan)
int RightBST(int _bstNodeIndex, int _numValuesPerNode, int _numOTNodes,
bool _bstRoot, __global __read_only int *_tsp) {
if (_bstRoot) {
return _bstNodeIndex + _numOTNodes*2;
} else {
return _tsp[_bstNodeIndex*_numValuesPerNode + 1] + _numOTNodes;
}
}
// Return child node index given a BST node, a time span and a timestep
// Updates timespan
int ChildNodeIndex(int _bstNodeIndex,
int *_timespanStart,
int *_timespanEnd,
int _timestep,
int _numValuesPerNode,
int _numOTNodes,
bool _bstRoot,
__global __read_only int *_tsp) {
// Choose left or right child
int middle = *_timespanStart + (*_timespanEnd - *_timespanStart)/2;
if (_timestep <= middle) {
// Left subtree
*_timespanEnd = middle;
return LeftBST(_bstNodeIndex, _numValuesPerNode, _numOTNodes,
_bstRoot, _tsp);
} else {
// Right subtree
*_timespanStart = middle+1;
return RightBST(_bstNodeIndex, _numValuesPerNode, _numOTNodes,
_bstRoot, _tsp);
}
}
// Return the brick index that a BST node represents
int BrickIndex(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return _tsp[_bstNodeIndex*_numValuesPerNode + 0];
}
// Checks if a BST node is a leaf ot not
bool IsBSTLeaf(int _bstNodeIndex, int _numValuesPerNode,
bool _bstRoot, __global __read_only int *_tsp) {
if (_bstRoot) return false;
return (_tsp[_bstNodeIndex*_numValuesPerNode + 1] == -1);
}
// Checks if an OT node is a leaf or not
bool IsOctreeLeaf(int _otNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
// CHILD_INDEX is at offset 1, and -1 represents leaf
return (_tsp[_otNodeIndex*_numValuesPerNode + 1] == -1);
}
// Return OT child index given current node and child number (0-7)
int OTChildIndex(int _otNodeIndex, int _numValuesPerNode,
int _child,
__global __read_only int *_tsp) {
int firstChild = _tsp[_otNodeIndex*_numValuesPerNode+1];
return firstChild + _child;
}
// Increment the count for a brick in the request list
void AddToList(int _brickIndex,
__global volatile int *_reqList) {
atomic_inc(&_reqList[_brickIndex]);
}
float TemporalError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return as_float(_tsp[_bstNodeIndex*_numValuesPerNode + 3]);
}
float SpatialError(int _bstNodeIndex, int _numValuesPerNode,
__global __read_only int *_tsp) {
return as_float(_tsp[_bstNodeIndex*_numValuesPerNode + 2]);
}
// Given an octree node index, traverse the corresponding BST tree and look
// for a useful brick.
bool TraverseBST(int _otNodeIndex,
int *_brickIndex,
int _timestep,
__constant struct TraversalConstants *_constants,
__global volatile int *_reqList,
__global __read_only int *_tsp) {
// Start at the root of the current BST
int bstNodeIndex = _otNodeIndex;
bool bstRoot = true;
int timespanStart = 0;
int timespanEnd = _constants->numTimesteps_;
// Rely on structure for termination
while (true) {
// Update brick index (regardless if we use it or not)
*_brickIndex = BrickIndex(bstNodeIndex,
_constants->numValuesPerNode_,
_tsp);
// If temporal error is ok
// TODO float and <= errors
if (TemporalError(bstNodeIndex, _constants->numValuesPerNode_,
_tsp) <= _constants->temporalTolerance_) {
// If the ot node is a leaf, we can't do any better spatially so we
// return the current brick
if (IsOctreeLeaf(_otNodeIndex, _constants->numValuesPerNode_, _tsp)) {
return true;
// All is well!
} else if (SpatialError(bstNodeIndex, _constants->numValuesPerNode_,
_tsp) <= _constants->spatialTolerance_) {
return true;
// If spatial failed and the BST node is a leaf
// The traversal will continue in the octree (we know that
// the octree node is not a leaf)
} else if (IsBSTLeaf(bstNodeIndex, _constants->numValuesPerNode_,
bstRoot, _tsp)) {
return false;
// Keep traversing BST
} else {
bstNodeIndex = ChildNodeIndex(bstNodeIndex,
&timespanStart,
&timespanEnd,
_timestep,
_constants->numValuesPerNode_,
_constants->numOTNodes_,
bstRoot,
_tsp);
}
// If temporal error is too big and the node is a leaf
// Return false to traverse OT
} else if (IsBSTLeaf(bstNodeIndex, _constants->numValuesPerNode_,
bstRoot, _tsp)) {
return false;
// If temporal error is too big and we can continue
} else {
bstNodeIndex = ChildNodeIndex(bstNodeIndex,
&timespanStart,
&timespanEnd,
_timestep,
_constants->numValuesPerNode_,
_constants->numOTNodes_,
bstRoot,
_tsp);
}
bstRoot = false;
}
}
// Given a point, a box mid value and an offset, return enclosing child
int EnclosingChild(float3 _P, float _boxMid, float3 _offset) {
if (_P.x < _boxMid+_offset.x) {
if (_P.y < _boxMid+_offset.y) {
if (_P.z < _boxMid+_offset.z) {
return 0;
} else {
return 4;
}
} else {
if (_P.z < _boxMid+_offset.z) {
return 2;
} else {
return 6;
}
}
} else {
if (_P.y < _boxMid+_offset.y) {
if (_P.z < _boxMid+_offset.z) {
return 1;
} else {
return 5;
}
} else {
if (_P.z < _boxMid+_offset.z) {
return 3;
} else {
return 7;
}
}
}
}
void UpdateOffset(float3 *_offset, float _boxDim, int _child) {
if (_child == 0) {
// do nothing
} else if (_child == 1) {
_offset->x += _boxDim;
} else if (_child == 2) {
_offset->y += _boxDim;
} else if (_child == 3) {
_offset->x += _boxDim;
_offset->y += _boxDim;
} else if (_child == 4) {
_offset->z += _boxDim;
} else if (_child == 5) {
_offset->x += _boxDim;
_offset->z += _boxDim;
} else if (_child == 6) {
_offset->y += _boxDim;
_offset->z += _boxDim;
} else if (_child == 7) {
*_offset += (float3)(_boxDim);
}
}
// Traverse one ray through the volume, build brick list
void TraverseOctree(float3 _rayO,
float3 _rayD,
float _maxDist,
__constant struct TraversalConstants *_constants,
__global volatile int *_reqList,
__global __read_only int *_tsp,
const int _timestep) {
// Choose a stepsize that guarantees that we don't miss any bricks
// TODO dynamic depending on brick dimensions
float stepsize = _constants->stepsize_;
float3 P = _rayO;
// Keep traversing until the sample point goes outside the unit cube
float traversed = 0.0f;
int max_iterations = 50;
int iterations = 0;
bool ok = stepsize > 0.0f && stepsize < fabs(_maxDist);
//while (traversed < _maxDist && iterations < max_iterations) {
while (traversed < _maxDist) {
// Reset traversal variables
float3 offset = (float3)(0.0f);
float boxDim = 1.0f;
int child;
// Init the octree node index to the root
int otNodeIndex = OctreeRootNodeIndex();
// Start traversing octree
// Rely on finding a leaf for loop termination
while (true) {
// See if the BST tree is good enough
int brickIndex = 0;
bool bstSuccess = TraverseBST(otNodeIndex,
&brickIndex,
_timestep,
_constants,
_reqList,
_tsp);
if (bstSuccess) {
// Add the found brick to brick list
AddToList(brickIndex, _reqList);
// We are now done with this node, so go to next
break;
// If the BST lookup failed but the octree node is a leaf,
// add the brick anyway (it is the BST leaf)
} else if (IsOctreeLeaf(otNodeIndex,
_constants->numValuesPerNode_, _tsp)) {
AddToList(brickIndex, _reqList);
// We are now done with this node, so go to next
break;
// If the BST lookup failed and we can traverse the octree,
// visit the child that encloses the point
} else {
// Next box dimension
boxDim = boxDim/2.0f;
// Current mid point
float boxMid = boxDim;
// Check which child encloses P
if (_constants->gridType_ == 0) { // Cartesian
child = EnclosingChild(P, boxMid, offset);
} else { // Spherical (==1)
child = EnclosingChild(CartesianToSpherical(P), boxMid, offset);
}
// Update offset
UpdateOffset(&offset, boxDim, child);
// Update node index to new node
//int oldIndex = otNodeIndex;
otNodeIndex = OTChildIndex(otNodeIndex, _constants->numValuesPerNode_,
child, _tsp);
}
} // while traversing
// Update
traversed = traversed + stepsize;
P = P + stepsize * _rayD;
} // while (traversed < maxDist)
}
__kernel void TSPTraversal(__read_only image2d_t _cubeFront,
__read_only image2d_t _cubeBack,
__constant struct TraversalConstants *_constants,
__global __read_only int *_tsp,
__global int *_reqList,
const int _timestep) {
// Kernel should be launched in 2D with one work item per pixel
int2 intCoords = (int2)(get_global_id(0), get_global_id(1));
// Sampler for color cube reading
const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE;
// Read from color cube textures
float4 cubeFrontColor = read_imagef(_cubeFront, sampler, intCoords);
if (length(cubeFrontColor.xyz) == 0.0) return;
float4 cubeBackColor = read_imagef(_cubeBack, sampler, intCoords);
// Figure out ray direction
float maxDist = length(cubeBackColor.xyz-cubeFrontColor.xyz);
float3 direction = normalize(cubeBackColor.xyz-cubeFrontColor.xyz);
// Traverse octree and fill the brick request list
TraverseOctree(cubeFrontColor.xyz, direction, maxDist,
_constants, _reqList, _tsp, _timestep);
return;
}

28
kernels/helpers/volume_helpers.cl
View File

@@ -1,28 +0,0 @@
float3 CartesianToSpherical(float3 _cartesian);
float intensityNormalizer(float intensity, float iMin, float iMax);
float3 CartesianToSpherical(float3 _cartesian) {
// Put cartesian in [-1..1] range first
_cartesian = (float3)(-1.0) + 2.0f* _cartesian;
float r = length(_cartesian);
float theta, phi;
if (r == 0.0) {
theta = phi = 0.0;
} else {
theta = acospi(_cartesian.z/r);
phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI);
}
r = r / native_sqrt(3.0f);
// Sampler ignores w component
return (float3)(r, theta, phi);
}
float intensityNormalizer(float intensity, float iMin, float iMax) {
float i = clamp(intensity, iMin, iMax);
i = (i - iMin) / (iMax - iMin);
return clamp(i, 0.0f, 1.0f);
}

113
kernels/helpers/volume_raycasting.cl
View File

@@ -1,113 +0,0 @@
#ifdef CL_VERSION_1_2
#define RC_TF_TYPE image1d_t
#define RC_TF_MAP(intensity) intensity
#else
#define RC_TF_TYPE image2d_t
#define RC_TF_MAP(intensity) (float2)(intensity,0.0f)
#endif
#define EARLY_RAY_TERMINATION_OPACITY 0.95
void raySetup(float3 first, float3 last, float3 dimension, float3* rayDirection, float* tIncr, float* tEnd);
bool earlyRayTermination(float4* color, float maxOpacity);
#define RC_DEFINE_TEXTUE_COORDINATES(coords) \
int2 coords = (int2)(get_global_id(0), get_global_id(1))
#define RC_DEFINE_VOLUME3D_DIMENSIONS(dimension, volume) \
float3 dimension; \
{ \
int4 idim = get_image_dim(volume); \
dimension = (float3)(idim.x,idim.y,idim.z); \
}
//#define RC_DEFINE_VARIABLES(dimension, intCoords)
/***
* Calculates the direction of the ray and returns the number
* of steps and the direction.
***/
void raySetup(float3 first, float3 last, float3 dimension, float3* rayDirection, float* tIncr, float* tEnd) {
float samplingRate_ = 1.0f;
*rayDirection = last - first;
*tEnd = length(*rayDirection);
*rayDirection = normalize(*rayDirection);
*tIncr = 1.0/(samplingRate_*length((*rayDirection)*dimension));
}
/***
* Applies early ray termination. The current opacity is compared to
* the maximum opacity. In case it is greater, the opacity is set to
* 1.0 and true is returned, otherwise false is returned.
***/
bool earlyRayTermination(float4* color, float maxOpacity) {
if ((*color).w >= maxOpacity) {
(*color).w = 1.0f;
return true;
} else {
return false;
}
}
#define RC_EARLY_RAY_TERMINATION(opacity, maxOpacity, finished) \
finished = earlyRayTermination(&opacity, maxOpacity)
#define RC_BEGIN_LOOP_FOR \
for (int loop=0; !finished && loop<RAYCASTING_LOOP_COUNT; ++loop) {
#define RC_END_LOOP_BRACES }
/***
* The beginning of a typical raycasting loop.
*/
#define RC_BEGIN_LOOP \
bool finished = false; \
float t = 0.0f; \
int RAYCASTING_LOOP_COUNT = tEnd / tIncr; \
RC_BEGIN_LOOP_FOR
/***
* The end of a typical raycasting loop. If adaptive sampling
* is used for rendering bricked volumes, t is increased by a
* multiple of tIncr, thereby skipping several samples.
*/
#ifdef ADAPTIVE_SAMPLING
#define RC_END_LOOP(result) \
RC_EARLY_RAY_TERMINATION(result, EARLY_RAY_TERMINATION_OPACITY, finished); \
t += (tIncr * float(numberOfSkippedSamples)); \
finished = finished || (t > tEnd); \
RC_END_LOOP_BRACES
#else
#define RC_END_LOOP(result) \
RC_EARLY_RAY_TERMINATION(result, EARLY_RAY_TERMINATION_OPACITY, finished); \
t += tIncr; \
finished = finished || (t > tEnd); \
RC_END_LOOP_BRACES
#endif
/**
* In order to keep the shaders as free as possible from dealing
* with bricking and adaptive sampling, these defines can be placed
* before and after the compositing function calls to enable adaptive
* sampling when bricking is used. For normal datasets these defines
* will have no impact at all.
*/
#ifdef ADAPTIVE_SAMPLING
#define RC_BEGIN_COMPOSITING \
for (int i=0; i<numberOfSkippedSamples; i++) {
#else
#define RC_BEGIN_COMPOSITING
#endif
#ifdef ADAPTIVE_SAMPLING
#define RC_END_COMPOSITING \
}
#else
#define RC_END_COMPOSITING
#endif

10
kernels/test.cl
View File

@@ -1,10 +0,0 @@
#ifndef OFFSET
#define OFFSET 0
#endif
__kernel void hello(__global int * out)
{
size_t tid = get_global_id(0);
out[tid] = tid + OFFSET;
}