A signal processing approach to fair surface design

SIGGRAPH 95 
A signal processing approach to fair surface design | Proceedings of the 22nd annual conference on Computer graphics and interactive techniques (acm.org)https://dl.acm.org/doi/10.1145/218380.218473

struct PrimAttribBlur : INode {
    void apply() override {
        auto prim = get_input<PrimitiveObject>("prim");
        auto prim_type = get_input2<std::string>("primType");

//        auto maskName = get_input2<std::string>("group");
//        if (!prim->verts.has_attr(maskName)) {
//            auto &_mask = prim->verts.add_attr<float>(maskName);
//            std::fill(_mask.begin(), _mask.end(), 1.0);
//        }
//        auto &mask = prim->verts.attr<float>(maskName);

        auto attr_name = get_input2<std::string>("attributes");
        auto attr_type = get_input2<std::string>("attributesType");

        auto useEdgeLength = get_input<NumericObject>("useEdgeLengthWeight")->get<bool>();

        auto iterations = get_input<NumericObject>("blurringIterations")->get<int>();

        auto mode = get_input2<std::string>("mode");
        auto mu = get_input<NumericObject>("stepSize")->get<float>();
        auto lambda = mu;
        if (mode == "VolumePreserving") {
            auto passband = get_input<NumericObject>("cutoffFrequency")->get<float>();
            if (passband < 1e-5) {
                float sqrt2_5 = 0.6324555320336758664;
                lambda = sqrt2_5;
                mu = -sqrt2_5;
            } else {
                // See:
                //   Gabriel Taubin. "A signal processing approach to fair surface design".
                //   SIGGRAPH 1995
                //
                // Let l be lambda, u be mu, and b be the passband frequency.
                // f(k) = (1-k*l)(1-k*u).  This function has maximum at 1/l + 1/u.
                // We want the maximum to occur at b, so we have the constraint
                //          1/l + 1/u = passband
                // so       u = l/(bl - 1)
                // We also want f(1) = -f(2).  This gives the constraint:
                //          u = (2-3l)/(3-5l)
                // Equating this equations gives:
                //          (3b-5)l^2 - 2bl + 2 = 0
                // Solve for l using the quadratic formula.  We want l>0.  Since b>0 and
                // 6b-10<0, subtract the square root of the discriminant to get a negative
                // numerator so the quotient is positive.
                auto discriminant = 4.0 * passband * passband - 24.0 * passband + 40.0;
                lambda = 2.0 * passband - sqrt(discriminant);
                lambda /= (6.0 * passband - 10.0);

                // Now solve for u.  Note that u is also the other root of the quadratic.
                mu = lambda / (passband * lambda - 1.0);
            }
        } else if (mode == "custom") {
            lambda = get_input<NumericObject>("oddStepSize")->get<float>();
            mu = get_input<NumericObject>("evenStepSize")->get<float>();
        }

        auto weightName = get_input2<std::string>("weightAttributes");
        if (!prim->verts.has_attr(weightName)) {
            auto &_weight = prim->verts.add_attr<float>(weightName);
            std::fill(_weight.begin(), _weight.end(), 1.0);
        }
        auto &weight = prim->verts.attr<float>(weightName);

        // 找临近点,假设最多 8 个临近点
        auto &neighbor_0 = prim->verts.add_attr<int>("_neighbor_0");
        auto &neighbor_1 = prim->verts.add_attr<int>("_neighbor_1");
        auto &neighbor_2 = prim->verts.add_attr<int>("_neighbor_2");
        auto &neighbor_3 = prim->verts.add_attr<int>("_neighbor_3");
        auto &neighbor_4 = prim->verts.add_attr<int>("_neighbor_4");
        auto &neighbor_5 = prim->verts.add_attr<int>("_neighbor_5");
        auto &neighbor_6 = prim->verts.add_attr<int>("_neighbor_6");
        auto &neighbor_7 = prim->verts.add_attr<int>("_neighbor_7");
        std::fill(neighbor_0.begin(), neighbor_0.end(), -1);
        std::fill(neighbor_1.begin(), neighbor_1.end(), -1);
        std::fill(neighbor_2.begin(), neighbor_2.end(), -1);
        std::fill(neighbor_3.begin(), neighbor_3.end(), -1);
        std::fill(neighbor_4.begin(), neighbor_4.end(), -1);
        std::fill(neighbor_5.begin(), neighbor_5.end(), -1);
        std::fill(neighbor_6.begin(), neighbor_6.end(), -1);
        std::fill(neighbor_7.begin(), neighbor_7.end(), -1);
        auto &edgeweight_0 = prim->verts.add_attr<float>("_edgeweight_0");
        auto &edgeweight_1 = prim->verts.add_attr<float>("_edgeweight_1");
        auto &edgeweight_2 = prim->verts.add_attr<float>("_edgeweight_2");
        auto &edgeweight_3 = prim->verts.add_attr<float>("_edgeweight_3");
        auto &edgeweight_4 = prim->verts.add_attr<float>("_edgeweight_4");
        auto &edgeweight_5 = prim->verts.add_attr<float>("_edgeweight_5");
        auto &edgeweight_6 = prim->verts.add_attr<float>("_edgeweight_6");
        auto &edgeweight_7 = prim->verts.add_attr<float>("_edgeweight_7");
        std::fill(edgeweight_0.begin(), edgeweight_0.end(), 0);
        std::fill(edgeweight_1.begin(), edgeweight_1.end(), 0);
        std::fill(edgeweight_2.begin(), edgeweight_2.end(), 0);
        std::fill(edgeweight_3.begin(), edgeweight_3.end(), 0);
        std::fill(edgeweight_4.begin(), edgeweight_4.end(), 0);
        std::fill(edgeweight_5.begin(), edgeweight_5.end(), 0);
        std::fill(edgeweight_6.begin(), edgeweight_6.end(), 0);
        std::fill(edgeweight_7.begin(), edgeweight_7.end(), 0);

        //========================================
//        LARGE_INTEGER t1_0,t2_0,tc_0;
//        LARGE_INTEGER t1_1,t2_1,tc_1;
//        LARGE_INTEGER t1_2,t2_2,tc_2;
        //========================================

        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//        QueryPerformanceFrequency(&tc_0);
//        QueryPerformanceCounter(&t1_0);

#pragma omp parallel for
        for (size_t point_idx = 0; point_idx < prim->verts.size(); point_idx++) {   // 遍历所有点,找它的邻居
            std::map<std::string, int> neighborVertID;
            std::map<std::string, float> neighborEdgeLength;
            for(int i = 0; i < 8; i++) {
                neighborVertID["neighbor_" + std::to_string(i)] = -1;
                neighborEdgeLength["edgeweight_" + std::to_string(i)] = 0;
            }

            int find_neighbor_count = 0;
            float edgeLengthSum = 0;
            volatile bool flag = false;

            if (prim_type == "line") {
#pragma omp parallel for shared(flag)
                for (size_t line_idx = 0; line_idx < prim->lines.size(); line_idx++) {
                    if(flag) continue;
                    if (prim->lines[line_idx][0] == point_idx) {
                        neighborVertID["neighbor_" + std::to_string(find_neighbor_count)] = prim->lines[line_idx][1];
                        if (useEdgeLength) {
                            float edgeLength = length(prim->verts[prim->lines[line_idx][1]] - prim->verts[point_idx]);
                            neighborEdgeLength["edgeweight_" + std::to_string(find_neighbor_count)] = edgeLength;
                            edgeLengthSum += edgeLength;
                        }
                        find_neighbor_count ++;
                    } else if (prim->lines[line_idx][1] == point_idx) {
                        neighborVertID["neighbor_" + std::to_string(find_neighbor_count)] = prim->lines[line_idx][0];
                        if (useEdgeLength) {
                            float edgeLength = length(prim->verts[prim->lines[line_idx][0]] - prim->verts[point_idx]);
                            neighborEdgeLength["edgeweight_" + std::to_string(find_neighbor_count)] = edgeLength;
                            edgeLengthSum += edgeLength;
                        }
                        find_neighbor_count++;
                    }
                    if (find_neighbor_count >= 7)
                        flag = true;
                }
            } else if (prim_type == "tri") {
                std::vector<int> pointNeighborSign(prim->verts.size());
                std::fill(pointNeighborSign.begin(), pointNeighborSign.end(), 0);

                //========================================
//                if (point_idx == 50000) {
//                    QueryPerformanceFrequency(&tc_1);
//                    QueryPerformanceCounter(&t1_1);
//                }
                //========================================

#pragma omp parallel for
                for (size_t tri_idx = 0; tri_idx < prim->tris.size(); tri_idx++) {
                    auto const &ind = prim->tris[tri_idx];
                    if (ind[0] == point_idx) {
                        pointNeighborSign[ind[1]] = 1;
                        pointNeighborSign[ind[2]] = 1;
                    } else if (ind[1] == point_idx) {
                        pointNeighborSign[ind[0]] = 1;
                        pointNeighborSign[ind[2]] = 1;
                    } else if (ind[2] == point_idx) {
                        pointNeighborSign[ind[0]] = 1;
                        pointNeighborSign[ind[1]] = 1;
                    }
                }

#pragma omp parallel for shared(flag)
                for (int i = 0; i < prim->verts.size(); i++) {
                    if(flag) continue;
                    if (pointNeighborSign[i]) {
                        neighborVertID["neighbor_" + std::to_string(find_neighbor_count)] = i;
                        if (useEdgeLength) {
                            float edgeLength = length(prim->verts[i] - prim->verts[point_idx]);
                            neighborEdgeLength["edgeweight_" + std::to_string(find_neighbor_count)] = edgeLength;
                            edgeLengthSum += edgeLength;
                        }
                        find_neighbor_count ++;
                    }
                    if (find_neighbor_count >= 7)
                        flag = true;
                }

                //========================================
//                if (point_idx == 50000) {
//                    QueryPerformanceCounter(&t2_1);
//                    double time_1 = (double)(t2_1.QuadPart - t1_1.QuadPart)/(double)tc_1.QuadPart;
//                    printf("time_1 = %f s\n", time_1);
//                }
                //========================================
            }

            neighbor_0[point_idx] = neighborVertID["neighbor_0"];
            neighbor_1[point_idx] = neighborVertID["neighbor_1"];
            neighbor_2[point_idx] = neighborVertID["neighbor_2"];
            neighbor_3[point_idx] = neighborVertID["neighbor_3"];
            neighbor_4[point_idx] = neighborVertID["neighbor_4"];
            neighbor_5[point_idx] = neighborVertID["neighbor_5"];
            neighbor_6[point_idx] = neighborVertID["neighbor_6"];
            neighbor_7[point_idx] = neighborVertID["neighbor_7"];

            if (useEdgeLength) {
                float min_length = (edgeLengthSum / (float)(find_neighbor_count)) * 0.001f;
                float sum = 0;

#pragma omp parallel for
                for (int i = 0; i < find_neighbor_count; i++)
                {
                    float length = neighborEdgeLength["edgeweight_" + std::to_string(i)];

                    if ( length > min_length )
                        neighborEdgeLength["edgeweight_" + std::to_string(i)] = 1.0 / length;
                    else    // 基本重合的点,不考虑其影响,权重打到 0
                        neighborEdgeLength["edgeweight_" + std::to_string(i)] = 0;

                    sum += neighborEdgeLength["edgeweight_" + std::to_string(i)];   // 累计总权重
                }
                if ( sum > 0 )
                {
#pragma omp parallel for
                    for (int i = 0; i < find_neighbor_count; ++i)
                    {
                        neighborEdgeLength["edgeweight_" + std::to_string(i)] /= sum;   // 权重归一化
                    }
                }

                edgeweight_0[point_idx] = neighborEdgeLength["edgeweight_0"];
                edgeweight_1[point_idx] = neighborEdgeLength["edgeweight_1"];
                edgeweight_2[point_idx] = neighborEdgeLength["edgeweight_2"];
                edgeweight_3[point_idx] = neighborEdgeLength["edgeweight_3"];
                edgeweight_4[point_idx] = neighborEdgeLength["edgeweight_4"];
                edgeweight_5[point_idx] = neighborEdgeLength["edgeweight_5"];
                edgeweight_6[point_idx] = neighborEdgeLength["edgeweight_6"];
                edgeweight_7[point_idx] = neighborEdgeLength["edgeweight_7"];
            }

        }

//        QueryPerformanceCounter(&t2_0);
//        double time_0 = (double)(t2_0.QuadPart - t1_0.QuadPart)/(double)tc_0.QuadPart;
//        printf("time_0 = %f s\n", time_0);
//        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//        QueryPerformanceFrequency(&tc_2);
//        QueryPerformanceCounter(&t1_2);

        // 平滑属性计算
        std::visit(
            [&](auto ty) {
              using T = decltype(ty);

              auto &data = prim->verts.attr<T>(attr_name);
              auto &data_temp = prim->verts.add_attr<T>("_data_temp");
              std::fill(data_temp.begin(), data_temp.end(), T(0));

              for (int loop = 0; loop < iterations; loop++) {
#pragma omp parallel for
                  // data => data_temp
                  for (size_t i = 0; i < prim->verts.size(); i++) {
                      std::vector<int> neighborIDs(8);
                      neighborIDs[0] = neighbor_0[i];
                      neighborIDs[1] = neighbor_1[i];
                      neighborIDs[2] = neighbor_2[i];
                      neighborIDs[3] = neighbor_3[i];
                      neighborIDs[4] = neighbor_4[i];
                      neighborIDs[5] = neighbor_5[i];
                      neighborIDs[6] = neighbor_6[i];
                      neighborIDs[7] = neighbor_7[i];
                      std::vector<T> neighborValues(8);
                      for (int i = 0; i < neighborIDs.size(); i++) {
                          if (neighborIDs[i] != -1)
                              neighborValues[i] = data[neighborIDs[i]];
                      }
                      std::vector<float> neighborEdgeWeights(8);
                      neighborEdgeWeights[0] = edgeweight_0[i];
                      neighborEdgeWeights[1] = edgeweight_1[i];
                      neighborEdgeWeights[2] = edgeweight_2[i];
                      neighborEdgeWeights[3] = edgeweight_3[i];
                      neighborEdgeWeights[4] = edgeweight_4[i];
                      neighborEdgeWeights[5] = edgeweight_5[i];
                      neighborEdgeWeights[6] = edgeweight_6[i];
                      neighborEdgeWeights[7] = edgeweight_7[i];
                      smooth(neighborIDs, neighborValues, neighborEdgeWeights, useEdgeLength, data[i], weight[i], lambda, data_temp[i]);
                  }
#pragma omp parallel for
                  // data_temp => data
                  for (size_t i = 0; i < prim->verts.size(); i++) {
                      std::vector<int> neighborIDs(8);
                      neighborIDs[0] = neighbor_0[i];
                      neighborIDs[1] = neighbor_1[i];
                      neighborIDs[2] = neighbor_2[i];
                      neighborIDs[3] = neighbor_3[i];
                      neighborIDs[4] = neighbor_4[i];
                      neighborIDs[5] = neighbor_5[i];
                      neighborIDs[6] = neighbor_6[i];
                      neighborIDs[7] = neighbor_7[i];
                      std::vector<T> neighborValues(8);
                      for(int i = 0; i < neighborIDs.size(); i++)
                      {
                          if (neighborIDs[i] != -1)
                              neighborValues[i] = data_temp[neighborIDs[i]];
                      }
                      std::vector<float> neighborEdgeWeights(8);
                      neighborEdgeWeights[0] = edgeweight_0[i];
                      neighborEdgeWeights[1] = edgeweight_1[i];
                      neighborEdgeWeights[2] = edgeweight_2[i];
                      neighborEdgeWeights[3] = edgeweight_3[i];
                      neighborEdgeWeights[4] = edgeweight_4[i];
                      neighborEdgeWeights[5] = edgeweight_5[i];
                      neighborEdgeWeights[6] = edgeweight_6[i];
                      neighborEdgeWeights[7] = edgeweight_7[i];
                      smooth(neighborIDs, neighborValues, neighborEdgeWeights, useEdgeLength, data_temp[i], weight[i], mu, data[i]);
                  }
              }
              prim->verts.erase_attr("_data_temp");

            },
            enum_variant<std::variant<float, vec3f>>(
                array_index({"float", "vec3f"}, attr_type)));

//        QueryPerformanceCounter(&t2_2);
//        double time_2 = (double)(t2_2.QuadPart - t1_2.QuadPart)/(double)tc_2.QuadPart;
//        printf("time_2 = %f s\n", time_2);
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        prim->verts.erase_attr("_neighbor_0");
        prim->verts.erase_attr("_neighbor_1");
        prim->verts.erase_attr("_neighbor_2");
        prim->verts.erase_attr("_neighbor_3");
        prim->verts.erase_attr("_neighbor_4");
        prim->verts.erase_attr("_neighbor_5");
        prim->verts.erase_attr("_neighbor_6");
        prim->verts.erase_attr("_neighbor_7");
        prim->verts.erase_attr("_edgeweight_0");
        prim->verts.erase_attr("_edgeweight_1");
        prim->verts.erase_attr("_edgeweight_2");
        prim->verts.erase_attr("_edgeweight_3");
        prim->verts.erase_attr("_edgeweight_4");
        prim->verts.erase_attr("_edgeweight_5");
        prim->verts.erase_attr("_edgeweight_6");
        prim->verts.erase_attr("_edgeweight_7");

        set_output("prim", std::move(prim));
    }
};
ZENDEFNODE(PrimAttribBlur,
           {/* inputs: */ {
                   "prim",
                   {"enum line tri", "primType", "tri"},
//                   {"string", "group", "mask"},
                   {"string", "attributes", "ratio"},
                   {"enum float vec3f ", "attributesType", "float"},
                   {"bool", "useEdgeLengthWeight", "false"},
                   {"int", "blurringIterations", "0"},
                   {"enum laplacian VolumePreserving custom", "mode", "laplacian"},
                   {"float", "stepSize", "0.683"},
                   {"float", "cutoffFrequency", "0.1"},
                   {"float", "evenStepSize", "0.5"},
                   {"float", "oddStepSize", "0.5"},
                   {"string", "weightAttributes", "weight"},
               }, /* outputs: */ {
                   "prim",
               }, /* params: */ {
               }, /* category: */ {
                   "primCurve",
               }});

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值