From f3a5b976adfcb0cda37f1fe005afc2f8d1e57a02 Mon Sep 17 00:00:00 2001 From: Alexander Bock Date: Sun, 1 May 2016 12:38:16 +0200 Subject: [PATCH] Remove old OpenCL kernels --- kernels/Raycaster.cl | 202 ---------- kernels/RaycasterBricks.cl | 297 -------------- kernels/RaycasterTSP.cl | 579 --------------------------- kernels/RaycasterTexture.cl | 124 ------ kernels/TSPTraversal.cl | 445 -------------------- kernels/helpers/volume_helpers.cl | 28 -- kernels/helpers/volume_raycasting.cl | 113 ------ kernels/test.cl | 10 - 8 files changed, 1798 deletions(-) delete mode 100644 kernels/Raycaster.cl delete mode 100644 kernels/RaycasterBricks.cl delete mode 100644 kernels/RaycasterTSP.cl delete mode 100644 kernels/RaycasterTexture.cl delete mode 100644 kernels/TSPTraversal.cl delete mode 100644 kernels/helpers/volume_helpers.cl delete mode 100644 kernels/helpers/volume_raycasting.cl delete mode 100644 kernels/test.cl diff --git a/kernels/Raycaster.cl b/kernels/Raycaster.cl deleted file mode 100644 index 98474ced33..0000000000 --- a/kernels/Raycaster.cl +++ /dev/null @@ -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); - -} - diff --git a/kernels/RaycasterBricks.cl b/kernels/RaycasterBricks.cl deleted file mode 100644 index 278008e32d..0000000000 --- a/kernels/RaycasterBricks.cl +++ /dev/null @@ -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); - -} - diff --git a/kernels/RaycasterTSP.cl b/kernels/RaycasterTSP.cl deleted file mode 100644 index de239b4c07..0000000000 --- a/kernels/RaycasterTSP.cl +++ /dev/null @@ -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, ×panStart, - ×panEnd, - _timestep, - _constants->numValuesPerNode_, - _constants->numOTNodes_, - bstRoot, _tsp); - } - - } else if (IsBSTLeaf(bstNodeIndex, _constants->numValuesPerNode_, - bstRoot, _tsp)) { - return false; - - } else { - - // Keep traversing - bstNodeIndex = ChildNodeIndex(bstNodeIndex, ×panStart, - ×panEnd, - _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; -} - diff --git a/kernels/RaycasterTexture.cl b/kernels/RaycasterTexture.cl deleted file mode 100644 index a4c1e6f422..0000000000 --- a/kernels/RaycasterTexture.cl +++ /dev/null @@ -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); - -} - diff --git a/kernels/TSPTraversal.cl b/kernels/TSPTraversal.cl deleted file mode 100644 index 0ea26357b4..0000000000 --- a/kernels/TSPTraversal.cl +++ /dev/null @@ -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, - ×panStart, - ×panEnd, - _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, - ×panStart, - ×panEnd, - _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; - -} diff --git a/kernels/helpers/volume_helpers.cl b/kernels/helpers/volume_helpers.cl deleted file mode 100644 index b3ba1d07ad..0000000000 --- a/kernels/helpers/volume_helpers.cl +++ /dev/null @@ -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); -} - - - diff --git a/kernels/helpers/volume_raycasting.cl b/kernels/helpers/volume_raycasting.cl deleted file mode 100644 index ae8402606f..0000000000 --- a/kernels/helpers/volume_raycasting.cl +++ /dev/null @@ -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 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