From 25b0582a34f3d202d6e54f901eeced5c915e4748 Mon Sep 17 00:00:00 2001 From: Cyrille Pierre Henri Favreau Date: Fri, 1 Dec 2023 19:43:47 +0100 Subject: [PATCH] Improved neurons representation using SDF displacement --- .../backend/science/common/SDFGeometries.cpp | 65 ++++--- .../science/common/ThreadSafeContainer.cpp | 15 +- .../science/common/ThreadSafeContainer.h | 2 + .../backend/science/morphologies/Neurons.cpp | 172 +++++++++++------- .../backend/science/morphologies/Neurons.h | 9 +- .../science/morphologies/NeuronsLoader.cpp | 14 +- .../neurons/BioExplorer_neurons.ipynb | 59 ++++-- .../utils/BioExplorer_sdf_testing.ipynb | 45 ++++- platform/core/engineapi/Model.h | 14 +- platform/engines/optix6/OptiXModel.cpp | 99 +++++----- platform/engines/optix6/OptiXModel.h | 1 + .../optix6/cuda/geometry/SDFGeometries.cu | 81 ++++++--- .../ospray/ispc/geometry/SDFGeometries.ispc | 79 ++------ 13 files changed, 390 insertions(+), 265 deletions(-) diff --git a/bioexplorer/backend/science/common/SDFGeometries.cpp b/bioexplorer/backend/science/common/SDFGeometries.cpp index 339662216..8a130e3d8 100644 --- a/bioexplorer/backend/science/common/SDFGeometries.cpp +++ b/bioexplorer/backend/science/common/SDFGeometries.cpp @@ -88,30 +88,51 @@ void SDFGeometries::addSDFDemo(Model& model) const Vector3f displacement{0.05f, 10.f, 0.f}; ThreadSafeContainer modelContainer(model, 0.f, Vector3d(), Quaterniond()); + if (rand() % 2 == 0) - for (size_t materialId = 0; materialId < 10; ++materialId) - { - const float x = materialId * 3.0f; - Neighbours neighbours; - neighbours.insert(modelContainer.addSphere(Vector3f(0.f + x, 0.f, 0.f), 0.5f, materialId, useSdf, NO_USER_DATA, - neighbours, displacement)); - neighbours.insert(modelContainer.addCone(Vector3f(-1.f + x, 0.f, 0.f), 0.25, Vector3d(0.f + x, 0.f, 0.f), 0.1f, - materialId, useSdf, NO_USER_DATA, neighbours, displacement)); - neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, 0.f, 0.f), 0.1, Vector3f(1.f + x, 0.f, 0.f), 0.25f, - materialId, useSdf, NO_USER_DATA, neighbours, displacement)); - neighbours.insert(modelContainer.addSphere(Vector3f(-0.5 + x, 0.f, 0.f), 0.25f, materialId, useSdf, - NO_USER_DATA, neighbours, displacement)); - neighbours.insert(modelContainer.addSphere(Vector3f(0.5 + x, 0.f, 0.f), 0.25f, materialId, useSdf, NO_USER_DATA, - neighbours, displacement)); - neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, 0.25, 0.f), 0.5f, Vector3f(0.f + x, 1.f, 0.f), 0.f, - materialId, useSdf, NO_USER_DATA, neighbours, displacement)); - neighbours.insert(modelContainer.addTorus(Vector3f(0.f + x, 0.f, 0.f), 1.5f, 0.5f, materialId, NO_USER_DATA, - neighbours, displacement)); - neighbours.insert(modelContainer.addCutSphere(Vector3f(0.f + x, 1.f, 0.f), 1.0f, 0.5f, materialId, NO_USER_DATA, + for (size_t materialId = 0; materialId < 10; ++materialId) + { + const float x = materialId * 3.0f; + Neighbours neighbours; + neighbours.insert(modelContainer.addSphere(Vector3f(0.f + x, 0.f, 0.f), 0.5f, materialId, useSdf, + NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(-1.f + x, 0.f, 0.f), 0.25, Vector3d(0.f + x, 0.f, 0.f), + 0.1f, materialId, useSdf, NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, 0.f, 0.f), 0.1, Vector3f(1.f + x, 0.f, 0.f), + 0.25f, materialId, useSdf, NO_USER_DATA, neighbours, + displacement)); + neighbours.insert(modelContainer.addSphere(Vector3f(-0.5 + x, 0.f, 0.f), 0.25f, materialId, useSdf, + NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addSphere(Vector3f(0.5 + x, 0.f, 0.f), 0.25f, materialId, useSdf, + NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, 0.25, 0.f), 0.5f, Vector3f(0.f + x, 1.f, 0.f), + 0.f, materialId, useSdf, NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addTorus(Vector3f(0.f + x, 0.f, 0.f), 1.5f, 0.5f, materialId, NO_USER_DATA, neighbours, displacement)); - neighbours.insert(modelContainer.addVesica(Vector3f(0.f + x, -1.f, 0.f), Vector3f(0.f + x, -1.5f, 0.f), 1.0f, - materialId, NO_USER_DATA, neighbours, displacement)); - } + neighbours.insert(modelContainer.addCutSphere(Vector3f(0.f + x, 1.f, 0.f), 1.0f, 0.5f, materialId, + NO_USER_DATA, neighbours, displacement)); + neighbours.insert(modelContainer.addVesica(Vector3f(0.f + x, -1.f, 0.f), Vector3f(0.f + x, -1.5f, 0.f), + 1.0f, materialId, NO_USER_DATA, neighbours, displacement)); + } + else + for (size_t materialId = 0; materialId < 10; ++materialId) + { + const float x = materialId * 3.0f; + Neighbours neighbours; + neighbours.insert(modelContainer.addSphere(Vector3f(0.f + x, 0.f, 0.f), 1.0f, materialId, useSdf, + NO_USER_DATA, {}, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(0.5f + x, 0.f, 0.f), 0.75f, Vector3f(2.f + x, 0.f, 0.f), + 0.f, materialId, useSdf, NO_USER_DATA, {}, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(-0.5f + x, 0.f, 0.f), 0.75f, Vector3f(-2.f + x, 0.f, 0.f), + 0.f, materialId, useSdf, NO_USER_DATA, {}, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, 0.5f, 0.f), 0.75f, Vector3f(0.f + x, 2.f, 0.f), + 0.f, materialId, useSdf, NO_USER_DATA, {}, displacement)); + neighbours.insert(modelContainer.addCone(Vector3f(0.f + x, -0.5f, 0.f), 0.75f, Vector3f(0.f + x, -2.f, 0.f), + 0.f, materialId, useSdf, NO_USER_DATA, {}, displacement)); + + for (const auto index : neighbours) + modelContainer.setSDFGeometryNeighbours(index, neighbours); + } modelContainer.commitToModel(); model.applyDefaultColormap(); } diff --git a/bioexplorer/backend/science/common/ThreadSafeContainer.cpp b/bioexplorer/backend/science/common/ThreadSafeContainer.cpp index 274c57a4c..2c8cc72e5 100644 --- a/bioexplorer/backend/science/common/ThreadSafeContainer.cpp +++ b/bioexplorer/backend/science/common/ThreadSafeContainer.cpp @@ -25,6 +25,8 @@ #include "Utils.h" +#include + #include #include @@ -199,11 +201,11 @@ uint64_t ThreadSafeContainer::_addCone(const size_t materialId, const Cone& cone return 0; // Only used by SDF geometry } -uint64_t ThreadSafeContainer::_addSDFGeometry(const size_t materialId, const SDFGeometry& geom, +uint64_t ThreadSafeContainer::_addSDFGeometry(const size_t materialId, const SDFGeometry& geometry, const std::set& neighbours) { const uint64_t geometryIndex = _sdfMorphologyData.geometries.size(); - _sdfMorphologyData.geometries.push_back(geom); + _sdfMorphologyData.geometries.push_back(geometry); _sdfMorphologyData.neighbours.push_back(neighbours); _sdfMorphologyData.materials.push_back(materialId); return geometryIndex; @@ -356,5 +358,14 @@ void ThreadSafeContainer::_commitStreamlinesToModel() } _streamlinesMap.clear(); } + +void ThreadSafeContainer::setSDFGeometryNeighbours(const uint64_t geometryIndex, const std::set& neighbours) +{ + if (geometryIndex >= _sdfMorphologyData.neighbours.size()) + PLUGIN_THROW("Invalid SDF geometry Id"); + auto n = neighbours; + n.erase(geometryIndex); + _sdfMorphologyData.neighbours[geometryIndex] = n; +} } // namespace common } // namespace bioexplorer diff --git a/bioexplorer/backend/science/common/ThreadSafeContainer.h b/bioexplorer/backend/science/common/ThreadSafeContainer.h index 6e8d1e360..c36d76c77 100644 --- a/bioexplorer/backend/science/common/ThreadSafeContainer.h +++ b/bioexplorer/backend/science/common/ThreadSafeContainer.h @@ -167,6 +167,8 @@ class ThreadSafeContainer */ void commitToModel(); + void setSDFGeometryNeighbours(const uint64_t geometryIndex, const std::set& neighbours); + MaterialSet& getMaterialIds() { return _materialIds; } private: diff --git a/bioexplorer/backend/science/morphologies/Neurons.cpp b/bioexplorer/backend/science/morphologies/Neurons.cpp index 1bc3ff8df..0776cb41a 100644 --- a/bioexplorer/backend/science/morphologies/Neurons.cpp +++ b/bioexplorer/backend/science/morphologies/Neurons.cpp @@ -60,7 +60,9 @@ namespace morphology { const uint64_t NB_MYELIN_FREE_SEGMENTS = 4; const double DEFAULT_ARROW_RADIUS_RATIO = 10.0; +const uint64_t DEFAULT_DEBUG_SYNAPSE_DENSITY_RATIO = 5; const double MAX_SOMA_RADIUS = 10.0; + const Vector2d DEFAULT_SIMULATION_VALUE_RANGE = {-80.0, -10.0}; std::map reportTypeAsString = {{ReportType::undefined, "undefined"}, @@ -315,6 +317,8 @@ void Neurons::_buildModel(const LoaderProgress& callback) container.commitToModel(); } model->applyDefaultColormap(); + if (_details.simulationReportId != -1) + setDefaultTransferFunction(*model, DEFAULT_SIMULATION_VALUE_RANGE); ModelMetadata metadata = {{"Number of Neurons", std::to_string(somas.size())}, {"Number of Spines", std::to_string(_nbSpines)}, @@ -411,30 +415,32 @@ SectionSynapseMap Neurons::_buildDebugSynapses(const uint64_t neuronId, const Se if (static_cast(section.second.type) == NeuronSectionType::axon) continue; - const auto& points = section.second.points; + // Process points according to representation + const auto points = _getProcessedSectionPoints(_details.morphologyRepresentation, section.second.points); double sectionLength = 0.0; doubles segmentEnds; for (size_t i = 0; i < points.size() - 1; ++i) { - const auto segmentLength = length(points[i + 1] - points[i]); - segmentEnds.push_back(sectionLength); + const double segmentLength = length(points[i + 1] - points[i]); sectionLength += segmentLength; + segmentEnds.push_back(sectionLength); } - const size_t potentialNumberOfSynapses = sectionLength / DEFAULT_SPINE_RADIUS + 1; - const size_t effectiveNumberOfSynapses = potentialNumberOfSynapses * (0.5 + 0.5 * rnd0()); + const size_t potentialNumberOfSynapses = + DEFAULT_DEBUG_SYNAPSE_DENSITY_RATIO * sectionLength / DEFAULT_SPINE_RADIUS + 1; + size_t effectiveNumberOfSynapses = potentialNumberOfSynapses * (0.5 + 0.5 * rnd0()); for (size_t i = 0; i < effectiveNumberOfSynapses; ++i) { const double distance = rnd0() * sectionLength; size_t segmentId = 0; - while (segmentEnds[segmentId] < distance && segmentId < segmentEnds.size()) + while (distance > segmentEnds[segmentId] && segmentId < segmentEnds.size()) ++segmentId; const auto preSynapticSectionId = section.first; const auto preSynapticSegmentId = segmentId; Synapse synapse; synapse.type = (rand() % 2 == 0 ? MorphologySynapseType::afferent : MorphologySynapseType::efferent); - synapse.preSynapticSegmentDistance = segmentEnds[segmentId] - distance; + synapse.preSynapticSegmentDistance = distance - (segmentId > 0 ? segmentEnds[segmentId - 1] : 0.f); synapse.postSynapticNeuronId = neuronId; synapse.postSynapticSectionId = 0; synapse.postSynapticSegmentId = 0; @@ -533,28 +539,6 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne } } - uint64_t somaGeometryIndex = 0; - if (_details.loadSomas) - { - if (_details.showMembrane) - { - const bool useSdf = andCheck(static_cast(_details.realismLevel), - static_cast(MorphologyRealismLevel::soma)); - somaGeometryIndex = - container.addSphere(somaPosition, somaRadius, somaMaterialId, useSdf, somaUserData, {}, - Vector3f(_getDisplacementValue(DisplacementElement::morphology_soma_strength), - _getDisplacementValue(DisplacementElement::morphology_soma_frequency), - 0.f)); - } - if (_details.generateInternals) - { - const bool useSdf = andCheck(static_cast(_details.realismLevel), - static_cast(MorphologyRealismLevel::internals)); - _addSomaInternals(container, baseMaterialId, somaPosition, somaRadius, mitochondriaDensity, useSdf, - _details.radiusMultiplier); - } - } - // Load synapses for all sections SectionSynapseMap synapses; switch (_details.synapsesType) @@ -572,11 +556,11 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne } // Sections (dendrites and axon) - uint64_t geometryIndex = 0; - Neighbours neighbours{somaGeometryIndex}; - + Neighbours somaNeighbours; + double correctedSomaRadius = 0.f; for (const auto& section : sections) { + Neighbours sectionNeighbours; const auto sectionType = static_cast(section.second.type); const auto& points = section.second.points; bool useSdf = andCheck(static_cast(_details.realismLevel), static_cast(sectionType)); @@ -601,6 +585,7 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne continue; } + uint64_t count = 0.0; // Sections connected to the soma if (_details.showMembrane && _details.loadSomas && section.second.parentId == SOMA_AS_PARENT) { @@ -614,8 +599,8 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne useSdf = andCheck(static_cast(_details.realismLevel), static_cast(MorphologyRealismLevel::soma)); - const float srcRadius = _getCorrectedRadius(somaRadius * 0.75f, _details.radiusMultiplier); - const float dstRadius = _getCorrectedRadius(point.w * 0.5f, _details.radiusMultiplier); + const double srcRadius = _getCorrectedRadius(somaRadius * 0.75, _details.radiusMultiplier); + const double dstRadius = _getCorrectedRadius(point.w * 0.5, _details.radiusMultiplier); const auto sectionType = static_cast(section.second.type); const bool loadSection = @@ -626,24 +611,68 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne if (!loadSection) continue; - const auto dst = + const Vector3d dst = _animatedPosition(Vector4d(somaPosition + somaRotation * Vector3d(point), dstRadius), neuronId); const Vector3f displacement = { Vector3f(_getDisplacementValue(DisplacementElement::morphology_soma_strength), _getDisplacementValue(DisplacementElement::morphology_soma_frequency), 0.f)}; - geometryIndex = container.addCone(somaPosition, srcRadius, dst, dstRadius, somaMaterialId, useSdf, - somaUserData, neighbours, displacement); - neighbours.insert(geometryIndex); - geometryIndex = - container.addSphere(dst, dstRadius, somaMaterialId, useSdf, somaUserData, neighbours, displacement); - neighbours.insert(geometryIndex); + + const Vector3d segmentDirection = normalize(lastPoint - firstPoint); + const double halfDistanceToSoma = length(Vector3d(point)) * 0.5f; + const Vector3d p1 = Vector3d(point) - halfDistanceToSoma * segmentDirection; + + const Vector3d p2 = p1 * dstRadius / somaRadius * 0.95; + const auto src = _animatedPosition(Vector4d(somaPosition + somaRotation * p2, srcRadius), neuronId); + + correctedSomaRadius = std::max(correctedSomaRadius, length(p2)) * 2.0; + const uint64_t geometryIndex = container.addCone(src, srcRadius, dst, dstRadius, somaMaterialId, useSdf, + somaUserData, {}, displacement); + somaNeighbours.insert(geometryIndex); + sectionNeighbours.insert(geometryIndex); + if (!useSdf) + container.addSphere(dst, dstRadius, somaMaterialId, useSdf, somaUserData); + ++count; + } + correctedSomaRadius = count == 0 ? somaRadius : correctedSomaRadius / count; + + float parentRadius = section.second.points[0].w; + if (sections.find(section.second.parentId) != sections.end()) + { + const auto& parentSection = sections[section.second.parentId]; + const auto& parentSectionPoints = parentSection.points; + parentRadius = parentSection.points[parentSection.points.size() - 1].w; } if (distanceToSoma <= _details.maxDistanceToSoma) - _addSection(container, neuronId, soma.morphologyId, section.first, section.second, geometryIndex, - somaPosition, somaRotation, somaRadius, baseMaterialId, mitochondriaDensity, somaUserData, - synapses, distanceToSoma, neighbours); + _addSection(container, neuronId, soma.morphologyId, section.first, section.second, somaPosition, + somaRotation, parentRadius, baseMaterialId, mitochondriaDensity, somaUserData, synapses, + distanceToSoma, sectionNeighbours); + } + + if (_details.loadSomas) + { + if (_details.showMembrane) + { + const bool useSdf = andCheck(static_cast(_details.realismLevel), + static_cast(MorphologyRealismLevel::soma)); + somaNeighbours.insert( + container.addSphere(somaPosition, correctedSomaRadius, somaMaterialId, useSdf, somaUserData, {}, + Vector3f(_getDisplacementValue(DisplacementElement::morphology_soma_strength), + _getDisplacementValue(DisplacementElement::morphology_soma_frequency), + 0.f))); + if (useSdf) + for (const auto index : somaNeighbours) + container.setSDFGeometryNeighbours(index, somaNeighbours); + } + if (_details.generateInternals) + { + const bool useSdf = andCheck(static_cast(_details.realismLevel), + static_cast(MorphologyRealismLevel::internals)); + _addSomaInternals(container, baseMaterialId, somaPosition, correctedSomaRadius, mitochondriaDensity, useSdf, + _details.radiusMultiplier); + } } + } // namespace morphology void Neurons::_addVaricosity(Vector4fs& points) @@ -727,11 +756,10 @@ void Neurons::_addArrow(ThreadSafeContainer& container, const uint64_t neuronId, } void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronId, const uint64_t morphologyId, - const uint64_t sectionId, const Section& section, const uint64_t somaGeometryIndex, - const Vector3d& somaPosition, const Quaterniond& somaRotation, const double somaRadius, - const size_t baseMaterialId, const double mitochondriaDensity, const uint64_t somaUserData, - const SectionSynapseMap& synapses, const double distanceToSoma, - const Neighbours& somaNeighbours) + const uint64_t sectionId, const Section& section, const Vector3d& somaPosition, + const Quaterniond& somaRotation, const double parentRadius, const size_t baseMaterialId, + const double mitochondriaDensity, const uint64_t somaUserData, + const SectionSynapseMap& synapses, const double distanceToSoma, const Neighbours& neighbours) { const auto& connector = DBConnector::getInstance(); const auto sectionType = static_cast(section.type); @@ -750,7 +778,9 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI } auto userData = NO_USER_DATA; - const auto& points = section.points; + auto points = section.points; + points[0].w = parentRadius; + size_t sectionMaterialId; switch (_details.morphologyColorScheme) { @@ -778,11 +808,6 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI // Section surface and volume double sectionLength = 0.0; double sectionVolume = 0.0; - uint64_t geometryIndex = 0; - Neighbours neighbours = somaNeighbours; - if (_details.morphologyColorScheme == MorphologyColorScheme::none) - neighbours.insert(somaGeometryIndex); - uint64_ts compartments; switch (_simulationReport.type) { @@ -810,7 +835,10 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI segmentSynapses = (*it).second; // Section points - for (uint64_t i = 0; i < localPoints.size() - 1; ++i) + Neighbours sectionNeighbours = neighbours; + uint64_t previousGeometryIndex = 0; + const uint64_t startIndex = (section.parentId == SOMA_AS_PARENT ? 1 : 0); + for (uint64_t i = startIndex; i < localPoints.size() - 1; ++i) { if (!compartments.empty()) { @@ -832,10 +860,8 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI if (_details.showMembrane) { - if (i > 0 && _details.morphologyRepresentation != MorphologyRepresentation::segment) - neighbours = {geometryIndex}; - - Vector3f displacement{srcRadius * _getDisplacementValue(DisplacementElement::morphology_section_strength), + Vector3f displacement{std::min(std::min(srcRadius, dstRadius) * 0.5f, + _getDisplacementValue(DisplacementElement::morphology_section_strength)), _getDisplacementValue(DisplacementElement::morphology_section_frequency), 0.f}; size_t materialId = _details.morphologyColorScheme == MorphologyColorScheme::distance_to_soma @@ -843,6 +869,7 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI : sectionMaterialId; + // Varicosity (axon only) if (addVaricosity && _details.morphologyColorScheme == MorphologyColorScheme::section_type) { if (i > middlePointIndex && i < middlePointIndex + 3) @@ -851,14 +878,16 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI displacement = Vector3f(2.f * srcRadius * _getDisplacementValue(DisplacementElement::morphology_section_strength), - _getDisplacementValue(DisplacementElement::morphology_section_frequency), 0.f); + _getDisplacementValue(DisplacementElement::morphology_section_frequency) / srcRadius, + 0.f); } if (i == middlePointIndex + 1 || i == middlePointIndex + 3) - neighbours = {}; + sectionNeighbours = {}; if (i == middlePointIndex + 1) _varicosities[neuronId].push_back(dst); } + // Synapses const auto it = segmentSynapses.find(i); if (it != segmentSynapses.end()) { @@ -887,9 +916,19 @@ void Neurons::_addSection(ThreadSafeContainer& container, const uint64_t neuronI if (!useSdf) container.addSphere(dst, dstRadius, materialId, useSdf, userData); - geometryIndex = container.addCone(src, srcRadius, dst, dstRadius, materialId, useSdf, userData, neighbours, - displacement); - neighbours.insert(geometryIndex); +#if 0 + const uint64_t geometryIndex = container.addCone(src, srcRadius, dst, dstRadius, materialId, useSdf, + userData, sectionNeighbours, displacement); + + if (i > 0) + container.setSDFGeometryNeighbours(previousGeometryIndex, {previousGeometryIndex}); +#else + const uint64_t geometryIndex = + container.addCone(src, srcRadius, dst, dstRadius, materialId, useSdf, userData, {}, displacement); + +#endif + previousGeometryIndex = geometryIndex; + sectionNeighbours = {geometryIndex}; // Stop if distance to soma in greater than the specified max value _maxDistanceToSoma = std::max(_maxDistanceToSoma, distanceToSoma + sectionLength); @@ -1174,7 +1213,6 @@ void Neurons::_attachSimulationReport(Model& model) "guids to the simulated ones"); auto handler = std::make_shared(_details.populationName, _details.simulationReportId); model.setSimulationHandler(handler); - setDefaultTransferFunction(model, DEFAULT_SIMULATION_VALUE_RANGE); break; } case ReportType::soma: @@ -1184,7 +1222,6 @@ void Neurons::_attachSimulationReport(Model& model) "to the simulated ones"); auto handler = std::make_shared(_details.populationName, _details.simulationReportId); model.setSimulationHandler(handler); - setDefaultTransferFunction(model, DEFAULT_SIMULATION_VALUE_RANGE); break; } case ReportType::compartment: @@ -1195,7 +1232,6 @@ void Neurons::_attachSimulationReport(Model& model) auto handler = std::make_shared(_details.populationName, _details.simulationReportId); model.setSimulationHandler(handler); - setDefaultTransferFunction(model, DEFAULT_SIMULATION_VALUE_RANGE); break; } } diff --git a/bioexplorer/backend/science/morphologies/Neurons.h b/bioexplorer/backend/science/morphologies/Neurons.h index 31514917b..d44ac1302 100644 --- a/bioexplorer/backend/science/morphologies/Neurons.h +++ b/bioexplorer/backend/science/morphologies/Neurons.h @@ -93,11 +93,10 @@ class Neurons : public Morphologies const double distanceToSoma); void _addSection(common::ThreadSafeContainer& container, const uint64_t neuronId, const uint64_t morphologyId, - const uint64_t sectionId, const Section& section, const uint64_t somaGeometryIndex, - const core::Vector3d& somaPosition, const core::Quaterniond& somaRotation, const double somaRadius, - const size_t baseMaterialId, const double mitochondriaDensity, const uint64_t somaUserData, - const SectionSynapseMap& synapses, const double distanceToSoma, - const common::Neighbours& somaNeighbours); + const uint64_t sectionId, const Section& section, const core::Vector3d& somaPosition, + const core::Quaterniond& somaRotation, const double somaRadius, const size_t baseMaterialId, + const double mitochondriaDensity, const uint64_t somaUserData, const SectionSynapseMap& synapses, + const double distanceToSoma, const common::Neighbours& somaNeighbours); void _addSpine(common::ThreadSafeContainer& container, const uint64_t neuronId, const uint64_t morphologyId, const uint64_t sectionId, const Synapse& synapse, const size_t baseMaterialId, diff --git a/bioexplorer/backend/science/morphologies/NeuronsLoader.cpp b/bioexplorer/backend/science/morphologies/NeuronsLoader.cpp index b9a610d78..24e7e334c 100644 --- a/bioexplorer/backend/science/morphologies/NeuronsLoader.cpp +++ b/bioexplorer/backend/science/morphologies/NeuronsLoader.cpp @@ -85,13 +85,25 @@ ModelDescriptorPtr NeuronsLoader::importFromStorage(const std::string& storage, props.getProperty(LOADER_PROPERTY_POPULATION_COLOR_SCHEME.name)); details.morphologyColorScheme = stringToEnum( props.getProperty(LOADER_PROPERTY_MORPHOLOGY_COLOR_SCHEME.name)); - details.realismLevel = 0; + details.realismLevel = static_cast(morphology::MorphologyRealismLevel::none); details.realismLevel += (props.getProperty(LOADER_PROPERTY_MORPHOLOGY_REALISM_LEVEL_SOMA.name) ? static_cast(morphology::MorphologyRealismLevel::soma) : 0); details.realismLevel += (props.getProperty(LOADER_PROPERTY_MORPHOLOGY_REALISM_LEVEL_DENDRITE.name) ? static_cast(morphology::MorphologyRealismLevel::dendrite) : 0); + details.realismLevel += (props.getProperty(LOADER_PROPERTY_NEURONS_REALISM_LEVEL_SPINE.name) + ? static_cast(morphology::MorphologyRealismLevel::spine) + : 0); + details.realismLevel += (props.getProperty(LOADER_PROPERTY_MORPHOLOGY_REALISM_LEVEL_AXON.name) + ? static_cast(morphology::MorphologyRealismLevel::axon) + : 0); + details.realismLevel += (props.getProperty(LOADER_PROPERTY_NEURONS_REALISM_LEVEL_EXTERNALS.name) + ? static_cast(morphology::MorphologyRealismLevel::externals) + : 0); + details.realismLevel += (props.getProperty(LOADER_PROPERTY_MORPHOLOGY_REALISM_LEVEL_INTERNALS.name) + ? static_cast(morphology::MorphologyRealismLevel::internals) + : 0); details.loadSomas = props.getProperty(LOADER_PROPERTY_MORPHOLOGY_LOAD_SOMA.name); details.loadAxon = props.getProperty(LOADER_PROPERTY_NEURONS_LOAD_AXON.name); details.loadBasalDendrites = props.getProperty(LOADER_PROPERTY_NEURONS_LOAD_BASAL_DENDRITES.name); diff --git a/bioexplorer/pythonsdk/notebooks/neurons/BioExplorer_neurons.ipynb b/bioexplorer/pythonsdk/notebooks/neurons/BioExplorer_neurons.ipynb index 77d1d94de..9286dd803 100644 --- a/bioexplorer/pythonsdk/notebooks/neurons/BioExplorer_neurons.ipynb +++ b/bioexplorer/pythonsdk/notebooks/neurons/BioExplorer_neurons.ipynb @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 413, "id": "697e36dd", "metadata": { "scrolled": false @@ -45,20 +45,31 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 414, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 414, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "core.set_geometry_parameters(\n", - " sdf_nb_march_iterations=32, sdf_epsilon=0.001,\n", - " sdf_blend_factor=0.8, sdf_blend_lerp_factor=0.1,\n", - " sdf_ray_marching_omega=1.0, sdf_distance=100.0\n", - ")" + " sdf_nb_march_iterations=64, sdf_epsilon=0.001,\n", + " sdf_blend_factor=0.2, sdf_blend_lerp_factor=0.1,\n", + " sdf_ray_marching_omega=0.4, sdf_distance=50.0\n", + ")\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 415, "id": "3007d98c", "metadata": {}, "outputs": [], @@ -73,15 +84,16 @@ " morphology_representation=be.MORPHOLOGY_REPRESENTATION_SEGMENT,\n", " morphology_color_scheme=be.MORPHOLOGY_COLOR_SCHEME_SECTION_TYPE,\n", " realism_level=be.MORPHOLOGY_REALISM_LEVEL_ALL,\n", + " # realism_level=be.MORPHOLOGY_REALISM_LEVEL_SOMA,\n", " synapses_type=be.NEURON_SYNAPSES_DEBUG,\n", - " load_axon=True,\n", + " load_axon=False,\n", " generate_varicosities=False,\n", " load_basal_dendrites=True, load_apical_dendrites=True,\n", - " generate_internals=False, generate_externals=True,\n", - " # sql_node_filter='guid=69',\n", - " sql_node_filter='guid=89',\n", + " # generate_internals=False, generate_externals=True,\n", + " sql_node_filter='guid=49',\n", + " # sql_node_filter='guid=89',\n", " displacement_params=NeuronDisplacementParams(\n", - " soma=Vector2(0.25, 2.0), section=Vector2(0.1, 1.0),\n", + " soma=Vector2(0.25, 2.0), section=Vector2(0.25, 2.0),\n", " nucleus=Vector2(0.02, 1.0), mitochondrion=Vector2(0.4, 1.0),\n", " myelin_steath=Vector2(0.15, 4.0), spine=Vector2(0.01, 50.0))\n", ")" @@ -89,7 +101,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 416, "id": "2c667b48", "metadata": {}, "outputs": [], @@ -104,7 +116,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 417, "id": "3e5b7003", "metadata": {}, "outputs": [], @@ -184,10 +196,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 418, "id": "26dc02eb", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 418, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "model_ids = be.get_model_ids()['ids']\n", "offset = model_ids[0]\n", @@ -197,7 +220,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 419, "id": "9fd33ec8", "metadata": {}, "outputs": [], diff --git a/bioexplorer/pythonsdk/notebooks/utils/BioExplorer_sdf_testing.ipynb b/bioexplorer/pythonsdk/notebooks/utils/BioExplorer_sdf_testing.ipynb index 8cb38e020..15d538a5e 100644 --- a/bioexplorer/pythonsdk/notebooks/utils/BioExplorer_sdf_testing.ipynb +++ b/bioexplorer/pythonsdk/notebooks/utils/BioExplorer_sdf_testing.ipynb @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": 643, + "execution_count": 95, "metadata": {}, "outputs": [], "source": [ @@ -20,20 +20,49 @@ "be = BioExplorer('localhost:5000')\n", "core = be.core_api()\n", "be.reset_scene()\n", - "core.set_geometry_parameters(\n", - " sdf_nb_march_iterations=32, sdf_epsilon=0.001,\n", - " sdf_blend_factor=0.4, sdf_blend_lerp_factor=0.2,\n", - " sdf_ray_marching_omega=1.0, sdf_distance=5.0\n", - ")\n", "be.add_sdf_demo()\n", "status = be.reset_camera()" ] }, { "cell_type": "code", - "execution_count": 644, + "execution_count": 96, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 96, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "core.set_geometry_parameters(\n", + " sdf_nb_march_iterations=64, sdf_epsilon=0.001,\n", + " sdf_blend_factor=0.2, sdf_blend_lerp_factor=0.1,\n", + " sdf_ray_marching_omega=0.4, sdf_distance=5.0\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 97, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Future exception was never retrieved\n", + "future: \n", + "websockets.exceptions.ConnectionClosedError: no close frame received or sent\n" + ] + } + ], "source": [ "status = core.set_renderer(\n", " current='advanced',\n", diff --git a/platform/core/engineapi/Model.h b/platform/core/engineapi/Model.h index 4a3e8548f..87da647bc 100644 --- a/platform/core/engineapi/Model.h +++ b/platform/core/engineapi/Model.h @@ -53,13 +53,13 @@ namespace core */ struct SDFGeometryData { - std::vector geometries; // A vector of SDFGeometry objects - std::map> geometryIndices; // A map where the keys are size_t values that represent - // the indices of geometries vector. The map values are - // vectors of uint64_t values that represent the indices of - // vertices that belong to the geometry. - std::vector> neighbours; // A vector of vectors of uint64_t values that represent the indices - // of neighbouring vertices for each vertex. + std::vector geometries; // A vector of SDFGeometry objects + std::map geometryIndices; // A map where the keys are size_t values that represent + // the indices of geometries vector. The map values are + // vectors of uint64_t values that represent the indices of + // vertices that belong to the geometry. + std::vector neighbours; // A vector of vectors of uint64_t values that represent the indices + // of neighbouring vertices for each vertex. std::vector neighboursFlat; // A vector of uint64_t values that represents the indices of neighbouring // vertices flattened into a single vector. diff --git a/platform/engines/optix6/OptiXModel.cpp b/platform/engines/optix6/OptiXModel.cpp index 9068c1e2d..426084015 100644 --- a/platform/engines/optix6/OptiXModel.cpp +++ b/platform/engines/optix6/OptiXModel.cpp @@ -34,6 +34,8 @@ #include +#include + using namespace optix; namespace core @@ -87,6 +89,7 @@ OptiXModel::~OptiXModel() { RT_DESTROY(sdfGeometriesBuffers.second.geometries_buffer); RT_DESTROY(sdfGeometriesBuffers.second.neighbours_buffer); + RT_DESTROY(sdfGeometriesBuffers.second.indices_buffer); } RT_DESTROY_MAP(_optixVolumes) @@ -292,72 +295,88 @@ uint64_t OptiXModel::_commitSDFGeometries() auto context = OptiXContext::get().getOptixContext(); auto& sdfGeometries = _geometries->_sdf; - context[CONTEXT_SDF_GEOMETRY_SIZE]->setUint(sizeof(SDFGeometry)); + const uint32_t sdfGeometrySize = sizeof(SDFGeometry); + context[CONTEXT_SDF_GEOMETRY_SIZE]->setUint(sdfGeometrySize); const uint64_t nbGeometries = sdfGeometries.geometries.size(); const uint64_t geometriesBufferSize = sizeof(SDFGeometry) * nbGeometries; uint64_ts geometryMaterials; - for (const auto& geometryIndices : sdfGeometries.geometryIndices) - geometryMaterials.push_back(geometryIndices.first); + for (const auto& geometryIndex : sdfGeometries.geometryIndices) + geometryMaterials.push_back(geometryIndex.first); + const uint64_t nbMaterials = geometryMaterials.size(); + uint64_t i = 0; -#pragma omp parallel for - for (i = 0; i < sdfGeometries.geometryIndices.size(); ++i) +#pragma omp parallel for private(memoryFootPrint) + for (i = 0; i < nbMaterials; ++i) { const auto materialId = geometryMaterials[i]; - const auto& indices = sdfGeometries.geometryIndices[materialId]; - _optixSdfGeometries[materialId] = OptiXContext::get().createGeometry(OptixGeometryType::sdfGeometry); + + // Indices of geometries for current material id + const auto& geometryIndices = sdfGeometries.geometryIndices[materialId]; // Create a local copy of SDF geometries attached to the material id std::map localGeometries; - for (const auto primIdx : indices) + + for (const auto geometryIndex : geometryIndices) { - localGeometries[primIdx] = &sdfGeometries.geometries[primIdx]; - const auto& neighbours = sdfGeometries.neighbours[primIdx]; - if (!neighbours.empty()) - { - localGeometries[primIdx]->numNeighbours = neighbours.size(); - // Not used by OptiX, every geometry holds it own vector of neighbours - localGeometries[primIdx]->neighboursIndex = 0; - for (const auto neighbour : neighbours) - { - localGeometries[neighbour] = &sdfGeometries.geometries[neighbour]; - localGeometries[neighbour]->numNeighbours = sdfGeometries.neighbours[neighbour].size(); - } - } + localGeometries[geometryIndex] = &sdfGeometries.geometries[geometryIndex]; + const auto& neighbours = sdfGeometries.neighbours[geometryIndex]; + for (const auto neighbour : neighbours) + localGeometries[neighbour] = &sdfGeometries.geometries[neighbour]; } - // Prepare an index between the local SDF geometries and their position in the localGeometries map + uint64_t index = 0; uint64_tm localGeometriesMapping; - uint64_t i = 0; - for (const auto localGeometry : localGeometries) + for (const auto& localGeometry : localGeometries) { - localGeometriesMapping[localGeometry.first] = i; - ++i; + localGeometriesMapping[localGeometry.first] = index; + ++index; } + uint64_ts localIndicesFlat; + localIndicesFlat.reserve(geometryIndices.size()); + for (const auto geometryIndex : geometryIndices) + localIndicesFlat.push_back(localGeometriesMapping[geometryIndex]); + + if (localIndicesFlat.empty()) + continue; + // Create a local flat representation of the SDF geometries neighbours uint64_ts localNeighboursFlat; - for (const auto primIdx : indices) + for (const auto geometryIndex : geometryIndices) { - const auto& neighbours = sdfGeometries.neighbours[primIdx]; - for (const auto neighbour : neighbours) + const auto& neighbours = sdfGeometries.neighbours[geometryIndex]; + localGeometries[geometryIndex]->neighboursIndex = localNeighboursFlat.size(); + localGeometries[geometryIndex]->numNeighbours = static_cast(neighbours.size()); + for (uint64_t i = 0; i < neighbours.size(); ++i) { - const auto it = localGeometriesMapping.find(neighbour); + const auto it = localGeometriesMapping.find(neighbours[i]); if (it == localGeometriesMapping.end()) CORE_THROW("Invalid neighbour index"); localNeighboursFlat.push_back((*it).second); } } - // Make sure we don't create an empty buffer in the case of no neighbours - if (localNeighboursFlat.empty()) - localNeighboursFlat.resize(1, 0); + // Convert map to vector in order to send the data to the device + std::vector localGeometriesFlat; + localGeometriesFlat.reserve(localGeometries.size()); + for (const auto geometry : localGeometries) + localGeometriesFlat.push_back(*geometry.second); + + // Prepare OptiX geometry and corresponding buffers + _optixSdfGeometries[materialId] = OptiXContext::get().createGeometry(OptixGeometryType::sdfGeometry); + _optixSdfGeometries[materialId]->setPrimitiveCount(localIndicesFlat.size()); - // Prepare buffers - _optixSdfGeometries[materialId]->setPrimitiveCount(indices.size()); + uint64_t bufferSize = localIndicesFlat.size(); +#pragma omp critical + setBuffer(RT_BUFFER_INPUT, RT_FORMAT_UNSIGNED_LONG_LONG, _sdfGeometriesBuffers[materialId].indices_buffer, + _optixSdfGeometries[materialId][OPTIX_GEOMETRY_PROPERTY_SDF_GEOMETRIES_INDICES], localIndicesFlat, + bufferSize); +#pragma omp critical + memoryFootPrint += bufferSize; - uint64_t bufferSize = localNeighboursFlat.size(); + bufferSize = localNeighboursFlat.size(); #pragma omp critical setBuffer(RT_BUFFER_INPUT, RT_FORMAT_UNSIGNED_LONG_LONG, _sdfGeometriesBuffers[materialId].neighbours_buffer, _optixSdfGeometries[materialId][OPTIX_GEOMETRY_PROPERTY_SDF_GEOMETRIES_NEIGHBOURS], @@ -365,14 +384,10 @@ uint64_t OptiXModel::_commitSDFGeometries() #pragma omp critical memoryFootPrint += bufferSize; - std::vector localGeometriesAsVector; - localGeometriesAsVector.reserve(localGeometries.size()); - for (const auto geometry : localGeometries) - localGeometriesAsVector.push_back(*geometry.second); - bufferSize = sizeof(SDFGeometry) * localGeometriesAsVector.size(); + bufferSize = sizeof(SDFGeometry) * localGeometriesFlat.size(); #pragma omp critical setBuffer(RT_BUFFER_INPUT, RT_FORMAT_BYTE, _sdfGeometriesBuffers[materialId].geometries_buffer, - _optixSdfGeometries[materialId][OPTIX_GEOMETRY_PROPERTY_SDF_GEOMETRIES], localGeometriesAsVector, + _optixSdfGeometries[materialId][OPTIX_GEOMETRY_PROPERTY_SDF_GEOMETRIES], localGeometriesFlat, bufferSize); #pragma omp critical memoryFootPrint += bufferSize; diff --git a/platform/engines/optix6/OptiXModel.h b/platform/engines/optix6/OptiXModel.h index 3efaed11b..516301ad0 100644 --- a/platform/engines/optix6/OptiXModel.h +++ b/platform/engines/optix6/OptiXModel.h @@ -103,6 +103,7 @@ class OptiXModel : public Model // SDF geometries struct OptiXSDFGeometryBuffers { + ::optix::Buffer indices_buffer{nullptr}; ::optix::Buffer geometries_buffer{nullptr}; ::optix::Buffer neighbours_buffer{nullptr}; }; diff --git a/platform/engines/optix6/cuda/geometry/SDFGeometries.cu b/platform/engines/optix6/cuda/geometry/SDFGeometries.cu index fa6760425..88b3f8f04 100644 --- a/platform/engines/optix6/cuda/geometry/SDFGeometries.cu +++ b/platform/engines/optix6/cuda/geometry/SDFGeometries.cu @@ -24,13 +24,16 @@ // Ray-cone intersection: based on Ching-Kuang Shene (Graphics Gems 5, p. 227-230) #include +#include // Global variables rtDeclareVariable(uint, sdf_geometry_size, , ); rtBuffer sdf_geometries_buffer; +rtBuffer sdf_geometries_indices_buffer; rtBuffer sdf_geometries_neighbours_buffer; #define SDF_NO_INTERSECTION -1.f +#define SDF_NORMAL_EPSILON_FACTOR 10.f enum SDFType : uint8_t { @@ -74,21 +77,36 @@ static __device__ inline float sign(const float x) return (x >= 0.f ? 1.f : -1.f); } -static __device__ inline float lerp(const float factor, const float a, const float b) -{ - return (1.f - factor) * a + factor * b; -} - // polynomial smooth min (k = 0.1); static __device__ inline float sminPoly(const float a, const float b, const float k) { const float h = ::optix::clamp(0.5f + 0.5f * (b - a) / k, 0.f, 1.f); - return lerp(h, b, a) - k * h * (1.f - h); + return ::optix::lerp(b, a, h) - k * h * (1.f - h); +} + +static __device__ inline float opSmoothUnion(float d1, float d2, float k) +{ + const float h = ::optix::clamp(0.5f + 0.5f * (d2 - d1) / k, 0.f, 1.f); + return mix(d2, d1, h) - k * h * (1.f - h); } -static __device__ inline float opDisplacement(const float3& p, const float a, const float b) +static __device__ inline float opSmoothSubtraction(const float d1, const float d2, const float k) { - return a * sin(b * p.x) * sin(b * p.y * 0.65f) * sin(b * p.z * 0.81f); + const float h = ::optix::clamp(0.5f - 0.5f * (d2 + d1) / k, 0.f, 1.f); + return mix(d2, -d1, h) + k * h * (1.f - h); +} + +static __device__ inline float opSmoothIntersection(const float d1, const float d2, const float k) +{ + const float h = ::optix::clamp(0.5f - 0.5f * (d2 - d1) / k, 0.f, 1.f); + return mix(d2, d1, h) + k * h * (1.f - h); +} + +static __device__ inline float opDisplacement(const float3& p, const float amplitude, const float frequency) +{ + return amplitude * + (0.7f * sin(frequency * p.x * 0.72f) * sin(frequency * p.y * 0.65f) * sin(frequency * p.z * 0.81f) + + 0.3f * cos(p.x * 2.12f) * cos(p.y * 2.23f) * cos(p.z * 2.41f)); } static __device__ inline float sdSphere(const float3& p, const float3& c, float r) @@ -150,7 +168,7 @@ static __device__ inline float sdCone(const float3& p, const float3 a, const flo float cbx = x - ra - f * rba; float cby = paba - f; - float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0; + float s = (cbx < 0.f && cay < 0.f) ? -1.f : 1.f; return s * sqrt(min(cax * cax + cay * cay * baba, cbx * cbx + cby * cby * baba)); } @@ -188,25 +206,33 @@ static __device__ inline float sdVesica(const float3& p, const float3 a, const f return ::optix::length(q - make_float3(h.x, h.y, 0.f)) - h.z; } -static __device__ inline bool intersectBox(const ::optix::Aabb& box, float& t0, float& t1) +static __device__ inline float reduce_min(const float4& a) +{ + return min(min(a.x, a.y), min(a.z, a.w)); +} + +static __device__ inline float reduce_max(const float4& a) { - const float3 a = (box.m_min - ray.origin) / ray.direction; - const float3 b = (box.m_max - ray.origin) / ray.direction; - const float3 near = fminf(a, b); - const float3 far = fmaxf(a, b); - t0 = fmaxf(near); - t1 = fminf(far); + return max(max(a.x, a.y), max(a.z, a.w)); +} +static __device__ inline bool intersectBox(const ::optix::Aabb& box, float& t0, float& t1) +{ + const float3 mins = (box.m_min - ray.origin) / ray.direction; + const float3 maxs = (box.m_max - ray.origin) / ray.direction; + t0 = reduce_max(make_float4(fminf(mins, maxs), ray.tmin)); + t1 = reduce_min(make_float4(fmaxf(mins, maxs), ray.tmax)); return (t0 <= t1); } static __device__ inline SDFGeometry* getPrimitive(const int primIdx) { - const uint64_t bufferIndex = primIdx * sdf_geometry_size; + const uint64_t geometryIndex = sdf_geometries_indices_buffer[primIdx]; + const uint64_t bufferIndex = geometryIndex * sdf_geometry_size; return (SDFGeometry*)&sdf_geometries_buffer[bufferIndex]; } -static __device__ inline uint64_t getNeighbourIndex(const uint8_t index) +static __device__ inline uint64_t getNeighbourIndex(const uint64_t index) { return sdf_geometries_neighbours_buffer[index]; } @@ -278,13 +304,13 @@ static __device__ inline float sdfDistance(const float3& position, const SDFGeom const float r0 = max(primitive->r0, primitive->r1); for (uint8_t i = 0; i < primitive->numNeighbours; ++i) { - const uint64_t neighbourIndex = getNeighbourIndex(i); + const uint64_t neighbourIndex = getNeighbourIndex(primitive->neighboursIndex + i); const SDFGeometry* neighbourGeometry = getPrimitive(neighbourIndex); const float neighbourDistance = calcDistance(neighbourGeometry, position, processDisplacement); if (neighbourDistance < 0.f) continue; const float r1 = max(neighbourGeometry->r0, neighbourGeometry->r1); - const float blendFactor = lerp(geometrySdfBlendLerpFactor, min(r0, r1), max(r0, r1)); + const float blendFactor = ::optix::lerp(min(r0, r1), max(r0, r1), geometrySdfBlendLerpFactor); distance = sminPoly(neighbourDistance, distance, blendFactor * geometrySdfBlendFactor); } } @@ -300,10 +326,11 @@ static __device__ inline float3 computeNormal(const float3& position, const SDFG const float3 k1 = make_float3(-t, -t, t); const float3 k2 = make_float3(-t, t, -t); const float3 k3 = make_float3(t, t, t); - return ::optix::normalize(k0 * sdfDistance(position + geometrySdfEpsilon * k0, primitive, processDisplacement) + - k1 * sdfDistance(position + geometrySdfEpsilon * k1, primitive, processDisplacement) + - k2 * sdfDistance(position + geometrySdfEpsilon * k2, primitive, processDisplacement) + - k3 * sdfDistance(position + geometrySdfEpsilon * k3, primitive, processDisplacement)); + const float e = geometrySdfEpsilon * SDF_NORMAL_EPSILON_FACTOR; + return ::optix::normalize(k0 * sdfDistance(position + e * k0, primitive, processDisplacement) + + k1 * sdfDistance(position + e * k1, primitive, processDisplacement) + + k2 * sdfDistance(position + e * k2, primitive, processDisplacement) + + k3 * sdfDistance(position + e * k3, primitive, processDisplacement)); } static __device__ inline float rayMarching(const SDFGeometry* primitive, bool& processDisplacement) @@ -318,9 +345,9 @@ static __device__ inline float rayMarching(const SDFGeometry* primitive, bool& p const float pixel_radius = geometrySdfEpsilon; float omega = geometrySdfOmega; - float t = t0; float candidateError = 1e6f; - float tCandidate = t0; + float t = t0; + float tCandidate = t; float previousRadius = 0.f; float stepLength = 0.f; uint64_t stepCount = 0; @@ -333,7 +360,7 @@ static __device__ inline float rayMarching(const SDFGeometry* primitive, bool& p { const float3 p = ray.origin + ray.direction * t; const float3 tp = rtTransformPoint(RT_OBJECT_TO_WORLD, p); - processDisplacement = (/*ray.flags == RAY_FLAG_PRIMARY && */ ::optix::length(tp - eye) < geometrySdfDistance); + processDisplacement = (/*prd.depth == 0 && */ ::optix::length(tp - eye) < geometrySdfDistance); float signed_radius = sdfSign * sdfDistance(p, primitive, processDisplacement); float radius = abs(signed_radius); diff --git a/platform/engines/ospray/ispc/geometry/SDFGeometries.ispc b/platform/engines/ospray/ispc/geometry/SDFGeometries.ispc index 3b5c8a5e4..be28a2074 100644 --- a/platform/engines/ospray/ispc/geometry/SDFGeometries.ispc +++ b/platform/engines/ospray/ispc/geometry/SDFGeometries.ispc @@ -35,6 +35,8 @@ #include "RayMarching.isph" +#define SDF_NO_INTERSECTION -1.f + enum SDFType { sdf_sphere = 0, @@ -47,8 +49,6 @@ enum SDFType sdf_vesica = 7 }; -///////////////////////////////////////////////////////////////////////////// - inline vec3f mod(const vec3f& v, const int m) { return make_vec3f(v.x - m * floor(v.x / m), v.y - m * floor(v.y / m), v.z - m * floor(v.z / m)); @@ -71,9 +71,11 @@ inline float smootherStep(const float x) return x * x * x * (x * (x * 6 - 15) + 10); } -inline float opDisplacement(const vec3f& p, const float a, const float b) +inline float opDisplacement(const vec3f& p, const float amplitude, const float frequency) { - return a * sin(b * p.x) * sin(b * p.y * 0.6) * sin(b * p.z * 0.3); + return amplitude * + (0.7f * sin(frequency * p.x * 0.72f) * sin(frequency * p.y * 0.65f) * sin(frequency * p.z * 0.81f) + + 0.3f * cos(p.x * 2.12f) * cos(p.y * 2.23f) * cos(p.z * 2.41f)); } inline float sdSphere(const vec3f& p, const vec3f c, float r) @@ -126,33 +128,6 @@ inline float mix(const float x, const float y, const float a) return x * (1.f - a) + y * a; } -#if 0 -inline float sdConePill(const vec3f& p, const vec3f p0, const vec3f p1, const float r0, const float r1, - const bool useSigmoid) -{ - const vec3f v = p1 - p0; - const vec3f w = p - p0; - - // distance to p0 along cone axis - const float c1 = dot(w, v); - if (c1 <= 0.f) // p0 Sphere - return length(p - p0) - r0; - - // cone length - const float c2 = dot(v, v); - if (c2 <= c1) // p1 Sphere - return length(p - p1) - r1; - - const float b = c1 / c2; - const vec3f Pb = p0 + b * v; - - const float thicknessAt = lerp(b, r0, r1); - const float thickness = useSigmoid ? mix(r0, r1, smootherStep(b)) : thicknessAt; - - return length(p - Pb) - thickness; -} -#else - inline float sign(const float x) { return (x >= 0.f ? 1.f : -1.f); @@ -185,7 +160,6 @@ inline float sdConePill(const vec3f& p, const vec3f a, const vec3f b, const floa return sqrt(x2 + y2) * il2 - r1; return (sqrt(x2 * a2 * il2) + y * rr) * il2 - r1; } -#endif inline float sdCone(const vec3f& p, const vec3f a, const vec3f b, float ra, float rb) { @@ -205,7 +179,7 @@ inline float sdCone(const vec3f& p, const vec3f a, const vec3f b, float ra, floa float cbx = x - ra - f * rba; float cby = paba - f; - float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0; + float s = (cbx < 0.f && cay < 0.f) ? -1.f : 1.f; return s * sqrt(min(cax * cax + cay * cay * baba, cbx * cbx + cby * cby * baba)); } @@ -242,8 +216,6 @@ inline float sdVesica(const vec3f& p, const vec3f a, const vec3f b, float w) return length(q - make_vec3f(h.x, h.y, 0.f)) - h.z; } -////////////////////////////////////////////////////////////////////// - // NOTE: This layout must match exactly the 'SDFGeometry' struct in // 'SDFGeometry.h' struct SDFGeometry @@ -276,18 +248,12 @@ struct SDFGeometries uniform float distance; }; -////////////////////////////////////////////////////////////////////// - typedef uniform SDFGeometries* uniform Geo_ptr; typedef uniform SDFGeometry* uniform Prim_ptr; -////////////////////////////////////////////////////////////////////// - DEFINE_SAFE_INCREMENT(SDFGeometry); DEFINE_SAFE_INCREMENT(uint64); -////////////////////////////////////////////////////////////////////// - uniform uint64 primToIdx(const uniform SDFGeometries& geometry, const uniform uint64 primID) { return *safeIncrement(geometry.useSafeIncrement, geometry.geometryRefs, primID); @@ -314,8 +280,6 @@ const uniform SDFGeometry* varying getPrimitiveVarying(const uniform SDFGeometri return safeIncrement(geometry.useSafeIncrement, geometry.geometries, idx); } -////////////////////////////////////////////////////////////////////// - inline float calcDistance(const uniform SDFGeometry& primitive, const vec3f& p, const bool processDisplacement) { const float displacement = (processDisplacement && primitive.userParams.x > 0.f) @@ -340,11 +304,9 @@ inline float calcDistance(const uniform SDFGeometry& primitive, const vec3f& p, case sdf_vesica: return displacement + sdVesica(p, primitive.p0, primitive.p1, primitive.r0); } - return -1.0; + return SDF_NO_INTERSECTION; } -////////////////////////////////////////////////////////////////////// - uniform box3fa sdfBounds(uDataPtr_t geo_, uDataPtr_t prim_) { const Prim_ptr prim = (const Prim_ptr)prim_; @@ -365,36 +327,31 @@ uniform box3fa sdfBounds(uDataPtr_t geo_, uDataPtr_t prim_) } } -////////////////////////////////////////////////////////////////////// - float sdfDistance(const vec3f& p, uDataPtr_t geo_, uDataPtr_t prim_, const SDFParams& params, const bool processDisplacement) { const Geo_ptr geo = (const Geo_ptr)geo_; const Prim_ptr prim = (const Prim_ptr)prim_; - float d = calcDistance(*prim, p, processDisplacement); - + float distance = calcDistance(*prim, p, processDisplacement); if (processDisplacement && prim->numNeighbours > 0) { const float r0 = max(prim->r0, prim->r1); - for (uniform int i = 0; i < prim->numNeighbours; i++) + for (uniform uint8 i = 0; i < prim->numNeighbours; ++i) { const uniform uint64 index = getNeighbourIdx(*geo, prim->neighboursIndex, i); const uniform SDFGeometry& neighbour = *getPrimitive(*geo, index); - const float dOther = calcDistance(neighbour, p, processDisplacement); - if (dOther < 0.f) + const float neighbourDistance = calcDistance(neighbour, p, processDisplacement); + if (neighbourDistance < 0.f) continue; const float r1 = max(neighbour.r0, neighbour.r1); const float blendFactor = lerp(params.blendLerpFactor, min(r0, r1), max(r0, r1)); - d = sminPoly(dOther, d, blendFactor * params.blendFactor); + distance = sminPoly(neighbourDistance, distance, blendFactor * params.blendFactor); } } - return d; + return distance; } -////////////////////////////////////////////////////////////////////// - void SDFGeometries_bounds(const RTCBoundsFunctionArguments* uniform args) { const Geo_ptr geo = (Geo_ptr)args->geometryUserPtr; @@ -404,8 +361,6 @@ void SDFGeometries_bounds(const RTCBoundsFunctionArguments* uniform args) *((box3fa * uniform) args->bounds_o) = sdfBounds((uDataPtr_t)geo, (uDataPtr_t)prim); } -////////////////////////////////////////////////////////////////////// - unmasked void SDFGeometries_intersect(const RTCIntersectFunctionNArguments* uniform args) { const Geo_ptr geo = (Geo_ptr)args->geometryUserPtr; @@ -443,8 +398,6 @@ unmasked void SDFGeometries_intersect(const RTCIntersectFunctionNArguments* unif } } -////////////////////////////////////////////////////////////////////// - static void SDFGeometries_postIntersect(uniform Geometry* uniform geometry, uniform Model* uniform model, varying DifferentialGeometry& dg, const varying Ray& ray, uniform int64 flags) { @@ -469,8 +422,6 @@ static void SDFGeometries_postIntersect(uniform Geometry* uniform geometry, unif dg.Ns = Ns; } -////////////////////////////////////////////////////////////////////// - export void* uniform SDFGeometries_create(void* uniform cppEquivalent) { Geo_ptr geom = uniform new uniform SDFGeometries; @@ -478,8 +429,6 @@ export void* uniform SDFGeometries_create(void* uniform cppEquivalent) return geom; } -////////////////////////////////////////////////////////////////////// - export void SDFGeometriesGeometry_set(void* uniform _self, void* uniform _model, void* uniform data, int uniform numPrimitives, void* uniform neighbours, int uniform numNeighbours, void* uniform geometries, float uniform epsilon, int uniform nbMarchIterations,