Merge pull request #295 from favreau/master

Added Streamlines geometry to OptiX 6 engine

README.md

@@ -252,7 +252,9 @@ In the context of brain simulation, detecting touches between neurons is a essen
 
 When a ray hits a geometry, a random secondary ray is sent in a direction belonging to an hemisphere defined by the normal to the surface. If that secondary ray hits another geometry, the distance between the initial hit and the new intersection is computed, and the corresponding color is assigned to the pixel. By default, red is for short distances (including touches), and green for longer ones. The notion of short and long is defined in the settings of the renderer.
 
-### White matter
+### White matter and tractography
+
+#### White matter
 
 White matter is composed of bundles, which connect various grey matter areas (the locations of nerve cell bodies) of the brain to each other, and carry nerve impulses between neurons. Myelin acts as an insulator, which allows electrical signals to jump, rather than coursing through the axon, increasing the speed of transmission of all nerve signals.
 
@@ -277,6 +279,30 @@ The white matter feature contributed to generating images for the [null model of
       </a>
 </div>
 
+##### Tractography
+
+Tractography is a neuroimaging technique used in the field of neuroscience and medical imaging to visualize and study the white matter pathways in the human brain. White matter consists of bundles of nerve fibers (axons) that connect different regions of the brain and enable communication between them. Tractography helps researchers and clinicians map and understand these complex neural pathways.
+
+The basic principle of tractography involves tracking the diffusion of water molecules in brain tissue. This technique is often referred to as diffusion-weighted magnetic resonance imaging (DW-MRI) or diffusion tensor imaging (DTI). Here's how it works:
+
+1. Diffusion of Water Molecules: Water molecules in brain tissue naturally diffuse in a preferred direction along the axonal fibers in white matter. This diffusion is hindered by cellular structures, membranes, and myelin sheaths surrounding axons.
+
+2. Imaging Process: During a DW-MRI scan, multiple images are acquired with different gradients of magnetic fields, which allows for the measurement of water diffusion in multiple directions within each voxel (3D pixel) of the brain.
+
+3. Data Analysis: The data from DW-MRI is processed to calculate a diffusion tensor, which provides information about the direction and magnitude of water diffusion in each voxel. From this information, researchers can infer the orientation of axonal fibers in that region.
+
+4. Tractography Reconstruction: Tractography algorithms use the diffusion tensor information to reconstruct and visualize the trajectories of white matter pathways throughout the brain. This results in colorful images and 3D representations that show the connections between different brain regions.
+
+Tractography has numerous applications in neuroscience and clinical practice. Researchers use it to study brain connectivity, understand the organization of neural circuits, and investigate neurological disorders. In the clinical setting, tractography can help surgeons plan brain surgeries, evaluate the extent of brain lesions, and assess the impact of brain injuries or diseases on white matter pathways. It has also been used to study conditions like Alzheimer's disease, multiple sclerosis, and stroke.
+
+Different variations of tractography algorithms have been developed to improve the accuracy of tracking neural pathways, such as deterministic tractography and probabilistic tractography. Overall, tractography is a valuable tool for advancing our understanding of brain structure and function.
+
+![___](./bioexplorer/pythonsdk/doc/source/images/tractography.png)
+
+A [Python notebook example](./bioexplorer/pythonsdk/notebooks/tractography/BioExplorer_ISMRM_2015_Tracto_challenge_ground_truth_bundles_TRK_v2.ipynb) demonstrates how to download streamlines from the [Tractography Challenge ISMRM 2015](https://zenodo.org/record/572345) dataset and visualize them in the _BBBE_.
+
+References:
+* [Tractography Challenge ISMRM 2015](https://zenodo.org/record/572345)
 
 ### Enzyme reactions
 

bioexplorer/backend/science/common/SDFGeometries.cpp

@@ -26,6 +26,7 @@
 #include <science/common/UniqueId.h>
 #include <science/common/Utils.h>
 
+#include <platform/core/common/utils/Utils.h>
 #include <platform/core/engineapi/Material.h>
 #include <platform/core/engineapi/Model.h>
 

bioexplorer/backend/science/common/Utils.cpp

@@ -418,25 +418,6 @@ bool rayBoxIntersection(const Vector3d& origin, const Vector3d& direction, const
     return (tmin < t1 && tmax > t0);
 }
 
-Vector4f getBezierPoint(const Vector4fs& controlPoints, const double t)
-{
-    if (t < 0.0 || t > 1.0)
-        PLUGIN_THROW("Invalid value with t=" + std::to_string(t) + ". Must be between 0 and 1");
-    const uint64_t nbControlPoints = controlPoints.size();
-    // 3D points
-    Vector3fs points;
-    points.reserve(nbControlPoints);
-    for (const auto& controlPoint : controlPoints)
-        points.push_back({controlPoint.x, controlPoint.y, controlPoint.z});
-    for (int64_t i = nbControlPoints - 1; i >= 0; --i)
-        for (uint64_t j = 0; j < i; ++j)
-            points[j] += t * (points[j + 1] - points[j]);
-
-    // Radius
-    const double radius = controlPoints[floor(t * double(nbControlPoints))].w;
-    return Vector4f(points[0].x, points[0].y, points[0].z, radius);
-}
-
 double sphereVolume(const double radius)
 {
     return 4.0 * M_PI * pow(radius, 3) / 3.0;

bioexplorer/backend/science/common/Utils.h

@@ -194,17 +194,6 @@ Transformation combineTransformations(const Transformations& transformations);
 bool rayBoxIntersection(const Vector3d& origin, const Vector3d& direction, const Boxd& box, const double t0,
                         const double t1, double& t);
 
-/**
- * @brief Get the Bezier Point from a curve defined by the provided control
- * points
- *
- * @param controlPoints Curve control points with radius
- * @param t The t in the function for a curve can be thought of as describing
- * how far B(t) is from first to last control point.
- * @return Vector3f
- */
-Vector4f getBezierPoint(const Vector4fs& controlPoints, const double t);
-
 // Volumes
 double sphereVolume(const double radius);
 double cylinderVolume(const double height, const double radius);

bioexplorer/pythonsdk/notebooks/tractography/BioExplorer_ISMRM_2015_Tracto_challenge_ground_truth_bundles_TRK_v2.ipynb

@@ -0,0 +1,145 @@
+{
+    "cells": [
+        {
+            "cell_type": "markdown",
+            "metadata": {},
+            "source": [
+                "# Blue Brain BioExplorer - Tractography Challenge ISMRM 2015\n",
+                "\n",
+                "https://zenodo.org/record/572345\n",
+                "\n",
+                "![](../bioexplorer_white_matter_banner.png)"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "metadata": {},
+            "source": [
+                "### Prerequeries"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "!pip install dipy"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "metadata": {},
+            "source": [
+                "### Connect to BioExplorer backend"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "from bioexplorer import BioExplorer\n",
+                "\n",
+                "be = BioExplorer('localhost:5000')\n",
+                "core = be.core_api()\n",
+                "be.reset_scene()"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "metadata": {},
+            "source": [
+                "### Dataset"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "import os, fnmatch\n",
+                "\n",
+                "home_folder = '<TO BE DEFINED>'\n",
+                "folder = os.path.join(home_folder, 'ISMRM_2015_Tracto_challenge_ground_truth_bundles_TRK_v2')\n",
+                "files = fnmatch.filter(os.listdir(folder), '*.trk')\n",
+                "files.sort()"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "def add_streamlines(name, streams, radius):\n",
+                "    ''' Build streamlines'''\n",
+                "    count = 0\n",
+                "    indices = list()\n",
+                "    vertices = list()\n",
+                "    colors = list()\n",
+                "    for stream in streams:\n",
+                "        l = len(stream)\n",
+                "        for coordinates in stream:\n",
+                "            for coordinate in coordinates:\n",
+                "                vertices.append(float(coordinate))\n",
+                "            vertices.append(radius)\n",
+                "        count = count + l\n",
+                "        indices.append(count)\n",
+                "    '''Send streamlines to the BioExplorer'''\n",
+                "    return be.add_streamlines(name, indices, vertices, colors)"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "from dipy.io.streamline import load_trk\n",
+                "for file in files:\n",
+                "    streams = load_trk(os.path.join(folder, file), reference='same', bbox_valid_check=False)\n",
+                "    streamlines = streams.streamlines\n",
+                "    add_streamlines(file, streamlines, radius=0.1)"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "status = be.reset_camera()\n",
+                "status = core.set_renderer(current='advanced')"
+            ]
+        }
+    ],
+    "metadata": {
+        "kernelspec": {
+            "display_name": "Python 3.8.10 64-bit ('env': venv)",
+            "language": "python",
+            "name": "python3"
+        },
+        "language_info": {
+            "codemirror_mode": {
+                "name": "ipython",
+                "version": 3
+            },
+            "file_extension": ".py",
+            "mimetype": "text/x-python",
+            "name": "python",
+            "nbconvert_exporter": "python",
+            "pygments_lexer": "ipython3",
+            "version": "3.10.12"
+        },
+        "vscode": {
+            "interpreter": {
+                "hash": "ea9a5fa46eb6bad2806a8ea1d08e15bb1e255a2d4320b81e765591579963c56b"
+            }
+        }
+    },
+    "nbformat": 4,
+    "nbformat_minor": 4
+}

extensions/medicalimaging/dicom/module/cuda/renderer/DICOM.cu

@@ -32,15 +32,21 @@ using namespace optix;
 
 rtDeclareVariable(float, surfaceOffset, , );
 
-static __device__ void dicomShade()
+static __device__ void dicomShade(const bool textured)
 {
     float3 result = make_float3(0.f);
 
     float4 voxelColor = make_float4(Kd, luminance(Ko));
-    if (volume_map != 0)
+    if (textured)
     {
-        const float voxelValue = rtTex3D<float>(volume_map, texcoord3d.x, texcoord3d.y, texcoord3d.z);
-        voxelColor = calcTransferFunctionColor(transfer_function_map, value_range, voxelValue);
+        if (volume_map != 0)
+        {
+            const float voxelValue = rtTex3D<float>(volume_map, texcoord3d.x, texcoord3d.y, texcoord3d.z);
+            voxelColor = calcTransferFunctionColor(transfer_function_map, value_range, voxelValue);
+        }
+        else
+            voxelColor = make_float4(make_float3(optix::rtTex2D<float4>(albedoMetallic_map, texcoord.x, texcoord.y)),
+                                     luminance(Ko));
     }
 
     const float3 hit_point = ray.origin + t_hit * ray.direction;
@@ -262,12 +268,12 @@ RT_PROGRAM void any_hit_shadow()
 
 RT_PROGRAM void closest_hit_radiance()
 {
-    dicomShade();
+    dicomShade(false);
 }
 
 RT_PROGRAM void closest_hit_radiance_textured()
 {
-    dicomShade();
+    dicomShade(true);
 }
 
 RT_PROGRAM void exception()

platform/core/common/utils/Utils.cpp

@@ -70,4 +70,23 @@ std::string extractExtension(const std::string& filename)
     return extension;
 }
 
+Vector4f getBezierPoint(const Vector4fs& controlPoints, const double t)
+{
+    if (t < 0.0 || t > 1.0)
+        CORE_THROW("Invalid value with t=" + std::to_string(t) + ". Must be between 0 and 1");
+    const uint64_t nbControlPoints = controlPoints.size();
+    // 3D points
+    Vector3fs points;
+    points.reserve(nbControlPoints);
+    for (const auto& controlPoint : controlPoints)
+        points.push_back({controlPoint.x, controlPoint.y, controlPoint.z});
+    for (int64_t i = nbControlPoints - 1; i >= 0; --i)
+        for (uint64_t j = 0; j < i; ++j)
+            points[j] += t * (points[j + 1] - points[j]);
+
+    // Radius
+    const double radius = controlPoints[floor(t * double(nbControlPoints))].w;
+    return Vector4f(points[0].x, points[0].y, points[0].z, radius);
+}
+
 } // namespace core

platform/core/common/utils/Utils.h

@@ -51,4 +51,15 @@ inline std::array<T, M> toArray(const glm::vec<M, T>& input)
     return output;
 }
 
+/**
+ * @brief Get the Bezier Point from a curve defined by the provided control
+ * points
+ *
+ * @param controlPoints Curve control points with radius
+ * @param t The t in the function for a curve can be thought of as describing
+ * how far B(t) is from first to last control point.
+ * @return Vector3f
+ */
+Vector4f getBezierPoint(const Vector4fs& controlPoints, const double t);
+
 } // namespace core

platform/engines/optix6/CMakeLists.txt

@@ -48,6 +48,7 @@ set(${NAME}_CU
   cuda/geometry/Spheres.cu
   cuda/geometry/TriangleMesh.cu
   cuda/geometry/Volumes.cu
+  cuda/geometry/Streamlines.cu
   cuda/renderer/Basic.cu
   cuda/renderer/Advanced.cu
 )

platform/engines/optix6/OptiXCommonStructs.h

@@ -31,6 +31,7 @@ const size_t OPTIX_ENTRY_POINT_COUNT = 1;
 const size_t OPTIX_MAX_TRACE_DEPTH = 31;
 
 const float EPSILON = 1e-2f;
+const size_t MAX_TEXTURE_SIZE = 16384;
 
 struct BasicLight
 {

platform/engines/optix6/OptiXContext.cpp

@@ -32,6 +32,7 @@
 #include <platform/engines/optix6/OptiX6Engine_generated_Cones.cu.ptx.h>
 #include <platform/engines/optix6/OptiX6Engine_generated_Cylinders.cu.ptx.h>
 #include <platform/engines/optix6/OptiX6Engine_generated_Spheres.cu.ptx.h>
+#include <platform/engines/optix6/OptiX6Engine_generated_Streamlines.cu.ptx.h>
 #include <platform/engines/optix6/OptiX6Engine_generated_TriangleMesh.cu.ptx.h>
 #include <platform/engines/optix6/OptiX6Engine_generated_Volumes.cu.ptx.h>
 
@@ -45,6 +46,7 @@ const std::string CUDA_CYLINDERS = OptiX6Engine_generated_Cylinders_cu_ptx;
 const std::string CUDA_CONES = OptiX6Engine_generated_Cones_cu_ptx;
 const std::string CUDA_TRIANGLES_MESH = OptiX6Engine_generated_TriangleMesh_cu_ptx;
 const std::string CUDA_VOLUMES = OptiX6Engine_generated_Volumes_cu_ptx;
+const std::string CUDA_STREAMLINES = OptiX6Engine_generated_Streamlines_cu_ptx;
 
 template <typename T>
 T white();
@@ -367,6 +369,11 @@ void OptiXContext::_initialize()
     _intersects[OptixGeometryType::volume] =
         _optixContext->createProgramFromPTXString(CUDA_VOLUMES, CUDA_FUNC_INTERSECTION);
 
+    _bounds[OptixGeometryType::streamline] =
+        _optixContext->createProgramFromPTXString(CUDA_STREAMLINES, CUDA_FUNC_BOUNDS);
+    _intersects[OptixGeometryType::streamline] =
+        _optixContext->createProgramFromPTXString(CUDA_STREAMLINES, CUDA_FUNC_INTERSECTION);
+
     // Exceptions
     _optixContext[CONTEXT_EXCEPTION_BAD_COLOR]->setFloat(1.0f, 0.0f, 0.0f, 1.f);
 

platform/engines/optix6/OptiXContext.h

@@ -157,7 +157,8 @@ enum class OptixGeometryType
     cone,
     cylinder,
     triangleMesh,
-    volume
+    volume,
+    streamline
 };
 
 struct OptixShaderProgram