Varied node radius for CurveNetwork.

include/polyscope/curve_network.h

@@ -124,6 +124,15 @@ class CurveNetwork : public QuantityStructure<CurveNetwork> {
   CurveNetwork* setColor(glm::vec3 newVal);
   glm::vec3 getColor();
 
+
+  // === Set radius from a scalar quantity
+  // effect is multiplicative with pointRadius
+  // negative values are always clamped to 0
+  // if autoScale==true, values are rescaled such that the largest has size pointRadius
+  void setNodeRadiusQuantity(CurveNetworkNodeScalarQuantity* quantity, bool autoScale = true);
+  void setNodeRadiusQuantity(std::string name, bool autoScale = true);
+  void clearNodeRadiusQuantity();
+
   // set the radius of the points
   CurveNetwork* setRadius(float newVal, bool isRelative = true);
   float getRadius();
@@ -168,6 +177,11 @@ class CurveNetwork : public QuantityStructure<CurveNetwork> {
   CurveNetworkNodeVectorQuantity* addNodeVectorQuantityImpl(std::string name, const std::vector<glm::vec3>& vectors, VectorType vectorType);
   CurveNetworkEdgeVectorQuantity* addEdgeVectorQuantityImpl(std::string name, const std::vector<glm::vec3>& vectors, VectorType vectorType);
   // clang-format on
+
+  // Manage varying node, edge size
+  std::string nodeRadiusQuantityName = ""; // e empty string means none
+  bool nodeRadiusQuantityAutoscale = true;
+  std::vector<double> resolveNodeRadiusQuantity(); // helper
 };
 
 

include/polyscope/render/opengl/shaders/cylinder_shaders.h

@@ -18,6 +18,7 @@ extern const ShaderReplacementRule CYLINDER_PROPAGATE_COLOR;
 extern const ShaderReplacementRule CYLINDER_PROPAGATE_BLEND_COLOR;
 extern const ShaderReplacementRule CYLINDER_PROPAGATE_PICK;
 extern const ShaderReplacementRule CYLINDER_CULLPOS_FROM_MID;
+extern const ShaderReplacementRule CYLINDER_VARIABLE_SIZE;
 
 
 } // namespace backend_openGL3_glfw

src/curve_network.cpp

@@ -56,7 +56,11 @@ void CurveNetwork::setCurveNetworkNodeUniforms(render::ShaderProgram& p) {
   glm::mat4 Pinv = glm::inverse(P);
   p.setUniform("u_invProjMatrix", glm::value_ptr(Pinv));
   p.setUniform("u_viewport", render::engine->getCurrentViewport());
-  p.setUniform("u_pointRadius", getRadius());
+  if (nodeRadiusQuantityName != "" && nodeRadiusQuantityAutoscale) {
+    p.setUniform("u_pointRadius", 1.); // u_pointRadius in sphere shader
+  } else {
+    p.setUniform("u_pointRadius", getRadius());
+  }
 }
 
 void CurveNetwork::setCurveNetworkEdgeUniforms(render::ShaderProgram& p) {
@@ -124,13 +128,23 @@ void CurveNetwork::drawPick() {
 
 std::vector<std::string> CurveNetwork::addCurveNetworkNodeRules(std::vector<std::string> initRules) {
   initRules = addStructureRules(initRules);
+
+  if (nodeRadiusQuantityName != "") {
+    initRules.push_back("SPHERE_VARIABLE_SIZE");
+  }
   if (wantsCullPosition()) {
     initRules.push_back("SPHERE_CULLPOS_FROM_CENTER");
   }
   return initRules;
 }
 std::vector<std::string> CurveNetwork::addCurveNetworkEdgeRules(std::vector<std::string> initRules) {
   initRules = addStructureRules(initRules);
+
+  // use node radius to blend cylinder radius
+  if (nodeRadiusQuantityName != "") {
+    initRules.push_back("CYLINDER_VARIABLE_SIZE");
+  }
+
   if (wantsCullPosition()) {
     initRules.push_back("CYLINDER_CULLPOS_FROM_MID");
   }
@@ -223,6 +237,12 @@ void CurveNetwork::preparePick() {
 
 void CurveNetwork::fillNodeGeometryBuffers(render::ShaderProgram& program) {
   program.setAttribute("a_position", nodes);
+
+  if (nodeRadiusQuantityName != "") {
+    // Resolve the quantity
+    std::vector<double> nodeRadiusQuantityVals = resolveNodeRadiusQuantity();
+    program.setAttribute("a_pointRadius", nodeRadiusQuantityVals);
+  }
 }
 
 void CurveNetwork::fillEdgeGeometryBuffers(render::ShaderProgram& program) {
@@ -239,6 +259,26 @@ void CurveNetwork::fillEdgeGeometryBuffers(render::ShaderProgram& program) {
   }
   program.setAttribute("a_position_tail", posTail);
   program.setAttribute("a_position_tip", posTip);
+
+
+  // blend the cylinder's tip/tail radii from the node radii
+  if (nodeRadiusQuantityName != "") {
+    // Resolve the quantity
+    std::vector<double> nodeRadiusQuantityVals = resolveNodeRadiusQuantity();
+    std::vector<double> cylinderTipRadiusQuantityVals(nEdges());
+    std::vector<double> cylinderTailRadiusQuantityVals(nEdges());
+
+    for (size_t iE = 0; iE < nEdges(); iE++) {
+      auto& edge = edges[iE];
+      size_t eTip = std::get<0>(edge);
+      size_t eTail = std::get<1>(edge);
+      cylinderTipRadiusQuantityVals[iE] = nodeRadiusQuantityVals[eTip];
+      cylinderTailRadiusQuantityVals[iE] = nodeRadiusQuantityVals[eTail];
+    }
+
+    program.setAttribute("a_tipRadius", cylinderTipRadiusQuantityVals);
+    program.setAttribute("a_tailRadius", cylinderTailRadiusQuantityVals);
+  }
 }
 
 void CurveNetwork::refresh() {
@@ -322,7 +362,6 @@ void CurveNetwork::buildEdgePickUI(size_t edgeInd) {
   ImGui::Indent(-20.);
 }
 
-
 void CurveNetwork::buildCustomUI() {
   ImGui::Text("nodes: %lld  edges: %lld", static_cast<long long int>(nNodes()), static_cast<long long int>(nEdges()));
   if (ImGui::ColorEdit3("Color", &color.get()[0], ImGuiColorEditFlags_NoInputs)) {
@@ -371,6 +410,24 @@ CurveNetwork* CurveNetwork::setColor(glm::vec3 newVal) {
 }
 glm::vec3 CurveNetwork::getColor() { return color.get(); }
 
+void CurveNetwork::setNodeRadiusQuantity(CurveNetworkNodeScalarQuantity* quantity, bool autoScale) {
+  setNodeRadiusQuantity(quantity->name, autoScale);
+}
+
+void CurveNetwork::setNodeRadiusQuantity(std::string name, bool autoScale) {
+  nodeRadiusQuantityName = name;
+  nodeRadiusQuantityAutoscale = autoScale;
+
+  resolveNodeRadiusQuantity(); // do it once, just so we fail fast if it doesn't exist
+
+  refresh(); // TODO this is a bit overkill
+}
+
+void CurveNetwork::clearNodeRadiusQuantity() {
+  nodeRadiusQuantityName = "";
+  refresh();
+};
+
 CurveNetwork* CurveNetwork::setRadius(float newVal, bool isRelative) {
   radius = ScaledValue<float>(newVal, isRelative);
   polyscope::requestRedraw();
@@ -445,5 +502,42 @@ CurveNetworkEdgeVectorQuantity* CurveNetwork::addEdgeVectorQuantityImpl(std::str
   return q;
 }
 
+std::vector<double> CurveNetwork::resolveNodeRadiusQuantity() {
+  CurveNetworkScalarQuantity* sizeScalarQ = nullptr;
+  CurveNetworkQuantity* sizeQ = getQuantity(nodeRadiusQuantityName);
+  if (sizeQ != nullptr) {
+    sizeScalarQ = dynamic_cast<CurveNetworkScalarQuantity*>(sizeQ);
+    if (sizeScalarQ == nullptr) {
+      polyscope::error("Cannot populate point size from quantity [" + name + "], it is not a scalar quantity");
+    }
+  } else {
+    polyscope::error("Cannot populate point size from quantity [" + name + "], it does not exist");
+  }
+
+  std::vector<double> sizes;
+  if (sizeScalarQ == nullptr || sizeScalarQ->values.size() != nNodes()) {
+    // we failed to resolve above; populate with dummy data so we can continue processing
+    std::vector<double> ones(nNodes(), 1.);
+    sizes = ones;
+    polyscope::error("quantity # != node #");
+  } else {
+    sizes = sizeScalarQ->values;
+  }
+
+  // clamp to nonnegative and autoscale (if requested)
+  double max = 0;
+  for (double& x : sizes) {
+    if (!(x > 0)) x = 0; // ensure all nonnegative
+    max = std::fmax(max, x);
+  }
+  if (max == 0) max = 1e-6;
+  if (nodeRadiusQuantityAutoscale) {
+    for (double& x : sizes) {
+      x /= max;
+    }
+  }
+
+  return sizes;
+}
 
 } // namespace polyscope

src/curve_network_scalar_quantity.cpp

@@ -79,8 +79,8 @@ void CurveNetworkNodeScalarQuantity::createProgram() {
       "RAYCAST_CYLINDER", addScalarRules(parent.addCurveNetworkEdgeRules({"CYLINDER_PROPAGATE_BLEND_VALUE"})));
 
   // Fill geometry buffers
-  parent.fillEdgeGeometryBuffers(*edgeProgram);
   parent.fillNodeGeometryBuffers(*nodeProgram);
+  parent.fillEdgeGeometryBuffers(*edgeProgram);
 
   { // Fill node color buffers
     nodeProgram->setAttribute("a_value", values);

src/render/mock_opengl/mock_gl_engine.cpp

@@ -1456,6 +1456,7 @@ void MockGLEngine::populateDefaultShadersAndRules() {
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_BLEND_COLOR", CYLINDER_PROPAGATE_BLEND_COLOR});
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_PICK", CYLINDER_PROPAGATE_PICK});
   registeredShaderRules.insert({"CYLINDER_CULLPOS_FROM_MID", CYLINDER_CULLPOS_FROM_MID});
+  registeredShaderRules.insert({"CYLINDER_VARIABLE_SIZE", CYLINDER_VARIABLE_SIZE});
 
   // marching tets things
   registeredShaderRules.insert({"SLICE_TETS_BASECOLOR_SHADE", SLICE_TETS_BASECOLOR_SHADE});

src/render/opengl/gl_engine.cpp

@@ -2144,6 +2144,7 @@ void GLEngine::populateDefaultShadersAndRules() {
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_BLEND_COLOR", CYLINDER_PROPAGATE_BLEND_COLOR});
   registeredShaderRules.insert({"CYLINDER_PROPAGATE_PICK", CYLINDER_PROPAGATE_PICK});
   registeredShaderRules.insert({"CYLINDER_CULLPOS_FROM_MID", CYLINDER_CULLPOS_FROM_MID});
+  registeredShaderRules.insert({"CYLINDER_VARIABLE_SIZE", CYLINDER_VARIABLE_SIZE});
 
   // marching tets things
   registeredShaderRules.insert({"SLICE_TETS_BASECOLOR_SHADE", SLICE_TETS_BASECOLOR_SHADE});

src/render/opengl/shaders/common.cpp

@@ -197,7 +197,7 @@ bool rayDiskIntersection(vec3 rayStart, vec3 rayDir, vec3 planePos, vec3 planeDi
   return true;
 }
 
-bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylTip, float cylRad, out float tHit, out vec3 pHit, out vec3 nHit) {
+bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylTip, float cylTipRad, float cylTailRad, out float tHit, out vec3 pHit, out vec3 nHit) {
 
     rayDir = normalize(rayDir);
     float cylLen = max(length(cylTip - cylTail), 1e-6);
@@ -209,7 +209,7 @@ bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylT
     vec3 pVec = o - dot(o, cylDir)*cylDir;
     float a = length2(qVec);
     float b = 2.0 * dot(qVec, pVec);
-    float c = length2(pVec) - cylRad*cylRad;
+    float c = length2(pVec) - cylTipRad*cylTailRad; // not sure about here. What are we doing in this section? Need comment for the root-solving.
     float disc = b*b - 4*a*c;
     if(disc < 0){
       tHit = LARGE_FLOAT();
@@ -240,7 +240,7 @@ bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylT
     float tHitTail;
     vec3 pHitTail;
     vec3 nHitTail;
-    rayDiskIntersection(rayStart, rayDir, cylTail, -cylDir, cylRad, tHitTail, pHitTail, nHitTail);
+    rayDiskIntersection(rayStart, rayDir, cylTail, -cylDir, cylTipRad, tHitTail, pHitTail, nHitTail);
     if(tHitTail < tHit) {
       tHit = tHitTail;
       pHit = pHitTail;
@@ -251,7 +251,7 @@ bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylT
     float tHitTip;
     vec3 pHitTip;
     vec3 nHitTip;
-    rayDiskIntersection(rayStart, rayDir, cylTip, cylDir, cylRad, tHitTip, pHitTip, nHitTip);
+    rayDiskIntersection(rayStart, rayDir, cylTip, cylDir, cylTailRad, tHitTip, pHitTip, nHitTip);
     if(tHitTip < tHit) {
       tHit = tHitTip;
       pHit = pHitTip;

src/render/opengl/shaders/cylinder_shaders.cpp

@@ -4,7 +4,7 @@
 
 namespace polyscope {
 namespace render {
-namespace backend_openGL3_glfw { 
+namespace backend_openGL3_glfw {
 
 // clang-format off
 
@@ -79,6 +79,9 @@ R"(
         void buildTangentBasis(vec3 unitNormal, out vec3 basisX, out vec3 basisY);
 
         void main() {
+            float tipRadius = u_radius;
+            float tailRadius = u_radius;
+            ${ CYLINDER_SET_RADIUS_GEOM }$
 
             // Build an orthogonal basis
             vec3 tailViewVal = gl_in[0].gl_Position.xyz / gl_in[0].gl_Position.w;
@@ -89,17 +92,19 @@ R"(
             // Compute corners of cube
             vec4 tailProj = u_projMatrix * gl_in[0].gl_Position;
             vec4 tipProj = u_projMatrix * position_tip[0];
-            vec4 dx = u_projMatrix * vec4(basisX * u_radius, 0.);
-            vec4 dy = u_projMatrix * vec4(basisY * u_radius, 0.);
-
-            vec4 p1 = tailProj - dx - dy;
-            vec4 p2 = tailProj + dx - dy;
-            vec4 p3 = tailProj - dx + dy;
-            vec4 p4 = tailProj + dx + dy;
-            vec4 p5 = tipProj - dx - dy;
-            vec4 p6 = tipProj + dx - dy;
-            vec4 p7 = tipProj - dx + dy;
-            vec4 p8 = tipProj + dx + dy;
+            vec4 dxTip = u_projMatrix * vec4(basisX * tipRadius, 0.);
+            vec4 dyTip = u_projMatrix * vec4(basisY * tipRadius, 0.);
+            vec4 dxTail = u_projMatrix * vec4(basisX * tailRadius, 0.);
+            vec4 dyTail = u_projMatrix * vec4(basisY * tailRadius, 0.);
+
+            vec4 p1 = tailProj - dxTail - dyTail;
+            vec4 p2 = tailProj + dxTail - dyTail;
+            vec4 p3 = tailProj - dxTail + dyTail;
+            vec4 p4 = tailProj + dxTail + dyTail;
+            vec4 p5 = tipProj - dxTip - dyTip;
+            vec4 p6 = tipProj + dxTip - dyTip;
+            vec4 p7 = tipProj - dxTip + dyTip;
+            vec4 p8 = tipProj + dxTip + dyTip;
 
             // Other data to emit   
 
@@ -158,7 +163,7 @@ R"(
 
         float LARGE_FLOAT();
         vec3 fragmentViewPosition(vec4 viewport, vec2 depthRange, mat4 invProjMat, vec4 fragCoord);
-        bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylTip, float cylRad, out float tHit, out vec3 pHit, out vec3 nHit);
+        bool rayCylinderIntersection(vec3 rayStart, vec3 rayDir, vec3 cylTail, vec3 cylTip, float cylTipRad, float cylTailRad, out float tHit, out vec3 pHit, out vec3 nHit);
         float fragDepthFromView(mat4 projMat, vec2 depthRange, vec3 viewPoint);
 
         ${ FRAG_DECLARATIONS }$
@@ -169,11 +174,16 @@ R"(
            vec2 depthRange = vec2(gl_DepthRange.near, gl_DepthRange.far);
            vec3 viewRay = fragmentViewPosition(u_viewport, depthRange, u_invProjMatrix, gl_FragCoord);
 
+
+           float tipRadius = u_radius;
+           float tailRadius = u_radius;
+           ${ CYLINDER_SET_RADIUS_FRAG }$
+
            // Raycast to the cylinder
            float tHit;
            vec3 pHit;
            vec3 nHit;
-           rayCylinderIntersection(vec3(0., 0., 0), viewRay, tailView, tipView, u_radius, tHit, pHit, nHit);
+           rayCylinderIntersection(vec3(0., 0., 0), viewRay, tailView, tipView, tipRadius, tailRadius, tHit, pHit, nHit);
            if(tHit >= LARGE_FLOAT()) {
               discard;
            }
@@ -421,6 +431,49 @@ const ShaderReplacementRule CYLINDER_PROPAGATE_PICK (
     /* textures */ {}
 );
 
+const ShaderReplacementRule CYLINDER_VARIABLE_SIZE (
+    /* rule name */ "CYLINDER_VARIABLE_SIZE",
+    { /* replacement sources */
+      {"VERT_DECLARATIONS", R"(
+          in float a_tipRadius;
+          in float a_tailRadius;
+          out float a_tipRadiusToGeom;
+          out float a_tailRadiusToGeom;
+        )"},
+      {"VERT_ASSIGNMENTS", R"(
+          a_tipRadiusToGeom = a_tipRadius;
+          a_tailRadiusToGeom = a_tailRadius;
+        )"},
+      {"GEOM_DECLARATIONS", R"(
+          in float a_tipRadiusToGeom[];
+          in float a_tailRadiusToGeom[];
+          out float a_tipRadiusToFrag;
+          out float a_tailRadiusToFrag;
+        )"},
+      {"GEOM_PER_EMIT", R"(
+          a_tipRadiusToFrag = a_tipRadiusToGeom[0]; 
+          a_tailRadiusToFrag = a_tailRadiusToGeom[0]; 
+        )"},
+      {"FRAG_DECLARATIONS", R"(
+          in float a_tipRadiusToFrag;
+          in float a_tailRadiusToFrag;
+        )"},
+      {"CYLINDER_SET_RADIUS_GEOM", R"(
+          tipRadius *= a_tipRadiusToGeom[0];
+          tailRadius *= a_tailRadiusToGeom[0];
+        )"},
+      {"CYLINDER_SET_RADIUS_FRAG", R"(
+          tipRadius *= a_tipRadiusToFrag;
+          tailRadius *= a_tailRadiusToFrag;
+        )"},
+    },
+    /* uniforms */ {},
+    /* attributes */ {
+      {"a_tipRadius", DataType::Float},
+      {"a_tailRadius", DataType::Float},
+    },
+    /* textures */ {}
+);
 
 // clang-format on