Sample ray grid traversal (#776)

* Camera: add parameter to allow text centering
* Camera: Improve offset calculation using spatial X value
* Camera: move message center logic to update function
* Add: sampleRayGridTraversal function
* Fix: grid traversal out of bounds error
* RAIVolume: fix grid traversal by setting voxel sizes to 1
* SamplingExamples: update diagrams and diagram readability
* Change sampling log to debug
* DefaultSpatial: add optional argument to ignore children in AABB intersection
* RAIVolumeSamplingExample: fix diagram display, use sampleRayGridTraversal
* RAIVolume: potentially fix sampling over/underflows
* Fix: sampleRay and sampleRayGridTraversal both correctly sample anisotropy now
* RAIVolumeSamplingExample: change used sampling method back to original ray sampling
* BufferedVolume: fix sampleRay functionality, fix sampling beyond limit issue
* VolumeSamplingExample: change info to debug logs, move diagram lower for better visibility
* RAIVolume: clamp UV space in sample() to prevent index out of range errors
* RAIVolume: simplify bounding box scaling calculation
* RAIVolume: better handle getVoxelScale() case of null dimensions
* RAIVolume: remove duplicate try-catch clause
* OrientedBoundingBox: add isInside method to check if a point lies inside the bounding box
* RAIVolume: fix typo in docs

src/main/kotlin/graphics/scenery/OrientedBoundingBox.kt

@@ -94,6 +94,14 @@ open class OrientedBoundingBox(val n: Node, val min: Vector3f, val max: Vector3f
         }
     }
 
+    /**
+     * Checks whether a [point] is inside the given [OrientedBoundingBox] or not.
+     */
+    fun isInside(point: Vector3f): Boolean {
+        return (point.x > min.x && point.y > min.y && point.z > min.z &&
+                point.x < max.x && point.y < max.y && point.z < max.z)
+    }
+
     /**
      * Returns the hash code of this [OrientedBoundingBox], taking [min] and [max] into consideration.
      */

src/main/kotlin/graphics/scenery/attribute/spatial/DefaultSpatial.kt

@@ -285,9 +285,9 @@ open class DefaultSpatial(@Transient protected var node: Node = DefaultNode()) :
      *
      * Code adapted from [zachamarz](http://gamedev.stackexchange.com/a/18459).
      */
-    override fun intersectAABB(origin: Vector3f, dir: Vector3f): MaybeIntersects {
-        val bbmin = node.getMaximumBoundingBox().min.xyzw()
-        val bbmax = node.getMaximumBoundingBox().max.xyzw()
+    override fun intersectAABB(origin: Vector3f, dir: Vector3f, ignoreChildren: Boolean): MaybeIntersects {
+        val bbmin = if (ignoreChildren) node.boundingBox!!.min.xyzw() else node.getMaximumBoundingBox().min.xyzw()
+        val bbmax = if (ignoreChildren) node.boundingBox!!.max.xyzw() else node.getMaximumBoundingBox().max.xyzw()
 
         val min = world.transform(bbmin)
         val max = world.transform(bbmax)

src/main/kotlin/graphics/scenery/attribute/spatial/Spatial.kt

@@ -68,7 +68,7 @@ interface Spatial: RealLocalizable, RealPositionable, Networkable {
      */
     fun worldRotation(): Quaternionf
 
-    fun intersectAABB(origin: Vector3f, dir: Vector3f): MaybeIntersects
+    fun intersectAABB(origin: Vector3f, dir: Vector3f, ignoreChildren: Boolean = false): MaybeIntersects
 
     /**
      * Returns the [Node]'s world position

src/main/kotlin/graphics/scenery/volumes/BufferedVolume.kt

@@ -165,19 +165,21 @@ class BufferedVolume(val ds: VolumeDataSource.RAISource<*>, options: VolumeViewe
         }
 
         val absoluteCoords = Vector3f(uv.x() * dimensions.x(), uv.y() * dimensions.y(), uv.z() * dimensions.z())
-//        val index: Int = (floor(gt.dimensions.x() * gt.dimensions.y() * absoluteCoords.z()).toInt()
-//            + floor(gt.dimensions.x() * absoluteCoords.y()).toInt()
-//            + floor(absoluteCoords.x()).toInt())
         val absoluteCoordsD = Vector3f(floor(absoluteCoords.x()), floor(absoluteCoords.y()), floor(absoluteCoords.z()))
         val diff = absoluteCoords - absoluteCoordsD
 
-        fun toIndex(absoluteCoords: Vector3f): Int = (
-            absoluteCoords.x().roundToInt().dec()
-                + (dimensions.x() * absoluteCoords.y()).roundToInt().dec()
-                + (dimensions.x() * dimensions.y() * absoluteCoords.z()).roundToInt().dec()
-            )
+        // Turn the absolute volume coordinates into a sampling index,
+        // clamping the value to 1 below the maximum dimension to prevent overflow conditions.
+        fun clampedToIndex(absoluteCoords: Vector3f): Int {
+            val x = absoluteCoords.x().coerceIn(0f, dimensions.x() - 1f)
+            val y = absoluteCoords.y().coerceIn(0f, dimensions.y() - 1f)
+            val z = absoluteCoords.z().coerceIn(0f, dimensions.z() - 1f)
+            return (x.roundToInt() +
+                (dimensions.x() * y).roundToInt() +
+                (dimensions.x() * dimensions.y() * z).roundToInt())
+        }
 
-        val index = toIndex(absoluteCoordsD)
+        val index = clampedToIndex(absoluteCoordsD)
 
         val contents = texture.contents.duplicate().order(ByteOrder.LITTLE_ENDIAN)
 
@@ -186,7 +188,6 @@ class BufferedVolume(val ds: VolumeDataSource.RAISource<*>, options: VolumeViewe
             return 0.0f
         }
 
-
         fun density(index: Int): Float {
             if (index * bpp >= contents.limit()) {
                 logger.warn("Sampling beyond limit")
@@ -207,38 +208,43 @@ class BufferedVolume(val ds: VolumeDataSource.RAISource<*>, options: VolumeViewe
             val transferRangeMax = range.second
 
             val final = transferFunction.evaluate(s / transferRangeMax)
-            logger.info("Sample at $index is $s, final is $final $transferRangeMax")
+            logger.debug("Sample at $index is $s, final is $final $transferRangeMax")
             return final
         }
 
         return if (interpolate) {
             val offset = 1.0f
-
-            val d00 = lerp(diff.x(), density(index), density(toIndex(absoluteCoordsD + Vector3f(offset, 0.0f, 0.0f))))
+            // first step: interpolate between the four parallel sides of the voxel
+            val d00 = lerp(
+                diff.x(),
+                density(index),
+                density(clampedToIndex(absoluteCoordsD + Vector3f(offset, 0.0f, 0.0f)))
+            )
             val d10 = lerp(
                 diff.x(),
-                density(toIndex(absoluteCoordsD + Vector3f(0.0f, offset, 0.0f))),
-                density(toIndex(absoluteCoordsD + Vector3f(offset, offset, 0.0f)))
+                density(clampedToIndex(absoluteCoordsD + Vector3f(0.0f, offset, 0.0f))),
+                density(clampedToIndex(absoluteCoordsD + Vector3f(offset, offset, 0.0f)))
             )
             val d01 = lerp(
                 diff.x(),
-                density(toIndex(absoluteCoordsD + Vector3f(0.0f, 0.0f, offset))),
-                density(toIndex(absoluteCoordsD + Vector3f(offset, 0.0f, offset)))
+                density(clampedToIndex(absoluteCoordsD + Vector3f(0.0f, 0.0f, offset))),
+                density(clampedToIndex(absoluteCoordsD + Vector3f(offset, 0.0f, offset)))
             )
             val d11 = lerp(
                 diff.x(),
-                density(toIndex(absoluteCoordsD + Vector3f(0.0f, offset, offset))),
-                density(toIndex(absoluteCoordsD + Vector3f(offset, offset, offset)))
+                density(clampedToIndex(absoluteCoordsD + Vector3f(0.0f, offset, offset))),
+                density(clampedToIndex(absoluteCoordsD + Vector3f(offset, offset, offset)))
             )
+            // then interpolate between the two calculated lines
             val d0 = lerp(diff.y(), d00, d10)
             val d1 = lerp(diff.y(), d01, d11)
+            // lastly, interpolate along the single remaining line
             lerp(diff.z(), d0, d1)
         } else {
             density(index)
         }
     }
 
-
     private fun lerp(t: Float, v0: Float, v1: Float): Float {
         return (1.0f - t) * v0 + t * v1
     }
@@ -250,9 +256,10 @@ class BufferedVolume(val ds: VolumeDataSource.RAISource<*>, options: VolumeViewe
      * Returns the list of samples (which might include `null` values in case a sample failed),
      * as well as the delta used along the ray, or null if the start/end coordinates are invalid.
      */
-
     override fun sampleRay(start: Vector3f, end: Vector3f): Pair<List<Float?>, Vector3f>? {
         val dimensions = Vector3f(getDimensions())
+        val rayStart = start / dimensions
+        val rayEnd = end / dimensions
         if (start.x() < 0.0f || start.x() > 1.0f || start.y() < 0.0f || start.y() > 1.0f || start.z() < 0.0f || start.z() > 1.0f) {
             logger.debug("Invalid UV coords for ray start: {} -- will clamp values to [0.0, 1.0].", start)
         }
@@ -262,15 +269,15 @@ class BufferedVolume(val ds: VolumeDataSource.RAISource<*>, options: VolumeViewe
         }
 
         val startClamped = Vector3f(
-                min(max(start.x(), 0.0f), 1.0f),
-                min(max(start.y(), 0.0f), 1.0f),
-                min(max(start.z(), 0.0f), 1.0f)
+                min(max(rayStart.x(), 0.0f), 1.0f),
+                min(max(rayStart.y(), 0.0f), 1.0f),
+                min(max(rayStart.z(), 0.0f), 1.0f)
         )
 
         val endClamped = Vector3f(
-                min(max(end.x(), 0.0f), 1.0f),
-                min(max(end.y(), 0.0f), 1.0f),
-                min(max(end.z(), 0.0f), 1.0f)
+                min(max(rayEnd.x(), 0.0f), 1.0f),
+                min(max(rayEnd.y(), 0.0f), 1.0f),
+                min(max(rayEnd.z(), 0.0f), 1.0f)
         )
 
         val direction = (endClamped - startClamped)

src/main/kotlin/graphics/scenery/volumes/RAIVolume.kt

@@ -9,12 +9,15 @@ import graphics.scenery.Origin
 import graphics.scenery.utils.extensions.minus
 import graphics.scenery.utils.extensions.plus
 import graphics.scenery.utils.extensions.times
+import graphics.scenery.utils.extensions.toFloatArray
+import net.imglib2.realtransform.AffineTransform3D
 import net.imglib2.type.numeric.NumericType
 import net.imglib2.type.numeric.integer.*
 import net.imglib2.type.numeric.real.FloatType
 import org.joml.Matrix4f
 import org.joml.Vector3f
 import org.joml.Vector3i
+import kotlin.math.abs
 import kotlin.math.floor
 import kotlin.math.roundToInt
 
@@ -44,9 +47,11 @@ class RAIVolume(@Transient val ds: VolumeDataSource, options: VolumeViewerOption
     }
 
     override fun generateBoundingBox(): OrientedBoundingBox {
-        return OrientedBoundingBox(this,
-            Vector3f(-0.0f, -0.0f, -0.0f),
-            Vector3f(getDimensions()))
+        val scale = getVoxelScale() ?: Vector3f(0f,0f,0f)
+        return OrientedBoundingBox(
+            this, Vector3f(-0.0f, -0.0f, -0.0f),
+            Vector3f(getDimensions()) * scale
+        )
     }
 
     override fun localScale(): Vector3f {
@@ -76,6 +81,22 @@ class RAIVolume(@Transient val ds: VolumeDataSource, options: VolumeViewerOption
         }
     }
 
+    /** Return the dimensions of a voxel as [Vector3f]. Isotropic volumes will return Vector3f(1.0f),
+     * anisotropic volumes will return their scaling factors. */
+    fun getVoxelScale(): Vector3f? {
+        val d = firstSource()?.spimSource?.voxelDimensions
+        if (d != null) {
+            val dims = d.dimensionsAsDoubleArray()
+            return Vector3f(
+                dims[0].toFloat(),
+                dims[1].toFloat(),
+                dims[2].toFloat()
+            )
+        } else {
+            return null
+        }
+    }
+
     override fun createSpatial(): VolumeSpatial {
         return RAIVolumeSpatial(this)
     }
@@ -115,8 +136,10 @@ class RAIVolume(@Transient val ds: VolumeDataSource, options: VolumeViewerOption
         val d = getDimensions()
         val dimensions = Vector3f(d.x.toFloat(), d.y.toFloat(), d.z.toFloat())
 
-        val start = rayStart
-        val end = rayEnd
+        val voxelDims = getVoxelScale() ?: return null
+
+        val start = rayStart / (dimensions * voxelDims )
+        val end = rayEnd / (dimensions * voxelDims )
 
         if (start.x !in 0.0f..1.0f || start.y !in 0.0f..1.0f || start.z !in 0.0f..1.0f) {
             logger.debug("Invalid UV coords for ray start: {} -- will clamp values to [0.0, 1.0].", start)
@@ -158,13 +181,114 @@ class RAIVolume(@Transient val ds: VolumeDataSource, options: VolumeViewerOption
         else -> 1.0f
     }
 
+
+    override fun sampleRayGridTraversal(rayStart: Vector3f, rayEnd: Vector3f): Pair<List<Float?>, List<Vector3f?>> {
+        val d = getDimensions()
+        val dimensions = Vector3f(d.x.toFloat(), d.y.toFloat(), d.z.toFloat())
+        val voxelDimArray = firstSource()!!.spimSource.voxelDimensions.dimensionsAsDoubleArray()
+        // this contains the anisotropic scaling factors
+        val voxelDims = Vector3f(
+            voxelDimArray[0].toFloat(),
+            voxelDimArray[1].toFloat(),
+            voxelDimArray[2].toFloat()
+        )
+        // for the traversal we assume that every voxel has identical dimensions
+        val voxelSize = Vector3f(1f)
+
+        // the start and end points handed over by the AABB test assume isotropic scaling,
+        // so we need to scale it down to the actual dimensions (start and end points would exceed the domain otherwise)
+        val start = rayStart / voxelDims
+        val end = rayEnd / voxelDims
+
+        // clamp it just in case
+        val startClamped = Vector3f(
+            start.x.coerceIn(0.0f, dimensions.x),
+            start.y.coerceIn(0.0f, dimensions.y),
+            start.z.coerceIn(0.0f, dimensions.z),
+        )
+
+        val endClamped = Vector3f(
+            end.x.coerceIn(0.0f, dimensions.x),
+            end.y.coerceIn(0.0f, dimensions.y),
+            end.z.coerceIn(0.0f, dimensions.z),
+        )
+
+        val ray = endClamped - startClamped
+        val rayDir = Vector3f(ray).normalize()
+        val rayLength = ray.length()
+
+        // determine the initial grid direction
+        val stepX = if (rayDir.x >= 0) 1 else -1
+        val stepY = if (rayDir.y >= 0) 1 else -1
+        val stepZ = if (rayDir.z >= 0) 1 else -1
+
+        val tDeltaX = abs(voxelSize.x / rayDir.x)
+        val tDeltaY = abs(voxelSize.y / rayDir.y)
+        val tDeltaZ = abs(voxelSize.z / rayDir.z)
+
+        // this is where it all started
+        val voxelPos = Vector3f(
+            floor(startClamped.x),
+            floor(startClamped.y),
+            floor(startClamped.z)
+        )
+        logger.debug("starting with voxel pos $voxelPos")
+        val tMax = Vector3f(
+            if (stepX > 0) ((voxelPos.x + 1) * voxelSize.x - rayStart.x) / rayDir.x
+            else (voxelPos.x * voxelSize.x - rayStart.x) / rayDir.x,
+            if (stepY > 0) ((voxelPos.y + 1) * voxelSize.y - rayStart.y) / rayDir.y
+            else (voxelPos.y * voxelSize.y - rayStart.y) / rayDir.y,
+            if (stepZ > 0) ((voxelPos.z + 1) * voxelSize.z - rayStart.z) / rayDir.z
+            else (voxelPos.z * voxelSize.z - rayStart.z) / rayDir.z
+        )
+
+        val samplesList = mutableListOf<Float?>()
+        val samplesPosList = mutableListOf<Vector3f>()
+
+        // Start traversing the grid, with t being the already traversed length
+        var t = 0f
+        while (true) {
+            val currentPos = startClamped + rayDir * t
+            // For sampling, we need coordinates in UV space, so we normalize
+            val sampleValue = sample(Vector3f(currentPos).div(dimensions), false)
+            samplesList.add(sampleValue)
+            samplesPosList.add(currentPos)
+
+            // traversal has finished if the accumulated length exceeds the total ray length
+            if ((currentPos - rayStart).length() > rayLength) break
+
+            // decide the voxel direction to travel to next
+            if (tMax.x <= tMax.y && tMax.x <= tMax.z) {
+                t = tMax.x
+                voxelPos.x += stepX
+                tMax.x += tDeltaX
+            }
+            if (tMax.y <= tMax.z && tMax.y <= tMax.x) {
+                t = tMax.y
+                voxelPos.y += stepY
+                tMax.y += tDeltaY
+            }
+            if (tMax.z <= tMax.x && tMax.z <= tMax.y) {
+                t = tMax.z
+                voxelPos.z += stepZ
+                tMax.z += tDeltaZ
+            }
+        }
+
+        return Pair(samplesList, samplesPosList)
+    }
+
     /**
     This sample function is not finished yet, transferRangeMax function should be improved to fit different data type
      **/
     override fun sample(uv: Vector3f, interpolate: Boolean): Float? {
-         val d = getDimensions()
-
-        val absoluteCoords = Vector3f(uv.x() * d.x(), uv.y() * d.y(), uv.z() * d.z())
+        val d = getDimensions()
+        val clampedUV = Vector3f(
+            uv.x.coerceIn(0f, 1f),
+            uv.y.coerceIn(0f, 1f),
+            uv.z.coerceIn(0f, 1f)
+        )
+        val absoluteCoords = Vector3f(clampedUV.x() * d.x(), clampedUV.y() * d.y(), clampedUV.z() * d.z())
         val absoluteCoordsD = Vector3i(floor(absoluteCoords.x()).toInt(), floor(absoluteCoords.y()).toInt(), floor(absoluteCoords.z()).toInt())
 
         val r = when(ds) {
@@ -183,12 +307,14 @@ class RAIVolume(@Transient val ds: VolumeDataSource, options: VolumeViewerOption
         r.setPosition(absoluteCoordsD.z(),2)
 
         val value = r.get()
+        val finalresult: Any
 
-        val finalresult = when(value) {
+        finalresult = when (value) {
             is UnsignedShortType -> value.realFloat
             else -> throw java.lang.IllegalStateException("Can't determine density for ${value?.javaClass} data")
         }
 
+
         val transferRangeMax = when(ds)
         {
             is VolumeDataSource.RAISource<*> -> ds.converterSetups.firstOrNull()?.displayRangeMax?.toFloat()?:ds.type.maxValue()

src/main/kotlin/graphics/scenery/volumes/Volume.kt

@@ -636,6 +636,12 @@ open class Volume(
         return null
     }
 
+    /** This method traverses the grid voxel to voxel, using the algorithm by Amanatides and Woo.
+     * It returns a pair of two lists with the volume samples as [Float] and positions as [Vector3f] each. */
+    open fun sampleRayGridTraversal(rayStart: Vector3f, rayEnd: Vector3f): Pair<List<Float?>, List<Vector3f?>>? {
+        return null
+    }
+
     /**
      * Returns the volume's physical (voxel) dimensions.
      */