Skip to content

Commit

Permalink
Add generalized winding number computation for meshToVolume conversion
Browse files Browse the repository at this point in the history
  • Loading branch information
@Alex23087
Alex23087 committed Mar 3, 2024
1 parent ff57414 commit 4e7029d
Show file tree
Hide file tree
Showing 4 changed files with 190 additions and 43 deletions.
5 changes: 5 additions & 0 deletions apps/openvdb/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,17 +11,22 @@ find_package(EIGEN3 REQUIRED NO_MODULE)

set(VCGLIB_DIR "../.." CACHE INTERNAL "Path to vcglib")
set(EIGEN_DIR = "${VCGLIB_DIR}/eigenlib" CACHE INTERNAL "")
set(WINDING_DIR ${VCGLIB_DIR}/../WindingNumber CACHE INTERNAL "")
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${ECM_FIND_MODULE_DIR})

file(GLOB WNFILES ${WINDING_DIR}/*.cpp)

add_executable(${PROJECT_NAME} OpenVDBRemeshing.cpp
${VCGLIB_DIR}/wrap/ply/plylib.cpp
${WNFILES}
)

target_include_directories(${PROJECT_NAME} PRIVATE ${VCGLIB_DIR})
target_include_directories(${PROJECT_NAME} PRIVATE ${TBB_INCLUDE_DIRS})
target_include_directories(${PROJECT_NAME} PRIVATE ${Boost_INCLUDE_DIRS})
target_include_directories(${PROJECT_NAME} PRIVATE ${OpenVDB_INCLUDE_DIR})
target_include_directories(${PROJECT_NAME} PRIVATE ${EIGEN_DIR})
target_include_directories(${PROJECT_NAME} PRIVATE ${WINDING_DIR})


target_link_libraries(${PROJECT_NAME} PRIVATE
Expand Down
32 changes: 20 additions & 12 deletions apps/openvdb/OpenVDBRemeshing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004-2012 \/)\/ *
* Copyright(C) 2004-2024 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
Expand Down Expand Up @@ -40,12 +40,12 @@ class MyFace;
class MyVertex;
struct MyUsedTypes : public UsedTypes<
Use<MyVertex> ::AsVertexType,
Use<MyFace> ::AsFaceType>{};
Use<MyFace> ::AsFaceType
>{};

class MyVertex : public Vertex<MyUsedTypes, vertex::Coord3f, vertex::BitFlags >{};
class MyFace : public Face<MyUsedTypes, face::VertexRef, vertex::BitFlags > {};
class MyMesh : public tri::TriMesh< vector<MyVertex>, vector<MyFace>
> {};
class MyVertex : public Vertex<MyUsedTypes, vertex::Coord3f, vertex::BitFlags>{};
class MyFace : public Face<MyUsedTypes, face::VertexRef, vertex::BitFlags> {};
class MyMesh : public tri::TriMesh< vector<MyVertex>, vector<MyFace>> {};

int main( int argc, char **argv )
{
Expand All @@ -64,14 +64,18 @@ int main( int argc, char **argv )
}

assert(original.VN()>0);

// OpenVDB remeshing parameters
double targetLenPerc=0.2d;
double isovalue=0.0d;
double adaptivity=0.0d;
double targetLenPerc = 0.2;
double isovalue = 0.0;
double adaptivity = 0.0;
bool useLevelSet = false; // Use level set instead of volume
if(argc>=3) targetLenPerc = atof(argv[2]);
if(argc>=4) isovalue = atof(argv[3]);
if(argc>=5) adaptivity = atof(argv[4]);
double voxelSize = targetLenPerc*(original.bbox.Diag()/100.0d);
if(argc>=6) useLevelSet = atoi(argv[5]);

double voxelSize = targetLenPerc * (original.bbox.Diag() / 100.0);

// Mesh cleaning
tri::Clean<MyMesh>::RemoveUnreferencedVertex(original);
Expand All @@ -83,9 +87,13 @@ int main( int argc, char **argv )

printf(" Input mesh %8i v %8i f\n",original.VN(),original.FN());

adapter.loadMesh(original);
adapter.setMesh(&original);
adapter.setVoxelSize(voxelSize);
adapter.meshToVolume();
if(useLevelSet){
adapter.meshToLevelSet();
}else{
adapter.meshToVolume();
}

// OpenVDB volume to mesh
adapter.setIsovalue(isovalue);
Expand Down
147 changes: 116 additions & 31 deletions wrap/openvdb/OpenVDBAdapter.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,36 @@
#ifndef OPENVDB_VCGLIB_ADAPTER_H
#define OPENVDB_VCGLIB_ADAPTER_H
#ifndef OPENVDB_VCGLIB_ADAPTER_CPP
#define OPENVDB_VCGLIB_ADAPTER_CPP

#include <vcg/complex/complex.h>
#include <vcg/complex/allocate.h>
#include <openvdb/openvdb.h>
#include <openvdb/tools/MeshToVolume.h>
#include <openvdb/tools/VolumeToMesh.h>
#include <openvdb/util/Util.h>

double isovalue = 0.01;
double adaptivity = 1;
double voxelSize = 0.005;
#include "WindingNumber.cpp"

namespace vcg{
namespace tri{

template<class TRI_MESH_TYPE>
class OpenVDBAdapter {
private:
typedef TRI_MESH_TYPE MeshType;
typedef typename MeshType::FaceType FaceType;
typedef typename FaceType::VertexType VertexType;
typedef typename VertexType::ScalarType ScalarType;
using GridType = openvdb::Grid<openvdb::v11_0::tree::Tree<openvdb::v11_0::tree::RootNode<openvdb::v11_0::tree::InternalNode<openvdb::v11_0::tree::InternalNode<openvdb::v11_0::tree::LeafNode<ScalarType, 3U>, 4U>, 5U>>>>;
class MeshTypeDataAdapter;


private:
std::vector<openvdb::v11_0::Vec3s> mVtx;
std::vector<openvdb::v11_0::math::Vec3<ScalarType>> mVtx;
std::vector<openvdb::v11_0::Vec3I> mTri;
std::vector<openvdb::v11_0::Vec4I> mQuad;
openvdb::FloatGrid::Ptr grid;
MeshType *m = nullptr;
typename GridType::Ptr grid;
WindingNumber<MeshType> windingNumber;
MeshTypeDataAdapter meshDataAdapter;

//MARK: Parameters
double isovalue = 0;
Expand All @@ -36,6 +40,7 @@ class OpenVDBAdapter {
public:
OpenVDBAdapter(){
openvdb::initialize();
meshDataAdapter = MeshTypeDataAdapter();
}

void setIsovalue(double isovalue){
Expand All @@ -50,63 +55,143 @@ class OpenVDBAdapter {
this->voxelSize = voxelSize;
}

void loadMesh(MeshType& m){
clearVectors();

for (size_t i = 0; i < m.vert.size(); i++){
mVtx.emplace_back(m.vert[i].P()[0], m.vert[i].P()[1], m.vert[i].P()[2]);
}
for (size_t i = 0; i < m.face.size(); i++){
mTri.emplace_back(Index(m, m.face[i].V(0)), Index(m, m.face[i].V(1)), Index(m, m.face[i].V(2)));
}
void setMesh(MeshType* m){
this->m = m;
}

void meshToVolume(){
assert(voxelSize > 0);

openvdb::math::Transform::Ptr xform = openvdb::math::Transform::createLinearTransform(voxelSize);

meshDataAdapter.setMesh(m);
meshDataAdapter.setTransform(xform);

windingNumber.init(*m);
auto interiorTest = [xform, &windingNumber = windingNumber](const openvdb::Coord &coord) -> bool
{
auto worldCoord = xform->indexToWorld(coord);
auto coordV = std::vector<ScalarType>{
static_cast<ScalarType>(worldCoord.x()),
static_cast<ScalarType>(worldCoord.y()),
static_cast<ScalarType>(worldCoord.z())
};
auto wn = windingNumber.computeWindingNumber(coordV);
return fabs(wn) >= 0.5 ? true : false;
};

float outerBand = (isovalue > 0 ? isovalue/voxelSize : 0) + 0.5f;
float innerBand = (isovalue < 0 ? (-isovalue)/voxelSize : 0) + 0.5f;

grid = openvdb::tools::meshToVolume<GridType>(meshDataAdapter, *xform, outerBand, innerBand, 0, nullptr, interiorTest, openvdb::tools::EVAL_EVERY_TILE);
}

void meshToLevelSet(){
copyMeshVectors();

assert(mVtx.size() > 0);
assert(mTri.size() > 0);
assert(voxelSize > 0);

openvdb::math::Transform::Ptr xform = openvdb::math::Transform::createLinearTransform(voxelSize);
grid = openvdb::tools::meshToLevelSet<openvdb::FloatGrid>(*xform, mVtx, mTri, 200.0f);

// The band half-width is dynamically computed as abs(isovalue/voxelSize) + 1.0f, which should be the smallest value to avoid artifacts.
grid = openvdb::tools::meshToLevelSet<GridType>(*xform, mVtx, mTri, std::fabs(isovalue/voxelSize) + 1.0f);
}

void volumeToMesh(MeshType& m){
void volumeToMesh(MeshType& outMesh){
assert(grid);

openvdb::tools::volumeToMesh(*grid, mVtx, mTri, mQuad, isovalue, adaptivity);

m.Clear();
auto vi = Allocator<MeshType>::AddVertices(m, mVtx.size());
auto fi = Allocator<MeshType>::AddFaces(m, mTri.size() + mQuad.size() * 2);
outMesh.Clear();
auto vi = Allocator<MeshType>::AddVertices(outMesh, mVtx.size());
auto fi = Allocator<MeshType>::AddFaces(outMesh, mTri.size() + mQuad.size() * 2);
for (size_t i = 0; i < mVtx.size(); vi++, i++){
vi->P()[0] = (ScalarType) mVtx[i][0];
vi->P()[1] = (ScalarType) mVtx[i][1];
vi->P()[2] = (ScalarType) mVtx[i][2];
}
for (size_t i = 0; i < mTri.size(); fi++, i++){
fi->V(0) = &m.vert[mTri[i][2]];
fi->V(1) = &m.vert[mTri[i][1]];
fi->V(2) = &m.vert[mTri[i][0]];
fi->V(0) = &outMesh.vert[mTri[i][2]];
fi->V(1) = &outMesh.vert[mTri[i][1]];
fi->V(2) = &outMesh.vert[mTri[i][0]];
}
for (size_t i = 0; i < mQuad.size(); i++){
fi->V(0) = &m.vert[mQuad[i][2]];
fi->V(1) = &m.vert[mQuad[i][1]];
fi->V(2) = &m.vert[mQuad[i][0]];
fi->V(0) = &outMesh.vert[mQuad[i][2]];
fi->V(1) = &outMesh.vert[mQuad[i][1]];
fi->V(2) = &outMesh.vert[mQuad[i][0]];
fi++;
fi->V(0) = &m.vert[mQuad[i][3]];
fi->V(1) = &m.vert[mQuad[i][2]];
fi->V(2) = &m.vert[mQuad[i][0]];
fi->V(0) = &outMesh.vert[mQuad[i][3]];
fi->V(1) = &outMesh.vert[mQuad[i][2]];
fi->V(2) = &outMesh.vert[mQuad[i][0]];
fi++;
}
}

private:
void copyMeshVectors(){
clearVectors();

for (size_t i = 0; i < m->vert.size(); i++){
mVtx.emplace_back(m->vert[i].P()[0], m->vert[i].P()[1], m->vert[i].P()[2]);
}
for (size_t i = 0; i < m->face.size(); i++){
mTri.emplace_back(Index(*m, m->face[i].V(0)), Index(*m, m->face[i].V(1)), Index(*m, m->face[i].V(2)));
}
}

void clearVectors(){
mVtx.clear();
mTri.clear();
mQuad.clear();
}
};



template <class TRI_MESH_TYPE>
class OpenVDBAdapter<TRI_MESH_TYPE>::MeshTypeDataAdapter{
private:
MeshType *m = nullptr;
openvdb::math::Transform::Ptr xform;

public:
MeshTypeDataAdapter(){}

MeshType* getMesh(){
return m;
}

void setMesh(MeshType* mesh){
m = mesh;
}

void setTransform(openvdb::math::Transform::Ptr xform){
this->xform = xform;
}

// Total number of polygons
size_t polygonCount() const {
return m->face.size();
}

// Total number of points
size_t pointCount() const {
return m->vert.size();
}

// Vertex count for polygon n
size_t vertexCount(size_t n) const {
return 3;
}

// Return position pos in local grid index space for polygon n and vertex v
void getIndexSpacePoint(size_t n, size_t v, openvdb::Vec3d& pos) const {
pos = openvdb::Vec3d(m->face[n].V(v)->P()[0], m->face[n].V(v)->P()[1], m->face[n].V(v)->P()[2]);
pos = xform->worldToIndex(pos);
}
};
}
}

Expand Down
49 changes: 49 additions & 0 deletions wrap/openvdb/WindingNumber.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
#ifndef WINDING_NUMBER_CPP
#define WINDING_NUMBER_CPP

#include <UT_SolidAngle.h>
#include <vcg/complex/complex.h>

template<class TRI_MESH_TYPE>
class WindingNumber {
typedef TRI_MESH_TYPE MeshType;
typedef typename MeshType::FaceType FaceType;
typedef typename FaceType::VertexType VertexType;
typedef typename VertexType::ScalarType ScalarType;

private:
HDK_Sample::UT_SolidAngle<float,float> solid_angle; // Only works with float

public:
void init(MeshType &m, int order = 2)
{
// Initialize vector of vertex positions
std::vector<HDK_Sample::UT_Vector3T<float>> U(m.vert.size());
for(int i = 0; i < m.vert.size(); i++)
{
for(int j = 0;j<3;j++)
{
U[i][j] = m.vert[i].P()[j];
}
}

// Initialize vector of triangle indices
std::vector<int> mTriIndex(m.face.size() * 3);
for(int i = 0; i < m.face.size(); i++)
{
for(int j = 0;j<3;j++)
{
mTriIndex[i*3+j] = vcg::tri::Index(m, m.face[i].V(j));
}
}

solid_angle.init(m.face.size(), mTriIndex.data(), m.vert.size(), &U[0], order);
}

ScalarType computeWindingNumber(std::vector<ScalarType> coordV, double accuracy_scale = 2.0){
auto pt = HDK_Sample::UT_Vector3T<float>(coordV.data());
return solid_angle.computeSolidAngle(pt, accuracy_scale) / (4.0 * M_PI);
}
};

#endif

0 comments on commit 4e7029d

Please sign in to comment.