logo标签列表

CGAL 3D函数插值

3dinterpolationcgal
3
在CGAL手册中，它在这里说：链接：
散点数据插值解决了以下问题：给定一组数据点上函数的测量值，任务是在任意查询点上插值该函数。更正式地说，设P={p1,…,pn}为R2或R3中的n个点集合，Φ为定义在P凸包上的标量函数。我们假设函数值在P点已知，即对于每个pi∈P，我们关联zi=Φ(pi)。
然后它继续展示了例子，但我发现最接近3D插值的例子是Interpolation/surface_neighbor_coordinates_3.cpp。据我所知，这只能找到一个点的相邻顶点。
是否有随后的步骤可以采取来找到3D中任意查询点处的插值函数值？（特别是在一个球体上？）
自七月中旬以来，您找到了解决方案吗？ - lrineau
基本上，我的解决方案现在是：进行3D Delaunay三角剖分，将凸包转换为多面体，制作该多面体的树形结构，然后使用树形结构中的“closest_point_and_primitive”查找来查询任何点。之后，我使用该点周围的3个邻居计算自己的球面三角形面积，并根据这些子区域占据的完整球面区域的百分比加权三个数据值。3D插值。简单易行... - kdottiemo
对于 Stack Overflow 的知识库，您能否发布一个自我回答并接受它？ - lrineau
2个回答
2
0

将2D示例转换为3D确实不是一件容易的事情；我会尝试这样做

#include <stdexcept>  // std::logic_error
#include <iterator>   // back_inserter
#include <algorithm>  // copy_n
#include <functional> // function<R(A)> 
#include <map>        // map, pair
#include <cmath>      // abs()
#include <string>     // string, stoul  
#include <limits>     // numeric_limits<size_t>::max
#include <iostream>
#include <vector>


const long  MIN_SAMPLE_SIZE     = 10;
const long  DEFAULT_SAMPLE_SIZE = 1000;
const long  MAX_SAMPLE_SIZE     = 100000;

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/natural_neighbor_coordinates_3.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/surface_neighbor_coordinates_3.h>
///@TODO check if using Delaunay tesselation could make it faster
//#include <CGAL/Delaunay_triangulation_3.h>

using Kernel             = CGAL::Exact_predicates_inexact_constructions_kernel;
using real_t             = Kernel::FT;
using point3_t           = Kernel::Point_3;
using vec3_t             = Kernel::Vector_3;
using value_map_t        = std::map<point3_t, real_t, Kernel::Less_xyz_3>;
using value_access_t     = CGAL::Data_access<value_map_t>;
using distance_vec_t     = std::vector<std::pair<point3_t, real_t>>;


using ScalarR3Field      = std::function<real_t(point3_t)>;

const ScalarR3Field maxN = [] (const point3_t& p) {
    std::vector<real_t> v {
        std::abs(p.x()),
        std::abs(p.y()),
        std::abs(p.z())
    };
    return *std::max_element(v.begin(),v.end());
    //return std::max(v.begin(), v.end());
};

const ScalarR3Field l1N = [] (const point3_t& p) {
    return std::abs(p.x()) +
        std::abs(p.y()) +
        std::abs(p.z());
};


std::map<std::string, ScalarR3Field> TESTS {
    {"maximum norm", maxN},
    {"l_1 norm",  l1N}
};

int main(int argc, char* argv[]) {
    auto sample_size = DEFAULT_SAMPLE_SIZE; 

    if (argc > 1) {
        try {
            auto val = std::stol(argv[1]);
            if (val < MIN_SAMPLE_SIZE) {
                throw std::logic_error("sample size suppplied is too small");
            } else if (val > MAX_SAMPLE_SIZE) {
                throw std::logic_error("sample size suppplied is too large");
            }
            sample_size = val;
        } catch (const std::logic_error& e) {
            std::cerr << e.what() << ";"
                      << " using default sample size of "
                      << (sample_size = DEFAULT_SAMPLE_SIZE) 
                      << "." << std::endl;  
        } // try/catch for argv[1] conversion 
    } // if (argc > 1) 

    std::vector<point3_t> points;
    points.reserve(sample_size);

    CGAL::Random_points_on_sphere_3<point3_t> gen (1.0);
    std::copy_n(gen, sample_size, std::back_inserter(points));

    auto point  = *(++gen);
    auto vec    = point - CGAL::ORIGIN;

    std::cout << std::endl;

    for (auto& test : TESTS) {
        value_map_t F;
        auto& f = test.second;

        for (auto& p: points)
            F.emplace(std::make_pair(p, f(p)));

        distance_vec_t coords;

        auto triple =  CGAL::surface_neighbor_coordinates_3(
                           points.begin(), points.end(),
                           point, vec,
                           std::back_inserter(coords),
                           Kernel()
                       );

        auto intp   =  CGAL::linear_interpolation(
                                coords.begin(), coords.end(),
                                triple.second,
                                value_access_t(F)
                       );

        auto nomn   = f(point); 

        std::cout   << test.first << "([" << point << "])"
                    << "\tnominal: "      << nomn
                    << "\tinterpolated: " << intp 
                    << "\t\tΔ=" << (intp-nomn)*real_t(100.0)/nomn 
                    << "%\n";
    } // for (auto& test : TESTS)

    std::cout << std::endl;
    return EXIT_SUCCESS;
} // int main

方式；看起来并不太遥远：

./cgal_test 10000

l_1 norm([-0.689748 -0.39429 -0.607274])    nominal: 1.69131    interpolated: 1.69079       Δ=-0.0309701%
maximum norm([-0.689748 -0.39429 -0.607274])    nominal: 0.689748   interpolated: 0.689535      Δ=-0.0309701%
非常感谢您的建议！我会看一下的，之前因为时间紧迫不得不放弃这个问题，但现在会重新关注它。 - kdottiemo
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接