easy-iso

Iso-surface extraction from volume data
git clone https://git.0xfab.ch/easy-iso.git
Log | Files | Refs | Submodules | README | LICENSE

commit 52ff767858d5e8a9769bb10a339184fb6dcc8e29
parent f4fba5b8f9d1569d1ab2f952310a1d98d222445a
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Wed, 30 Nov 2016 08:17:56 +0100

added template parameter for interpolator type

Diffstat:
Msrc/IsoExtractor.h | 111+++++++++++++++++++++++++++++++++++++++++++------------------------------------
1 file changed, 61 insertions(+), 50 deletions(-)

diff --git a/src/IsoExtractor.h b/src/IsoExtractor.h @@ -364,7 +364,7 @@ static const int gTetrahedronTriangles[16][7] = }; -template <typename Tinterp> +template <typename Tkernel, template<typename> class Tinterpolator=Interpolator> class IsoExtractor { public: @@ -378,10 +378,6 @@ public: void extract(const Real isoval, const std::string fout) { - const Real x0 = m_parser("xmin").asDouble(0.0); - const Real y0 = m_parser("ymin").asDouble(0.0); - const Real z0 = m_parser("zmin").asDouble(0.0); - const Real cscale = m_parser("cube_scale").asDouble(1.0); assert(cscale > 0.0); const int Nx = static_cast<Real>(m_interp.getNx()) / cscale; @@ -392,45 +388,15 @@ public: assert(Nz > 1); const Real ch = m_interp.getExtent() / (std::max(Nx, std::max(Ny, Nz)) - m_interp.isNodal()); - m_triList.clear(); - - if (m_parser("march").asString("cube") == "cube") - { - for (int k = 0; k < Nz; ++k) - for (int j = 0; j < Ny; ++j) - for (int i = 0; i < Nx; ++i) - { - const Real x = i*ch; - const Real y = j*ch; - const Real z = k*ch; - _marchCube(x, y, z, ch, isoval); - } - } - else if (m_parser("march").asString("cube") == "tetrahedron") - { - for (int k = 0; k < Nz; ++k) - for (int j = 0; j < Ny; ++j) - for (int i = 0; i < Nx; ++i) - { - const Real x = i*ch; - const Real y = j*ch; - const Real z = k*ch; - _marchTetrahedron(x, y, z, ch, isoval); - } - } - else - { - std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << m_parser("march").asString() << std::endl; - abort(); - } + _extractIso(isoval, Nx, Ny, Nz, ch); _writePLY(fout, m_triList); } -private: +protected: ArgumentParser& m_parser; - Interpolator<Tinterp> m_interp; + Tinterpolator<Tkernel> m_interp; const Real m_x0; const Real m_y0; const Real m_z0; @@ -445,8 +411,9 @@ private: return (t - v0)/dv; } - inline void _marchCube(const Real x, const Real y, const Real z, const Real h, const Real target); - inline void _marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target); + inline void _extractIso(const Real isoval, const int Nx, const int Ny, const int Nz, const Real ch, const Real Ox=0.0, const Real Oy=0.0, const Real Oz=0.0, const bool verbose=true); + inline void _marchCube(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz); + inline void _marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz); inline void _evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target); void _writePLY(const std::string basename, const std::vector<Triangle>& tList) @@ -488,7 +455,7 @@ private: } else { - ply.open(basename+".ply", std::ios::out | std::ios::app | std::ios::binary); + ply.open(basename+".ply", std::ios::out | std::ios::app); _writePLY_ascii(ply, tList); } ply.close(); @@ -528,8 +495,47 @@ private: }; -template <typename Tinterp> -void IsoExtractor<Tinterp>::_marchCube(const Real x, const Real y, const Real z, const Real h, const Real target) +template <typename Tkernel, template<typename> class Tinterpolator> +void IsoExtractor<Tkernel,Tinterpolator>::_extractIso(const Real isoval, const int Nx, const int Ny, const int Nz, const Real ch, const Real Ox, const Real Oy, const Real Oz, const bool verbose) +{ + const int speculate = m_parser("speculate_triangle").asInt(10000); + m_triList.clear(); + m_triList.reserve(speculate); + if (m_parser("march").asString("cube") == "cube") + { + for (int k = 0; k < Nz; ++k) + for (int j = 0; j < Ny; ++j) + for (int i = 0; i < Nx; ++i) + { + const Real x = i*ch; + const Real y = j*ch; + const Real z = k*ch; + _marchCube(x, y, z, ch, isoval, Ox, Oy, Oz); + } + } + else if (m_parser("march").asString("cube") == "tetrahedron") + { + for (int k = 0; k < Nz; ++k) + for (int j = 0; j < Ny; ++j) + for (int i = 0; i < Nx; ++i) + { + const Real x = i*ch; + const Real y = j*ch; + const Real z = k*ch; + _marchTetrahedron(x, y, z, ch, isoval, Ox, Oy, Oz); + } + } + else + { + if (verbose) + std::cerr << "ERROR: IsoExtractor: Unknown algorithm " << m_parser("march").asString() << std::endl; + abort(); + } +} + + +template <typename Tkernel, template<typename> class Tinterpolator> +void IsoExtractor<Tkernel,Tinterpolator>::_marchCube(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz) { Real valCubeVerts[8]; Real coordEdgeVert[12][3]; @@ -562,9 +568,9 @@ void IsoExtractor<Tinterp>::_marchCube(const Real x, const Real y, const Real z, { // linear interpolation const Real frac = _getEdgeFrac(valCubeVerts[gEdgeConnection[i][0]], valCubeVerts[gEdgeConnection[i][1]], target); - coordEdgeVert[i][0] = x + h*(gVertexOffset[gEdgeConnection[i][0]][0] + frac*gEdgeDirection[i][0]); - coordEdgeVert[i][1] = y + h*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]); - coordEdgeVert[i][2] = z + h*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]); + coordEdgeVert[i][0] = x + h*(gVertexOffset[gEdgeConnection[i][0]][0] + frac*gEdgeDirection[i][0]) + m_x0 + Ox; + coordEdgeVert[i][1] = y + h*(gVertexOffset[gEdgeConnection[i][0]][1] + frac*gEdgeDirection[i][1]) + m_y0 + Oy; + coordEdgeVert[i][2] = z + h*(gVertexOffset[gEdgeConnection[i][0]][2] + frac*gEdgeDirection[i][2]) + m_z0 + Oz; } } @@ -586,8 +592,8 @@ void IsoExtractor<Tinterp>::_marchCube(const Real x, const Real y, const Real z, } -template <typename Tinterp> -void IsoExtractor<Tinterp>::_marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target) +template <typename Tkernel, template<typename> class Tinterpolator> +void IsoExtractor<Tkernel,Tinterpolator>::_marchTetrahedron(const Real x, const Real y, const Real z, const Real h, const Real target, const Real Ox, const Real Oy, const Real Oz) { Real valCubeVerts[8]; Real coordCubeVerts[8][3]; @@ -604,10 +610,15 @@ void IsoExtractor<Tinterp>::_marchTetrahedron(const Real x, const Real y, const // compute cube vertex values for (int i = 0; i < 8; ++i) + { valCubeVerts[i] = m_interp( coordCubeVerts[i][0], coordCubeVerts[i][1], coordCubeVerts[i][2]); + coordCubeVerts[i][0] += m_x0 + Ox; + coordCubeVerts[i][1] += m_y0 + Oy; + coordCubeVerts[i][2] += m_z0 + Oz; + } for (int i = 0; i < 6; ++i) { @@ -624,8 +635,8 @@ void IsoExtractor<Tinterp>::_marchTetrahedron(const Real x, const Real y, const } -template <typename Tinterp> -void IsoExtractor<Tinterp>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target) +template <typename Tkernel, template<typename> class Tinterpolator> +void IsoExtractor<Tkernel,Tinterpolator>::_evalTetrahedron(const Real coords[4][3], const Real val[4], const Real target) { Real coordEdgeVert[6][3];