You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
141 lines
6.4 KiB
141 lines
6.4 KiB
// This file is part of libigl, a simple c++ geometry processing library.
|
|
//
|
|
// Copyright (C) 2021 Alec Jacobson <alecjacobson@gmail.com>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla Public License
|
|
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
|
|
// obtain one at http://mozilla.org/MPL/2.0/.
|
|
#include "marching_cubes.h"
|
|
#include "march_cube.h"
|
|
|
|
// Adapted from public domain code at
|
|
// http://paulbourke.net/geometry/polygonise/marchingsource.cpp
|
|
|
|
#include <unordered_map>
|
|
#include <iostream>
|
|
#include <cstdint>
|
|
|
|
template <typename DerivedS, typename DerivedGV, typename DerivedV, typename DerivedF>
|
|
IGL_INLINE void igl::marching_cubes(
|
|
const Eigen::MatrixBase<DerivedS> &S,
|
|
const Eigen::MatrixBase<DerivedGV> &GV,
|
|
const unsigned nx,
|
|
const unsigned ny,
|
|
const unsigned nz,
|
|
const typename DerivedS::Scalar isovalue,
|
|
Eigen::PlainObjectBase<DerivedV> &V,
|
|
Eigen::PlainObjectBase<DerivedF> &F)
|
|
{
|
|
typedef typename DerivedS::Scalar Scalar;
|
|
typedef unsigned Index;
|
|
// use same order as a2fVertexOffset
|
|
const unsigned ioffset[8] = {0,1,1+nx,nx,nx*ny,1+nx*ny,1+nx+nx*ny,nx+nx*ny};
|
|
|
|
|
|
std::unordered_map<std::int64_t,int> E2V;
|
|
V.resize(std::pow(nx*ny*nz,2./3.),3);
|
|
F.resize(std::pow(nx*ny*nz,2./3.),3);
|
|
Index n = 0;
|
|
Index m = 0;
|
|
|
|
const auto xyz2i = [&nx,&ny]
|
|
(const int & x, const int & y, const int & z)->unsigned
|
|
{
|
|
return x+nx*(y+ny*(z));
|
|
};
|
|
const auto cube =
|
|
[
|
|
&GV,&S,&V,&n,&F,&m,&isovalue,
|
|
&E2V,&xyz2i,&ioffset
|
|
]
|
|
(const int x, const int y, const int z)
|
|
{
|
|
const unsigned i = xyz2i(x,y,z);
|
|
|
|
//Make a local copy of the values at the cube's corners
|
|
Eigen::Matrix<Scalar,8,1> cS;
|
|
Eigen::Matrix<Index,8,1> cI;
|
|
//Find which vertices are inside of the surface and which are outside
|
|
for(int c = 0; c < 8; c++)
|
|
{
|
|
const unsigned ic = i + ioffset[c];
|
|
cI(c) = ic;
|
|
cS(c) = S(ic);
|
|
}
|
|
|
|
march_cube(GV,cS,cI,isovalue,V,n,F,m,E2V);
|
|
|
|
};
|
|
|
|
// march over all cubes (loop order chosen to match memory)
|
|
//
|
|
// Should be possible to parallelize safely if threads are "well separated".
|
|
// Like red-black Gauss Seidel. Probably each thread need's their own E2V,V,F,
|
|
// and then merge at the end. Annoying part are the edges lying on the
|
|
// interface between chunks.
|
|
for(int z=0;z<nz-1;z++)
|
|
{
|
|
for(int y=0;y<ny-1;y++)
|
|
{
|
|
for(int x=0;x<nx-1;x++)
|
|
{
|
|
cube(x,y,z);
|
|
}
|
|
}
|
|
}
|
|
V.conservativeResize(n,3);
|
|
F.conservativeResize(m,3);
|
|
}
|
|
|
|
template <
|
|
typename DerivedS,
|
|
typename DerivedGV,
|
|
typename DerivedGI,
|
|
typename DerivedV,
|
|
typename DerivedF>
|
|
IGL_INLINE void igl::marching_cubes(
|
|
const Eigen::MatrixBase<DerivedS> & S,
|
|
const Eigen::MatrixBase<DerivedGV> & GV,
|
|
const Eigen::MatrixBase<DerivedGI> & GI,
|
|
const typename DerivedS::Scalar isovalue,
|
|
Eigen::PlainObjectBase<DerivedV> &V,
|
|
Eigen::PlainObjectBase<DerivedF> &F)
|
|
{
|
|
typedef Eigen::Index Index;
|
|
typedef typename DerivedV::Scalar Scalar;
|
|
|
|
std::unordered_map<std::int64_t,int> E2V;
|
|
V.resize(4*GV.rows(),3);
|
|
F.resize(4*GV.rows(),3);
|
|
Index n = 0;
|
|
Index m = 0;
|
|
|
|
// march over cubes
|
|
|
|
//Make a local copy of the values at the cube's corners
|
|
Eigen::Matrix<Scalar, 8, 1> cS;
|
|
Eigen::Matrix<Index, 8, 1> cI;
|
|
for(Index c = 0;c<GI.rows();c++)
|
|
{
|
|
for(int v = 0; v < 8; v++)
|
|
{
|
|
cI(v) = GI(c,v);
|
|
cS(v) = S(GI(c,v));
|
|
}
|
|
march_cube(GV,cS,cI,isovalue,V,n,F,m,E2V);
|
|
}
|
|
V.conservativeResize(n,3);
|
|
F.conservativeResize(m,3);
|
|
}
|
|
|
|
#ifdef IGL_STATIC_LIBRARY
|
|
// Explicit template instantiation
|
|
// generated by autoexplicit.sh
|
|
template void igl::marching_cubes<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::Matrix<float, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
|
|
// generated by autoexplicit.sh
|
|
template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
|
|
// generated by autoexplicit.sh
|
|
template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, unsigned int, unsigned int, unsigned int, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
|
|
template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
|
|
template void igl::marching_cubes<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 8, 0, -1, 8>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 8, 0, -1, 8> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
|
|
#endif
|
|
|