Fix different shading between CPU and GPU subdivision

Reuse the same vertex normals calculation as for the GPU code, by
weighing each vertex normals by the angle of the edges incident to the
vertex on the face.

Additionally, remove limit normals, as the CPU code does not use them
either, and would also cause different shading issues when limit surface
is used.

Fixes T95242: shade smooth artifacts with edge crease and limit surface
Fixes T94919: subdivision, different shading between CPU and GPU

4 files changed, 69 insertions, 20 deletions

diff --git a/source/blender/draw/intern/shaders/common_subdiv_lib.glsl b/source/blender/draw/intern/shaders/common_subdiv_lib.glsl
index 9dd86c35ee4..a620905622a 100644
--- a/source/blender/draw/intern/shaders/common_subdiv_lib.glsl
+++ b/source/blender/draw/intern/shaders/common_subdiv_lib.glsl
@@ -140,6 +140,13 @@ void set_vertex_nor(inout PosNorLoop vertex_data, vec3 nor)
   set_vertex_nor(vertex_data, nor, 0);
 }
 
+void add_newell_cross_v3_v3v3(inout vec3 n, vec3 v_prev, vec3 v_curr)
+{
+  n[0] += (v_prev[1] - v_curr[1]) * (v_prev[2] + v_curr[2]);
+  n[1] += (v_prev[2] - v_curr[2]) * (v_prev[0] + v_curr[0]);
+  n[2] += (v_prev[0] - v_curr[0]) * (v_prev[1] + v_curr[1]);
+}
+
 #define ORIGINDEX_NONE -1
 
 #ifdef SUBDIV_POLYGON_OFFSET
diff --git a/source/blender/draw/intern/shaders/common_subdiv_normals_accumulate_comp.glsl b/source/blender/draw/intern/shaders/common_subdiv_normals_accumulate_comp.glsl
index 575090472b1..a6d6c205615 100644
--- a/source/blender/draw/intern/shaders/common_subdiv_normals_accumulate_comp.glsl
+++ b/source/blender/draw/intern/shaders/common_subdiv_normals_accumulate_comp.glsl
@@ -16,11 +16,33 @@ layout(std430, binding = 2) readonly buffer faceAdjacencyLists
   uint face_adjacency_lists[];
 };
 
-layout(std430, binding = 3) writeonly buffer vertexNormals
+layout(std430, binding = 3) readonly buffer vertexLoopMap
+{
+  uint vert_loop_map[];
+};
+
+layout(std430, binding = 4) writeonly buffer vertexNormals
 {
   vec3 normals[];
 };
 
+void find_prev_and_next_vertex_on_face(
+    uint face_index, uint vertex_index, out uint curr, out uint next, out uint prev)
+{
+  uint start_loop_index = face_index * 4;
+
+  for (uint i = 0; i < 4; i++) {
+    uint subdiv_vert_index = vert_loop_map[start_loop_index + i];
+
+    if (subdiv_vert_index == vertex_index) {
+      curr = i;
+      next = (i + 1) % 4;
+      prev = (i + 4 - 1) % 4;
+      break;
+    }
+  }
+}
+
 void main()
 {
   uint vertex_index = get_global_invocation_index();
@@ -39,18 +61,37 @@ void main()
     uint adjacent_face = face_adjacency_lists[first_adjacent_face_offset + i];
     uint start_loop_index = adjacent_face * 4;
 
-    /* Compute face normal. */
-    vec3 adjacent_verts[3];
-    for (uint j = 0; j < 3; j++) {
-      adjacent_verts[j] = get_vertex_pos(pos_nor[start_loop_index + j]);
+    /* Compute the face normal using Newell's method. */
+    vec3 verts[4];
+    for (uint j = 0; j < 4; j++) {
+      verts[j] = get_vertex_pos(pos_nor[start_loop_index + j]);
     }
 
-    vec3 face_normal = normalize(
-        cross(adjacent_verts[1] - adjacent_verts[0], adjacent_verts[2] - adjacent_verts[0]));
-    accumulated_normal += face_normal;
+    vec3 face_normal = vec3(0.0);
+    add_newell_cross_v3_v3v3(face_normal, verts[0], verts[1]);
+    add_newell_cross_v3_v3v3(face_normal, verts[1], verts[2]);
+    add_newell_cross_v3_v3v3(face_normal, verts[2], verts[3]);
+    add_newell_cross_v3_v3v3(face_normal, verts[3], verts[0]);
+
+    /* Accumulate angle weighted normal. */
+    uint curr_vert = 0;
+    uint next_vert = 0;
+    uint prev_vert = 0;
+    find_prev_and_next_vertex_on_face(
+        adjacent_face, vertex_index, curr_vert, next_vert, prev_vert);
+
+    vec3 curr_co = verts[curr_vert];
+    vec3 prev_co = verts[next_vert];
+    vec3 next_co = verts[prev_vert];
+
+    vec3 edvec_prev = normalize(prev_co - curr_co);
+    vec3 edvec_next = normalize(curr_co - next_co);
+
+    float fac = acos(-dot(edvec_prev, edvec_next));
+
+    accumulated_normal += face_normal * fac;
   }
 
-  float weight = 1.0 / float(number_of_adjacent_faces);
   vec3 normal = normalize(accumulated_normal);
   normals[vertex_index] = normal;
 }
diff --git a/source/blender/draw/intern/shaders/common_subdiv_patch_evaluation_comp.glsl b/source/blender/draw/intern/shaders/common_subdiv_patch_evaluation_comp.glsl
index 5dd7decf663..26410eb2d7f 100644
--- a/source/blender/draw/intern/shaders/common_subdiv_patch_evaluation_comp.glsl
+++ b/source/blender/draw/intern/shaders/common_subdiv_patch_evaluation_comp.glsl
@@ -394,12 +394,8 @@ void main()
 
     evaluate_patches_limits(patch_co.patch_index, uv.x, uv.y, pos, du, dv);
 
-#  if defined(LIMIT_NORMALS)
-    vec3 nor = normalize(cross(du, dv));
-#  else
     /* This will be computed later. */
     vec3 nor = vec3(0.0);
-#  endif
 
     int origindex = input_vert_origindex[loop_index];
     uint flag = 0;
diff --git a/source/blender/draw/intern/shaders/common_subdiv_vbo_lnor_comp.glsl b/source/blender/draw/intern/shaders/common_subdiv_vbo_lnor_comp.glsl
index 41a8df3cf82..8b4a5dd14f3 100644
--- a/source/blender/draw/intern/shaders/common_subdiv_vbo_lnor_comp.glsl
+++ b/source/blender/draw/intern/shaders/common_subdiv_vbo_lnor_comp.glsl
@@ -38,13 +38,18 @@ void main()
     }
   }
   else {
-    /* Face is flat shaded, compute flat face normal from an inscribed triangle. */
-    vec3 verts[3];
-    for (int i = 0; i < 3; i++) {
-      verts[i] = get_vertex_pos(pos_nor[start_loop_index + i]);
-    }
-
-    vec3 face_normal = normalize(cross(verts[1] - verts[0], verts[2] - verts[0]));
+    vec3 v0 = get_vertex_pos(pos_nor[start_loop_index + 0]);
+    vec3 v1 = get_vertex_pos(pos_nor[start_loop_index + 1]);
+    vec3 v2 = get_vertex_pos(pos_nor[start_loop_index + 2]);
+    vec3 v3 = get_vertex_pos(pos_nor[start_loop_index + 3]);
+
+    vec3 face_normal = vec3(0.0);
+    add_newell_cross_v3_v3v3(face_normal, v0, v1);
+    add_newell_cross_v3_v3v3(face_normal, v1, v2);
+    add_newell_cross_v3_v3v3(face_normal, v2, v3);
+    add_newell_cross_v3_v3v3(face_normal, v3, v0);
+
+    face_normal = normalize(face_normal);
     for (int i = 0; i < 4; i++) {
       output_lnor[start_loop_index + i] = face_normal;
     }