Build of somas as spheres

Author: Cyrille Pierre Henri Favreau

bioexplorer/backend/science/morphologies/Astrocytes.cpp

@@ -203,14 +203,15 @@ void Astrocytes::_buildModel(const LoaderProgress& callback, const doubles& radi
  const bool useSdf = andCheck(static_cast<uint32_t>(_details.realismLevel),
  static_cast<uint32_t>(MorphologyRealismLevel::soma));
 
- if (!_spheresRepresentation.enabled)
- {
+ if (_spheresRepresentation.enabled)
+ _addSomaAsSpheres(somaId, somaMaterialId, sections, somaPosition, Quaterniond(), somaRadius,
+ NO_USER_DATA, _details.radiusMultiplier, container);
+ else
  somaGeometryIndex = container.addSphere(
  somaPosition, somaRadius, somaMaterialId, useSdf, NO_USER_DATA, {},
  Vector3f(somaRadius * _getDisplacementValue(DisplacementElement::morphology_soma_strength),
  somaRadius * _getDisplacementValue(DisplacementElement::morphology_soma_frequency),
  0.f));
- }
 
  if (_details.generateInternals)
  {
@@ -259,7 +260,7 @@ void Astrocytes::_buildModel(const LoaderProgress& callback, const doubles& radi
  const bool useSdf = andCheck(static_cast<uint32_t>(_details.realismLevel),
  static_cast<uint32_t>(MorphologyRealismLevel::dendrite));
  uint64_t geometryIndex = 0;
- if (section.second.parentId == SOMA_AS_PARENT)
+ if (section.second.parentId == SOMA_AS_PARENT && !_spheresRepresentation.enabled)
  {
  // Section connected to the soma
  const auto& point = points[0];

bioexplorer/backend/science/morphologies/Morphologies.cpp

@@ -41,6 +41,76 @@ size_t Morphologies::_getNbMitochondrionSegments() const
  return 2 + rand() % 5;
 }
 
+double Morphologies::_addSomaAsSpheres(const uint64_t neuronId, const size_t somaMaterialId, const SectionMap& sections,
+ const Vector3d& somaPosition, const Quaterniond& somaRotation,
+ const double somaRadius, const uint64_t somaUserData,
+ const double radiusMultiplier, ThreadSafeContainer& container)
+{
+ Vector3ds sectionRoots;
+
+ double minRadius = std::numeric_limits<double>::max();
+ double maxRadius = std::numeric_limits<double>::min();
+ double sectionRootRadius = std::numeric_limits<double>::max();
+ Vector3d baryCenter;
+ for (const auto& section : sections)
+ if (section.second.parentId == SOMA_AS_PARENT)
+ {
+ const auto& points = section.second.points;
+ const uint64_t index = std::min(size_t(0), points.size());
+ const Vector3d p = points[index];
+ sectionRoots.push_back(p);
+ const double l = length(p);
+ minRadius = std::min(l, minRadius);
+ maxRadius = std::max(l, maxRadius);
+ sectionRootRadius = std::min(sectionRootRadius, static_cast<double>(points[index].w * 0.5));
+ baryCenter += p;
+ }
+ baryCenter /= sectionRoots.size();
+ // sectionRootRadius = sectionRootRadius / sectionRoots.size() * radiusMultiplier;
+
+ const double radius = (minRadius + maxRadius) / 2.0;
+ const double somaSurface = 4.0 * glm::pi<double>() * pow(radius, 2.0);
+ const double sphereSurfaceOnSoma = glm::pi<double>() * pow(sectionRootRadius, 2.0);
+ const uint64_t nbSpheres = somaSurface / sphereSurfaceOnSoma;
+
+ Vector3ds spheres(nbSpheres);
+
+ const double goldenRatio = (1.0 + std::sqrt(5.0)) / 2.0;
+ const double angleIncrement = 2.0 * M_PI * goldenRatio;
+ for (uint64_t i = 0; i < nbSpheres; ++i)
+ {
+ const double t = static_cast<double>(i) / static_cast<double>(nbSpheres - 1);
+ const double phi = std::acos(1.0 - 2.0 * t);
+ const double theta = angleIncrement * static_cast<double>(i);
+ const double r = minRadius;
+ spheres[i] = Vector3d(r * sin(phi) * cos(theta), r * sin(phi) * sin(theta), r * cos(phi));
+ }
+
+ // Smooth soma according to section root
+ for (auto& sphere : spheres)
+ {
+#if 1
+ const double smoothingFactor = 1.0;
+ double r = minRadius * smoothingFactor;
+ for (const auto& sr : sectionRoots)
+ {
+ const auto dir = sr - (baryCenter + sphere);
+ const double angle = dot(normalize(sr), normalize(baryCenter + sphere));
+ if (angle >= 0.0)
+ r += length(dir) * -angle * smoothingFactor;
+ }
+#endif
+
+ const auto src =
+ _animatedPosition(Vector4d(somaPosition + somaRotation * (baryCenter + normalize(sphere) * r * 0.5),
+ sectionRootRadius),
+ neuronId);
+ container.addSphere(src, sectionRootRadius, somaMaterialId, false, somaUserData);
+ }
+
+ return somaRadius;
+}
+
 void Morphologies::_addSomaInternals(ThreadSafeContainer& container, const size_t baseMaterialId,
  const Vector3d& somaPosition, const double somaRadius,
  const double mitochondriaDensity, const bool useSdf, const double radiusMultiplier)

bioexplorer/backend/science/morphologies/Morphologies.h

@@ -56,6 +56,11 @@ class Morphologies : public common::SDFGeometries
 protected:
  size_t _getNbMitochondrionSegments() const;
 
+ double _addSomaAsSpheres(const uint64_t neuronId, const size_t somaMaterialId, const SectionMap& sections,
+ const core::Vector3d& somaPosition, const core::Quaterniond& somaRotation,
+ const double somaRadius, const uint64_t somaUserData, const double radiusMultiplier,
+ common::ThreadSafeContainer& container);
+
  void _addSomaInternals(common::ThreadSafeContainer& container, const size_t materialId,
  const core::Vector3d& somaPosition, const double somaRadius,
  const double mitochondriaDensity, const bool useSdf, const double radiusMultiplier);

bioexplorer/backend/science/morphologies/Neurons.cpp

@@ -555,69 +555,6 @@ double Neurons::_addSoma(const uint64_t neuronId, const size_t somaMaterialId, c
  return count == 0 ? somaRadius : correctedSomaRadius / count;
 }
 
-double Neurons::_addSomaAsSpheres(const uint64_t neuronId, const size_t somaMaterialId, const SectionMap& sections,
- const Vector3d& somaPosition, const Quaterniond& somaRotation,
- const double somaRadius, const uint64_t somaUserData, ThreadSafeContainer& container)
-{
- // Paramètres
- uint64_t nbPoints = 2000;
- double small_sphere_radius = 0.5;
-
- // Vecteurs pour stocker les coordonnées
- Vector3ds points(nbPoints);
- Vector3ds stretched_points;
-
- double minRadius = std::numeric_limits<double>::max();
- double maxRadius = std::numeric_limits<double>::min();
- Vector3d baryCenter;
- for (const auto& section : sections)
- {
- if (section.second.parentId == SOMA_AS_PARENT)
- {
- const auto& points = section.second.points;
- const uint64_t index = std::min(size_t(1), points.size());
- const Vector3d p = points[index];
- stretched_points.push_back(p);
- const double l = length(p);
- minRadius = std::min(l, minRadius);
- maxRadius = std::max(l, maxRadius);
- baryCenter += p;
- }
- }
- baryCenter /= stretched_points.size();
-
- const auto radius = (minRadius + maxRadius) / 2.0;
-
- // Génération des points aléatoires sur une sphère
- for (uint64_t i = 0; i < nbPoints; ++i)
- {
- const double theta = glm::linearRand(0.0, 2.0 * glm::pi<double>());
- const double phi = acos(glm::linearRand(-1.0, 1.0));
- points[i] = baryCenter + Vector3d(minRadius * sin(theta) * cos(phi), minRadius * sin(theta) * sin(phi),
- minRadius * cos(theta));
- }
-
- for (auto& point : points)
- {
- auto p = Vector3d();
- for (const auto& sp : stretched_points)
- {
- const auto dir = sp - point;
- const double angle = dot(normalize(sp), normalize(point));
- if (angle > 0.7)
- p += dir * angle;
- }
- p /= stretched_points.size();
- point += 3.0 * p;
-
- const auto src =
- _animatedPosition(Vector4d(somaPosition + somaRotation * point, small_sphere_radius), neuronId);
- container.addSphere(src, small_sphere_radius, somaMaterialId, false, somaUserData);
- }
-
- return somaRadius;
-}
-
 void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t neuronId, const NeuronSoma& soma,
  const uint64_t neuronIndex, const float* voltages)
 {
@@ -733,7 +670,7 @@ void Neurons::_buildMorphology(ThreadSafeContainer& container, const uint64_t ne
  double correctedSomaRadius = 0.f;
  if (_spheresRepresentation.enabled)
  correctedSomaRadius = _addSomaAsSpheres(neuronId, somaMaterialId, sections, somaPosition, somaRotation,
- somaRadius, somaUserData, container);
+ somaRadius, somaUserData, _details.radiusMultiplier, container);
 
  // Sections (dendrites and axon)
  Neighbours somaNeighbours;

bioexplorer/backend/science/morphologies/Neurons.h

@@ -78,11 +78,6 @@ class Neurons : public Morphologies
  const uint64_t somaUserData, const double voltageScaling, common::ThreadSafeContainer& container,
  common::Neighbours& somaMeighbours, common::Neighbours& sectionNeighbours);
 
- double _addSomaAsSpheres(const uint64_t neuronId, const size_t somaMaterialId, const SectionMap& sections,
- const core::Vector3d& somaPosition, const core::Quaterniond& somaRotation,
- const double somaRadius, const uint64_t somaUserData,
- common::ThreadSafeContainer& container);
-
  void _buildOrientations(common::ThreadSafeContainer& container, const NeuronSomaMap& somas);
 
  void _buildMorphology(common::ThreadSafeContainer& container, const uint64_t neuronId, const NeuronSoma& soma,