kannonier8/blender
1
0
Fork 0
forked from bartvdbraak/blender
Code Pull Requests Activity
12c08ceee3
blender/extern/mantaflow/preprocessed/levelset.cpp
Sebastián Barschkis 311031ecd0 Cleanup: Use nullptr everywhere in fluid code 
Switched from NULL to nullptr.
2020-11-06 12:06:05 +01:00

1076 lines
30 KiB
C++
Raw Blame History

// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).
/******************************************************************************
*
* MantaFlow fluid solver framework
* Copyright 2011 Tobias Pfaff, Nils Thuerey
*
* This program is free software, distributed under the terms of the
* Apache License, Version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Levelset
*
******************************************************************************/
#include "levelset.h"
#include "fastmarch.h"
#include "kernel.h"
#include "mcubes.h"
#include "mesh.h"
#include <stack>
using namespace std;
namespace Manta {
//************************************************************************
// Helper functions and kernels for marching
static const int FlagInited = FastMarch<FmHeapEntryOut, +1>::FlagInited;
// neighbor lookup vectors
static const Vec3i neighbors[6] = {Vec3i(-1, 0, 0),
Vec3i(1, 0, 0),
Vec3i(0, -1, 0),
Vec3i(0, 1, 0),
Vec3i(0, 0, -1),
Vec3i(0, 0, 1)};
struct InitFmIn : public KernelBase {
InitFmIn(const FlagGrid &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
bool ignoreWalls,
int obstacleType)
: KernelBase(&flags, 1),
flags(flags),
fmFlags(fmFlags),
phi(phi),
ignoreWalls(ignoreWalls),
obstacleType(obstacleType)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
bool ignoreWalls,
int obstacleType) const
{
const IndexInt idx = flags.index(i, j, k);
const Real v = phi[idx];
if (ignoreWalls) {
if (v >= 0. && ((flags[idx] & obstacleType) == 0))
fmFlags[idx] = FlagInited;
else
fmFlags[idx] = 0;
}
else {
if (v >= 0)
fmFlags[idx] = FlagInited;
else
fmFlags[idx] = 0;
}
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline Grid<int> &getArg1()
{
return fmFlags;
}
typedef Grid<int> type1;
inline Grid<Real> &getArg2()
{
return phi;
}
typedef Grid<Real> type2;
inline bool &getArg3()
{
return ignoreWalls;
}
typedef bool type3;
inline int &getArg4()
{
return obstacleType;
}
typedef int type4;
void runMessage()
{
debMsg("Executing kernel InitFmIn ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
const int _maxX = maxX;
const int _maxY = maxY;
if (maxZ > 1) {
for (int k = __r.begin(); k != (int)__r.end(); k++)
for (int j = 1; j < _maxY; j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
}
else {
const int k = 0;
for (int j = __r.begin(); j != (int)__r.end(); j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
}
}
void run()
{
if (maxZ > 1)
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
else
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
}
const FlagGrid &flags;
Grid<int> &fmFlags;
Grid<Real> &phi;
bool ignoreWalls;
int obstacleType;
};
struct InitFmOut : public KernelBase {
InitFmOut(const FlagGrid &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
bool ignoreWalls,
int obstacleType)
: KernelBase(&flags, 1),
flags(flags),
fmFlags(fmFlags),
phi(phi),
ignoreWalls(ignoreWalls),
obstacleType(obstacleType)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
bool ignoreWalls,
int obstacleType) const
{
const IndexInt idx = flags.index(i, j, k);
const Real v = phi[idx];
if (ignoreWalls) {
fmFlags[idx] = (v < 0) ? FlagInited : 0;
if ((flags[idx] & obstacleType) != 0) {
fmFlags[idx] = 0;
phi[idx] = 0;
}
}
else {
fmFlags[idx] = (v < 0) ? FlagInited : 0;
}
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline Grid<int> &getArg1()
{
return fmFlags;
}
typedef Grid<int> type1;
inline Grid<Real> &getArg2()
{
return phi;
}
typedef Grid<Real> type2;
inline bool &getArg3()
{
return ignoreWalls;
}
typedef bool type3;
inline int &getArg4()
{
return obstacleType;
}
typedef int type4;
void runMessage()
{
debMsg("Executing kernel InitFmOut ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
const int _maxX = maxX;
const int _maxY = maxY;
if (maxZ > 1) {
for (int k = __r.begin(); k != (int)__r.end(); k++)
for (int j = 1; j < _maxY; j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
}
else {
const int k = 0;
for (int j = __r.begin(); j != (int)__r.end(); j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
}
}
void run()
{
if (maxZ > 1)
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
else
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
}
const FlagGrid &flags;
Grid<int> &fmFlags;
Grid<Real> &phi;
bool ignoreWalls;
int obstacleType;
};
struct SetUninitialized : public KernelBase {
SetUninitialized(const Grid<int> &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
const Real val,
int ignoreWalls,
int obstacleType)
: KernelBase(&flags, 1),
flags(flags),
fmFlags(fmFlags),
phi(phi),
val(val),
ignoreWalls(ignoreWalls),
obstacleType(obstacleType)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const Grid<int> &flags,
Grid<int> &fmFlags,
Grid<Real> &phi,
const Real val,
int ignoreWalls,
int obstacleType) const
{
if (ignoreWalls) {
if ((fmFlags(i, j, k) != FlagInited) && ((flags(i, j, k) & obstacleType) == 0)) {
phi(i, j, k) = val;
}
}
else {
if ((fmFlags(i, j, k) != FlagInited))
phi(i, j, k) = val;
}
}
inline const Grid<int> &getArg0()
{
return flags;
}
typedef Grid<int> type0;
inline Grid<int> &getArg1()
{
return fmFlags;
}
typedef Grid<int> type1;
inline Grid<Real> &getArg2()
{
return phi;
}
typedef Grid<Real> type2;
inline const Real &getArg3()
{
return val;
}
typedef Real type3;
inline int &getArg4()
{
return ignoreWalls;
}
typedef int type4;
inline int &getArg5()
{
return obstacleType;
}
typedef int type5;
void runMessage()
{
debMsg("Executing kernel SetUninitialized ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
const int _maxX = maxX;
const int _maxY = maxY;
if (maxZ > 1) {
for (int k = __r.begin(); k != (int)__r.end(); k++)
for (int j = 1; j < _maxY; j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, val, ignoreWalls, obstacleType);
}
else {
const int k = 0;
for (int j = __r.begin(); j != (int)__r.end(); j++)
for (int i = 1; i < _maxX; i++)
op(i, j, k, flags, fmFlags, phi, val, ignoreWalls, obstacleType);
}
}
void run()
{
if (maxZ > 1)
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
else
tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
}
const Grid<int> &flags;
Grid<int> &fmFlags;
Grid<Real> &phi;
const Real val;
int ignoreWalls;
int obstacleType;
};
template<bool inward>
inline bool isAtInterface(const Grid<int> &fmFlags, Grid<Real> &phi, const Vec3i &p)
{
// check for interface
int max = phi.is3D() ? 6 : 4;
for (int nb = 0; nb < max; nb++) {
const Vec3i pn(p + neighbors[nb]);
if (!fmFlags.isInBounds(pn))
continue;
if (fmFlags(pn) != FlagInited)
continue;
if ((inward && phi(pn) >= 0.) || (!inward && phi(pn) < 0.))
return true;
}
return false;
}
//************************************************************************
// Levelset class def
LevelsetGrid::LevelsetGrid(FluidSolver *parent, bool show) : Grid<Real>(parent, show)
{
mType = (GridType)(TypeLevelset | TypeReal);
}
LevelsetGrid::LevelsetGrid(FluidSolver *parent, Real *data, bool show)
: Grid<Real>(parent, data, show)
{
mType = (GridType)(TypeLevelset | TypeReal);
}
Real LevelsetGrid::invalidTimeValue()
{
return FastMarch<FmHeapEntryOut, 1>::InvalidTime();
}
//! Kernel: perform levelset union
struct KnJoin : public KernelBase {
KnJoin(Grid<Real> &a, const Grid<Real> &b) : KernelBase(&a, 0), a(a), b(b)
{
runMessage();
run();
}
inline void op(IndexInt idx, Grid<Real> &a, const Grid<Real> &b) const
{
a[idx] = min(a[idx], b[idx]);
}
inline Grid<Real> &getArg0()
{
return a;
}
typedef Grid<Real> type0;
inline const Grid<Real> &getArg1()
{
return b;
}
typedef Grid<Real> type1;
void runMessage()
{
debMsg("Executing kernel KnJoin ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
op(idx, a, b);
}
void run()
{
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
}
Grid<Real> &a;
const Grid<Real> &b;
};
void LevelsetGrid::join(const LevelsetGrid &o)
{
KnJoin(*this, o);
}
//! subtract b, note does not preserve SDF!
struct KnSubtract : public KernelBase {
KnSubtract(Grid<Real> &a, const Grid<Real> &b, const FlagGrid *flags, int subtractType)
: KernelBase(&a, 0), a(a), b(b), flags(flags), subtractType(subtractType)
{
runMessage();
run();
}
inline void op(IndexInt idx,
Grid<Real> &a,
const Grid<Real> &b,
const FlagGrid *flags,
int subtractType) const
{
if (flags && ((*flags)(idx)&subtractType) == 0)
return;
if (b[idx] < 0.)
a[idx] = b[idx] * -1.;
}
inline Grid<Real> &getArg0()
{
return a;
}
typedef Grid<Real> type0;
inline const Grid<Real> &getArg1()
{
return b;
}
typedef Grid<Real> type1;
inline const FlagGrid *getArg2()
{
return flags;
}
typedef FlagGrid type2;
inline int &getArg3()
{
return subtractType;
}
typedef int type3;
void runMessage()
{
debMsg("Executing kernel KnSubtract ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
op(idx, a, b, flags, subtractType);
}
void run()
{
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
}
Grid<Real> &a;
const Grid<Real> &b;
const FlagGrid *flags;
int subtractType;
};
void LevelsetGrid::subtract(const LevelsetGrid &o, const FlagGrid *flags, const int subtractType)
{
KnSubtract(*this, o, flags, subtractType);
}
//! re-init levelset and extrapolate velocities (in & out)
// note - uses flags to identify border (could also be done based on ls values)
static void doReinitMarch(Grid<Real> &phi,
const FlagGrid &flags,
Real maxTime,
MACGrid *velTransport,
bool ignoreWalls,
bool correctOuterLayer,
int obstacleType)
{
const int dim = (phi.is3D() ? 3 : 2);
Grid<int> fmFlags(phi.getParent());
FastMarch<FmHeapEntryIn, -1> marchIn(flags, fmFlags, phi, maxTime, nullptr);
// march inside
InitFmIn(flags, fmFlags, phi, ignoreWalls, obstacleType);
FOR_IJK_BND(flags, 1)
{
if (fmFlags(i, j, k) == FlagInited)
continue;
if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
continue;
const Vec3i p(i, j, k);
if (isAtInterface<true>(fmFlags, phi, p)) {
// set value
fmFlags(p) = FlagInited;
// add neighbors that are not at the interface
for (int nb = 0; nb < 2 * dim; nb++) {
const Vec3i pn(p + neighbors[nb]); // index always valid due to bnd=1
if (ignoreWalls && ((flags.get(pn) & obstacleType) != 0))
continue;
// check neighbors of neighbor
if (phi(pn) < 0. && !isAtInterface<true>(fmFlags, phi, pn)) {
marchIn.addToList(pn, p);
}
}
}
}
marchIn.performMarching();
// done with inwards marching
// now march out...
// set un initialized regions
SetUninitialized(flags, fmFlags, phi, -maxTime - 1., ignoreWalls, obstacleType);
InitFmOut(flags, fmFlags, phi, ignoreWalls, obstacleType);
FastMarch<FmHeapEntryOut, +1> marchOut(flags, fmFlags, phi, maxTime, velTransport);
// by default, correctOuterLayer is on
if (correctOuterLayer) {
// normal version, inwards march is done, now add all outside values (0..2] to list
// note, this might move the interface a bit! but keeps a nice signed distance field...
FOR_IJK_BND(flags, 1)
{
if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
continue;
const Vec3i p(i, j, k);
// check nbs
for (int nb = 0; nb < 2 * dim; nb++) {
const Vec3i pn(p + neighbors[nb]); // index always valid due to bnd=1
if (fmFlags(pn) != FlagInited)
continue;
if (ignoreWalls && ((flags.get(pn) & obstacleType)) != 0)
continue;
const Real nbPhi = phi(pn);
// only add nodes near interface, not e.g. outer boundary vs. invalid region
if (nbPhi < 0 && nbPhi >= -2)
marchOut.addToList(p, pn);
}
}
}
else {
// alternative version, keep interface, do not distort outer cells
// add all ouside values, but not those at the IF layer
FOR_IJK_BND(flags, 1)
{
if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
continue;
// only look at ouside values
const Vec3i p(i, j, k);
if (phi(p) < 0)
continue;
if (isAtInterface<false>(fmFlags, phi, p)) {
// now add all non, interface neighbors
fmFlags(p) = FlagInited;
// add neighbors that are not at the interface
for (int nb = 0; nb < 2 * dim; nb++) {
const Vec3i pn(p + neighbors[nb]); // index always valid due to bnd=1
if (ignoreWalls && ((flags.get(pn) & obstacleType) != 0))
continue;
// check neighbors of neighbor
if (phi(pn) > 0. && !isAtInterface<false>(fmFlags, phi, pn)) {
marchOut.addToList(pn, p);
}
}
}
}
}
marchOut.performMarching();
// set un initialized regions
SetUninitialized(flags, fmFlags, phi, +maxTime + 1., ignoreWalls, obstacleType);
}
//! call for levelset grids & external real grids
void LevelsetGrid::reinitMarching(const FlagGrid &flags,
Real maxTime,
MACGrid *velTransport,
bool ignoreWalls,
bool correctOuterLayer,
int obstacleType)
{
doReinitMarch(*this, flags, maxTime, velTransport, ignoreWalls, correctOuterLayer, obstacleType);
}
void LevelsetGrid::initFromFlags(const FlagGrid &flags, bool ignoreWalls)
{
FOR_IDX(*this)
{
if (flags.isFluid(idx) || (ignoreWalls && flags.isObstacle(idx)))
mData[idx] = -0.5;
else
mData[idx] = 0.5;
}
}
/* Helper variables that are used in flood-fill functions. */
static const int ID_UNKNOWN = 0;
static const int ID_VISITED = 1;
/* Fills all cells in the target grid that have not been marked during a flood-fill. */
struct KnFillApply : public KernelBase {
KnFillApply(Grid<Real> &target,
Grid<int> &visited,
const Real value,
const int boundaryWidth,
const bool outside)
: KernelBase(&target, boundaryWidth),
target(target),
visited(visited),
value(value),
boundaryWidth(boundaryWidth),
outside(outside)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
Grid<Real> &target,
Grid<int> &visited,
const Real value,
const int boundaryWidth,
const bool outside) const
{
if (visited(i, j, k) == ID_VISITED)
return;
if (outside && target(i, j, k) < 0)
return;
if (!outside && target(i, j, k) >= 0)
return;
/* Actual flood-fill override. */
target(i, j, k) = value;
}
inline Grid<Real> &getArg0()
{
return target;
}
typedef Grid<Real> type0;
inline Grid<int> &getArg1()
{
return visited;
}
typedef Grid<int> type1;
inline const Real &getArg2()
{
return value;
}
typedef Real type2;
inline const int &getArg3()
{
return boundaryWidth;
}
typedef int type3;
inline const bool &getArg4()
{
return outside;
}
typedef bool type4;
void runMessage()
{
debMsg("Executing kernel KnFillApply ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
const int _maxX = maxX;
const int _maxY = maxY;
if (maxZ > 1) {
for (int k = __r.begin(); k != (int)__r.end(); k++)
for (int j = boundaryWidth; j < _maxY; j++)
for (int i = boundaryWidth; i < _maxX; i++)
op(i, j, k, target, visited, value, boundaryWidth, outside);
}
else {
const int k = 0;
for (int j = __r.begin(); j != (int)__r.end(); j++)
for (int i = boundaryWidth; i < _maxX; i++)
op(i, j, k, target, visited, value, boundaryWidth, outside);
}
}
void run()
{
if (maxZ > 1)
tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
else
tbb::parallel_for(tbb::blocked_range<IndexInt>(boundaryWidth, maxY), *this);
}
Grid<Real> &target;
Grid<int> &visited;
const Real value;
const int boundaryWidth;
const bool outside;
};
/* Basic flood fill implementation used to fill inside / outside areas of levelset.
* Calling this function will ensure that there are no fluid cells inside obstacles.
* I.e. starting from walls, cells will be tagged in flood-fill fashion, stopping at 0 borders.
* All remaining cells will be filled with the fill value. Outside mode inverts search behavior. */
void LevelsetGrid::floodFill(const Real value, const bool outside, const int boundaryWidth)
{
/* Sanity check: Filling mode and filling value need to "match". */
if (outside) {
assertMsg(value < 0, "Cannot fill outside with (positive) value " << value);
}
else {
assertMsg(value >= 0, "Cannot fill inside with (negative) value " << value);
}
Grid<Real> levelsetCopy(this->getParent());
Grid<int> visited(this->getParent());
std::stack<Vec3i> todoPos;
const int maxNeighbors = this->is3D() ? 6 : 4;
const Vec3i maxSize(this->getSize() - 1);
Vec3i bnd(2 * boundaryWidth);
if (!this->is3D())
bnd.z = 0;
const int cellCntNoBnd = (this->getSizeX() - bnd.x) * (this->getSizeY() - bnd.y) *
(this->getSizeZ() - bnd.z);
/* Initialize temporary helper grids. */
levelsetCopy.copyFrom(*this);
visited.setConst(ID_UNKNOWN);
FOR_IJK_BND(visited, boundaryWidth)
{
/* Skip inside / outside cells depending on search mode. */
if (outside && levelsetCopy(i, j, k) < 0)
continue;
if (!outside && levelsetCopy(i, j, k) >= 0)
continue;
/* Skip cell if it already has been visited. */
if (visited(i, j, k) == ID_VISITED)
continue;
Vec3i c(i, j, k);
bool isWallCell = (c.x - boundaryWidth == 0 || c.x == maxSize.x - boundaryWidth);
isWallCell |= (c.y - boundaryWidth == 0 || c.y == maxSize.y - boundaryWidth);
if (this->is3D())
isWallCell |= (c.z - boundaryWidth == 0 || c.z == maxSize.z - boundaryWidth);
/* Only start searching from borders. */
if (!isWallCell)
continue;
/* Start flood-fill loop by initializing todo stack with current cell. */
todoPos.push(c);
visited(c) = ID_VISITED;
while (!todoPos.empty()) {
c = todoPos.top();
todoPos.pop();
/* Add all neighbor cells to search stack. */
for (int nb = 0; nb < maxNeighbors; nb++) {
const Vec3i neigh(c + neighbors[nb]);
if (!visited.isInBounds(neigh, boundaryWidth))
continue;
/* Skip inside / outside area depening on what we search for. */
if (outside && levelsetCopy(neigh) < 0)
continue;
if (!outside && levelsetCopy(neigh) >= 0)
continue;
/* Skip neighbor if it already has been visited. */
if (visited(neigh) == ID_VISITED)
continue;
assertMsg(visited(neigh) == ID_UNKNOWN,
"Cell must be of type 'unknown' at this point in the loop");
todoPos.push(neigh);
visited(neigh) = ID_VISITED;
}
assertMsg(todoPos.size() <= cellCntNoBnd,
"Flood-fill todo stack cannot be greater than domain cell count - "
<< todoPos.size() << " vs " << cellCntNoBnd);
}
}
KnFillApply(*this, visited, value, boundaryWidth, outside);
}
/* Deprecated: Use floodFill() function instead. */
void LevelsetGrid::fillHoles(int maxDepth, int boundaryWidth)
{
debMsg("Deprecated - do not use fillHoles() ... use floodFill() instead", 1);
Real curVal, i1, i2, j1, j2, k1, k2;
Vec3i c, cTmp;
std::stack<Vec3i> undoPos;
std::stack<Real> undoVal;
std::stack<Vec3i> todoPos;
FOR_IJK_BND(*this, boundaryWidth)
{
curVal = mData[index(i, j, k)];
i1 = mData[index(i - 1, j, k)];
i2 = mData[index(i + 1, j, k)];
j1 = mData[index(i, j - 1, k)];
j2 = mData[index(i, j + 1, k)];
k1 = mData[index(i, j, k - 1)];
k2 = mData[index(i, j, k + 1)];
/* Skip cells inside and cells outside with no inside neighbours early */
if (curVal < 0.)
continue;
if (curVal > 0. && i1 > 0. && i2 > 0. && j1 > 0. && j2 > 0. && k1 > 0. && k2 > 0.)
continue;
/* Cell at c is positive (outside) and has at least one negative (inside) neighbour cell */
c = Vec3i(i, j, k);
/* Current cell is outside and has inside neighbour(s) */
undoPos.push(c);
undoVal.push(curVal);
todoPos.push(c);
/* Enforce negative cell - if search depth gets exceeded this will be reverted to the original
* value */
mData[index(c.x, c.y, c.z)] = -0.5;
while (!todoPos.empty()) {
todoPos.pop();
/* Add neighbouring positive (inside) cells to stacks and set negavtive cell value */
if (c.x > 0 && mData[index(c.x - 1, c.y, c.z)] > 0.) {
cTmp = Vec3i(c.x - 1, c.y, c.z);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
if (c.y > 0 && mData[index(c.x, c.y - 1, c.z)] > 0.) {
cTmp = Vec3i(c.x, c.y - 1, c.z);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
if (c.z > 0 && mData[index(c.x, c.y, c.z - 1)] > 0.) {
cTmp = Vec3i(c.x, c.y, c.z - 1);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
if (c.x < (*this).getSizeX() - 1 && mData[index(c.x + 1, c.y, c.z)] > 0.) {
cTmp = Vec3i(c.x + 1, c.y, c.z);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
if (c.y < (*this).getSizeY() - 1 && mData[index(c.x, c.y + 1, c.z)] > 0.) {
cTmp = Vec3i(c.x, c.y + 1, c.z);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
if (c.z < (*this).getSizeZ() - 1 && mData[index(c.x, c.y, c.z + 1)] > 0.) {
cTmp = Vec3i(c.x, c.y, c.z + 1);
undoPos.push(cTmp);
undoVal.push(mData[index(cTmp)]);
todoPos.push(cTmp);
mData[index(cTmp)] = -0.5;
}
/* Restore original value in cells if undo needed ie once cell undo count exceeds given limit
*/
if (undoPos.size() > maxDepth) {
/* Clear todo stack */
while (!todoPos.empty()) {
todoPos.pop();
}
/* Clear undo stack and revert value */
while (!undoPos.empty()) {
c = undoPos.top();
curVal = undoVal.top();
undoPos.pop();
undoVal.pop();
mData[index(c.x, c.y, c.z)] = curVal;
}
break;
}
/* Ensure that undo stack is cleared at the end if no more items in todo stack left */
if (todoPos.empty()) {
while (!undoPos.empty()) {
undoPos.pop();
}
while (!undoVal.empty()) {
undoVal.pop();
}
}
/* Pop value for next while iteration */
else {
c = todoPos.top();
}
}
}
}
//! run marching cubes to create a mesh for the 0-levelset
void LevelsetGrid::createMesh(Mesh &mesh)
{
assertMsg(is3D(), "Only 3D grids supported so far");
mesh.clear();
const Real invalidTime = invalidTimeValue();
const Real isoValue = 1e-4;
// create some temp grids
Grid<int> edgeVX(mParent);
Grid<int> edgeVY(mParent);
Grid<int> edgeVZ(mParent);
for (int k = 0; k < mSize.z - 1; k++)
for (int j = 0; j < mSize.y - 1; j++)
for (int i = 0; i < mSize.x - 1; i++) {
Real value[8] = {get(i, j, k),
get(i + 1, j, k),
get(i + 1, j + 1, k),
get(i, j + 1, k),
get(i, j, k + 1),
get(i + 1, j, k + 1),
get(i + 1, j + 1, k + 1),
get(i, j + 1, k + 1)};
// build lookup index, check for invalid times
bool skip = false;
int cubeIdx = 0;
for (int l = 0; l < 8; l++) {
value[l] *= -1;
if (-value[l] <= invalidTime)
skip = true;
if (value[l] < isoValue)
cubeIdx |= 1 << l;
}
if (skip || (mcEdgeTable[cubeIdx] == 0))
continue;
// where to look up if this point already exists
int triIndices[12];
int *eVert[12] = {&edgeVX(i, j, k),
&edgeVY(i + 1, j, k),
&edgeVX(i, j + 1, k),
&edgeVY(i, j, k),
&edgeVX(i, j, k + 1),
&edgeVY(i + 1, j, k + 1),
&edgeVX(i, j + 1, k + 1),
&edgeVY(i, j, k + 1),
&edgeVZ(i, j, k),
&edgeVZ(i + 1, j, k),
&edgeVZ(i + 1, j + 1, k),
&edgeVZ(i, j + 1, k)};
const Vec3 pos[9] = {Vec3(i, j, k),
Vec3(i + 1, j, k),
Vec3(i + 1, j + 1, k),
Vec3(i, j + 1, k),
Vec3(i, j, k + 1),
Vec3(i + 1, j, k + 1),
Vec3(i + 1, j + 1, k + 1),
Vec3(i, j + 1, k + 1)};
for (int e = 0; e < 12; e++) {
if (mcEdgeTable[cubeIdx] & (1 << e)) {
// vertex already calculated ?
if (*eVert[e] == 0) {
// interpolate edge
const int e1 = mcEdges[e * 2];
const int e2 = mcEdges[e * 2 + 1];
const Vec3 p1 = pos[e1]; // scalar field pos 1
const Vec3 p2 = pos[e2]; // scalar field pos 2
const float valp1 = value[e1]; // scalar field val 1
const float valp2 = value[e2]; // scalar field val 2
const float mu = (isoValue - valp1) / (valp2 - valp1);
// init isolevel vertex
Node vertex;
vertex.pos = p1 + (p2 - p1) * mu + Vec3(Real(0.5));
vertex.normal = getNormalized(
getGradient(
*this, i + cubieOffsetX[e1], j + cubieOffsetY[e1], k + cubieOffsetZ[e1]) *
(1.0 - mu) +
getGradient(
*this, i + cubieOffsetX[e2], j + cubieOffsetY[e2], k + cubieOffsetZ[e2]) *
(mu));
triIndices[e] = mesh.addNode(vertex) + 1;
// store vertex
*eVert[e] = triIndices[e];
}
else {
// retrieve from vert array
triIndices[e] = *eVert[e];
}
}
}
// Create the triangles...
for (int e = 0; mcTriTable[cubeIdx][e] != -1; e += 3) {
mesh.addTri(Triangle(triIndices[mcTriTable[cubeIdx][e + 0]] - 1,
triIndices[mcTriTable[cubeIdx][e + 1]] - 1,
triIndices[mcTriTable[cubeIdx][e + 2]] - 1));
}
}
// mesh.rebuildCorners();
// mesh.rebuildLookup();
// Update mdata fields
mesh.updateDataFields();
}
} // namespace Manta
View Git Blame Copy Permalink