【毕业设计】基于程序化生成和音频检测的生态仿真与3D内容生成系统----程序化生成地形算法设计
2 程序化生成地形算法设计
Input:
Output:
2.1 地形的生成
程序化生成地形算法是一种在计算机中生成地形的方法，通常用于游戏开发和虚拟现实应用。下面是几种常见的程序化生成地形算法：
Diamond-Square Algorithm（钻石-正方形算法） 该算法通过随机值填充网格的四个角落，然后计算中间点的高度值，不断重复直到整个网格被填充完毕。该算法的优点是简单易用，可以生成具有自然纹理的地形。
Perlin Noise Algorithm（柏林噪声算法） 该算法利用随机数生成器生成一系列噪声函数，并通过叠加和缩放这些函数来生成地形。该算法可以生成平滑的连续地形，但需要更复杂的数学计算。
Fractal Brownian Motion Algorithm（分形布朗运动算法） 该算法与Perlin Noise Algorithm类似，但使用了分形函数来生成更复杂的纹理。该算法可以生成非常逼真的山脉和岩石等地形，但需要更高的计算成本。
Voronoi Diagram Algorithm（沃罗诺伊图算法） 该算法将地形分割成多个区域，每个区域具有随机高度和形状。这种算法通常用于生成岛屿和湖泊等地形。
以上算法都有其优缺点，本论文综合了上述四种基本算法，并自主添加了各种处理方案程序化生成地形。
2.1.1 初始化高度场
在3D计算机图形中，高度场是指一个二维网格，每个点上都有一个高度值。高度场通常用于模拟地形、水面、云等自然现象。创建低精度版的高度场可以在一定程度上提高计算速度和效率，同时也可以减小文件大小。
使用Make2DGridPrimitive节点来创建一个平面。
struct Make2DGridPrimitive : INode {virtual void apply() override {size_t nx = get_input("nx")->get();nx = std::max(nx, (size_t)1);size_t ny = has_input("ny") ?makepositive(get_input("ny")->get()) : 0;if (!ny) ny = nx;float dx = 1.f / std::max(nx - 1, (size_t)1);float dy = 1.f / std::max(ny - 1, (size_t)1);vec3f ax = has_input("sizeX") ?get_input("sizeX")->get(): vec3f(1, 0, 0);vec3f ay = has_input("sizeY") ?get_input("sizeY")->get(): vec3f(0, 1, 0);vec3f o = has_input("origin") ?get_input("origin")->get() : vec3f(0);if (has_input("scale")) {auto obj = get_input("scale");auto scale = obj->is() ? obj->get() : obj->get();ax *= scale;ay *= scale;}auto dir = get_param("Direction");if(dir == "YZ"){ax = zeno::vec3f(0,ax[0],0);ay = zeno::vec3f(0, 0, ay[1]);}if(dir == "XZ"){ay = zeno::vec3f(0,0,ay[1]);}if (get_param("isCentered"))o -= (ax + ay) / 2;ax *= dx; ay *= dy;auto prim = std::make_shared();prim->resize(nx * ny);auto &pos = prim->add_attr("pos");
#pragma omp parallel for collapse(2)for (intptr_t y = 0; y < ny; y++) for (intptr_t x = 0; x < nx; x++) {intptr_t index = y * nx + x;vec3f p = o + x * ax + y * ay;size_t i = x + y * nx;pos[i] = p;}if (get_param("hasFaces")) {prim->tris.resize((nx - 1) * (ny - 1) * 2);
#pragma omp parallel for collapse(2)for (intptr_t y = 0; y < ny-1; y++) for (intptr_t x = 0; x < nx-1; x++) {intptr_t index = y * (nx - 1) + x;prim->tris[index * 2][2] = y * nx + x;prim->tris[index * 2][1] = y * nx + x + 1;prim->tris[index * 2][0] = (y + 1) * nx + x + 1;prim->tris[index * 2 + 1][2] = (y + 1) * nx + x + 1;prim->tris[index * 2 + 1][1] = (y + 1) * nx + x;prim->tris[index * 2 + 1][0] = y * nx + x;}}prim->userData().set("nx", std::make_shared((int)nx));//zhxxprim->userData().set("ny", std::make_shared((int)ny));//zhxxset_output("prim", std::move(prim));}
};ZENDEFNODE(Make2DGridPrimitive,{ /* inputs: */ {{"int", "nx", "2"},{"int", "ny", "0"},{"vec3f", "sizeX", "1,0,0"},{"vec3f", "sizeY", "0,1,0"},{"float", "scale", "1"},{"vec3f", "origin", "0,0,0"},}, /* outputs: */ {{"PrimitiveObject", "prim"},}, /* params: */ {{"enum XZ XY YZ", "Direction", "XZ"}, // zhxxhappy{"bool", "isCentered", "0"},{"bool", "hasFaces", "1"},}, /* category: */ {"primitive",}}); 2.1.2 原始几何体投影到高度场
用户输入基础几何体，用来生成基础的地形。如山川、山谷、 海岛、平原等基础地形的体素化数据。通过将输入的数据，转化为基础地形的体素数据。在本模块中，用户先输入几何体数据将原始几何体投影到高度场可以通过将每个几何体的顶点映射到高度场上的对应位置来实现。这样做的过程是将原始几何体的顶点坐标(x,y,z)投影到高度场上的二维坐标(x',y')，并将高度场上的高度值赋给投影点的z坐标。
这里直接将原始几何体顶点的z坐标替换为高度场上的高度值，得到投影后的几何体。
步骤如下：
创建一个圆锥
利用ParticlesTwoWrangle将圆锥的y值赋给平面的height属性
将高度场可视化出来是这样的
2.1.3 创建低精度版的高度场
取出原始平面数据，重新布置平面顶点，这里将原始平面x轴和z轴上的点分别除以了16，得到由7442个三角面构成的平面。
原始：1996002个三角面
低精度版：7442个三角面
2.1.4 把造型采样到低精度高度场，并做平滑处理
2.1.4.1体素化：
体素化（Voxelization）是将三维模型转换为由三维像素组成的体素网格的过程。在体素网格中，每个像素代表空间中的一个体积元素，类似于二维图像中的像素。通过将三维模型转换为体素网格，可以使计算机更容易处理和分析模型，例如进行碰撞检测、流体模拟等操作。体素化通常在计算机图形学、虚拟现实、计算机辅助设计等领域中使用。
本模块使用了一种基于体素的地形可视化算法，具体步骤如下：
-
计算整个场景空间的 AABB 包围盒，得到场景的最小外接矩形。
-
基于体素分辨率 α 细分包围盒，将场景分割成一系列等间距的体素单元。
-
对于包围盒顶部的每一个体素单元 V，在其底部向下发射一条射线 Ray。
-
当射线与几何体的三角形 Tri 第一次相交时，记录其相交的距离 dismin。
-
根据 dismin 和 α 计算出相交三角形在 Z 轴方向上所覆盖的体素单元 Ve，将其标记为已经被覆盖。
-
从上向下处理所有体素单元，标记已经被覆盖的体素单元。
-
标记完所有被覆盖的体素单元后，整个地形就被体素化了。
这种算法的优点是可以高效地处理大规模的地形数据，同时也能够提供高质量的地形渲染效果。缺点是需要进行大量的计算，同时对于一些特殊形状的地形可能需要进行额外的处理。
(4条消息) 【开发日志】2022.09.09 ZENO viewportwidge---raycast---sphere&AABB Box（zry还没把一些基础交互加上去，暂时不能用OBB）_EndlessDaydream的博客-CSDN博客
https://blog.csdn.net/Angelloveyatou/article/details/126724054
（bmin返回的是顶点）
#include
#include
#include
#include
#include
#include
#include
#include
#include namespace zeno {ZENO_API std::pair primBoundingBox(PrimitiveObject *prim) {if (!prim->verts.size())return {{0, 0, 0}, {0, 0, 0}};return parallel_reduce_minmax(prim->verts.begin(), prim->verts.end());
}namespace {struct PrimBoundingBox : INode {virtual void apply() override {auto prim = get_input("prim");auto extraBound = get_input2("extraBound");auto [bmin, bmax] = primBoundingBox(prim.get());if (extraBound != 0) {bmin -= extraBound;bmax += extraBound;}auto center = (bmin + bmax) / 2;auto radius = (bmax - bmin) / 2;auto diameter = bmax - bmin;set_output2("bmin", bmin);set_output2("bmax", bmax);set_output2("center", center);set_output2("radius", radius);set_output2("diameter", diameter);}
};ZENO_DEFNODE(PrimBoundingBox)({{{"PrimitiveObject", "prim"},{"float", "extraBound", "0"},},{{"vec3f", "bmin"},{"vec3f", "bmax"},{"vec3f", "center"},{"vec3f", "radius"},{"vec3f", "diameter"},},{},{"primitive"},
});
将标记后的体素单元，从bmin接口 传入Grid2DSample结点，得到整个地形的体素信息。
（高度场可视化）：
struct Grid2DSample : zeno::INode {virtual void apply() override {auto nx = get_input("nx")->get();auto ny = get_input("ny")->get();auto bmin = get_input2("bmin");auto grid = get_input("grid");auto attrT = get_param("attrT");auto channel = get_input("channel")->get();auto sampleby = get_input("sampleBy")->get();auto h = get_input("h")->get();if(grid->has_attr(channel) && grid->has_attr(sampleby)){if(attrT == "float"){sample2D(grid->attr(sampleby), grid->attr(channel), nx, ny, h, bmin);}else if(attrT == "vec3f"){sample2D(grid->attr(sampleby), grid->attr(channel), nx, ny, h, bmin);}}set_output("prim", std::move(grid));}
};
ZENDEFNODE(Grid2DSample, {{{"PrimitiveObject", "grid"},{"int", "nx","1"},{"int", "ny", "1"},{"float", "h", "1"},{"vec3f","bmin", "0,0,0" },{"string", "channel", "pos"},{"string", "sampleBy", "pos"}},{{"PrimitiveObject", "prim"}},{{"enum vec3 float", "attrT", "float"},},{"zenofx"},});
} struct ParticlesWrangle : zeno::INode {virtual void apply() override {auto prim = get_input("prim");auto code = get_input("zfxCode")->get();// BEGIN张心欣快乐自动加@INDif (auto pos = code.find("@IND"); pos != code.npos && (code.size() <= pos + 4 || !(isalnum(code[pos + 4]) || strchr("_@$", code[pos + 4]))) && (pos == 0 || !(isalnum(code[pos - 1]) || strchr("_@$", code[pos - 1])))) {auto &indatt = prim->verts.add_attr("IND");for (size_t i = 0; i < indatt.size(); i++) indatt[i] = float(i);}// END张心欣快乐自动加@INDzfx::Options opts(zfx::Options::for_x64);opts.detect_new_symbols = true;prim->foreach_attr([&] (auto const &key, auto const &attr) {int dim = ([] (auto const &v) {using T = std::decay_t;if constexpr (std::is_same_v) return 3;else if constexpr (std::is_same_v) return 1;else return 0;})(attr);dbg_printf("define symbol: @%s dim %d\n", key.c_str(), dim);opts.define_symbol('@' + key, dim);});auto params = has_input("params") ?get_input("params") :std::make_shared();{// BEGIN心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动有$F$DT$T做参数auto const &gs = *this->getGlobalState();params->lut["PI"] = objectFromLiterial((float)(std::atan(1.f) * 4));params->lut["F"] = objectFromLiterial((float)gs.frameid);params->lut["DT"] = objectFromLiterial(gs.frame_time);params->lut["T"] = objectFromLiterial(gs.frame_time * gs.frameid + gs.frame_time_elapsed);// END心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动有$F$DT$T做参数// BEGIN心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动引用portal做参数for (auto const &[key, ref]: getThisGraph()->portalIns) {if (auto i = code.find('$' + key); i != std::string::npos) {i = i + key.size() + 1;if (code.size() <= i || !std::isalnum(code[i])) {if (params->lut.count(key)) continue;dbg_printf("ref portal %s\n", key.c_str());auto res = getThisGraph()->callTempNode("PortalOut",{{"name:", objectFromLiterial(key)}}).at("port");params->lut[key] = std::move(res);}}}// END心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动引用portal做参数// BEGIN伺候心欣伺候懒得extract出变量了std::vector keys;for (auto const &[key, val]: params->lut) {keys.push_back(key);}for (auto const &key: keys) {if (!dynamic_cast(params->lut.at(key).get())) {dbg_printf("ignored non-numeric %s\n", key.c_str());params->lut.erase(key);}}// END伺候心欣伺候懒得extract出变量了}std::vector parvals;std::vector> parnames;for (auto const &[key_, par]: params->getLiterial()) {auto key = '$' + key_;auto dim = std::visit([&] (auto const &v) {using T = std::decay_t;if constexpr (std::is_convertible_v) {parvals.push_back(v[0]);parvals.push_back(v[1]);parvals.push_back(v[2]);parnames.emplace_back(key, 0);parnames.emplace_back(key, 1);parnames.emplace_back(key, 2);return 3;} else if constexpr (std::is_convertible_v) {parvals.push_back(v);parnames.emplace_back(key, 0);return 1;} else if constexpr (std::is_convertible_v) {parvals.push_back(v[0]);parvals.push_back(v[1]);parnames.emplace_back(key, 0);parnames.emplace_back(key, 1);return 2;} else {printf("invalid parameter type encountered: `%s`\n",typeid(T).name());return 0;}}, par);dbg_printf("define param: %s dim %d\n", key.c_str(), dim);opts.define_param(key, dim);//auto par = zeno::safe_any_cast(obj);}auto prog = compiler.compile(code, opts);auto exec = assembler.assemble(prog->assembly);for (auto const &[name, dim]: prog->newsyms) {dbg_printf("auto-defined new attribute: %s with dim %d\n",name.c_str(), dim);assert(name[0] == '@');auto key = name.substr(1);if (dim == 3) {prim->add_attr(key);} else if (dim == 1) {prim->add_attr(key);} else {err_printf("ERROR: bad attribute dimension for primitive: %d\n",dim);}}for (int i = 0; i < prog->params.size(); i++) {auto [name, dimid] = prog->params[i];dbg_printf("parameter %d: %s.%d\n", i, name.c_str(), dimid);assert(name[0] == '$');auto it = std::find(parnames.begin(),parnames.end(), std::pair{name, dimid});auto value = parvals.at(it - parnames.begin());dbg_printf("(valued %f)\n", value);exec->parameter(prog->param_id(name, dimid)) = value;}std::vector chs(prog->symbols.size());for (int i = 0; i < chs.size(); i++) {auto [name, dimid] = prog->symbols[i];dbg_printf("channel %d: %s.%d\n", i, name.c_str(), dimid);assert(name[0] == '@');Buffer iob;prim->attr_visit(name.substr(1),[&, dimid_ = dimid] (auto const &arr) {iob.base = (float *)arr.data() + dimid_;iob.count = arr.size();iob.stride = sizeof(arr[0]) / sizeof(float);});chs[i] = iob;}vectors_wrangle(exec, chs);set_output("prim", std::move(prim));}
};ZENDEFNODE(ParticlesWrangle, {{{"PrimitiveObject", "prim"},{"string", "zfxCode"}, {"DictObject:NumericObject", "params"}},{{"PrimitiveObject", "prim"}},{},{"zenofx"},
});struct ParticlesTwoWrangle : zeno::INode {virtual void apply() override {auto prim = get_input("prim");auto prim2 = get_input("prim2");auto code = get_input("zfxCode")->get();if (prim->size() != prim2->size()) {dbg_printf("prim and prim2 size mismatch (%d != %d), using minimal\n",prim->size(), prim2->size());}// BEGIN张心欣快乐自动加@INDif (auto pos = code.find("@IND"); pos != code.npos && (code.size() <= pos + 4 || !(isalpha(code[pos + 4]) || strchr("_@$", code[pos + 4]))) && (pos == 0 || !(isalpha(code[pos - 1]) || strchr("_@$", code[pos - 1])))) {auto &indatt = prim->verts.add_attr("IND");for (size_t i = 0; i < indatt.size(); i++) indatt[i] = float(i);}if (auto pos = code.find("@@IND"); pos != code.npos && (code.size() <= pos + 4 || !(isalpha(code[pos + 4]) || strchr("_@$", code[pos + 4]))) && (pos == 0 || !(isalpha(code[pos - 1]) || strchr("_@$", code[pos - 1])))) {auto &indatt = prim2->verts.add_attr("IND");for (size_t i = 0; i < indatt.size(); i++) indatt[i] = float(i);}// END张心欣快乐自动加@INDzfx::Options opts(zfx::Options::for_x64);opts.detect_new_symbols = true;prim->foreach_attr([&] (auto const &key, auto const &attr) {int dim = ([] (auto const &v) {using T = std::decay_t;if constexpr (std::is_same_v) return 3;else if constexpr (std::is_same_v) return 1;else return 0;})(attr);dbg_printf("define symbol: @%s dim %d\n", key.c_str(), dim);opts.define_symbol('@' + key, dim);});prim2->foreach_attr([&] (auto const &key, auto const &attr) {int dim = ([] (auto const &v) {using T = std::decay_t;if constexpr (std::is_same_v) return 3;else if constexpr (std::is_same_v) return 1;else return 0;})(attr);dbg_printf("define symbol: @@%s dim %d\n", key.c_str(), dim);opts.define_symbol("@@" + key, dim);});auto params = has_input("params") ?get_input("params") :std::make_shared();{// BEGIN心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动有$F$DT$T做参数auto const &gs = *this->getGlobalState();params->lut["PI"] = objectFromLiterial((float)(std::atan(1.f) * 4));params->lut["F"] = objectFromLiterial((float)gs.frameid);params->lut["DT"] = objectFromLiterial(gs.frame_time);params->lut["T"] = objectFromLiterial(gs.frame_time * gs.frameid + gs.frame_time_elapsed);// END心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动有$F$DT$T做参数// BEGIN心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动引用portal做参数for (auto const &[key, ref]: getThisGraph()->portalIns) {if (auto i = code.find('$' + key); i != std::string::npos) {i = i + key.size() + 1;if (code.size() <= i || !std::isalnum(code[i])) {if (params->lut.count(key)) continue;dbg_printf("ref portal %s\n", key.c_str());auto res = getThisGraph()->callTempNode("PortalOut",{{"name:", objectFromLiterial(key)}}).at("port");params->lut[key] = std::move(res);}}}// END心欣你也可以把这段代码加到其他wrangle节点去，这样这些wrangle也可以自动引用portal做参数// BEGIN伺候心欣伺候懒得extract出变量了std::vector keys;for (auto const &[key, val]: params->lut) {keys.push_back(key);}for (auto const &key: keys) {if (!dynamic_cast(params->lut.at(key).get())) {dbg_printf("ignored non-numeric %s\n", key.c_str());params->lut.erase(key);}}// END伺候心欣伺候懒得extract出变量了}std::vector parvals;std::vector> parnames;for (auto const &[key_, par]: params->getLiterial()) {auto key = '$' + key_;auto dim = std::visit([&] (auto const &v) {using T = std::decay_t;if constexpr (std::is_convertible_v) {parvals.push_back(v[0]);parvals.push_back(v[1]);parvals.push_back(v[2]);parnames.emplace_back(key, 0);parnames.emplace_back(key, 1);parnames.emplace_back(key, 2);return 3;} else if constexpr (std::is_convertible_v) {parvals.push_back(v);parnames.emplace_back(key, 0);return 1;} else if constexpr (std::is_convertible_v) {parvals.push_back(v[0]);parvals.push_back(v[1]);parnames.emplace_back(key, 0);parnames.emplace_back(key, 1);return 2;} else {printf("invalid parameter type encountered: `%s`\n",typeid(T).name());return 0;}}, par);dbg_printf("define param: %s dim %d\n", key.c_str(), dim);opts.define_param(key, dim);//auto par = zeno::safe_any_cast(obj);}auto prog = compiler.compile(code, opts);auto exec = assembler.assemble(prog->assembly);for (auto const &[name, dim]: prog->newsyms) {dbg_printf("auto-defined new attribute: %s with dim %d\n",name.c_str(), dim);assert(name[0] == '@');if (name[1] == '@') {err_printf("ERROR: cannot define new attribute %s on prim2\n",name.c_str());}auto key = name.substr(1);if (dim == 3) {prim->add_attr(key);} else if (dim == 1) {prim->add_attr(key);} else {err_printf("ERROR: bad attribute dimension for primitive: %d\n",dim);}}for (int i = 0; i < prog->params.size(); i++) {auto [name, dimid] = prog->params[i];dbg_printf("parameter %d: %s.%d\n", i, name.c_str(), dimid);assert(name[0] == '$');auto it = std::find(parnames.begin(),parnames.end(), std::pair{name, dimid});auto value = parvals.at(it - parnames.begin());dbg_printf("(valued %f)\n", value);exec->parameter(prog->param_id(name, dimid)) = value;}std::vector chs(prog->symbols.size());for (int i = 0; i < chs.size(); i++) {auto [name, dimid] = prog->symbols[i];dbg_printf("channel %d: %s.%d\n", i, name.c_str(), dimid);assert(name[0] == '@');Buffer iob;zeno::PrimitiveObject *primPtr;if (name[1] == '@') {name = name.substr(2);primPtr = prim2.get();} else {name = name.substr(1);primPtr = prim.get();}primPtr->attr_visit(name,[&, dimid_ = dimid] (auto const &arr) {iob.base = (float *)arr.data() + dimid_;iob.count = arr.size();iob.stride = sizeof(arr[0]) / sizeof(float);});chs[i] = iob;}vectors_wrangle(exec, chs);set_output("prim", std::move(prim));}
};ZENDEFNODE(ParticlesTwoWrangle, {{{"PrimitiveObject", "prim"}, {"PrimitiveObject", "prim2"},{"string", "zfxCode"}, {"DictObject:NumericObject", "params"}},{{"PrimitiveObject", "prim"}},{},{"zenofx"},
});
平滑后（高度场可视化）：
2.1.5 平滑后的数据采样回高精度的高度场
平滑后回到高精度（高度场可视化）：
对低精度grid坐标周围的四个字段值进行采样，使用双线性插值。使用 OpenMP 并行化遍历坐标的两个循环。
template
T lerp(T a, T b, float c)
{return (1.0 - c) * a + c * b;
}
template
void sample2D(std::vector &coord, std::vector &field, int nx, int ny, float h, zeno::vec3f bmin)
{std::vector temp(field.size());
#pragma omp parallel forfor(size_t tidx=0;tidx(lerp(field[idx00], field[idx01], cx), lerp(field[idx10], field[idx11], cx), cy);}
#pragma omp parallel forfor(size_t tidx=0;tidx("nx")->get();auto ny = get_input("ny")->get();auto bmin = get_input2("bmin");auto grid = get_input("grid");auto attrT = get_param("attrT");auto channel = get_input("channel")->get();auto sampleby = get_input("sampleBy")->get();auto h = get_input("h")->get();if(grid->has_attr(channel) && grid->has_attr(sampleby)){if(attrT == "float"){sample2D(grid->attr(sampleby), grid->attr(channel), nx, ny, h, bmin);}else if(attrT == "vec3f"){sample2D(grid->attr(sampleby), grid->attr(channel), nx, ny, h, bmin);}}set_output("prim", std::move(grid));}
};
ZENDEFNODE(Grid2DSample, {{{"PrimitiveObject", "grid"},{"int", "nx","1"},{"int", "ny", "1"},{"float", "h", "1"},{"vec3f","bmin", "0,0,0" },{"string", "channel", "pos"},{"string", "sampleBy", "pos"}},{{"PrimitiveObject", "prim"}},{{"enum vec3 float", "attrT", "float"},},{"zenofx"},});
} 2.2 地形变形处理
对于基础地貌生成的地形，其表面还是太过光滑，缺少很多的细节信息。所以需要对地形进行变形处理。
常用的随机数生成器被称为伪随机数生成器（Pseudo-Random Number Generator，PRNG）。PRNG 是一种根据初始种子值生成随机数序列的算法，每次生成随机数时，PRNG 的内部状态都会发生变化，因此生成的随机数序列是不确定的、伪随机的。虽然 PRNG 可以在计算机程序中方便地生成随机数，但是对于一些应用场景，例如地形生成，PRNG 生成的随机数序列可能会导致不连续、不自然的效果。
Perlin Noise 是一种基于晶格的随机数生成器，它是由 Ken Perlin 在 1983 年所发明的，主要用于计算机图形学中生成自然、有机的纹理效果。Perlin Noise 可以产生连续的、自然的随机数序列，具有相当好的空间局部性，这些特性使其非常适合用于地形生成等需要高度连续性和自然性的应用场景。
Perlin Noise 生成的随机数在空间上具有连续性，也就是说，相邻的点之间的噪声值是相似的。这个特性使得 Perlin Noise 可以用于生成自然、有机的地形细节，例如山脉、峡谷、岩石等。在地形生成的过程中，Perlin Noise 会生成一组随机数，并将其映射到一个二维或三维的空间中。然后，可以使用这些随机数来调整地形的高度、形状、颜色等属性，从而使生成的地形更加逼真。特别适合用于生成连续、自然、逼真的地形和纹理效果。
我们采用柏林噪声生成的Curl Noise （卷曲噪声）和Hybrid Noise 分别对山脉和平地进行变形。
2.2.1 Curl Noise 变形处理
Curl Noise （卷曲噪声），也称为矢量噪声或涡度约束，是一种程序噪声，通常用于计算机图形和动画中，用于模拟流体动力学和自然现象，如烟雾、火焰和云。
Curl Noise 基于矢量场的卷曲，它测量空间中每个点的流体运动的旋转或循环。为了产生卷曲噪声，将生成一个随机向量场并计算其卷曲，以产生一个新的向量场，该向量场具有与真实流体相同的旋转和涡度的基本模式。由此产生的噪声的特点是其旋转、卷曲和扭曲图案，可用于创建视觉上有趣和动态的效果。
Curl Noise 可以使用各种算法来实现，例如柏林噪声、单纯形噪声和 Worley 噪声。它通常与其他噪声函数和技术（如湍流和平流）结合使用，以创建更复杂和逼真的模拟。
Curl Noise 在计算机图形学中具有广泛的应用，包括创建逼真的水、火、烟和云，以及程序动画、地形生成和视觉效果。
Curl Noise 变形后（山丘发生变化）
struct erode_noise_analytic_simplex_2d : INode {void apply() override {auto terrain = get_input("prim_2DGrid");auto attrName = get_param("attrName");if (!terrain->has_attr(attrName)) {terrain->add_attr(attrName);}auto& noise = terrain->verts.attr(attrName);auto posLikeAttrName = get_input("posLikeAttrName")->get();if (!terrain->verts.has_attr(posLikeAttrName)){zeno::log_error("no such data named '{}'.", posLikeAttrName);}auto& pos = terrain->verts.attr(posLikeAttrName);glm::vec3 ret{};// = glm::vec3(0, 0, 0);
//#pragma omp parallel forfor (int i = 0; i < terrain->verts.size(); i++){ret = sdnoise(glm::vec2(pos[i][0], pos[i][2]));noise[i] = vec3f(ret.x, ret.y, ret.z);}set_output("prim_2DGrid", get_input("prim_2DGrid"));}
};
ZENDEFNODE(erode_noise_analytic_simplex_2d,{ /* inputs: */ {"prim_2DGrid",{"string", "posLikeAttrName", "pos"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {{"string", "attrName", "analyticNoise"},}, /* category: */ {"erode",} });
glm::vec3 srdnoise(glm::vec2 pos, double rot) {// Offset y slightly to hide some rare artifactspos.y += 0.001;// Skew to hexagonal gridglm::vec2 uv = glm::vec2(pos.x + pos.y * 0.5, pos.y);glm::vec2 i0 = floor(uv);glm::vec2 f0 = fract(uv);// Traversal orderglm::vec2 i1 = (f0.x > f0.y) ? glm::vec2(1.0, 0.0) : glm::vec2(0.0, 1.0);// Unskewed grid points in (x,y) spaceglm::vec2 p0 = glm::vec2(i0.x - i0.y * 0.5, i0.y);glm::vec2 p1 = glm::vec2(p0.x + i1.x - i1.y * 0.5, p0.y + i1.y);glm::vec2 p2 = glm::vec2(p0.x + 0.5, p0.y + 1.0);// Integer grid point indices in (u,v) spacei1 = i0 + i1;glm::vec2 i2 = i0 + glm::vec2(1.0, 1.0);// Vectors in unskewed (x,y) coordinates from// each of the simplex corners to the evaluation pointglm::vec2 d0 = pos - p0;glm::vec2 d1 = pos - p1;glm::vec2 d2 = pos - p2;glm::vec3 x = glm::vec3(p0.x, p1.x, p2.x);glm::vec3 y = glm::vec3(p0.y, p1.y, p2.y);glm::vec3 iuw = x + glm::vec3(0.5, 0.5, 0.5) * y;glm::vec3 ivw = y;// Avoid precision issues in permutationiuw = mod289(iuw);ivw = mod289(ivw);// Create gradients from indicesglm::vec2 g0 = rgrad2(glm::vec2(iuw.x, ivw.x), rot);glm::vec2 g1 = rgrad2(glm::vec2(iuw.y, ivw.y), rot);glm::vec2 g2 = rgrad2(glm::vec2(iuw.z, ivw.z), rot);// Gradients dot vectors to corresponding corners// (The derivatives of this are simply the gradients)glm::vec3 w = glm::vec3(dot(g0, d0), dot(g1, d1), dot(g2, d2));// Radial weights from corners// 0.8 is the square of 2/sqrt(5), the distance from// a grid point to the nearest simplex boundaryglm::vec3 t = glm::vec3(0.8, 0.8, 0.8) - glm::vec3(dot(d0, d0), dot(d1, d1), dot(d2, d2));// Partial derivatives for analytical gradient computationglm::vec3 dtdx = glm::vec3(-2.0, -2.0, -2.0) * glm::vec3(d0.x, d1.x, d2.x);glm::vec3 dtdy = glm::vec3(-2.0, -2.0, -2.0) * glm::vec3(d0.y, d1.y, d2.y);// Set influence of each surflet to zero outside radius sqrt(0.8)if (t.x < 0.0) {dtdx.x = 0.0;dtdy.x = 0.0;t.x = 0.0;}if (t.y < 0.0) {dtdx.y = 0.0;dtdy.y = 0.0;t.y = 0.0;}if (t.z < 0.0) {dtdx.z = 0.0;dtdy.z = 0.0;t.z = 0.0;}// Fourth power of t (and third power for derivative)glm::vec3 t2 = t * t;glm::vec3 t4 = t2 * t2;glm::vec3 t3 = t2 * t;// Final noise value is:// sum of ((radial weights) times (gradient dot vector from corner))float n = dot(t4, w);// Final analytical derivative (gradient of a sum of scalar products)glm::vec2 dt0 = glm::vec2(dtdx.x, dtdy.x) * glm::vec2(4.0, 4.0) * t3.x;glm::vec2 dn0 = t4.x * g0 + dt0 * w.x;glm::vec2 dt1 = glm::vec2(dtdx.y, dtdy.y) * glm::vec2(4.0, 4.0) * t3.y;glm::vec2 dn1 = t4.y * g1 + dt1 * w.y;glm::vec2 dt2 = glm::vec2(dtdx.z, dtdy.z) * glm::vec2(4.0, 4.0) * t3.z;glm::vec2 dn2 = t4.z * g2 + dt2 * w.z;return glm::vec3(11.0, 11.0, 11.0) * glm::vec3(n, dn0 + dn1 + dn2);
}//
// 2-D non-tiling simplex noise with fixed gradients and analytical derivative.
// This function is implemented as a wrapper to "srdnoise",
// at the minimal cost of three extra additions.
//
glm::vec3 sdnoise(glm::vec2 pos) {return srdnoise(pos, 0.0);
}
2.2.2 整体 Hybrid Noise 变形
Hybrid Noise 混合噪声是一种程序噪声，它结合了两种或多种不同的噪声函数，以创建更复杂和多样化的模式。它通常用于计算机图形学和游戏开发中，以生成纹理、地形和其他视觉元素。
创建混合噪声的最常见方法是使用数学公式（例如加法、减法、乘法或除法）将两个或多个噪声函数混合在一起。每个噪点函数都会为最终结果贡献自己独特的特征，从而实现广泛的效果和纹理。
例如，一种常见的Hybrid Noise混合噪声技术是将Perlin Noise柏林噪声和Worley噪声组合在一起。柏林噪声以其平滑、流动的模式而闻名，而 Worley 噪声会产生更不规则的细胞状图案。通过将这两个噪点函数混合在一起，您可以创建一种既具有柏林噪点的平滑度又具有 Worley 噪点复杂性的纹理。
混合噪声的另一种方法是在不同尺度或频率下使用不同的噪声函数，从而创建细节和复杂性不同的多层纹理。此技术通常用于地形生成，以创建逼真的山脉、山谷和其他自然要素。
混合噪声可以成为在计算机图形中创建视觉上有趣且多样化的纹理和图案的强大工具。通过组合不同的噪点函数和技术，您可以创建几乎无限范围的效果和纹理，仅受您的创造力和想象力的限制。
整体 Hybrid Noise 变形后：
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Worley Noise
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
glm::vec3 noise_random3(glm::vec3 p) {glm::vec3 val = sin(glm::vec3(dot(p, glm::vec3(127.1, 311.7, 74.7)),dot(p, glm::vec3(269.5, 183.3, 246.1)),dot(p, glm::vec3(113.5, 271.9, 124.6))));val *= 43758.5453123;return fract(val);
}float noise_mydistance(glm::vec3 a, glm::vec3 b, int t) {if (t == 0) {return length(a - b);}else if (t == 1) {float xx = abs(a.x - b.x);float yy = abs(a.y - b.y);float zz = abs(a.z - b.z);return max(max(xx, yy), zz);}else {float xx = abs(a.x - b.x);float yy = abs(a.y - b.y);float zz = abs(a.z - b.z);return xx + yy + zz;}
}float noise_WorleyNoise3(float px, float py, float pz, int fType, int distType, float offsetX, float offsetY, float offsetZ) {glm::vec3 pos = glm::vec3(px, py, pz);glm::vec3 offset = glm::vec3(offsetX, offsetY, offsetZ);glm::vec3 i_pos = floor(pos);glm::vec3 f_pos = fract(pos);float f1 = 9e9;float f2 = f1;for (int z = -1; z <= 1; z++) {for (int y = -1; y <= 1; y++) {for (int x = -1; x <= 1; x++) {glm::vec3 neighbor = glm::vec3(float(x), float(y), float(z));glm::vec3 point = noise_random3(i_pos + neighbor);point = (float)0.5 + (float)0.5 * sin(offset + (float)6.2831 * point);glm::vec3 featurePoint = neighbor + point;float dist = noise_mydistance(featurePoint, f_pos, distType);if (dist < f1) {f2 = f1; f1 = dist;}else if (dist < f2) {f2 = dist;}}}}if (fType == 0) {return f1;}else {return f2 - f1;}
}struct erode_noise_worley : INode {void apply() override {auto terrain = get_input("prim_2DGrid");vec3f offset;if (!has_input("seed")) {std::mt19937 gen(std::random_device{}());std::uniform_real_distribution unif(0.f, 1.f);offset = vec3f(unif(gen), unif(gen), unif(gen));}else {offset = get_input("seed")->get();}int fType = 0;auto fTypeStr = get_input2("fType");// if (fTypeStr == "F1" ) fType = 0;if (fTypeStr == "F2-F1") fType = 1;int distType = 0;auto distTypeStr = get_input2("distType");// if (distTypeStr == "Euclidean") distType = 0;if (distTypeStr == "Chebyshev") distType = 1;if (distTypeStr == "Manhattan") distType = 2;auto attrName = get_param("attrName");auto attrType = get_param("attrType");auto& pos = terrain->verts;if (!terrain->has_attr(attrName)) {if (attrType == "float3") terrain->add_attr(attrName);else if (attrType == "float") terrain->add_attr(attrName);}terrain->attr_visit(attrName, [&](auto& arr) {
#pragma omp parallel forfor (int i = 0; i < arr.size(); i++){if constexpr (is_decay_same_v){float x = noise_WorleyNoise3(pos[i][0], pos[i][1], pos[i][2], fType, distType, offset[0], offset[1], offset[2]);float y = noise_WorleyNoise3(pos[i][1], pos[i][2], pos[i][0], fType, distType, offset[0], offset[1], offset[2]);float z = noise_WorleyNoise3(pos[i][2], pos[i][0], pos[i][1], fType, distType, offset[0], offset[1], offset[2]);arr[i] = vec3f(x, y, z);}else{arr[i] = noise_WorleyNoise3(pos[i][0], pos[i][1], pos[i][2], fType, distType, offset[0], offset[1], offset[2]);}}});set_output("prim_2DGrid", get_input("prim_2DGrid"));}
};
ZENDEFNODE(erode_noise_worley,{ /* inputs: */ {"prim_2DGrid","seed",{"enum Euclidean Chebyshev Manhattan", "distType", "Euclidean"},{"enum F1 F2-F1", "fType", "F1"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {{"string", "attrName", "noise"},{"enum float float3", "attrType", "float"},}, /* category: */ {"erode",} }); 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Perlin Noise
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const int noise_permutation[] = {151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
};float noise_fade(float t) {// Fade function as defined by Ken Perlin. This eases coordinate values// so that they will ease towards integral values. This ends up smoothing// the final output.return t * t * t * (t * (t * 6 - 15) + 10); // 6t^5 - 15t^4 + 10t^3
}int noise_inc(int num) {return num + 1;
}float noise_grad(int hash, float x, float y, float z) {switch (hash & 0xF) {case 0x0: return x + y;case 0x1: return -x + y;case 0x2: return x - y;case 0x3: return -x - y;case 0x4: return x + z;case 0x5: return -x + z;case 0x6: return x - z;case 0x7: return -x - z;case 0x8: return y + z;case 0x9: return -y + z;case 0xA: return y - z;case 0xB: return -y - z;case 0xC: return y + x;case 0xD: return -y + z;case 0xE: return y - x;case 0xF: return -y - z;default: return 0;}
}float noise_perlin(float x, float y, float z)
{x = fract(x / 256.f) * 256.f;y = fract(y / 256.f) * 256.f;z = fract(z / 256.f) * 256.f;int xi = (int)x & 255; // Calculate the "unit cube" that the point asked will be located inint yi = (int)y & 255; // The left bound is ( |_x_|,|_y_|,|_z_| ) and the right bound is thatint zi = (int)z & 255; // plus 1. Next we calculate the location (from 0.0 to 1.0) in that cube.float xf = x - (int)x;float yf = y - (int)y;float zf = z - (int)z;float u = noise_fade(xf);float v = noise_fade(yf);float w = noise_fade(zf);int aaa = noise_permutation[noise_permutation[noise_permutation[xi] + yi] + zi];int aba = noise_permutation[noise_permutation[noise_permutation[xi] + noise_inc(yi)] + zi];int aab = noise_permutation[noise_permutation[noise_permutation[xi] + yi] + noise_inc(zi)];int abb = noise_permutation[noise_permutation[noise_permutation[xi] + noise_inc(yi)] + noise_inc(zi)];int baa = noise_permutation[noise_permutation[noise_permutation[noise_inc(xi)] + yi] + zi];int bba = noise_permutation[noise_permutation[noise_permutation[noise_inc(xi)] + noise_inc(yi)] + zi];int bab = noise_permutation[noise_permutation[noise_permutation[noise_inc(xi)] + yi] + noise_inc(zi)];int bbb = noise_permutation[noise_permutation[noise_permutation[noise_inc(xi)] + noise_inc(yi)] + noise_inc(zi)];float x1 = mix(noise_grad(aaa, xf, yf, zf),noise_grad(baa, xf - 1, yf, zf),u);float x2 = mix(noise_grad(aba, xf, yf - 1, zf),noise_grad(bba, xf - 1, yf - 1, zf),u);float y1 = mix(x1, x2, v);x1 = mix(noise_grad(aab, xf, yf, zf - 1),noise_grad(bab, xf - 1, yf, zf - 1),u);x2 = mix(noise_grad(abb, xf, yf - 1, zf - 1),noise_grad(bbb, xf - 1, yf - 1, zf - 1),u);float y2 = mix(x1, x2, v);return mix(y1, y2, w);
}struct erode_noise_perlin : INode {void apply() override {auto terrain = get_input("prim_2DGrid");auto attrName = get_param("attrName");auto attrType = get_param("attrType");if (!terrain->has_attr(attrName)) {if (attrType == "float3") terrain->add_attr(attrName);else if (attrType == "float") terrain->add_attr(attrName);}auto vec3fAttrName = get_input("vec3fAttrName")->get();if (!terrain->verts.has_attr(vec3fAttrName)){zeno::log_error("no such data named '{}'.", vec3fAttrName);}auto& vec3fAttr = terrain->verts.attr(vec3fAttrName);terrain->attr_visit(attrName, [&](auto& arr) {
#pragma omp parallel forfor (int i = 0; i < arr.size(); i++){if constexpr (is_decay_same_v){float x = noise_perlin(vec3fAttr[i][0], vec3fAttr[i][1], vec3fAttr[i][2]);float y = noise_perlin(vec3fAttr[i][1], vec3fAttr[i][2], vec3fAttr[i][0]);float z = noise_perlin(vec3fAttr[i][2], vec3fAttr[i][0], vec3fAttr[i][1]);arr[i] = vec3f(x, y, z);}else{arr[i] = noise_perlin(vec3fAttr[i][0], vec3fAttr[i][1], vec3fAttr[i][2]);}}});set_output("prim_2DGrid", get_input("prim_2DGrid"));}
};
ZENDEFNODE(erode_noise_perlin,{ /* inputs: */ {"prim_2DGrid",{"string", "vec3fAttrName", "pos"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {{"string", "attrName", "noise"},{"enum float float3", "attrType", "float"},}, /* category: */ {"erode",} }); 2.2.3 初始化地形数据，并标记边界，为侵蚀处理做谁备
2.3 地形侵蚀
2.3.1 地表稳固
可以模拟地面的物理特性，增加摩擦力以确保地面的稳定性。
这里采用在山表面挖小洞洞的粗暴方法增加地表稳固性。
地表稳固后：
2.3.2 热侵蚀
struct erode_tumble_material_v0 : INode {void apply() override {// 地面网格标准处理过程// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");// 获取网格大小，目前只支持方格auto& pos = terrain->verts;float cellSize = std::abs(pos[0][0] - pos[1][0]);// 用于调试和可视化auto visualEnable = get_input("visualEnable")->get();// if (visualEnable) {if (!terrain->verts.has_attr("clr")){auto& _clr = terrain->verts.add_attr("clr");std::fill(_clr.begin(), _clr.end(), vec3f(1.0, 1.0, 1.0));}auto& attr_color = terrain->verts.attr("clr");if (!terrain->verts.has_attr("debug")){auto& _debug = terrain->verts.add_attr("debug");std::fill(_debug.begin(), _debug.end(), 0);}auto& attr_debug = terrain->verts.attr("debug");// }///auto gridbias = get_input("gridbias")->get();auto cut_angle = get_input("cut_angle")->get();auto global_erosionrate = get_input("global_erosionrate")->get();auto erosionrate = get_input("erosionrate")->get();auto erodability = get_input("erodability")->get();auto removalrate = get_input("removalrate")->get();auto maxdepth = get_input("maxdepth")->get();///std::uniform_real_distribution distr(0.0, 1.0); // 设置随机分布auto seed = get_input("seed")->get();auto iterations = get_input("iterations")->get(); // 外部迭代总次数 10auto iter = get_input("iter")->get(); // 外部迭代当前次数 1~10auto i = get_input("i")->get(); // 内部迭代当前次数 0~7auto openborder = get_input("openborder")->get(); // 获取边界标记auto perm = get_input("perm")->get2();auto p_dirs = get_input("p_dirs")->get2();auto x_dirs = get_input("x_dirs")->get2();// 计算用的临时属性，必须要有if (!terrain->verts.has_attr("_height") ||!terrain->verts.has_attr("_debris") ||!terrain->verts.has_attr("_temp_height") ||!terrain->verts.has_attr("_temp_debris")){// height 和 debris 数据要从外面读取，所以属性中要有 height 和 debriszeno::log_error("Node [erode_tumble_material_v0], no such data layer named '{}' or '{}' or '{}' or '{}'.","_height", "_debris", "_temp_height", "_temp_debris");}auto& _height = terrain->verts.attr("_height"); // 计算用的临时属性auto& _debris = terrain->verts.attr("_debris");auto& _temp_height = terrain->verts.attr("_temp_height"); // 备份用的临时属性auto& _temp_debris = terrain->verts.attr("_temp_debris");// 计算// 新的，确定的，随机方向，依据上次的计算结果进行计算
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int iterseed = iter * 134775813;int color = perm[i];// randomized color order，6 种网格随机取半模式int is_red = ((id_z & 1) == 1) && (color == 1);int is_green = ((id_x & 1) == 1) && (color == 2);int is_blue = ((id_z & 1) == 0) && (color == 3);int is_yellow = ((id_x & 1) == 0) && (color == 4);int is_x_turn_x = ((id_x & 1) == 1) && ((color == 5) || (color == 6));int is_x_turn_y = ((id_x & 1) == 0) && ((color == 7) || (color == 8));// randomized direction，其实只有 4 种模式int dxs[] = { 0, p_dirs[0], 0, p_dirs[0], x_dirs[0], x_dirs[1], x_dirs[0], x_dirs[1] };int dzs[] = { p_dirs[1], 0, p_dirs[1], 0, x_dirs[0],-x_dirs[1], x_dirs[0],-x_dirs[1] };if (is_red || is_green || is_blue || is_yellow || is_x_turn_x || is_x_turn_y){int idx = Pos2Idx(id_x, id_z, nx);int dx = dxs[color - 1];int dz = dzs[color - 1];int bound_x = nx;int bound_z = nz;int clamp_x = bound_x - 1;int clamp_z = bound_z - 1;// 读取上次计算的结果float i_debris = _temp_debris[idx];float i_height = _temp_height[idx];// 移除 邻格 被边界 clamp 的格子int samplex = clamp(id_x + dx, 0, clamp_x);int samplez = clamp(id_z + dz, 0, clamp_z);int validsource = (samplex == id_x + dx) && (samplez == id_z + dz);// If we have closed borders, pretend a valid source to create// a streak conditionif (validsource){// 移除被标记为边界的格子validsource = validsource || !openborder;// 邻格 的索引号int j_idx = Pos2Idx(samplex, samplez, nx);// 邻格 的 height 和 debrisfloat j_debris = validsource ? _temp_debris[j_idx] : 0.0f;float j_height = _temp_height[j_idx];// 邻格 跟 本格 比，高的是 中格，另一个是 邻格int cidx = 0; // 中格的 id_xint cidz = 0; // 中格的 id_yfloat c_height = 0.0f;float c_debris = 0.0f;float n_debris = 0.0f;int c_idx = 0; // 中格的 idxint n_idx = 0; // 邻格的 idxint dx_check = 0; // 中格 指向 邻格 的方向int dz_check = 0;float h_diff = 0.0f; // 高度差，>=0if ((j_height - i_height) > 0.0f) // TODO: What to do when validsource is FALSE?{// look at j's neighbourscidx = samplex;cidz = samplez;c_height = j_height;c_debris = j_debris;n_debris = i_debris;c_idx = j_idx;n_idx = idx;dx_check = -dx;dz_check = -dz;h_diff = j_height - i_height;}else{// look at i's neighbourscidx = id_x;cidz = id_z;c_height = i_height;c_debris = i_debris;n_debris = j_debris;c_idx = idx;n_idx = j_idx;dx_check = dx;dz_check = dz;h_diff = i_height - j_height;}float max_diff = 0.0f;float dir_prob = 0.0f;for (int tmp_dz = -1; tmp_dz <= 1; tmp_dz++){for (int tmp_dx = -1; tmp_dx <= 1; tmp_dx++){if (!tmp_dx && !tmp_dz)continue;int tmp_samplex = clamp(cidx + tmp_dx, 0, clamp_x);int tmp_samplez = clamp(cidz + tmp_dz, 0, clamp_z);int tmp_validsource = (tmp_samplex == (cidx + tmp_dx)) && (tmp_samplez == (cidz + tmp_dz));// If we have closed borders, pretend a valid source to create// a streak condition// TODO: what is streak condition?tmp_validsource = tmp_validsource || !openborder;int tmp_j_idx = Pos2Idx(tmp_samplex, tmp_samplez, nx);float n_height = _temp_height[tmp_j_idx];float tmp_diff = n_height - (c_height);float _gridbias = clamp(gridbias, -1.0f, 1.0f);if (tmp_dx && tmp_dz)tmp_diff *= clamp(1.0f - _gridbias, 0.0f, 1.0f) / 1.4142136f;else // !tmp_dx || !tmp_dztmp_diff *= clamp(1.0f + _gridbias, 0.0f, 1.0f);if (tmp_diff <= 0.0f){if ((dx_check == tmp_dx) && (dz_check == tmp_dz))dir_prob = tmp_diff;if (tmp_diff < max_diff)max_diff = tmp_diff;}}}if (max_diff > 0.001f || max_diff < -0.001f)dir_prob = dir_prob / max_diff;int cond = 0;if (dir_prob >= 1.0f)cond = 1;else{// Making sure all drops are movingdir_prob = dir_prob * dir_prob * dir_prob * dir_prob;// Get the lower 32 bits and divide by max intunsigned int cutoff = (unsigned int)(dir_prob * 4294967295.0); // 0 ~ 1 映射到 0 ~ 4294967295.0unsigned int randval = erode_random(seed, (idx + nx * nz) * 8 + color + iterseed);cond = randval < cutoff;}if (cond){float abs_h_diff = h_diff < 0.0f ? -h_diff : h_diff;// Note: Used the neighbour mask value.float _cut_angle = clamp(cut_angle, 0.0f, 90.0f);float delta_x = cellSize * (dx && dz ? 1.4142136f : 1.0f); // 用于计算斜率的底边长度float height_removed = _cut_angle < 90.0f ? tan(_cut_angle * M_PI / 180) * delta_x : 1e10f;float height_diff = abs_h_diff - height_removed;if (height_diff < 0.0f)height_diff = 0.0f;/*float rock_adj = height_diff / (abs_h_diff + 0.000001);rock_adj = clamp(rock_adj, 0.0f, 1.0f);if (rock_adj > 0.2){rock_adj = (rock_adj - 0.25) * 4;rock_adj = clamp(rock_adj, 0.0f, 1.0f);}else{rock_adj = 0.0f;}*/float prob = ((n_debris + c_debris) != 0.0f) ? clamp((height_diff / (n_debris + c_debris)), 0.0f, 1.0f) : 1.0f;// Get the lower 32 bits and divide by max intunsigned int cutoff = (unsigned int)(prob * 4294967295.0);unsigned int randval = erode_random(seed * 3.14, (idx + nx * nz) * 8 + color + iterseed);int do_erode = randval < cutoff;float height_removal_amt = do_erode * clamp(global_erosionrate * erosionrate * erodability, 0.0f, height_diff);//float height_removal_amt = do_erode * clamp(global_erosionrate * erosionrate * erodability, 0.0f, height_diff) * (1 - rock_adj);_height[c_idx] -= height_removal_amt;float bedrock_density = 1.0f - (removalrate);if (bedrock_density > 0.0f){float newdebris = bedrock_density * height_removal_amt;if (n_debris + newdebris > maxdepth){float rollback = n_debris + newdebris - maxdepth;rollback = min(rollback, newdebris);// return the excess debris as sediment to the higher elevation cell_height[c_idx] += rollback / bedrock_density;newdebris -= rollback;}_debris[c_idx] += newdebris;}}}}}}set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_tumble_material_v0,{ /* inputs: */ {"prim_2DGrid",{"ListObject", "perm"},{"ListObject", "p_dirs"},{"ListObject", "x_dirs"},{"float", "seed", "9676.79"},{"int", "iterations", "0"},{"int", "iter", "0"},{"int", "i", "0"},{"int", "openborder", "0"},{"float", "gridbias", "0.0"},{"int", "visualEnable", "0"},{"float", "cut_angle", "35"},{"float", "global_erosionrate", "1.0"},{"float", "erosionrate", "0.03"},{"float", "erodability", "0.4"},{"float", "removalrate", "0.7"},{"float", "maxdepth", "5.0"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {}, /* category: */ {"erode",} });
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2.3.3 崩塌
struct erode_slump_b2 : INode {void apply() override {// 地面网格标准处理过程// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");// 获取网格大小，目前只支持方格auto& pos = terrain->verts;float cellSize = std::abs(pos[0][0] - pos[1][0]);// 用于调试和可视化auto visualEnable = get_input("visualEnable")->get();// if (visualEnable) {if (!terrain->verts.has_attr("clr")){auto& _clr = terrain->verts.add_attr("clr");std::fill(_clr.begin(), _clr.end(), vec3f(1.0, 1.0, 1.0));}auto& attr_color = terrain->verts.attr("clr");if (!terrain->verts.has_attr("debug")){auto& _debug = terrain->verts.add_attr("debug");std::fill(_debug.begin(), _debug.end(), 0);}auto& attr_debug = terrain->verts.attr("debug");// }// 初始化数据层// height 和 debris 只能从外部获取，不能从内部创建，因为本节点要被嵌入循环中// 初始化 height 和 debris 的过程 应该 在此节点的外部auto heightLayerName = get_input("heightLayerName")->get();auto materialLayerName = get_input("materialLayerName")->get();auto stabilitymaskLayerName = get_input("stabilitymaskLayerName")->get();if (!terrain->verts.has_attr(heightLayerName) ||!terrain->verts.has_attr(materialLayerName) ||!terrain->verts.has_attr(stabilitymaskLayerName)){// height 和 debris 数据要从外面读取，所以属性中要有 height 和 debriszeno::log_error("no such data layer named '{}' or '{}' or '{}'.", heightLayerName, materialLayerName, stabilitymaskLayerName);}auto& height = terrain->verts.attr(heightLayerName); // 读取外部数据auto& debris = terrain->verts.attr(materialLayerName);auto& stabilitymask = terrain->verts.attr(stabilitymaskLayerName);// 创建临时属性，将外部数据拷贝到临时属性，我们将使用临时属性进行计算auto& _material = terrain->verts.add_attr("_material"); // 计算用的临时属性auto& _temp_material = terrain->verts.add_attr("_temp_material"); // 备份用的临时属性#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);_material[idx] = debris[idx]; // 外部数据拷贝到临时属性_temp_material[idx] = 0; // 正式计算前会把 _debris 的数据保存在这里}}// 获取计算所需参数std::uniform_real_distribution distr(0.0, 1.0); // 设置随机分布auto seed = get_input("seed")->get();auto iterations = get_input("iterations")->get(); // 获取迭代次数auto openborder = get_input("openborder")->get(); // 获取边界标记auto repose_angle = get_input("repose_angle")->get();auto gridbias = get_input("gridbias")->get();auto quant_amt = get_input("quant_amt")->get();auto flow_rate = get_input("flow_rate")->get();// 计算for (int iter = 1; iter <= iterations; iter++){// 准备随机数组，每次迭代都都会有变化，用于网格随机取半，以及产生随机方向int perm[] = { 1, 2, 3, 4, 5, 6, 7, 8 };for (int i = 0; i < 8; i++){vec2f vec;// std::mt19937 mt(i * iterations * iter + iter); // 梅森旋转算法std::mt19937 mt(iterations * iter * 8 * i + i); // 梅森旋转算法vec[0] = distr(mt);vec[1] = distr(mt);int idx1 = floor(vec[0] * 8);int idx2 = floor(vec[1] * 8);idx1 = idx1 == 8 ? 7 : idx1;idx2 = idx2 == 8 ? 7 : idx2;int temp = perm[idx1];perm[idx1] = perm[idx2];perm[idx2] = temp;}int p_dirs[] = { -1, -1 };for (int i = 0; i < 2; i++){// std::mt19937 mt(i * iterations * iter * 20 + iter);std::mt19937 mt(iterations * iter * 2 * i + i);float rand_val = distr(mt);if (rand_val > 0.5){p_dirs[i] = 1;}else{p_dirs[i] = -1;}}int x_dirs[] = { -1, -1 };for (int i = 0; i < 2; i++){// std::mt19937 mt(i * iterations * iter * 30 + iter);std::mt19937 mt(iterations * iter * 2 * i * 10 + i);float rand_val = distr(mt);if (rand_val > 0.5){x_dirs[i] = 1;}else{x_dirs[i] = -1;}}// 分别按 8 个随机方向，每个方向算一遍for (int i = 0; i < 8; i++){// 保存上次的计算结果
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);_temp_material[idx] = _material[idx];}}// 新的，确定的，随机方向，依据上次的计算结果进行计算
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int iterseed = iter * 134775813;int color = perm[i];// randomized color order，6 种网格随机取半模式int is_red = ((id_z & 1) == 1) && (color == 1);int is_green = ((id_x & 1) == 1) && (color == 2);int is_blue = ((id_z & 1) == 0) && (color == 3);int is_yellow = ((id_x & 1) == 0) && (color == 4);int is_x_turn_x = ((id_x & 1) == 1) && ((color == 5) || (color == 6));int is_x_turn_y = ((id_x & 1) == 0) && ((color == 7) || (color == 8));// randomized direction，其实只有 4 种模式int dxs[] = { 0, p_dirs[0], 0, p_dirs[0], x_dirs[0], x_dirs[1], x_dirs[0], x_dirs[1] };int dzs[] = { p_dirs[1], 0, p_dirs[1], 0, x_dirs[0],-x_dirs[1], x_dirs[0],-x_dirs[1] };if (is_red || is_green || is_blue || is_yellow || is_x_turn_x || is_x_turn_y){int idx = Pos2Idx(id_x, id_z, nx);int dx = dxs[color - 1];int dz = dzs[color - 1];int bound_x = nx;int bound_z = nz;int clamp_x = bound_x - 1;int clamp_z = bound_z - 1;//if (idx == 249494)//{//printf(" 1-----> flow_rate = %f\n", flow_rate);//}//float flow_rate = clamp(flow_rate, 0.0f, 1.0f); // 严重错误 哈哈哈flow_rate = clamp(flow_rate, 0.0f, 1.0f);//if (idx == 249494)//{//printf(" 2-----> flow_rate = %f\n", flow_rate);//}// 读取上次计算的结果float i_material = _temp_material[idx]; // 数据来自上次的计算结果 **************************float i_height = height[idx]; // height 数据只读// 移除 邻格 被边界 clamp 的格子int samplex = clamp(id_x + dx, 0, clamp_x);int samplez = clamp(id_z + dz, 0, clamp_z);int validsource = (samplex == id_x + dx) && (samplez == id_z + dz);// If we have closed borders, pretend a valid source to create// a streak conditionif (validsource){int same_node = !validsource; // 恒等于 0 ？？？ 干嘛？ 备份，因为后面 validsource 被修改了 ？？？// 移除被标记为边界的格子validsource = validsource || !openborder;// 邻格 的索引号int j_idx = Pos2Idx(samplex, samplez, nx);// 邻格 的 height 和 debrisfloat j_material = validsource ? _temp_material[j_idx] : 0.0f;float j_height = height[j_idx];// The maximum slope at which we stop slumpingfloat _repose_angle = repose_angle;_repose_angle = clamp(_repose_angle, 0.0f, 90.0f);float delta_x = cellSize * (dx && dz ? 1.4142136f : 1.0f); // 用于计算斜率的底边长度// repose_angle 对应的高度差，停止崩塌的高度差float static_diff = _repose_angle < 90.0f ? tan(_repose_angle * M_PI / 180.0) * delta_x : 1e10f;// 包含 height 和 debris 的高度差，注意这里是 邻格 - 本格float m_diff = (j_height + j_material) - (i_height + i_material);// 邻格 跟 本格 比，高的是 中格，另一个是 邻格int cidx = 0; // 中格的 id_xint cidz = 0; // 中格的 id_yfloat c_height = 0.0f; // 中格float c_material = 0.0f; // 中格float n_material = 0.0f; // 邻格int c_idx = 0; // 中格的 idxint n_idx = 0; // 邻格的 idxint dx_check = 0; // 中格 指向 邻格 的方向int dz_check = 0;// 如果邻格比本格高，邻格->中格，本格->邻格// 高的是 中格if (m_diff > 0.0f){// look at j's neighbourscidx = samplex;cidz = samplez;c_height = j_height;c_material = j_material;n_material = i_material;c_idx = j_idx;n_idx = idx;dx_check = -dx;dz_check = -dz;}else{// look at i's neighbourscidx = id_x;cidz = id_z;c_height = i_height;c_material = i_material;n_material = j_material;c_idx = idx;n_idx = j_idx;dx_check = dx;dz_check = dz;}float sum_diffs[] = { 0.0f, 0.0f };float dir_probs[] = { 0.0f, 0.0f };float dir_prob = 0.0f;for (int diff_idx = 0; diff_idx < 2; diff_idx++){for (int tmp_dz = -1; tmp_dz <= 1; tmp_dz++){for (int tmp_dx = -1; tmp_dx <= 1; tmp_dx++){if (!tmp_dx && !tmp_dz)continue;int tmp_samplex = clamp(cidx + tmp_dx, 0, clamp_x);int tmp_samplez = clamp(cidz + tmp_dz, 0, clamp_z);int tmp_validsource = (tmp_samplex == (cidx + tmp_dx)) && (tmp_samplez == (cidz + tmp_dz));// If we have closed borders, pretend a valid source to create// a streak condition// TODO: what is streak condition?tmp_validsource = tmp_validsource || !openborder;int tmp_j_idx = Pos2Idx(tmp_samplex, tmp_samplez, nx);// 中格周围的邻格 碎屑 的高度float n_material = tmp_validsource ? _temp_material[tmp_j_idx] : 0.0f;// 中格周围邻格 地面 高度float n_height = height[tmp_j_idx];// 中格周围的邻格 地面 高度 - 中格 地面 高度float tmp_h_diff = n_height - (c_height);// 中格周围的邻格 高度 - 中格 高度float tmp_m_diff = (n_height + n_material) - (c_height + c_material);// 地面高度差 : 总高度差float tmp_diff = diff_idx == 0 ? tmp_h_diff : tmp_m_diff;float _gridbias = gridbias;_gridbias = clamp(_gridbias, -1.0f, 1.0f);// 修正高度差if (tmp_dx && tmp_dz)tmp_diff *= clamp(1.0f - _gridbias, 0.0f, 1.0f) / 1.4142136f;else // !tmp_dx || !tmp_dztmp_diff *= clamp(1.0f + _gridbias, 0.0f, 1.0f);// diff_idx = 1 的时候，前面比较过格子的总高度差，此时// 如果周边格子 不比我高，因为前面有过交换，所以至少有一个格子满足这个要求// diff_idx = 0 的时候，下面的条件不一定能满足。格子的地面有可能是最低的if (tmp_diff <= 0.0f) // 只统计比我低的邻格，所以 高度差 的说法改为 深度差{// 指定方向上，中格(我) 与 邻格 的深度差// dir_probs[0] 可能此时 >0 不会进来，此时 dir_probs[0] 保持默认值 0if ((dx_check == tmp_dx) && (dz_check == tmp_dz))dir_probs[diff_idx] = tmp_diff;// 按格子总高度计算的时候，记录 tmp_diff 最深的深度，作为 dir_probif (diff_idx && dir_prob > tmp_diff){dir_prob = tmp_diff;}// 记录比 中格 低的邻格的深度和sum_diffs[diff_idx] += tmp_diff;}}}if (diff_idx && (dir_prob > 0.001f || dir_prob < -0.001f)){// 按 (地面高度差+碎屑高度差)来计算时，流动概率 = 指定方向上的深度差 / 最大深度差dir_prob = dir_probs[diff_idx] / dir_prob;}// 另一种计算方法：指定方向上的流动概率 = 指定方向上的深度差 / 所有比我低的邻格的深度差之和// 这种概率显然比上一种方法的计算结果要 低// diff_idx == 1 时，深度差 以 (地面高度差+碎屑高度差) 来计算时 // diff_idx == 0 时，深度差 以 (地面高度差) 来计算时，可能不存在，不过已经取默认值为 0 了if (sum_diffs[diff_idx] > 0.001f || sum_diffs[diff_idx] < -0.001f)dir_probs[diff_idx] = dir_probs[diff_idx] / sum_diffs[diff_idx];}// 最多可供流失的高度差float movable_mat = (m_diff < 0.0f) ? -m_diff : m_diff;float stability_val = 0.0f;/// 这里要非常注意 ！！！！！！！！！！！！！/// 大串联的时候，这里是要有输入的 ！！！！！！！！！！！！！！！！！stability_val = clamp(stabilitymask[c_idx], 0.0f, 1.0f);if (stability_val > 0.01f){//if (idx == 249494)//{//printf("-=WB1=- iter = %i, i = %i, movable_mat = %f, stability_val = %f, c_material = %f\n",// iter, i, movable_mat, stability_val, c_material);//}// movement is slowed down according to the stability mask and not the repose angle// 只要有一点点遮罩，流失量至少减半，不过默认没有遮罩movable_mat = clamp(movable_mat * (1.0f - stability_val) * 0.5f, 0.0f, c_material);}else{//if (idx == 249494)//{//printf("-=WB2=- iter = %i, i = %i, movable_mat = %f, static_diff = %f, c_material = %f\n",// iter, i, movable_mat, static_diff, c_material);//}// 流失量根据 static_diff 修正，static_diff 是 repose angle 对应的高度差// 问题是，repose_angle 默认为 0，但可流失量仍然减半了。。。movable_mat = clamp((movable_mat - static_diff) * 0.5f, 0.0f, c_material);}// 以 height + debris 来计算float l_rat = dir_probs[1];// TODO: What is a good limit here?// 让水流继续保持足够的水量if (quant_amt > 0.001) // 默认 = 1.0movable_mat = clamp(quant_amt * ceil((movable_mat * l_rat) / quant_amt), 0.0f, c_material);elsemovable_mat *= l_rat; // 乘上概率，这样随着水量快速减少，水流很快就消失了float diff = (m_diff > 0.0f) ? movable_mat : -movable_mat;//if (idx == 249494)//{//printf("diff = %f, m_diff = %f, movable_mat = %f\n", diff, m_diff, movable_mat);//}int cond = 0;if (dir_prob >= 1.0f)cond = 1;else{// Making sure all drops are movingdir_prob = dir_prob * dir_prob * dir_prob * dir_prob;unsigned int cutoff = (unsigned int)(dir_prob * 4294967295.0); // 0 ~ 1 映射到 0 ~ 4294967295.0unsigned int randval = erode_random(seed, (idx + nx * nz) * 8 + color + iterseed);cond = randval < cutoff;}// 不参与计算的格子，或者没有流动概率的格子if (!cond || same_node)diff = 0.0f;//if (idx == 249494)//{//printf("flow_rate = %f, diff = %f, movable_mat = %f\n", flow_rate, diff, movable_mat);//}// TODO: Check if this one should be here or before quantizationdiff *= flow_rate; // 1.0float abs_diff = (diff < 0.0f) ? -diff : diff;//if (idx == 249494)//{//printf(" flow_rate = %f, diff = %f, abs_diff = %f\n", flow_rate, diff, abs_diff);//}// Update the material level// 中格失去碎屑_material[c_idx] = c_material - abs_diff;// 邻格得到碎屑_material[n_idx] = n_material + abs_diff;}}}}}}// 将计算结果返回给外部数据，并删除临时属性#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);debris[idx] = _material[idx]; // 计算结果返回给外部数据if (visualEnable){float coef = min(1, (debris[idx] / 1.0));attr_color[idx] = (1 - coef) * attr_color[idx] + coef * vec3f(0.8, 0.6, 0.4);}}}terrain->verts.erase_attr("_material");terrain->verts.erase_attr("_temp_material");set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_slump_b2,{ /* inputs: */ {"prim_2DGrid",{"string", "heightLayerName", "height"},{"string", "materialLayerName", "debris"},{"string", "stabilitymaskLayerName", "_stability"},{"int", "iterations", "10"},{"int", "openborder", "0"},{"float", "gridbias", "0.0"},{"float", "seed", "15231.3"},{"float", "repose_angle", "15.0"},{"float", "quant_amt", "0.25"},{"float", "flow_rate", "1.0"},{"int", "visualEnable", "0"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {}, /* category: */ {"deprecated",} });
// 上面的（崩塌）节点 erode_slump_b2 可以废弃了，由下面的
// erode_tumble_material_v2
// 节点代替
struct erode_tumble_material_v2 : INode {void apply() override {// 地面网格标准处理过程// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");// 获取网格大小，目前只支持方格auto& pos = terrain->verts;float cellSize = std::abs(pos[0][0] - pos[1][0]);// 用于调试和可视化auto visualEnable = get_input("visualEnable")->get();// if (visualEnable) {if (!terrain->verts.has_attr("clr")){auto& _clr = terrain->verts.add_attr("clr");std::fill(_clr.begin(), _clr.end(), vec3f(1.0, 1.0, 1.0));}auto& attr_color = terrain->verts.attr("clr");if (!terrain->verts.has_attr("debug")){auto& _debug = terrain->verts.add_attr("debug");std::fill(_debug.begin(), _debug.end(), 0);}auto& attr_debug = terrain->verts.attr("debug");// }///auto gridbias = get_input("gridbias")->get();auto repose_angle = get_input("repose_angle")->get();auto quant_amt = get_input("quant_amt")->get();auto flow_rate = get_input("flow_rate")->get();///std::uniform_real_distribution distr(0.0, 1.0);auto seed = get_input("seed")->get();auto iterations = get_input("iterations")->get();auto iter = get_input("iter")->get();auto i = get_input("i")->get();auto openborder = get_input("openborder")->get();auto perm = get_input("perm")->get2();auto p_dirs = get_input("p_dirs")->get2();auto x_dirs = get_input("x_dirs")->get2();// 计算用的临时属性，必须要有if (!terrain->verts.has_attr("height") ||!terrain->verts.has_attr("_stability") ||!terrain->verts.has_attr("_material") ||!terrain->verts.has_attr("_temp_material")){// height 和 debris 数据要从外面读取，所以属性中要有 height 和 debriszeno::log_error("Node [erode_tumble_material_v2], no such data layer named '{}' or '{}' or '{}' or '{}'.","height", "_stability", "_material", "_temp_material");}auto& height = terrain->verts.add_attr("height");auto& stabilitymask = terrain->verts.add_attr("_stability");auto& _material = terrain->verts.add_attr("_material");auto& _temp_material = terrain->verts.add_attr("_temp_material");// 计算// 新的，确定的，随机方向，依据上次的计算结果进行计算
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int iterseed = iter * 134775813;int color = perm[i];// randomized color order，6 种网格随机取半模式int is_red = ((id_z & 1) == 1) && (color == 1);int is_green = ((id_x & 1) == 1) && (color == 2);int is_blue = ((id_z & 1) == 0) && (color == 3);int is_yellow = ((id_x & 1) == 0) && (color == 4);int is_x_turn_x = ((id_x & 1) == 1) && ((color == 5) || (color == 6));int is_x_turn_y = ((id_x & 1) == 0) && ((color == 7) || (color == 8));// randomized direction，其实只有 4 种模式int dxs[] = { 0, p_dirs[0], 0, p_dirs[0], x_dirs[0], x_dirs[1], x_dirs[0], x_dirs[1] };int dzs[] = { p_dirs[1], 0, p_dirs[1], 0, x_dirs[0],-x_dirs[1], x_dirs[0],-x_dirs[1] };if (is_red || is_green || is_blue || is_yellow || is_x_turn_x || is_x_turn_y){int idx = Pos2Idx(id_x, id_z, nx);int dx = dxs[color - 1];int dz = dzs[color - 1];int bound_x = nx;int bound_z = nz;int clamp_x = bound_x - 1;int clamp_z = bound_z - 1;//if (idx == 249494)//{//printf(" 1-----> flow_rate = %f\n", flow_rate);//}//float flow_rate = clamp(flow_rate, 0.0f, 1.0f); // 严重错误 哈哈哈flow_rate = clamp(flow_rate, 0.0f, 1.0f);//if (idx == 249494)//{//printf(" 2-----> flow_rate = %f\n", flow_rate);//}// 读取上次计算的结果float i_material = _temp_material[idx]; // 数据来自上次的计算结果 **************************float i_height = height[idx]; // height 数据只读// 移除 邻格 被边界 clamp 的格子int samplex = clamp(id_x + dx, 0, clamp_x);int samplez = clamp(id_z + dz, 0, clamp_z);int validsource = (samplex == id_x + dx) && (samplez == id_z + dz);// If we have closed borders, pretend a valid source to create// a streak conditionif (validsource){int same_node = !validsource; // 恒等于 0 ？？？ 干嘛？ 备份，因为后面 validsource 被修改了 ？？？// 移除被标记为边界的格子validsource = validsource || !openborder;// 邻格 的索引号int j_idx = Pos2Idx(samplex, samplez, nx);// 邻格 的 height 和 debrisfloat j_material = validsource ? _temp_material[j_idx] : 0.0f;float j_height = height[j_idx];// The maximum slope at which we stop slumpingfloat _repose_angle = repose_angle;_repose_angle = clamp(_repose_angle, 0.0f, 90.0f);float delta_x = cellSize * (dx && dz ? 1.4142136f : 1.0f); // 用于计算斜率的底边长度// repose_angle 对应的高度差，停止崩塌的高度差float static_diff = _repose_angle < 90.0f ? tan(_repose_angle * M_PI / 180.0) * delta_x : 1e10f;// 包含 height 和 debris 的高度差，注意这里是 邻格 - 本格float m_diff = (j_height + j_material) - (i_height + i_material);// 邻格 跟 本格 比，高的是 中格，另一个是 邻格int cidx = 0; // 中格的 id_xint cidz = 0; // 中格的 id_yfloat c_height = 0.0f; // 中格float c_material = 0.0f; // 中格float n_material = 0.0f; // 邻格int c_idx = 0; // 中格的 idxint n_idx = 0; // 邻格的 idxint dx_check = 0; // 中格 指向 邻格 的方向int dz_check = 0;// 如果邻格比本格高，邻格->中格，本格->邻格// 高的是 中格if (m_diff > 0.0f){// look at j's neighbourscidx = samplex;cidz = samplez;c_height = j_height;c_material = j_material;n_material = i_material;c_idx = j_idx;n_idx = idx;dx_check = -dx;dz_check = -dz;}else{// look at i's neighbourscidx = id_x;cidz = id_z;c_height = i_height;c_material = i_material;n_material = j_material;c_idx = idx;n_idx = j_idx;dx_check = dx;dz_check = dz;}float sum_diffs[] = { 0.0f, 0.0f };float dir_probs[] = { 0.0f, 0.0f };float dir_prob = 0.0f;for (int diff_idx = 0; diff_idx < 2; diff_idx++){for (int tmp_dz = -1; tmp_dz <= 1; tmp_dz++){for (int tmp_dx = -1; tmp_dx <= 1; tmp_dx++){if (!tmp_dx && !tmp_dz)continue;int tmp_samplex = clamp(cidx + tmp_dx, 0, clamp_x);int tmp_samplez = clamp(cidz + tmp_dz, 0, clamp_z);int tmp_validsource = (tmp_samplex == (cidx + tmp_dx)) && (tmp_samplez == (cidz + tmp_dz));// If we have closed borders, pretend a valid source to create// a streak condition// TODO: what is streak condition?tmp_validsource = tmp_validsource || !openborder;int tmp_j_idx = Pos2Idx(tmp_samplex, tmp_samplez, nx);// 中格周围的邻格 碎屑 的高度float n_material = tmp_validsource ? _temp_material[tmp_j_idx] : 0.0f;// 中格周围邻格 地面 高度float n_height = height[tmp_j_idx];// 中格周围的邻格 地面 高度 - 中格 地面 高度float tmp_h_diff = n_height - (c_height);// 中格周围的邻格 高度 - 中格 高度float tmp_m_diff = (n_height + n_material) - (c_height + c_material);// 地面高度差 : 总高度差float tmp_diff = diff_idx == 0 ? tmp_h_diff : tmp_m_diff;float _gridbias = gridbias;_gridbias = clamp(_gridbias, -1.0f, 1.0f);// 修正高度差if (tmp_dx && tmp_dz)tmp_diff *= clamp(1.0f - _gridbias, 0.0f, 1.0f) / 1.4142136f;else // !tmp_dx || !tmp_dztmp_diff *= clamp(1.0f + _gridbias, 0.0f, 1.0f);// diff_idx = 1 的时候，前面比较过格子的总高度差，此时// 如果周边格子 不比我高，因为前面有过交换，所以至少有一个格子满足这个要求// diff_idx = 0 的时候，下面的条件不一定能满足。格子的地面有可能是最低的if (tmp_diff <= 0.0f) // 只统计比我低的邻格，所以 高度差 的说法改为 深度差{// 指定方向上，中格(我) 与 邻格 的深度差// dir_probs[0] 可能此时 >0 不会进来，此时 dir_probs[0] 保持默认值 0if ((dx_check == tmp_dx) && (dz_check == tmp_dz))dir_probs[diff_idx] = tmp_diff;// 按格子总高度计算的时候，记录 tmp_diff 最深的深度，作为 dir_probif (diff_idx && dir_prob > tmp_diff){dir_prob = tmp_diff;}// 记录比 中格 低的邻格的深度和sum_diffs[diff_idx] += tmp_diff;}}}if (diff_idx && (dir_prob > 0.001f || dir_prob < -0.001f)){// 按 (地面高度差+碎屑高度差)来计算时，流动概率 = 指定方向上的深度差 / 最大深度差dir_prob = dir_probs[diff_idx] / dir_prob;}// 另一种计算方法：指定方向上的流动概率 = 指定方向上的深度差 / 所有比我低的邻格的深度差之和// 这种概率显然比上一种方法的计算结果要 低// diff_idx == 1 时，深度差 以 (地面高度差+碎屑高度差) 来计算时 // diff_idx == 0 时，深度差 以 (地面高度差) 来计算时，可能不存在，不过已经取默认值为 0 了if (sum_diffs[diff_idx] > 0.001f || sum_diffs[diff_idx] < -0.001f)dir_probs[diff_idx] = dir_probs[diff_idx] / sum_diffs[diff_idx];}// 最多可供流失的高度差float movable_mat = (m_diff < 0.0f) ? -m_diff : m_diff;float stability_val = 0.0f;/// 这里要非常注意 ！！！！！！！！！！！！！/// 大串联的时候，这里是要有输入的 ！！！！！！！！！！！！！！！！！stability_val = clamp(stabilitymask[c_idx], 0.0f, 1.0f);if (stability_val > 0.01f){//if (idx == 249494)//{//printf("-=WB1=- iter = %i, i = %i, movable_mat = %f, stability_val = %f, c_material = %f\n",// iter, i, movable_mat, stability_val, c_material);//}// movement is slowed down according to the stability mask and not the repose angle// 只要有一点点遮罩，流失量至少减半，不过默认没有遮罩movable_mat = clamp(movable_mat * (1.0f - stability_val) * 0.5f, 0.0f, c_material);}else{//if (idx == 249494)//{//printf("-=WB2=- iter = %i, i = %i, movable_mat = %f, static_diff = %f, c_material = %f\n",// iter, i, movable_mat, static_diff, c_material);//}// 流失量根据 static_diff 修正，static_diff 是 repose angle 对应的高度差// 问题是，repose_angle 默认为 0，但可流失量仍然减半了。。。movable_mat = clamp((movable_mat - static_diff) * 0.5f, 0.0f, c_material);}// 以 height + debris 来计算float l_rat = dir_probs[1];// TODO: What is a good limit here?// 让水流继续保持足够的水量if (quant_amt > 0.001) // 默认 = 1.0movable_mat = clamp(quant_amt * ceil((movable_mat * l_rat) / quant_amt), 0.0f, c_material);elsemovable_mat *= l_rat; // 乘上概率，这样随着水量快速减少，水流很快就消失了float diff = (m_diff > 0.0f) ? movable_mat : -movable_mat;//if (idx == 249494)//{//printf("diff = %f, m_diff = %f, movable_mat = %f\n", diff, m_diff, movable_mat);//}int cond = 0;if (dir_prob >= 1.0f)cond = 1;else{// Making sure all drops are movingdir_prob = dir_prob * dir_prob * dir_prob * dir_prob;unsigned int cutoff = (unsigned int)(dir_prob * 4294967295.0); // 0 ~ 1 映射到 0 ~ 4294967295.0unsigned int randval = erode_random(seed, (idx + nx * nz) * 8 + color + iterseed);cond = randval < cutoff;}// 不参与计算的格子，或者没有流动概率的格子if (!cond || same_node)diff = 0.0f;//if (idx == 249494)//{//printf("flow_rate = %f, diff = %f, movable_mat = %f\n", flow_rate, diff, movable_mat);//}// TODO: Check if this one should be here or before quantizationdiff *= flow_rate; // 1.0float abs_diff = (diff < 0.0f) ? -diff : diff;//if (idx == 249494)//{//printf(" flow_rate = %f, diff = %f, abs_diff = %f\n", flow_rate, diff, abs_diff);//}// Update the material level// 中格失去碎屑_material[c_idx] = c_material - abs_diff;// 邻格得到碎屑_material[n_idx] = n_material + abs_diff;}}}}set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_tumble_material_v2,{ /* inputs: */ {"prim_2DGrid",{"ListObject", "perm"},{"ListObject", "p_dirs"},{"ListObject", "x_dirs"},{"float", "seed", "15231.3"},{"int", "iterations", "0"},{"int", "iter", "0"},{"int", "i", "0"},{"int", "openborder", "0"},{"float", "gridbias", "0.0"},{"int", "visualEnable", "0"},// 崩塌流淌相关{"float", "repose_angle", "15.0"},{"float", "quant_amt", "0.25"},{"float", "flow_rate", "1.0"},}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {}, /* category: */ {"erode",} });
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2.3.4 降水
struct erode_value2cond : INode {void apply() override {// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");///auto value = get_input("value")->get();auto seed = get_input("seed")->get();if (!terrain->verts.has_attr("cond")){terrain->verts.add_attr("cond");}auto& attr_cond = terrain->verts.attr("cond");std::fill(attr_cond.begin(), attr_cond.end(), 0.0);#pragma omp parallel forfor (int z = 0; z < nz; z++){
#pragma omp parallel forfor (int x = 0; x < nx; x++){int idx = Pos2Idx(x, z, nx);if (value >= 1.0f){attr_cond[idx] = 1;}else{value = clamp(value, 0, 1);unsigned int cutoff = (unsigned int)(value * 4294967295.0);unsigned int randval = erode_random(seed, idx + nx * nz);attr_cond[idx] = randval < cutoff;}}}set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_value2cond,{ /* inputs: */ {"prim_2DGrid",{ "float", "value", "1.0" }, // 0.0 ~ 1.0{ "float", "seed", "0.0" },}, /* outputs: */{"prim_2DGrid",}, /* params: */{}, /* category: */{"erode",} }); 2.3.5 崩塌 + 水侵蚀
struct erode_slump_b4 : INode {void apply() override {// 地面网格标准处理过程// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");// 获取网格大小，目前只支持方格auto& pos = terrain->verts;float cellSize = std::abs(pos[0][0] - pos[1][0]);// 用于调试和可视化auto visualEnable = get_input("visualEnable")->get();// if (visualEnable) {if (!terrain->verts.has_attr("clr")){auto& _clr = terrain->verts.add_attr("clr");std::fill(_clr.begin(), _clr.end(), vec3f(1.0, 1.0, 1.0));}auto& attr_color = terrain->verts.attr("clr");if (!terrain->verts.has_attr("debug")){auto& _debug = terrain->verts.add_attr("debug");std::fill(_debug.begin(), _debug.end(), 0);}auto& attr_debug = terrain->verts.attr("debug");// }// 初始化数据层// height 和 debris 只能从外部获取，不能从内部创建，因为本节点要被嵌入循环中// 初始化 height 和 debris 的过程 应该 在此节点的外部auto heightLayerName = get_input("heightLayerName")->get(); // heightauto materialLayerName = get_input("materialLayerName")->get(); // waterauto sedimentLayerName = get_input("sedimentLayerName")->get(); // sedimentauto debrisLayerName = get_input("debrisLayerName")->get(); // debrisif (!terrain->verts.has_attr(heightLayerName) ||!terrain->verts.has_attr(materialLayerName) ||!terrain->verts.has_attr(sedimentLayerName) ||!terrain->verts.has_attr(debrisLayerName)){// 需要从外部读取的数据zeno::log_error("no such data layer named '{}' or '{}' or '{}' or '{}'.",heightLayerName, materialLayerName, sedimentLayerName, debrisLayerName);}auto& height = terrain->verts.attr(heightLayerName); // 读取外部数据auto& material = terrain->verts.attr(materialLayerName);auto& sediment = terrain->verts.attr(sedimentLayerName);auto& debris = terrain->verts.attr(debrisLayerName);// 创建临时属性，将外部数据拷贝到临时属性，我们将使用临时属性进行计算auto& _height = terrain->verts.add_attr("_height");auto& _temp_height = terrain->verts.add_attr("_temp_height");auto& _material = terrain->verts.add_attr("_material");auto& _temp_material = terrain->verts.add_attr("_temp_material");auto& _debris = terrain->verts.add_attr("_debris");auto& _temp_debris = terrain->verts.add_attr("_temp_debris");auto& _sediment = terrain->verts.add_attr("_sediment");#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);_height[idx] = height[idx]; // 用于获取外部数据_temp_height[idx] = 0; // 用于存放数据备份_material[idx] = material[idx];_temp_material[idx] = 0;_debris[idx] = debris[idx];_temp_debris[idx] = 0;_sediment[idx] = sediment[idx];}}// 获取计算所需参数std::uniform_real_distribution distr(0.0, 1.0); // 设置随机分布auto openborder = get_input("openborder")->get(); // 获取边界标记// 侵蚀主参数auto global_erosionrate = get_input("global_erosionrate")->get(); // 1 全局侵蚀率auto erodability = get_input("erodability")->get(); // 1.0 侵蚀能力auto erosionrate = get_input("erosionrate")->get(); // 0.4 侵蚀率auto bank_angle = get_input("bank_angle")->get(); // 70.0 河堤侵蚀角度auto seed = get_input("seed")->get(); // 12.34// 高级参数auto removalrate = get_input("removalrate")->get(); // 0.0 风化率/水吸收率auto max_debris_depth = get_input("max_debris_depth")->get();// 5 碎屑最大深度auto gridbias = get_input("gridbias")->get(); // 0.0// 侵蚀能力调整auto max_erodability_iteration = get_input("max_erodability_iteration")->get(); // 5auto initial_erodability_factor = get_input("initial_erodability_factor")->get(); // 0.5auto slope_contribution_factor = get_input("slope_contribution_factor")->get(); // 0.8// 河床参数auto bed_erosionrate_factor = get_input("bed_erosionrate_factor")->get(); // 1 河床侵蚀率因子auto depositionrate = get_input("depositionrate")->get(); // 0.01 沉积率auto sedimentcap = get_input("sedimentcap")->get(); // 10.0 高度差转变为沉积物的比率 / 泥沙容量，每单位流动水可携带的泥沙量// 河堤参数auto bank_erosionrate_factor = get_input("bank_erosionrate_factor")->get(); // 1.0 河堤侵蚀率因子auto max_bank_bed_ratio = get_input("max_bank_bed_ratio")->get(); // 0.5 The maximum of bank to bed water column height ratio// 高于这个比值的河岸将不会在侵蚀中被视为河岸，会停止侵蚀// 河流控制auto quant_amt = get_input("quant_amt")->get(); // 0.05 流量维持率，越高流量越稳定auto iterations = get_input("iterations")->get(); // 流淌的总迭代次数// 计算// 测量计算时间 开始clock_t start, finish;printf("~~~~~~~~~~~~~~~~~~~~\n");printf("start ... ...\n");start = clock();//#pragma omp parallel for // 会出错for (int iter = 1; iter <= iterations; iter++){// 准备随机数组，每次迭代都都会有变化，用于网格随机取半，以及产生随机方向int perm[] = { 1, 2, 3, 4, 5, 6, 7, 8 };//#pragma omp parallel for // 会导致每次结果不一样for (int i = 0; i < 8; i++){vec2f vec;// std::mt19937 mt(i * iterations * iter + iter); // 梅森旋转算法std::mt19937 mt(iterations * iter * 8 * i + i); // 梅森旋转算法vec[0] = distr(mt);vec[1] = distr(mt);int idx1 = floor(vec[0] * 8);int idx2 = floor(vec[1] * 8);idx1 = idx1 == 8 ? 7 : idx1;idx2 = idx2 == 8 ? 7 : idx2;int temp = perm[idx1];perm[idx1] = perm[idx2];perm[idx2] = temp;}int p_dirs[] = { -1, -1 };for (int i = 0; i < 2; i++){// std::mt19937 mt(i * iterations * iter * 20 + iter);std::mt19937 mt(iterations * iter * 2 * i + i);float rand_val = distr(mt);if (rand_val > 0.5){p_dirs[i] = 1;}else{p_dirs[i] = -1;}}int x_dirs[] = { -1, -1 };for (int i = 0; i < 2; i++){// std::mt19937 mt(i * iterations * iter * 30 + iter);std::mt19937 mt(iterations * iter * 2 * i * 10 + i);float rand_val = distr(mt);if (rand_val > 0.5){x_dirs[i] = 1;}else{x_dirs[i] = -1;}}// 分别按 8 个随机方向，每个方向算一遍，其实只有 4 个方向模式
//#pragma omp parallel for // 会出错for (int i = 0; i < 8; i++){// 保存上次的计算结果
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);_temp_height[idx] = _height[idx];_temp_material[idx] = _material[idx];_temp_debris[idx] = _debris[idx];}}// 新的，确定的，随机方向，依据上次的计算结果进行计算
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int iterseed = iter * 134775813;int color = perm[i];// randomized color order，6 种网格随机取半模式int is_red = ((id_z & 1) == 1) && (color == 1);int is_green = ((id_x & 1) == 1) && (color == 2);int is_blue = ((id_z & 1) == 0) && (color == 3);int is_yellow = ((id_x & 1) == 0) && (color == 4);int is_x_turn_x = ((id_x & 1) == 1) && ((color == 5) || (color == 6));int is_x_turn_y = ((id_x & 1) == 0) && ((color == 7) || (color == 8));// randomized direction，其实只有 4 种模式int dxs[] = { 0, p_dirs[0], 0, p_dirs[0], x_dirs[0], x_dirs[1], x_dirs[0], x_dirs[1] };int dzs[] = { p_dirs[1], 0, p_dirs[1], 0, x_dirs[0],-x_dirs[1], x_dirs[0],-x_dirs[1] };if (is_red || is_green || is_blue || is_yellow || is_x_turn_x || is_x_turn_y){int idx = Pos2Idx(id_x, id_z, nx);int dx = dxs[color - 1];int dz = dzs[color - 1];int bound_x = nx;int bound_z = nz;int clamp_x = bound_x - 1;int clamp_z = bound_z - 1;// 读取上次计算的结果float i_height = _temp_height[idx];float i_material = _temp_material[idx];float i_debris = _temp_debris[idx];float i_sediment = _sediment[idx];// 移除 邻格 被边界 clamp 的格子int samplex = clamp(id_x + dx, 0, clamp_x);int samplez = clamp(id_z + dz, 0, clamp_z);int validsource = (samplex == id_x + dx) && (samplez == id_z + dz);// If we have closed borders, pretend a valid source to create// a streak conditionif (validsource){// 移除被标记为边界的格子 ？？？validsource = validsource || !openborder;// 邻格 的索引号int j_idx = Pos2Idx(samplex, samplez, nx);// 邻格 的数据float j_height = _temp_height[j_idx]; // height 的值一定是有的float j_material = validsource ? _temp_material[j_idx] : 0.0f; // 无效的格子不会被计算，所以可能没有值float j_debris = validsource ? _temp_debris[j_idx] : 0.0f;float j_sediment = validsource ? _sediment[j_idx] : 0.0f;// 包含 height，debris，water 的高度差，注意这里是 邻格 - 本格float m_diff = (j_height + j_debris + j_material) - (i_height + i_debris + i_material);float delta_x = cellSize * (dx && dz ? 1.4142136f : 1.0f); // 用于计算斜率的底边长度// 邻格 跟 本格 比，高的是 中格，另一个是 邻格int cidx = 0; // 中格的 id_xint cidz = 0; // 中格的 id_zfloat c_height = 0.0f; // 中格float c_material = 0.0f; // 中格float n_material = 0.0f; // 邻格float c_sediment = 0.0f;float n_sediment = 0.0f;float c_debris = 0.0f;float n_debris = 0.0f;float h_diff = 0.0f;int c_idx = 0; // 中格的 idxint n_idx = 0; // 邻格的 idxint dx_check = 0; // 中格 指向 邻格 的方向int dz_check = 0;int is_mh_diff_same_sign = 0;// 如果邻格比本格高，邻格->中格，本格->邻格// 高的是 中格if (m_diff > 0.0f){// look at j's neighbourscidx = samplex;cidz = samplez;c_height = j_height;c_material = j_material;n_material = i_material;c_sediment = j_sediment;n_sediment = i_sediment;c_debris = j_debris;n_debris = i_debris;c_idx = j_idx;n_idx = idx;dx_check = -dx;dz_check = -dz;h_diff = j_height + j_debris - (i_height + i_debris);is_mh_diff_same_sign = (h_diff * m_diff) > 0.0f;}else{// look at i's neighbourscidx = id_x;cidz = id_z;c_height = i_height;c_material = i_material;n_material = j_material;c_sediment = i_sediment;n_sediment = j_sediment;c_debris = i_debris;n_debris = j_debris;c_idx = idx;n_idx = j_idx;dx_check = dx;dz_check = dz;h_diff = i_height + i_debris - (j_height + j_debris);is_mh_diff_same_sign = (h_diff * m_diff) > 0.0f;}h_diff = (h_diff < 0.0f) ? -h_diff : h_diff;float sum_diffs[] = { 0.0f, 0.0f };float dir_probs[] = { 0.0f, 0.0f };float dir_prob = 0.0f;for (int diff_idx = 0; diff_idx < 2; diff_idx++){for (int tmp_dz = -1; tmp_dz <= 1; tmp_dz++){for (int tmp_dx = -1; tmp_dx <= 1; tmp_dx++){if (!tmp_dx && !tmp_dz)continue;int tmp_samplex = clamp(cidx + tmp_dx, 0, clamp_x);int tmp_samplez = clamp(cidz + tmp_dz, 0, clamp_z);int tmp_validsource = (tmp_samplex == (cidx + tmp_dx)) && (tmp_samplez == (cidz + tmp_dz));// If we have closed borders, pretend a valid source to create// a streak condition// TODO: what is streak condition?tmp_validsource = tmp_validsource || !openborder;int tmp_j_idx = Pos2Idx(tmp_samplex, tmp_samplez, nx);// 中格周围的邻格 水，碎屑 的高度float tmp_n_material = tmp_validsource ? _temp_material[tmp_j_idx] : 0.0f;float tmp_n_debris = tmp_validsource ? _temp_debris[tmp_j_idx] : 0.0f;// 中格周围的邻格 地面 的高度float n_height = _temp_height[tmp_j_idx];// 中格周围的邻格 无水高度 - 中格 无水高度float tmp_h_diff = n_height + tmp_n_debris - (c_height + c_debris);// 中格周围的邻格 带水高度 - 中格 带水高度float tmp_m_diff = (n_height + tmp_n_debris + tmp_n_material) - (c_height + c_debris + c_material);float tmp_diff = diff_idx == 0 ? tmp_h_diff : tmp_m_diff;float _gridbias = gridbias;_gridbias = clamp(_gridbias, -1.0f, 1.0f);// 修正高度差if (tmp_dx && tmp_dz)tmp_diff *= clamp(1.0f - _gridbias, 0.0f, 1.0f) / 1.4142136f;else // !tmp_dx || !tmp_dztmp_diff *= clamp(1.0f + _gridbias, 0.0f, 1.0f);// diff_idx = 1 的时候，前面比较过格子的总高度差，此时// 如果周边格子 不比我高，因为前面有过交换，所以至少有一个格子满足这个要求// diff_idx = 0 的时候，下面的条件不一定能满足。格子的地面有可能是最低的if (tmp_diff <= 0.0f) // 只统计比我低的邻格，所以 高度差 的说法改为 深度差{// 指定方向上，中格(我) 与 邻格 的深度差// dir_probs[0] 可能此时 >0 不会进来，此时 dir_probs[0] 保持默认值 0if ((dx_check == tmp_dx) && (dz_check == tmp_dz))dir_probs[diff_idx] = tmp_diff;// 按格子总高度计算的时候，记录 tmp_diff 最深的深度，作为 dir_probif (diff_idx && (tmp_diff < dir_prob))dir_prob = tmp_diff;// 记录比 中格 低的邻格的深度和sum_diffs[diff_idx] += tmp_diff;}}}if (diff_idx && (dir_prob > 0.001f || dir_prob < -0.001f)){// 按 (地面高度差+碎屑高度差)来计算时，流动概率 = 指定方向上的深度差 / 最大深度差dir_prob = dir_probs[diff_idx] / dir_prob;}else // add by wangbo{dir_prob = 0.0f;}// 另一种计算方法：指定方向上的流动概率 = 指定方向上的深度差 / 所有比我低的邻格的深度差之和// 这种概率显然比上一种方法的计算结果要 低// diff_idx == 1 时，深度差 以 (地面高度差+碎屑高度差) 来计算时// diff_idx == 0 时，深度差 以 (地面高度差) 来计算时，可能不存在，不过已经取默认值为 0 了if (sum_diffs[diff_idx] > 0.001f || sum_diffs[diff_idx] < -0.001f){dir_probs[diff_idx] = dir_probs[diff_idx] / sum_diffs[diff_idx];}else // add by wangbo{dir_probs[diff_idx] = 0.0f;}}// 最多可供流失的高度差float movable_mat = (m_diff < 0.0f) ? -m_diff : m_diff;// 它首先会被clamp(0,c_material)，以保证有足够的材料被移动movable_mat = clamp(movable_mat * 0.5f, 0.0f, c_material);// 以 height + debris + water 来计算float l_rat = dir_probs[1];// TODO: What is a good limit here?// 让水流继续保持足够的水量if (quant_amt > 0.001) // 默认 = 1.0movable_mat = clamp(quant_amt * ceil((movable_mat * l_rat) / quant_amt), 0.0f, c_material);elsemovable_mat *= l_rat; // 乘上概率，这样随着水量快速减少，水流很快就消失了float diff = (m_diff > 0.0f) ? movable_mat : -movable_mat;int cond = 0;if (dir_prob >= 1.0f)cond = 1;else{// Making sure all drops are movingdir_prob = dir_prob * dir_prob * dir_prob * dir_prob;unsigned int cutoff = (unsigned int)(dir_prob * 4294967295.0); // 0 ~ 1 映射到 0 ~ 4294967295.0unsigned int randval = erode_random(seed, (idx + nx * nz) * 8 + color + iterseed);cond = randval < cutoff;}// 不参与计算的格子，或者没有流动概率的格子if (!cond)diff = 0.0f;/// 下面开始计算侵蚀，河床，河堤/// 通过 h_diff 计算 沉积条件，用于产生 河床 和 河堤float slope_cont = (delta_x > 0.0f) ? (h_diff / delta_x) : 0.0f; // 斜率=对边/临边// 沉积因子 = 1 / (1 + 斜率)，斜率大的地方沉积的倾向小float kd_factor = clamp((1 / (1 + (slope_contribution_factor * slope_cont))), 0.0f, 1.0f);// 当前迭代索引iter(1~50) / 最大侵蚀迭代次数(参数面板填写)float norm_iter = clamp(((float)iter / (float)max_erodability_iteration), 0.0f, 1.0f);// 侵蚀因子：地面(height+debris)斜率，斜率贡献因子，深度差/深度差之和，初始侵蚀因子，迭代递加侵蚀因子float ks_factor = clamp((1 - (slope_contribution_factor * exp(-slope_cont))) * sqrt(dir_probs[0]) *(initial_erodability_factor + ((1.0f - initial_erodability_factor) * sqrt(norm_iter))),0.0f, 1.0f);// 中格侵蚀率float c_ks = global_erosionrate * erosionrate * erodability * ks_factor;// 邻格沉积率float n_kd = depositionrate * kd_factor;n_kd = clamp(n_kd, 0.0f, 1.0f);// 类似通过风化率计算侵蚀产生的碎屑float _removalrate = removalrate;float bedrock_density = 1.0f - _removalrate;// 通过 m_diff 可能包含水面的高度差//Kc //Kd //Ks float abs_diff = (diff < 0.0f) ? -diff : diff; // 可能包含水面的高度差// sedimentcap：泥沙容量，每单位流动水可携带的泥沙量。// 容量越大，材料在开始沉积多余沉积物之前被侵蚀的时间就越长。float sediment_limit = sedimentcap * abs_diff; // 根据泥沙容量计算的水中泥沙上限float ent_check_diff = sediment_limit - c_sediment;// sediment_limit - c_sediment > 0，// 意味着水中可以携带的泥沙上限超过了 中格 的沉积物// 这会导致更大的侵蚀，倾向挖宽河床，这主要是一个侵蚀过程// sediment_limit - c_sediment < 0，// 意味着水中可以携带的泥沙的能力降低，// 倾向于水中的泥沙向地面沉积，这是主要是一个沉积过程// for current cellif (ent_check_diff > 0.0f) // sediment_limit > c_sediment{// 中格因为侵蚀而被溶解的物质float dissolve_amt = c_ks * bed_erosionrate_factor * abs_diff;// 优先溶解碎屑层，但是碎屑层最大的量也只有 c_debrisfloat dissolved_debris = min(c_debris, dissolve_amt);// 中格碎屑被溶解后，还剩下的部分_debris[c_idx] -= dissolved_debris; // 碎屑被侵蚀// 如果中格碎屑被溶完了，还不够，就开始溶解 height 层_height[c_idx] -= (dissolve_amt - dissolved_debris); // height 被侵蚀// 沉积，数据来自上一 frame 计算的结果，50{8{}} 循环内简单重复计算// 中格的沉积物被冲走一半_sediment[c_idx] -= c_sediment / 2; // 沉积物被侵蚀// 风化后仍有剩余if (bedrock_density > 0.0f){// 被冲走的那一半沉积物 + 溶解物的沉积，这些沉积会堆积到 邻格float newsediment = c_sediment / 2 + (dissolve_amt * bedrock_density);// 假设沉积物都会堆积到邻格，如果超过最大碎屑高度if (n_sediment + newsediment > max_debris_depth){// 回滚float rollback = n_sediment + newsediment - max_debris_depth;// 回滚量不可以超过刚刚计算的沉积高度rollback = min(rollback, newsediment);// return the excess sediment_height[c_idx] += rollback / bedrock_density; // 向上修正 height 高度newsediment -= rollback; // 向下修正 沉积高度}_sediment[n_idx] += newsediment; // 邻格沉积物增加}}else // sediment_limit <= c_sediment，这主要是一个沉积过程{float c_kd = depositionrate * kd_factor; // 计算沉积系数c_kd = clamp(c_kd, 0.0f, 1.0f);{// -ent_check_diff = 高度差产生的泥沙 - 能被水流携带走的泥沙// 这些过剩的泥沙会成为 碎屑 和 沉积物// 碎屑的定义：高度差侵蚀的直接结果// 沉积的定义：泥沙被被水搬运到邻格沉降的结果_debris[c_idx] += (c_kd * -ent_check_diff); // 中格碎屑增加_sediment[c_idx] = (1 - c_kd) * -ent_check_diff; // 剩下的变成中格的沉积物n_sediment += sediment_limit; // 被带走的沉积物到了邻格_debris[n_idx] += (n_kd * n_sediment); // 邻格的碎屑增加_sediment[n_idx] = (1 - n_kd) * n_sediment; // 剩下的变成邻格的沉积物}// 河岸 河床的侵蚀，碎屑转移过程，不涉及 sedimentint b_idx = 0; // 岸的位置索引int r_idx = 0; // 河的位置索引float b_material = 0.0f; // 岸的水高度float r_material = 0.0f; // 河的水高度float b_debris = 0.0f; // 岸的碎屑高度float r_debris = 0.0f; // 河的碎屑高度float r_sediment = 0.0f; // 河的沉积高度if (is_mh_diff_same_sign) // 中格的水高，地也高{b_idx = c_idx; // 中格地高，是岸r_idx = n_idx; // 邻格地低，是河b_material = c_material; // 岸的水高度r_material = n_material; // 河的水高度b_debris = c_debris; // 岸的碎屑高度r_debris = n_debris; // 河的碎屑高度r_sediment = n_sediment; // 河的沉积高度}else // 中格 水高，地低{b_idx = n_idx; // 邻格地高，是岸r_idx = c_idx; // 中格地低，是河b_material = n_material;r_material = c_material;b_debris = n_debris;r_debris = c_debris;r_sediment = c_sediment;}// 河中每单位水的侵蚀量float erosion_per_unit_water = global_erosionrate * erosionrate * bed_erosionrate_factor * erodability * ks_factor;// 河中有水 && 岸的水量/河的水量if (r_material != 0.0f &&(b_material / r_material) < max_bank_bed_ratio &&r_sediment > (erosion_per_unit_water * max_bank_bed_ratio)){// NOTE: Increase the bank erosion to get a certain// angle faster. This would make the river cuts less// deep.float height_to_erode = global_erosionrate * erosionrate * bank_erosionrate_factor * erodability * ks_factor;float _bank_angle = bank_angle;_bank_angle = clamp(_bank_angle, 0.0f, 90.0f);float safe_diff = _bank_angle < 90.0f ? tan(_bank_angle * M_PI / 180.0) * delta_x : 1e10f;float target_height_removal = (h_diff - safe_diff) < 0.0f ? 0.0f : h_diff - safe_diff;float dissolve_amt = clamp(height_to_erode, 0.0f, target_height_removal);float dissolved_debris = min(b_debris, dissolve_amt);_debris[b_idx] -= dissolved_debris; // 岸的碎屑被侵蚀float division = 1 / (1 + safe_diff);_height[b_idx] -= (dissolve_amt - dissolved_debris); // 岸的 height 被侵蚀if (bedrock_density > 0.0f) // 有沉积{float newdebris = (1 - division) * (dissolve_amt * bedrock_density);if (b_debris + newdebris > max_debris_depth){float rollback = b_debris + newdebris - max_debris_depth;rollback = min(rollback, newdebris);// return the excess debris_height[b_idx] += rollback / bedrock_density;newdebris -= rollback;}_debris[b_idx] += newdebris; // 河岸沉积newdebris = division * (dissolve_amt * bedrock_density);if (r_debris + newdebris > max_debris_depth){float rollback = r_debris + newdebris - max_debris_depth;rollback = min(rollback, newdebris);// return the excess debris_height[b_idx] += rollback / bedrock_density;newdebris -= rollback;}_debris[r_idx] += newdebris; // 河床沉积}}}// Update the material level// 水往低处流，这很简单，麻烦的是上面，注意这里的索引号，不是中格-邻格模式，而是本格-邻格模式_material[idx] = i_material + diff; // 本格更新水的高度_material[j_idx] = j_material - diff; // 邻格更新水的高度}}}}}}// 测量计算时间 结束finish = clock();double duration = ((double)finish - (double)start) / CLOCKS_PER_SEC;printf("... ... end!\n");printf("-= %f seconds =-\n", duration);// 将计算结果返回给外部数据，并删除临时属性#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int idx = Pos2Idx(id_x, id_z, nx);height[idx] = _height[idx]; // 计算结果返回给外部数据material[idx] = _material[idx];debris[idx] = _debris[idx];sediment[idx] = _sediment[idx];if (visualEnable){float coef = min(1, (material[idx] / 1.0));attr_color[idx] = (1 - coef) * attr_color[idx] + coef * vec3f(0.15, 0.45, 0.9);}}}terrain->verts.erase_attr("_sediment");terrain->verts.erase_attr("_height");terrain->verts.erase_attr("_temp_height");terrain->verts.erase_attr("_debris");terrain->verts.erase_attr("_temp_debris");terrain->verts.erase_attr("_material");terrain->verts.erase_attr("_temp_material");set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_slump_b4,{ /* inputs: */ {"prim_2DGrid",// 需要用到的属性/数据{"string", "heightLayerName", "height"},{"string", "materialLayerName", "water"},{"string", "sedimentLayerName", "sediment"},{"string", "debrisLayerName", "debris"},// 杂项{"int", "openborder", "0"}, // 获取边界标记{"int", "visualEnable", "0"}, // 开启可视化// 侵蚀主参数{"float", "global_erosionrate", "1.0"}, // 全局侵蚀率{"float", "erodability", "1.0"}, // 侵蚀能力{"float", "erosionrate", "0.4"}, // 侵蚀率{"float", "bank_angle", "70.0"}, // 河堤侵蚀角度{"float", "seed", "12.34"},// 高级参数{"float", "removalrate", "0.1"}, // 风化率/水吸收率{"float", "max_debris_depth", "5.0"}, // 碎屑最大深度{"float", "gridbias", "0.0"},// 侵蚀能力调整{"int", "max_erodability_iteration", "5"}, // 最大侵蚀能力迭代次数{"float", "initial_erodability_factor", "0.5"}, // 初始侵蚀能力因子{"float", "slope_contribution_factor", "0.8"}, // “地面斜率”对“侵蚀”和“沉积”的影响，“地面斜率大” -> 侵蚀因子大，沉积因子小// 河床参数{"float", "bed_erosionrate_factor", "1.0"}, // 河床侵蚀率因子{"float", "depositionrate", "0.01"}, // 沉积率{"float", "sedimentcap", "10.0"}, // 高度差转变为沉积物的比率 / 泥沙容量，每单位流动水可携带的泥沙量// 河堤参数{"float", "bank_erosionrate_factor", "1.0"}, // 河堤侵蚀率因子{"float", "max_bank_bed_ratio", "0.5"}, // 高于这个比值的河岸将不会在侵蚀中被视为河岸，会停止侵蚀// 河网控制{"float", "quant_amt", "0.05"}, // 流量维持率，越高河流流量越稳定{"int", "iterations", "40"}, // 流淌的总迭代次数}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {}, /* category: */ {"deprecated",} }); struct erode_tumble_material_v4 : INode {void apply() override {// 地面网格标准处理过程// 获取地形auto terrain = get_input("prim_2DGrid");// 获取用户数据，里面存有网格精度int nx, nz;auto& ud = terrain->userData();if ((!ud.has("nx")) || (!ud.has("nz"))){zeno::log_error("no such UserData named '{}' and '{}'.", "nx", "nz");}nx = ud.get2("nx");nz = ud.get2("nz");// 获取网格大小，目前只支持方格auto& pos = terrain->verts;float cellSize = std::abs(pos[0][0] - pos[1][0]);// 用于调试和可视化auto visualEnable = get_input("visualEnable")->get();// if (visualEnable) {if (!terrain->verts.has_attr("clr")){auto& _clr = terrain->verts.add_attr("clr");std::fill(_clr.begin(), _clr.end(), vec3f(1.0, 1.0, 1.0));}auto& attr_color = terrain->verts.attr("clr");if (!terrain->verts.has_attr("debug")){auto& _debug = terrain->verts.add_attr("debug");std::fill(_debug.begin(), _debug.end(), 0);}auto& attr_debug = terrain->verts.attr("debug");// }///// 侵蚀主参数auto global_erosionrate = get_input("global_erosionrate")->get(); // 1 全局侵蚀率auto erodability = get_input("erodability")->get(); // 1.0 侵蚀能力auto erosionrate = get_input("erosionrate")->get(); // 0.4 侵蚀率auto bank_angle = get_input("bank_angle")->get(); // 70.0 河堤侵蚀角度auto seed = get_input("seed")->get(); // 12.34// 高级参数auto removalrate = get_input("removalrate")->get(); // 0.0 风化率/水吸收率auto max_debris_depth = get_input("max_debris_depth")->get();// 5 碎屑最大深度auto gridbias = get_input("gridbias")->get(); // 0.0// 侵蚀能力调整auto max_erodability_iteration = get_input("max_erodability_iteration")->get(); // 5auto initial_erodability_factor = get_input("initial_erodability_factor")->get(); // 0.5auto slope_contribution_factor = get_input("slope_contribution_factor")->get(); // 0.8// 河床参数auto bed_erosionrate_factor = get_input("bed_erosionrate_factor")->get(); // 1 河床侵蚀率因子auto depositionrate = get_input("depositionrate")->get(); // 0.01 沉积率auto sedimentcap = get_input("sedimentcap")->get(); // 10.0 高度差转变为沉积物的比率 / 泥沙容量，每单位流动水可携带的泥沙量// 河堤参数auto bank_erosionrate_factor = get_input("bank_erosionrate_factor")->get(); // 1.0 河堤侵蚀率因子auto max_bank_bed_ratio = get_input("max_bank_bed_ratio")->get(); // 0.5 The maximum of bank to bed water column height ratio// 高于这个比值的河岸将不会在侵蚀中被视为河岸，会停止侵蚀// 河流控制auto quant_amt = get_input("quant_amt")->get(); // 0.05 流量维持率，越高流量越稳定auto iterations = get_input("iterations")->get(); // 流淌的总迭代次数///std::uniform_real_distribution distr(0.0, 1.0);// auto seed = get_input("seed")->get();// auto iterations = get_input("iterations")->get();auto iter = get_input("iter")->get();auto i = get_input("i")->get();auto openborder = get_input("openborder")->get();auto perm = get_input("perm")->get2();auto p_dirs = get_input("p_dirs")->get2();auto x_dirs = get_input("x_dirs")->get2();// 计算用的临时属性，必须要有if (!terrain->verts.has_attr("_height") ||!terrain->verts.has_attr("_temp_height") ||!terrain->verts.has_attr("_material") ||!terrain->verts.has_attr("_temp_material") ||!terrain->verts.has_attr("_debris") ||!terrain->verts.has_attr("_temp_debris") ||!terrain->verts.has_attr("_sediment")){// height 和 debris 数据要从外面读取，所以属性中要有 height 和 debriszeno::log_error("Node [erode_tumble_material_v4], no such data layer named '{}' or '{}' or '{}' or '{}' or '{}' or '{}' or '{}'.","_height", "_temp_height", "_material", "_temp_material", "_debris", "_temp_debris", "_sediment");}auto& _height = terrain->verts.add_attr("_height");auto& _temp_height = terrain->verts.add_attr("_temp_height");auto& _material = terrain->verts.add_attr("_material");auto& _temp_material = terrain->verts.add_attr("_temp_material");auto& _debris = terrain->verts.add_attr("_debris");auto& _temp_debris = terrain->verts.add_attr("_temp_debris");auto& _sediment = terrain->verts.add_attr("_sediment");// 计算// 新的，确定的，随机方向，依据上次的计算结果进行计算
#pragma omp parallel forfor (int id_z = 0; id_z < nz; id_z++){
#pragma omp parallel forfor (int id_x = 0; id_x < nx; id_x++){int iterseed = iter * 134775813;int color = perm[i];// randomized color order，6 种网格随机取半模式int is_red = ((id_z & 1) == 1) && (color == 1);int is_green = ((id_x & 1) == 1) && (color == 2);int is_blue = ((id_z & 1) == 0) && (color == 3);int is_yellow = ((id_x & 1) == 0) && (color == 4);int is_x_turn_x = ((id_x & 1) == 1) && ((color == 5) || (color == 6));int is_x_turn_y = ((id_x & 1) == 0) && ((color == 7) || (color == 8));// randomized direction，其实只有 4 种模式int dxs[] = { 0, p_dirs[0], 0, p_dirs[0], x_dirs[0], x_dirs[1], x_dirs[0], x_dirs[1] };int dzs[] = { p_dirs[1], 0, p_dirs[1], 0, x_dirs[0],-x_dirs[1], x_dirs[0],-x_dirs[1] };if (is_red || is_green || is_blue || is_yellow || is_x_turn_x || is_x_turn_y){int idx = Pos2Idx(id_x, id_z, nx);int dx = dxs[color - 1];int dz = dzs[color - 1];int bound_x = nx;int bound_z = nz;int clamp_x = bound_x - 1;int clamp_z = bound_z - 1;// 读取上次计算的结果float i_height = _temp_height[idx];float i_material = _temp_material[idx];float i_debris = _temp_debris[idx];float i_sediment = _sediment[idx];// 移除 邻格 被边界 clamp 的格子int samplex = clamp(id_x + dx, 0, clamp_x);int samplez = clamp(id_z + dz, 0, clamp_z);int validsource = (samplex == id_x + dx) && (samplez == id_z + dz);// If we have closed borders, pretend a valid source to create// a streak conditionif (validsource){// 移除被标记为边界的格子 ？？？validsource = validsource || !openborder;// 邻格 的索引号int j_idx = Pos2Idx(samplex, samplez, nx);// 邻格 的数据float j_height = _temp_height[j_idx]; // height 的值一定是有的float j_material = validsource ? _temp_material[j_idx] : 0.0f; // 无效的格子不会被计算，所以可能没有值float j_debris = validsource ? _temp_debris[j_idx] : 0.0f;float j_sediment = validsource ? _sediment[j_idx] : 0.0f;// 包含 height，debris，water 的高度差，注意这里是 邻格 - 本格float m_diff = (j_height + j_debris + j_material) - (i_height + i_debris + i_material);float delta_x = cellSize * (dx && dz ? 1.4142136f : 1.0f); // 用于计算斜率的底边长度// 邻格 跟 本格 比，高的是 中格，另一个是 邻格int cidx = 0; // 中格的 id_xint cidz = 0; // 中格的 id_zfloat c_height = 0.0f; // 中格float c_material = 0.0f; // 中格float n_material = 0.0f; // 邻格float c_sediment = 0.0f;float n_sediment = 0.0f;float c_debris = 0.0f;float n_debris = 0.0f;float h_diff = 0.0f;int c_idx = 0; // 中格的 idxint n_idx = 0; // 邻格的 idxint dx_check = 0; // 中格 指向 邻格 的方向int dz_check = 0;int is_mh_diff_same_sign = 0;// 如果邻格比本格高，邻格->中格，本格->邻格// 高的是 中格if (m_diff > 0.0f){// look at j's neighbourscidx = samplex;cidz = samplez;c_height = j_height;c_material = j_material;n_material = i_material;c_sediment = j_sediment;n_sediment = i_sediment;c_debris = j_debris;n_debris = i_debris;c_idx = j_idx;n_idx = idx;dx_check = -dx;dz_check = -dz;h_diff = j_height + j_debris - (i_height + i_debris);is_mh_diff_same_sign = (h_diff * m_diff) > 0.0f;}else{// look at i's neighbourscidx = id_x;cidz = id_z;c_height = i_height;c_material = i_material;n_material = j_material;c_sediment = i_sediment;n_sediment = j_sediment;c_debris = i_debris;n_debris = j_debris;c_idx = idx;n_idx = j_idx;dx_check = dx;dz_check = dz;h_diff = i_height + i_debris - (j_height + j_debris);is_mh_diff_same_sign = (h_diff * m_diff) > 0.0f;}h_diff = (h_diff < 0.0f) ? -h_diff : h_diff;float sum_diffs[] = { 0.0f, 0.0f };float dir_probs[] = { 0.0f, 0.0f };float dir_prob = 0.0f;for (int diff_idx = 0; diff_idx < 2; diff_idx++){for (int tmp_dz = -1; tmp_dz <= 1; tmp_dz++){for (int tmp_dx = -1; tmp_dx <= 1; tmp_dx++){if (!tmp_dx && !tmp_dz)continue;int tmp_samplex = clamp(cidx + tmp_dx, 0, clamp_x);int tmp_samplez = clamp(cidz + tmp_dz, 0, clamp_z);int tmp_validsource = (tmp_samplex == (cidx + tmp_dx)) && (tmp_samplez == (cidz + tmp_dz));// If we have closed borders, pretend a valid source to create// a streak condition// TODO: what is streak condition?tmp_validsource = tmp_validsource || !openborder;int tmp_j_idx = Pos2Idx(tmp_samplex, tmp_samplez, nx);// 中格周围的邻格 水，碎屑 的高度float tmp_n_material = tmp_validsource ? _temp_material[tmp_j_idx] : 0.0f;float tmp_n_debris = tmp_validsource ? _temp_debris[tmp_j_idx] : 0.0f;// 中格周围的邻格 地面 的高度float n_height = _temp_height[tmp_j_idx];// 中格周围的邻格 无水高度 - 中格 无水高度float tmp_h_diff = n_height + tmp_n_debris - (c_height + c_debris);// 中格周围的邻格 带水高度 - 中格 带水高度float tmp_m_diff = (n_height + tmp_n_debris + tmp_n_material) - (c_height + c_debris + c_material);float tmp_diff = diff_idx == 0 ? tmp_h_diff : tmp_m_diff;float _gridbias = gridbias;_gridbias = clamp(_gridbias, -1.0f, 1.0f);// 修正高度差if (tmp_dx && tmp_dz)tmp_diff *= clamp(1.0f - _gridbias, 0.0f, 1.0f) / 1.4142136f;else // !tmp_dx || !tmp_dztmp_diff *= clamp(1.0f + _gridbias, 0.0f, 1.0f);// diff_idx = 1 的时候，前面比较过格子的总高度差，此时// 如果周边格子 不比我高，因为前面有过交换，所以至少有一个格子满足这个要求// diff_idx = 0 的时候，下面的条件不一定能满足。格子的地面有可能是最低的if (tmp_diff <= 0.0f) // 只统计比我低的邻格，所以 高度差 的说法改为 深度差{// 指定方向上，中格(我) 与 邻格 的深度差// dir_probs[0] 可能此时 >0 不会进来，此时 dir_probs[0] 保持默认值 0if ((dx_check == tmp_dx) && (dz_check == tmp_dz))dir_probs[diff_idx] = tmp_diff;// 按格子总高度计算的时候，记录 tmp_diff 最深的深度，作为 dir_probif (diff_idx && (tmp_diff < dir_prob))dir_prob = tmp_diff;// 记录比 中格 低的邻格的深度和sum_diffs[diff_idx] += tmp_diff;}}}if (diff_idx && (dir_prob > 0.001f || dir_prob < -0.001f)){// 按 (地面高度差+碎屑高度差)来计算时，流动概率 = 指定方向上的深度差 / 最大深度差dir_prob = dir_probs[diff_idx] / dir_prob;}else // add by wangbo{dir_prob = 0.0f;}// 另一种计算方法：指定方向上的流动概率 = 指定方向上的深度差 / 所有比我低的邻格的深度差之和// 这种概率显然比上一种方法的计算结果要 低// diff_idx == 1 时，深度差 以 (地面高度差+碎屑高度差) 来计算时// diff_idx == 0 时，深度差 以 (地面高度差) 来计算时，可能不存在，不过已经取默认值为 0 了if (sum_diffs[diff_idx] > 0.001f || sum_diffs[diff_idx] < -0.001f){dir_probs[diff_idx] = dir_probs[diff_idx] / sum_diffs[diff_idx];}else // add by wangbo{dir_probs[diff_idx] = 0.0f;}}// 最多可供流失的高度差float movable_mat = (m_diff < 0.0f) ? -m_diff : m_diff;// 它首先会被clamp(0,c_material)，以保证有足够的材料被移动movable_mat = clamp(movable_mat * 0.5f, 0.0f, c_material);// 以 height + debris + water 来计算float l_rat = dir_probs[1];// TODO: What is a good limit here?// 让水流继续保持足够的水量if (quant_amt > 0.001) // 默认 = 1.0movable_mat = clamp(quant_amt * ceil((movable_mat * l_rat) / quant_amt), 0.0f, c_material);elsemovable_mat *= l_rat; // 乘上概率，这样随着水量快速减少，水流很快就消失了float diff = (m_diff > 0.0f) ? movable_mat : -movable_mat;int cond = 0;if (dir_prob >= 1.0f)cond = 1;else{// Making sure all drops are movingdir_prob = dir_prob * dir_prob * dir_prob * dir_prob;unsigned int cutoff = (unsigned int)(dir_prob * 4294967295.0); // 0 ~ 1 映射到 0 ~ 4294967295.0unsigned int randval = erode_random(seed, (idx + nx * nz) * 8 + color + iterseed);cond = randval < cutoff;}// 不参与计算的格子，或者没有流动概率的格子if (!cond)diff = 0.0f;/// 下面开始计算侵蚀，河床，河堤/// 通过 h_diff 计算 沉积条件，用于产生 河床 和 河堤float slope_cont = (delta_x > 0.0f) ? (h_diff / delta_x) : 0.0f; // 斜率=对边/临边// 沉积因子 = 1 / (1 + 斜率)，斜率大的地方沉积的倾向小float kd_factor = clamp((1 / (1 + (slope_contribution_factor * slope_cont))), 0.0f, 1.0f);// 当前迭代索引iter(1~50) / 最大侵蚀迭代次数(参数面板填写)float norm_iter = clamp(((float)iter / (float)max_erodability_iteration), 0.0f, 1.0f);// 侵蚀因子：地面(height+debris)斜率，斜率贡献因子，深度差/深度差之和，初始侵蚀因子，迭代递加侵蚀因子float ks_factor = clamp((1 - (slope_contribution_factor * exp(-slope_cont))) * sqrt(dir_probs[0]) *(initial_erodability_factor + ((1.0f - initial_erodability_factor) * sqrt(norm_iter))),0.0f, 1.0f);// 中格侵蚀率float c_ks = global_erosionrate * erosionrate * erodability * ks_factor;// 邻格沉积率float n_kd = depositionrate * kd_factor;n_kd = clamp(n_kd, 0.0f, 1.0f);// 类似通过风化率计算侵蚀产生的碎屑float _removalrate = removalrate;float bedrock_density = 1.0f - _removalrate;// 通过 m_diff 可能包含水面的高度差//Kc //Kd //Ks float abs_diff = (diff < 0.0f) ? -diff : diff; // 可能包含水面的高度差// sedimentcap：泥沙容量，每单位流动水可携带的泥沙量。// 容量越大，材料在开始沉积多余沉积物之前被侵蚀的时间就越长。float sediment_limit = sedimentcap * abs_diff; // 根据泥沙容量计算的水中泥沙上限float ent_check_diff = sediment_limit - c_sediment;// sediment_limit - c_sediment > 0，// 意味着水中可以携带的泥沙上限超过了 中格 的沉积物// 这会导致更大的侵蚀，倾向挖宽河床，这主要是一个侵蚀过程// sediment_limit - c_sediment < 0，// 意味着水中可以携带的泥沙的能力降低，// 倾向于水中的泥沙向地面沉积，这是主要是一个沉积过程// for current cellif (ent_check_diff > 0.0f) // sediment_limit > c_sediment{// 中格因为侵蚀而被溶解的物质float dissolve_amt = c_ks * bed_erosionrate_factor * abs_diff;// 优先溶解碎屑层，但是碎屑层最大的量也只有 c_debrisfloat dissolved_debris = min(c_debris, dissolve_amt);// 中格碎屑被溶解后，还剩下的部分_debris[c_idx] -= dissolved_debris; // 碎屑被侵蚀// 如果中格碎屑被溶完了，还不够，就开始溶解 height 层_height[c_idx] -= (dissolve_amt - dissolved_debris); // height 被侵蚀// 沉积，数据来自上一 frame 计算的结果，50{8{}} 循环内简单重复计算// 中格的沉积物被冲走一半_sediment[c_idx] -= c_sediment / 2; // 沉积物被侵蚀// 风化后仍有剩余if (bedrock_density > 0.0f){// 被冲走的那一半沉积物 + 溶解物的沉积，这些沉积会堆积到 邻格float newsediment = c_sediment / 2 + (dissolve_amt * bedrock_density);// 假设沉积物都会堆积到邻格，如果超过最大碎屑高度if (n_sediment + newsediment > max_debris_depth){// 回滚float rollback = n_sediment + newsediment - max_debris_depth;// 回滚量不可以超过刚刚计算的沉积高度rollback = min(rollback, newsediment);// return the excess sediment_height[c_idx] += rollback / bedrock_density; // 向上修正 height 高度newsediment -= rollback; // 向下修正 沉积高度}_sediment[n_idx] += newsediment; // 邻格沉积物增加}}else // sediment_limit <= c_sediment，这主要是一个沉积过程{float c_kd = depositionrate * kd_factor; // 计算沉积系数c_kd = clamp(c_kd, 0.0f, 1.0f);{// -ent_check_diff = 高度差产生的泥沙 - 能被水流携带走的泥沙// 这些过剩的泥沙会成为 碎屑 和 沉积物// 碎屑的定义：高度差侵蚀的直接结果// 沉积的定义：泥沙被被水搬运到邻格沉降的结果_debris[c_idx] += (c_kd * -ent_check_diff); // 中格碎屑增加_sediment[c_idx] = (1 - c_kd) * -ent_check_diff; // 剩下的变成中格的沉积物n_sediment += sediment_limit; // 被带走的沉积物到了邻格_debris[n_idx] += (n_kd * n_sediment); // 邻格的碎屑增加_sediment[n_idx] = (1 - n_kd) * n_sediment; // 剩下的变成邻格的沉积物}// 河岸 河床的侵蚀，碎屑转移过程，不涉及 sedimentint b_idx = 0; // 岸的位置索引int r_idx = 0; // 河的位置索引float b_material = 0.0f; // 岸的水高度float r_material = 0.0f; // 河的水高度float b_debris = 0.0f; // 岸的碎屑高度float r_debris = 0.0f; // 河的碎屑高度float r_sediment = 0.0f; // 河的沉积高度if (is_mh_diff_same_sign) // 中格的水高，地也高{b_idx = c_idx; // 中格地高，是岸r_idx = n_idx; // 邻格地低，是河b_material = c_material; // 岸的水高度r_material = n_material; // 河的水高度b_debris = c_debris; // 岸的碎屑高度r_debris = n_debris; // 河的碎屑高度r_sediment = n_sediment; // 河的沉积高度}else // 中格 水高，地低{b_idx = n_idx; // 邻格地高，是岸r_idx = c_idx; // 中格地低，是河b_material = n_material;r_material = c_material;b_debris = n_debris;r_debris = c_debris;r_sediment = c_sediment;}// 河中每单位水的侵蚀量float erosion_per_unit_water = global_erosionrate * erosionrate * bed_erosionrate_factor * erodability * ks_factor;// 河中有水 && 岸的水量/河的水量if (r_material != 0.0f &&(b_material / r_material) < max_bank_bed_ratio &&r_sediment > (erosion_per_unit_water * max_bank_bed_ratio)){// NOTE: Increase the bank erosion to get a certain// angle faster. This would make the river cuts less// deep.float height_to_erode = global_erosionrate * erosionrate * bank_erosionrate_factor * erodability * ks_factor;float _bank_angle = bank_angle;_bank_angle = clamp(_bank_angle, 0.0f, 90.0f);float safe_diff = _bank_angle < 90.0f ? tan(_bank_angle * M_PI / 180.0) * delta_x : 1e10f;float target_height_removal = (h_diff - safe_diff) < 0.0f ? 0.0f : h_diff - safe_diff;float dissolve_amt = clamp(height_to_erode, 0.0f, target_height_removal);float dissolved_debris = min(b_debris, dissolve_amt);_debris[b_idx] -= dissolved_debris; // 岸的碎屑被侵蚀float division = 1 / (1 + safe_diff);_height[b_idx] -= (dissolve_amt - dissolved_debris); // 岸的 height 被侵蚀if (bedrock_density > 0.0f) // 有沉积{float newdebris = (1 - division) * (dissolve_amt * bedrock_density);if (b_debris + newdebris > max_debris_depth){float rollback = b_debris + newdebris - max_debris_depth;rollback = min(rollback, newdebris);// return the excess debris_height[b_idx] += rollback / bedrock_density;newdebris -= rollback;}_debris[b_idx] += newdebris; // 河岸沉积newdebris = division * (dissolve_amt * bedrock_density);if (r_debris + newdebris > max_debris_depth){float rollback = r_debris + newdebris - max_debris_depth;rollback = min(rollback, newdebris);// return the excess debris_height[b_idx] += rollback / bedrock_density;newdebris -= rollback;}_debris[r_idx] += newdebris; // 河床沉积}}}// Update the material level// 水往低处流，这很简单，麻烦的是上面，注意这里的索引号，不是中格-邻格模式，而是本格-邻格模式_material[idx] = i_material + diff; // 本格更新水的高度_material[j_idx] = j_material - diff; // 邻格更新水的高度}}}}set_output("prim_2DGrid", std::move(terrain));}
};
ZENDEFNODE(erode_tumble_material_v4,{ /* inputs: */ {"prim_2DGrid",{"ListObject", "perm"},{"ListObject", "p_dirs"},{"ListObject", "x_dirs"},{"float", "seed", "12.34"},{"int", "iterations", "40"}, // 流淌的总迭代次数{"int", "iter", "0"},{"int", "i", "0"},{"int", "openborder", "0"},{"float", "gridbias", "0.0"},{"int", "visualEnable", "0"},// 侵蚀主参数{"float", "global_erosionrate", "1.0"}, // 全局侵蚀率{"float", "erodability", "1.0"}, // 侵蚀能力{"float", "erosionrate", "0.4"}, // 侵蚀率{"float", "bank_angle", "70.0"}, // 河堤侵蚀角度// 高级参数{"float", "removalrate", "0.1"}, // 风化率/水吸收率{"float", "max_debris_depth", "5.0"}, // 碎屑最大深度// 侵蚀能力调整{"int", "max_erodability_iteration", "5"}, // 最大侵蚀能力迭代次数{"float", "initial_erodability_factor", "0.5"}, // 初始侵蚀能力因子{"float", "slope_contribution_factor", "0.8"}, // “地面斜率”对“侵蚀”和“沉积”的影响，“地面斜率大” -> 侵蚀因子大，沉积因子小// 河床参数{"float", "bed_erosionrate_factor", "1.0"}, // 河床侵蚀率因子{"float", "depositionrate", "0.01"}, // 沉积率{"float", "sedimentcap", "10.0"}, // 高度差转变为沉积物的比率 / 泥沙容量，每单位流动水可携带的泥沙量// 河堤参数{"float", "bank_erosionrate_factor", "1.0"}, // 河堤侵蚀率因子{"float", "max_bank_bed_ratio", "0.5"}, // 高于这个比值的河岸将不会在侵蚀中被视为河岸，会停止侵蚀// 河网控制{"float", "quant_amt", "0.05"}, // 流量维持率，越高河流流量越稳定}, /* outputs: */ {"prim_2DGrid",}, /* params: */ {}, /* category: */ {"erode",} });
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2.3.6 蒸发
（参考降水）
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！