Spherical Harmonic static example

Author: lenixlobo

fury/actor.py

@@ -2498,3 +2498,66 @@ def texture_on_sphere(rgb, theta=60, phi=60, interpolate=True):
     earthActor.SetTexture(atext)
 
     return earthActor
+
+
+def sdf_sh(centers, directions=(1, 0, 0), colors=(255, 0, 0), scale=1):
+    """Create a SDF actor
+
+    Parameters
+    ----------
+    centers : ndarray, shape (N, 3)
+        SDF primitive positions
+    colors : ndarray (N,3) or (N, 4) or tuple (3,) or tuple (4,)
+        RGB or RGBA (for opacity) R, G, B and A should be at the range [0, 1]
+    directions : ndarray, shape (N, 3)
+        The orientation vector of the cube.
+    scale : float
+        The size of the cube
+
+
+    Returns
+    -------
+    vtkActor
+    """
+
+    verts, faces = fp.prim_box()
+
+    repeated = fp.repeat_primitive(verts, faces, centers=centers,
+                                   colors=colors, directions=directions,
+                                   scale=scale)
+
+    rep_verts, rep_faces, rep_colors, rep_centers = repeated
+    box_actor = get_actor_from_primitive(rep_verts, rep_faces, rep_colors)
+
+    vtk_center = numpy_support.numpy_to_vtk(rep_centers)
+    vtk_center.SetNumberOfComponents(3)
+    vtk_center.SetName("center")
+    box_actor.GetMapper().GetInput().GetPointData().AddArray(vtk_center)
+
+    vs_dec_code = fs.load("sdf_dec.vert")
+    vs_impl_code = fs.load("sdf_impl.vert")
+    fs_dec_code = fs.load("sdfsh_dec.frag")
+    fs_impl_code = fs.load("sdfsh_impl.frag")
+
+    mapper = box_actor.GetMapper()
+    mapper.MapDataArrayToVertexAttribute(
+        "center", "center", vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, -1)
+
+    mapper.AddShaderReplacement(
+        vtk.vtkShader.Vertex, "//VTK::ValuePass::Dec", True,
+        vs_dec_code, False)
+
+    mapper.AddShaderReplacement(
+        vtk.vtkShader.Vertex, "//VTK::ValuePass::Impl", True,
+        vs_impl_code, False)
+
+    mapper.AddShaderReplacement(
+        vtk.vtkShader.Fragment, "//VTK::ValuePass::Dec", True,
+        fs_dec_code, False)
+
+    mapper.AddShaderReplacement(
+        vtk.vtkShader.Fragment, "//VTK::Light::Impl", True,
+        fs_impl_code, False)
+
+    return box_actor
+

fury/shaders/sdf_dec.vert

@@ -0,0 +1,8 @@
+/* SDF vertex shader declaration */
+
+//VTK::ValuePass::Dec
+
+in vec3 center;
+out vec4 vertexMCVSOutput;
+
+out vec3 centerWCVSOutput;

fury/shaders/sdf_impl.vert

@@ -0,0 +1,6 @@
+/* SDF vertex shader implementation */
+
+//VTK::ValuePass::Impl
+
+vertexMCVSOutput = vertexMC;
+centerWCVSOutput = center;

fury/shaders/sdfsh_dec.frag

@@ -0,0 +1,219 @@
+/* SDF fragment shader declaration */
+
+
+in vec4 vertexMCVSOutput;
+in vec3 centerWCVSOutput;
+
+uniform mat4 MCVCMatrix;
+uniform mat4 MCWCMatrix;
+uniform mat3 WCVCMatrix;
+
+
+#define PI 3.1415926535898
+#define E 2.7182818284591
+
+//DIPY default basis 
+
+
+#define c01 2.82064703e-01
+#define c02 6.90495017e-02 
+#define c03 -7.97628948e-02
+#define c04 9.29992151e-02 
+#define c05 -1.98478828e-06 
+#define c06 2.22295059e-06 
+#define c07 2.21050672e-02 
+#define c08 -3.61304561e-02 
+#define c09 1.67251336e-02 
+#define c10 2.27264752e-02 
+#define c11 3.77708897e-02 
+#define c12 -3.55761267e-06 
+#define c13 -2.24799994e-05 
+#define c14 4.40539515e-06 
+#define c15 2.36901715e-06
+
+
+float P(int l, int m, float x)
+{
+    float pmm = 1.0;
+    if (m > 0)
+    {
+        float somx2 = sqrt((1.0-x)*(1.0+x));
+        float fact  = 1.0;
+        for (int i = 1; i <= m; i++)
+        {
+            pmm *= (-fact) * somx2;
+            fact += 2.0;
+        }
+    }
+    if (l == m) return pmm;
+    float pmmp1 = x * float(2*m + 1) * pmm;
+    if (l == m+1) return pmmp1;
+    float pll = 0.0;
+    for (int ll = m+2; ll <= l; ll++)
+    {
+        pll = (float(2*ll-1)*x*pmmp1-float(ll+m-1)*pmm) / float(ll - m);
+        pmm = pmmp1;
+        pmmp1 = pll;
+    }
+    return pll;
+}
+
+
+
+// Clenshaw Legendre normalized
+float Pgn(int l, int m, float x)
+{
+	float p0 = 0.;    
+	float p1 = 0.;
+	float p2 = 0.;
+    
+    for (int k = l; k >= 0; k--)
+	{
+		float k1 = float(k + 1);
+		float m1 = float(2 * m) + k1;
+        float m2 = float(2 * (m + k) + 1);
+        
+		p2 = p1;
+        p1 = p0;
+        
+		p0 = 0.;
+        if (l == m + k)
+            p0 += 1.;
+        
+        float u0 = sqrt(
+			(m2 * (m2 + 2.0)) /
+			(k1 * m1)
+		);
+        
+        float u1 = sqrt(
+            (k1 * m1 * (m2 + 4.0)) / 
+			((k1 + 1.0) * (m1 + 1.0) * m2)
+        );
+        p0 += p1 * u0 * x;
+        p0 += -u1 * p2;
+	}
+
+    for (int k = 1; k <= m; k++){
+		p0 *= sqrt(
+            (1.0 - 0.5/float(k)) * (1.0 - x) * (1.0 + x)
+        );
+	}
+    
+    p0 *= sqrt((0.5 * float(m) + 0.25)/PI);
+    return p0;
+}
+
+
+
+float gammaln(float x){
+	return log(x);
+}
+
+float SHH(in int l, in int m, in vec3 s )
+{
+	vec3 ns = normalize(s);
+	float val = P(l, m, ns.x);
+	val *= sqrt( (2*m + 1) / 4.0 / PI);
+	val *= exp(0.5 * (gammaln(l - m + 1) - gammaln(l + m + 1)));
+	val = val * exp(-1 * m * ns.y);
+	return val;
+}
+
+
+
+// Y_l_m(s), where l is the band and m the range in [-l..l] 
+float SH( in int l, in int m, in vec3 s ) 
+{
+	vec3 ns = s;
+    //vec3 ns = normalize(s);
+    
+    if (m < 0) {
+        float c = ns.x;
+        ns.x = ns.z;
+        ns.z = c;
+        m = -m;
+    }
+    
+    // spherical coordinates
+    float thetax = ns.y;
+    float phi = atan(ns.z, ns.x);
+    
+    float pl = Pgn(l, m, thetax);
+    
+    float r = pow(-1.0, float(m)) * cos(float(m) * phi) * pl;
+    
+    return r;
+}
+
+vec3 map( in vec3 position)
+{
+	vec3 p00 = position - centerWCVSOutput;
+	float r, d;
+	vec3 n, s, res;
+
+	d =  length(p00);
+	n = p00/d;
+
+
+	r =  c01 * SH(0, 0, n);
+	r  += c02 * SH(2, -2, n);
+	r  += c03 * SH(2, -1, n); 
+	r  += c04 * SH(2, 0, n);
+	r  += c05 * SH(2, 1, n);
+	r  += c06 * SH(2, 2, n);
+	r  += c07 * SH(4, -4, n);
+	r  += c08 * SH(4, -3, n);
+	r  += c09 * SH(4, -2, n);
+	r  += c10 * SH(4, -1, n);
+	r  += c11 * SH(4, 0, n);
+	r  += c12 * SH(4, 1, n);
+	r  += c13 * SH(4, 2, n);
+	r  += c14 * SH(4, 3, n);
+	r += c15 * SH(4, 4, n);
+
+	#define SHAPE (vec3(d-abs(r), sign(r),d))
+
+	s = SHAPE;
+	res = s;
+
+	return vec3(res.x, 0.5 + 0.5 * res.y, res.z);
+}
+
+vec3 calculateNormal( in vec3 pos )
+{
+    vec2 e = vec2(1.0,-1.0)*0.5773*0.0005;
+    return normalize( e.xyy*map( pos + e.xyy ).x + 
+                      e.yyx*map( pos + e.yyx ).x + 
+                      e.yxy*map( pos + e.yxy ).x + 
+                      e.xxx*map( pos + e.xxx ).x );
+
+}
+
+
+vec3 castRay(in vec3 ro, vec3 rd)
+{
+	vec3 res  = vec3(1e10, -1.0, 1.0);
+
+	float t = 0.0;
+	float h = 1.0;
+	vec2 m = vec2(-1.0);
+
+	for(int i=0;i<200;i++){
+
+		if(h<0.001) break;
+		
+		vec3 position = ro+t*rd;
+
+		vec3 res = map( position );
+		
+		h = res.x;
+		m = res.yz;
+		t += h*0.3;
+	}
+
+		if( t<res.x ) res = vec3(t, m);
+
+
+		return res;
+}
+

fury/shaders/sdfsh_impl.frag

@@ -0,0 +1,42 @@
+/* SDF fragment shader implementation */
+
+//VKT::Light::Impl
+
+vec3 point = vertexMCVSOutput.xyz;
+
+//ray origin
+vec4 ro = -MCVCMatrix[3] * MCVCMatrix;  // camera position in world space
+
+vec3 col = vertexColorVSOutput.rgb;
+
+//ray direction
+vec3 rd = normalize(point - ro.xyz);
+
+
+vec3 t = castRay(ro.xyz, rd);
+    
+if(t.y > -0.05)
+{
+
+	vec3 pos = ro.xyz + t.x*rd;
+   	vec3 normal = calculateNormal(pos);
+    vec3 ref = reflect( rd, normal );
+
+
+    float occ = clamp( 2.0*t.z, 0.0, 1.0 );
+    float sss = pow( clamp( 1.0 + dot(normal, rd), 0.0, 1.0 ), 1.0 );
+
+    // lights
+    vec3 lin  = 2.5*occ*vec3(1.0,1.00,1.00)*(0.6+0.4*normal.y);
+    lin += 1.0*sss*vec3(1.0,0.95,0.70)*occ;	
+
+    vec3 mater = 0.5*mix( vec3(1.0,0.6,0.15), vec3(0.2,0.4,0.5), t.y ); 	
+
+    fragOutput0 = vec4( mater * lin , 1.0);
+    //fragOutput0 = vec4( norm, 1.0);
+    	
+}
+else{
+    //fragOutput0 = vec4(0, 0, 0, 0.3);
+    discard;
+}

test_sh.py

@@ -0,0 +1,49 @@
+from fury import actor, window, ui
+import numpy as np
+import itertools
+
+
+tb = ui.TextBlock2D(font_size=20, color=(1, 0.5, 1))
+panel = ui.Panel2D(position=(400,400), size = (500, 150))
+panel.add_element(tb, (0.2, 0.5))
+
+
+scene = window.Scene()
+scene.background((1.0, 0.8, 0.8))
+# dirs = np.random.rand(7, 3)
+# colors = np.random.rand(7, 3) * 255
+# centers = np.array([[2, 0, 0], [0, 0, 0], [-2, 0, 0], [0, 2, 0], [0, -2, 0],
+#                     [0, 0, -2], [0, 0, 2]])
+centers = np.array([[0, 0, 0]])
+sdfactor = actor.sdf_sh(centers=centers)
+
+fpss = []
+
+counter = itertools.count()
+showm = window.ShowManager(scene, reset_camera=False )
+
+showm.initialize()
+
+
+scene.add(sdfactor)
+scene.add(panel)
+
+
+def timer_callback(_obj, _event):
+	cnt = next(counter)
+	showm.render()
+
+	if cnt % 1 == 0:
+		fps = scene.frame_rate
+		fpss.append(fps)
+		msg = "FPS: " + str(fps) 
+		tb.message = msg
+		
+	if(cnt==2000):
+		showm.exit()
+
+
+
+
+showm.add_timer_callback(True, 10, timer_callback)
+showm.start()