Cycles: Adaptive subdivision smooth UV interpolation

Using OpenSubdiv FVar interpolation.

Pull Request: https://projects.blender.org/blender/blender/pulls/135681

intern/cycles/blender/mesh.cpp

@@ -539,6 +539,8 @@ static void attr_create_subd_uv_map(Scene *scene,
           uv_attr = mesh->subd_attributes.add(uv_name, TypeFloat2, ATTR_ELEMENT_CORNER);
         }
 
+        uv_attr->flags |= ATTR_SUBDIVIDE_SMOOTH_FVAR;
+
         const blender::VArraySpan b_uv_map = *b_attributes.lookup<blender::float2>(
             uv_name.c_str(), blender::bke::AttrDomain::Corner);
         float2 *fdata = uv_attr->data_float2();
@@ -808,8 +810,7 @@ static void create_mesh(Scene *scene,
                         const array<Node *> &used_shaders,
                         const bool need_motion,
                         const float motion_scale,
-                        const bool subdivision = false,
-                        const bool subdivide_uvs = true)
+                        const bool subdivision = false)
 {
   const blender::Span<blender::float3> positions = b_mesh.vert_positions();
   const blender::OffsetIndices faces = b_mesh.faces();
@@ -1043,9 +1044,8 @@ static void create_subd_mesh(Scene *scene,
   BL::Object b_ob = b_ob_info.real_object;
 
   BL::SubsurfModifier subsurf_mod(b_ob.modifiers[b_ob.modifiers.length() - 1]);
-  const bool subdivide_uvs = subsurf_mod.uv_smooth() != BL::SubsurfModifier::uv_smooth_NONE;
 
-  create_mesh(scene, mesh, b_mesh, used_shaders, need_motion, motion_scale, true, subdivide_uvs);
+  create_mesh(scene, mesh, b_mesh, used_shaders, need_motion, motion_scale, true);
 
   const blender::VArraySpan creases = *b_mesh.attributes().lookup<float>(
       "crease_edge", blender::bke::AttrDomain::Edge);

intern/cycles/subd/interpolation.cpp

@@ -195,6 +195,100 @@ void SubdAttributeInterpolation::interp_attribute_vertex_linear(const Attribute
   }
 }
 
+#ifdef WITH_OPENSUBDIV
+template<typename T>
+void SubdAttributeInterpolation::interp_attribute_vertex_smooth(const Attribute &subd_attr,
+                                                                Attribute &mesh_attr,
+                                                                const int motion_step)
+{
+  // TODO: Avoid computing derivative weights when not needed
+  // TODO: overhead of FindPatch and EvaluateBasis with vertex position
+  const int num_refiner_verts = osd_data.refiner->GetNumVerticesTotal();
+  const int num_local_points = osd_data.patch_table->GetNumLocalPoints();
+  const int num_base_verts = mesh.get_num_subd_base_verts();
+
+  /* Refine attribute data to get patch coordinates. */
+  array<typename T::AccumType> refined_array(num_refiner_verts + num_local_points);
+
+  const typename T::Type *base_src = reinterpret_cast<const typename T::Type *>(subd_attr.data()) +
+                                     num_base_verts * motion_step;
+  typename T::AccumType *base_dst = refined_array.data();
+  for (int i = 0; i < num_base_verts; i++) {
+    base_dst[i] = T::read(base_src[i]);
+  }
+
+  Far::PrimvarRefiner primvar_refiner(*osd_data.refiner);
+  typename T::AccumType *src = refined_array.data();
+  for (int i = 0; i < osd_data.refiner->GetMaxLevel(); i++) {
+    typename T::AccumType *dest = src + osd_data.refiner->GetLevel(i).GetNumVertices();
+    primvar_refiner.Interpolate(
+        i + 1, (OsdValue<typename T::AccumType> *)src, (OsdValue<typename T::AccumType> *&)dest);
+    src = dest;
+  }
+
+  if (num_local_points) {
+    osd_data.patch_table->ComputeLocalPointValues(
+        (OsdValue<typename T::AccumType> *)refined_array.data(),
+        (OsdValue<typename T::AccumType> *)(refined_array.data() + num_refiner_verts));
+  }
+
+  /* Evaluate patches at limit. */
+  const size_t triangles_size = mesh.num_triangles();
+  const int *patch_index = mesh.subd_triangle_patch_index.data();
+  const float2 *patch_uv = mesh.subd_corner_patch_uv.data();
+  const typename T::AccumType *subd_data = refined_array.data();
+  typename T::Type *mesh_data = reinterpret_cast<typename T::Type *>(mesh_attr.data()) +
+                                mesh.get_verts().size() * motion_step;
+
+  /* Compute motion normals alongside positions. */
+  float3 *mesh_normal_data = nullptr;
+  if constexpr (std::is_same_v<typename T::Type, float3>) {
+    if (mesh_attr.std == ATTR_STD_MOTION_VERTEX_POSITION) {
+      Attribute *attr_normal = mesh.attributes.add(ATTR_STD_MOTION_VERTEX_NORMAL);
+      mesh_normal_data = attr_normal->data_float3() + mesh.get_verts().size() * motion_step;
+    }
+  }
+
+  for (size_t i = 0; i < triangles_size; i++) {
+    const int p = patch_index[i];
+    Mesh::Triangle triangle = mesh.get_triangle(i);
+
+    for (int j = 0; j < 3; j++) {
+      /* Compute patch weights. */
+      const float2 uv = patch_uv[(i * 3) + j];
+      const Far::PatchTable::PatchHandle &handle = *osd_data.patch_map->FindPatch(
+          p, (double)uv.x, (double)uv.y);
+
+      float p_weights[20], du_weights[20], dv_weights[20];
+      osd_data.patch_table->EvaluateBasis(handle, uv.x, uv.y, p_weights, du_weights, dv_weights);
+      Far::ConstIndexArray cv = osd_data.patch_table->GetPatchVertices(handle);
+
+      /* Compution position. */
+      typename T::AccumType value = subd_data[cv[0]] * p_weights[0];
+      for (int k = 1; k < cv.size(); k++) {
+        value += subd_data[cv[k]] * p_weights[k];
+      }
+      mesh_data[triangle.v[j]] = T::output(value);
+
+      /* Optionally compute normal. */
+      if constexpr (std::is_same_v<typename T::Type, float3>) {
+        if (mesh_normal_data) {
+          float3 du = zero_float3();
+          float3 dv = zero_float3();
+          for (int k = 0; k < cv.size(); k++) {
+            const float3 p = subd_data[cv[k]];
+            du += p * du_weights[k];
+            dv += p * dv_weights[k];
+          }
+          mesh_normal_data[triangle.v[j]] = safe_normalize_fallback(cross(du, dv),
+                                                                    make_float3(0.0f, 0.0f, 1.0f));
+        }
+      }
+    }
+  }
+}
+#endif
+
 template<typename T>
 void SubdAttributeInterpolation::interp_attribute_corner_linear(const Attribute &subd_attr,
                                                                 Attribute &mesh_attr)
@@ -259,6 +353,85 @@ void SubdAttributeInterpolation::interp_attribute_corner_linear(const Attribute
   }
 }
 
+#ifdef WITH_OPENSUBDIV
+template<typename T>
+void SubdAttributeInterpolation::interp_attribute_corner_smooth(Attribute &mesh_attr,
+                                                                const int channel,
+                                                                const vector<char> &merged_values)
+{
+  // TODO: Avoid computing derivative weights when not needed
+  const int num_refiner_fvars = osd_data.refiner->GetNumFVarValuesTotal(channel);
+  const int num_local_points = osd_data.patch_table->GetNumLocalPointsFaceVarying(channel);
+  const int num_base_fvars = osd_data.refiner->GetLevel(0).GetNumFVarValues(channel);
+
+  /* Refine attribute data to get patch coordinates. */
+  array<typename T::AccumType> refined_array(num_refiner_fvars + num_local_points);
+
+  const typename T::Type *base_src = reinterpret_cast<const typename T::Type *>(
+      merged_values.data());
+  typename T::AccumType *base_dst = refined_array.data();
+  for (int i = 0; i < num_base_fvars; i++) {
+    base_dst[i] = T::read(base_src[i]);
+  }
+
+  Far::PrimvarRefiner primvar_refiner(*osd_data.refiner);
+  typename T::AccumType *src = refined_array.data();
+  for (int i = 0; i < osd_data.refiner->GetMaxLevel(); i++) {
+    typename T::AccumType *dest = src + osd_data.refiner->GetLevel(i).GetNumFVarValues(channel);
+    primvar_refiner.InterpolateFaceVarying(i + 1,
+                                           (OsdValue<typename T::AccumType> *)src,
+                                           (OsdValue<typename T::AccumType> *&)dest,
+                                           channel);
+    src = dest;
+  }
+
+  if (num_local_points) {
+    osd_data.patch_table->ComputeLocalPointValuesFaceVarying(
+        (OsdValue<typename T::AccumType> *)refined_array.data(),
+        (OsdValue<typename T::AccumType> *)(refined_array.data() + num_refiner_fvars),
+        channel);
+  }
+
+  /* Evaluate patches at limit. */
+  const size_t triangles_size = mesh.num_triangles();
+  const int *patch_index = mesh.subd_triangle_patch_index.data();
+  const float2 *patch_uv = mesh.subd_corner_patch_uv.data();
+  const typename T::AccumType *subd_data = refined_array.data();
+  typename T::Type *mesh_data = reinterpret_cast<typename T::Type *>(mesh_attr.data());
+
+  for (size_t i = 0; i < triangles_size; i++) {
+    const int p = patch_index[i];
+
+    for (int j = 0; j < 3; j++) {
+      /* Compute patch weights. */
+      const float2 uv = patch_uv[(i * 3) + j];
+      const Far::PatchTable::PatchHandle &handle = *osd_data.patch_map->FindPatch(
+          p, (double)uv.x, (double)uv.y);
+
+      float p_weights[20], du_weights[20], dv_weights[20];
+      osd_data.patch_table->EvaluateBasisFaceVarying(handle,
+                                                     uv.x,
+                                                     uv.y,
+                                                     p_weights,
+                                                     du_weights,
+                                                     dv_weights,
+                                                     nullptr,
+                                                     nullptr,
+                                                     nullptr,
+                                                     channel);
+      Far::ConstIndexArray cv = osd_data.patch_table->GetPatchFVarValues(handle, channel);
+
+      /* Compution position. */
+      typename T::AccumType value = subd_data[cv[0]] * p_weights[0];
+      for (int k = 1; k < cv.size(); k++) {
+        value += subd_data[cv[k]] * p_weights[k];
+      }
+      mesh_data[(i * 3) + j] = T::output(value);
+    }
+  }
+}
+#endif
+
 template<typename T>
 void SubdAttributeInterpolation::interp_attribute_face(const Attribute &subd_attr,
                                                        Attribute &mesh_attr)
@@ -284,11 +457,42 @@ void SubdAttributeInterpolation::interp_attribute_type(const Attribute &subd_att
 {
   switch (subd_attr.element) {
     case ATTR_ELEMENT_VERTEX: {
-      interp_attribute_vertex_linear<T>(subd_attr, mesh_attr);
+#ifdef WITH_OPENSUBDIV
+      if (mesh.get_subdivision_type() == Mesh::SUBDIVISION_CATMULL_CLARK) {
+        /* Only smoothly interpolation known position-like attributes. */
+        switch (subd_attr.std) {
+          case ATTR_STD_GENERATED:
+          case ATTR_STD_POSITION_UNDEFORMED:
+          case ATTR_STD_POSITION_UNDISPLACED:
+            interp_attribute_vertex_smooth<T>(subd_attr, mesh_attr);
+            break;
+          default:
+            interp_attribute_vertex_linear<T>(subd_attr, mesh_attr);
+            break;
+        }
+      }
+      else
+#endif
+      {
+        interp_attribute_vertex_linear<T>(subd_attr, mesh_attr);
+      }
       break;
     }
     case ATTR_ELEMENT_CORNER:
     case ATTR_ELEMENT_CORNER_BYTE: {
+#ifdef WITH_OPENSUBDIV
+      if (osd_mesh.use_smooth_fvar(subd_attr)) {
+        for (const auto &merged_fvar : osd_mesh.merged_fvars) {
+          if (&merged_fvar.attr == &subd_attr) {
+            if constexpr (std::is_same_v<typename T::Type, float2>) {
+              interp_attribute_corner_smooth<T>(
+                  mesh_attr, merged_fvar.channel, merged_fvar.values);
+              return;
+            }
+          }
+        }
+      }
+#endif
       interp_attribute_corner_linear<T>(subd_attr, mesh_attr);
       break;
     }

intern/cycles/subd/interpolation.h

@@ -57,6 +57,17 @@ class SubdAttributeInterpolation {
   template<typename T>
   void interp_attribute_corner_linear(const Attribute &subd_attr, Attribute &mesh_attr);
 
+#ifdef WITH_OPENSUBDIV
+  template<typename T>
+  void interp_attribute_vertex_smooth(const Attribute &subd_attr,
+                                      Attribute &mesh_attr,
+                                      const int motion_step = 0);
+  template<typename T>
+  void interp_attribute_corner_smooth(Attribute &mesh_attr,
+                                      const int channel,
+                                      const vector<char> &merged_values);
+#endif
+
   template<typename T>
   void interp_attribute_face(const Attribute &subd_attr, Attribute &mesh_attr);
 

intern/cycles/subd/osd.cpp

@@ -126,10 +126,115 @@ bool TopologyRefinerFactory<OsdMesh>::assignComponentTags(TopologyRefiner &refin
   return true;
 }
 
-template<>
-bool TopologyRefinerFactory<OsdMesh>::assignFaceVaryingTopology(TopologyRefiner & /*refiner*/,
-                                                                OsdMesh const & /*osd_mesh*/)
+template<typename T>
+static void merge_smooth_fvar(const Mesh &mesh,
+                              const Attribute &subd_attr,
+                              OsdMesh::MergedFVar &merged_fvar,
+                              vector<int> &merged_next,
+                              vector<int> &merged_face_corners)
 {
+  const int num_base_verts = mesh.get_num_subd_base_verts();
+  const int num_base_faces = mesh.get_num_subd_faces();
+  const int *subd_face_corners = mesh.get_subd_face_corners().data();
+
+  const T *values = reinterpret_cast<const T *>(subd_attr.data());
+
+  merged_fvar.values.resize(num_base_verts * sizeof(T));
+
+  // Merge identical corner values with the same vertex. The first value is stored at the vertex
+  // index, and any different values are pushed backed onto the array. merged_next creates a
+  // linked list between all values for the same vertex.
+  const int state_uninitialized = 0;
+  const int state_end = -1;
+  merged_next.resize(num_base_verts, state_uninitialized);
+
+  for (int f = 0, i = 0; f < num_base_faces; f++) {
+    Mesh::SubdFace face = mesh.get_subd_face(f);
+
+    for (int corner = 0; corner < face.num_corners; corner++) {
+      int v = subd_face_corners[face.start_corner + corner];
+      const T value = values[i++];
+
+      if (merged_next[v] == state_uninitialized) {
+        // First corner to initialize vertex.
+        reinterpret_cast<T *>(merged_fvar.values.data())[v] = value;
+        merged_next[v] = state_end;
+        merged_face_corners.push_back(v);
+      }
+      else {
+        // Find vertex with matching value, following linked list per vertex.
+        int v_prev = v;
+        for (; v != state_end; v_prev = v, v = merged_next[v]) {
+          if (reinterpret_cast<T *>(merged_fvar.values.data())[v] == value) {
+            // Matching value found, reuse merged vertex.
+            merged_face_corners.push_back(v);
+            break;
+          }
+        }
+
+        if (v == state_end) {
+          // Non-matching value, add new merged vertex and add to linked list.
+          const int next = merged_next.size();
+          merged_fvar.values.resize((next + 1) * sizeof(T));
+          reinterpret_cast<T *>(merged_fvar.values.data())[next] = value;
+          merged_next.push_back(state_end);
+          merged_next[v_prev] = next;
+          merged_face_corners.push_back(next);
+        }
+      }
+    }
+  }
+}
+
+template<>
+bool TopologyRefinerFactory<OsdMesh>::assignFaceVaryingTopology(TopologyRefiner &refiner,
+                                                                OsdMesh const &osd_mesh)
+{
+  const Mesh &mesh = osd_mesh.mesh;
+  auto &merged_fvars = const_cast<OsdMesh &>(osd_mesh).merged_fvars;
+
+  for (const Attribute &subd_attr : mesh.subd_attributes.attributes) {
+    if (!osd_mesh.use_smooth_fvar(subd_attr)) {
+      continue;
+    }
+
+    // Created merged FVar, for use in subdivide_attribute_corner_smooth.
+    OsdMesh::MergedFVar merged_fvar{subd_attr};
+
+    vector<int> merged_next;
+    vector<int> merged_face_corners;
+
+    if (subd_attr.element == ATTR_ELEMENT_CORNER_BYTE) {
+      merge_smooth_fvar<uchar4>(mesh, subd_attr, merged_fvar, merged_next, merged_face_corners);
+    }
+    else if (Attribute::same_storage(subd_attr.type, TypeFloat)) {
+      merge_smooth_fvar<float>(mesh, subd_attr, merged_fvar, merged_next, merged_face_corners);
+    }
+    else if (Attribute::same_storage(subd_attr.type, TypeFloat2)) {
+      merge_smooth_fvar<float2>(mesh, subd_attr, merged_fvar, merged_next, merged_face_corners);
+    }
+    else if (Attribute::same_storage(subd_attr.type, TypeVector)) {
+      merge_smooth_fvar<float3>(mesh, subd_attr, merged_fvar, merged_next, merged_face_corners);
+    }
+    else if (Attribute::same_storage(subd_attr.type, TypeFloat4)) {
+      merge_smooth_fvar<float4>(mesh, subd_attr, merged_fvar, merged_next, merged_face_corners);
+    }
+
+    // Create FVar channel and topology for OpenUSD.
+    merged_fvar.channel = createBaseFVarChannel(refiner, merged_next.size());
+
+    const int num_base_faces = mesh.get_num_subd_faces();
+    for (int f = 0, i = 0; f < num_base_faces; f++) {
+      Far::IndexArray dst_face_uvs = getBaseFaceFVarValues(refiner, f, merged_fvar.channel);
+      const int num_corners = dst_face_uvs.size();
+      for (int corner = 0; corner < num_corners; corner++) {
+        dst_face_uvs[corner] = merged_face_corners[i++];
+      }
+    }
+
+    merged_fvars.push_back(std::move(merged_fvar));
+  }
+
   return true;
 }
 
@@ -145,6 +250,25 @@ void TopologyRefinerFactory<OsdMesh>::reportInvalidTopology(TopologyError /*err_
 
 CCL_NAMESPACE_BEGIN
 
+/* OsdMesh */
+
+bool OsdMesh::use_smooth_fvar(const Attribute &attr) const
+{
+  return attr.element == ATTR_ELEMENT_CORNER &&
+         (attr.std == ATTR_STD_UV || (attr.flags & ATTR_SUBDIVIDE_SMOOTH_FVAR));
+}
+
+bool OsdMesh::use_smooth_fvar() const
+{
+  for (const Attribute &attr : mesh.subd_attributes.attributes) {
+    if (use_smooth_fvar(attr)) {
+      return true;
+    }
+  }
+
+  return false;
+}
+
 /* OsdData */
 
 void OsdData::build(OsdMesh &osd_mesh)
@@ -153,18 +277,27 @@ void OsdData::build(OsdMesh &osd_mesh)
 
   Sdc::Options options;
   options.SetVtxBoundaryInterpolation(Sdc::Options::VTX_BOUNDARY_EDGE_ONLY);
+  options.SetFVarLinearInterpolation(Sdc::Options::FVAR_LINEAR_CORNERS_PLUS1);
 
   /* create refiner */
   refiner.reset(Far::TopologyRefinerFactory<OsdMesh>::Create(
       osd_mesh, Far::TopologyRefinerFactory<OsdMesh>::Options(type, options)));
 
   /* adaptive refinement */
+  const bool has_fvar = osd_mesh.use_smooth_fvar();
   const int max_isolation = 3;  // TODO: get from Blender
-  refiner->RefineAdaptive(Far::TopologyRefiner::AdaptiveOptions(max_isolation));
+
+  Far::TopologyRefiner::AdaptiveOptions adaptive_options(max_isolation);
+  adaptive_options.considerFVarChannels = has_fvar;
+  adaptive_options.useInfSharpPatch = true;
+  refiner->RefineAdaptive(adaptive_options);
 
   /* create patch table */
   Far::PatchTableFactory::Options patch_options;
   patch_options.endCapType = Far::PatchTableFactory::Options::ENDCAP_GREGORY_BASIS;
+  patch_options.generateFVarTables = has_fvar;
+  patch_options.generateFVarLegacyLinearPatches = false;
+  patch_options.useInfSharpPatch = true;
 
   patch_table.reset(Far::PatchTableFactory::Create(*refiner, patch_options));
 

intern/cycles/subd/osd.h

@@ -48,9 +48,22 @@ template<typename T> struct OsdValue {
 
 class OsdMesh {
  public:
+  /* Face-varying attribute that requires merging of corners with the same value, typically a UV
+   * map. The resulting topology after merging is stored in a topology refiner fvar channel. The
+   * merged attribute values are stored here, in a generic buffer used for different data types. */
+  struct MergedFVar {
+    const Attribute &attr;
+    int channel = -1;
+    vector<char> values;
+  };
+
   Mesh &mesh;
+  vector<MergedFVar> merged_fvars;
 
   explicit OsdMesh(Mesh &mesh) : mesh(mesh) {}
+
+  bool use_smooth_fvar(const Attribute &attr) const;
+  bool use_smooth_fvar() const;
 };
 
 /* OpenSubdiv refiner and patch data structures. */