findFeaturePoint for sharp-edge extraction?

[phaedon]

Hi all,
I am looking at the findFeaturePoint function which references the Kobbelt (2001) paper on the topic:
https://www.openvdb.org/documentation/doxygen/VolumeToMesh_8h.html

Am not quite sure how to incorporate this into my project. Roughly speaking, I am using meshToVolume to generate a voxel grid, doing some custom processing, and then volumeToMesh to reconstruct the surface.

I am uncertain both for:
(a) How to use this API. For example: Do I need to iterate over each grid cell and use the original input triangle mesh to manually search for all the triangles which intersect the grid cell, extract the barycenter (for example, which might be outside the cell region) and the normal, and pass them all into the function? If so is there already a built-in function for identifying triangle-voxel intersections?

(b) Once I have the feature points for each relevant cell, where do I add these extra samples? Do I put them in a second FloatGrid? And do I need to implement the extended Marching Cubes from scratch, or is this already provided somewhere in the library?

I'd appreciate any guidance and especially any example code that could help clarify this. I couldn't find any unit tests for findFeaturePoint, but even if there were, it wouldn't necessarily answer the second question.

[phaedon]

Updating my question here:

What I've figured out so far
Looking at some of the code in this example:
https://github.com/AcademySoftwareFoundation/openvdb/blob/master/openvdb_houdini/openvdb_houdini/GeometryUtil.cc
I was able to make a bit of progress. It appears that:

• MeshToVoxelEdgeData is the key to getting intersection points and primitives, and most critically, the input point_list when initialising this object with the convert() method must be in index space.
• Always work in index space; the output of openvdb::tools::findFeaturePoint is also in index space
• [I think?] don't bother calling findFeaturePoint if all of the normals are pointing in the same direction (but I'm not sure why the Houdini library doesn't seem to do this check). Is my intuition correct that I should be filtering out all of these superfluous "feature points", which are already captured by the isosurface in the signed grid?

My remaining questions:

• I am still missing all 8 corners, and the edges incident to the upper right of my simple cube test model. Am wondering if I need to make some constant adjustment for the notion of voxels-as-point-samples versus voxels-as-tiny-bounding-box-cubes (in Houdini, it looks like they are represented as unit voxel cubes with an extra padding of 0.3).
• I can't get the ray tracing to work (i.e. the branch that says if (!cell.isInside(feature_point))) -- and perhaps this is why my corners are missing, because those feature points are necessarily outside the voxel of interest?
• Once the above two issues are fixed, I'm still not sure where to put these feature points on my grid as per my original question.