tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

点の集合を包含する球

点の集合を包含する球のアルゴリズムを参照した。

#include "structure.hpp"
#include "mesh.hpp"

// https://www.ipsj.or.jp/07editj/promenade/4309.pdf
Vector find_center(const Mesh& mesh)
{
    auto norm = [](const Vector& pos, const Vertex& vert) -> real_t
    {
        const real_t dx{vert.x() - pos.x};
        const real_t dy{vert.y() - pos.y};
        const real_t dz{vert.z() - pos.z};
        return dx * dx + dy * dy + dz * dz;
    };

    const Vertex start{mesh.vertex[0]};
    // vector to represent currently estimated sphere center
    Vector p{start.x(), start.y(), start.z()};

    real_t move{0.5};
    while (move > 1e-8)
    {
        for (int t = 0; t < 50; ++t)
        {
            // index to farthest vertex from p
            std::size_t k{0};
            real_t max{0};
            for (std::size_t i = 0; i < mesh.vertex.size(); ++i)
            {
                if (norm(p, mesh.vertex[i]) > max)
                {
                    max = norm(p, mesh.vertex[i]);
                    k = i;
                }
            }
            p.x += (mesh.vertex[k].x() - p.x) * move;
            p.y += (mesh.vertex[k].y() - p.y) * move;
            p.z += (mesh.vertex[k].z() - p.z) * move;
        }
        move *= 0.5;
    }

    return p;
}
#define BOOST_TEST_MAIN
#include <boost/test/included/unit_test.hpp>
#include <boost/test/tools/detail/tolerance_manip.hpp>
#include <iostream>

#include "find_center.hpp"
#include "mesh.hpp"
#include "structure.hpp"

using namespace cnthd;

BOOST_AUTO_TEST_CASE(test_find_center)
{
    {
        // Mesh mesh = read_obj("./cnthd/test/data/cube_small.obj");
        Mesh mesh = read_obj("./test/data/cube_small.obj");
        // now center of the mesh is (0, 0, 0)
        mesh.move(Vector{-0.5, -0.5, -0.5});

        Vector center{find_center(mesh)};
        std::cout << center << std::endl;

        BOOST_TEST(center.x == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.y == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.z == 0., boost::test_tools::tolerance(1e-8));
    }

    {
        // Mesh mesh = read_obj("./cnthd/test/data/sphere.obj");
        Mesh mesh = read_obj("./test/data/sphere.obj");

        Vector center{find_center(mesh)};
        std::cout << center << std::endl;

        BOOST_TEST(center.x == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.y == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.z == 0., boost::test_tools::tolerance(1e-8));
    }

    {
        // Mesh mesh = read_obj("./cnthd/test/data/sphere.obj");
        Mesh mesh = read_obj("./test/data/DDGSpring2016/sphere_large.obj");

        Vector center{find_center(mesh)};
        std::cout << center << std::endl;

        BOOST_TEST(center.x == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.y == 0., boost::test_tools::tolerance(1e-8));
        BOOST_TEST(center.z == 0., boost::test_tools::tolerance(1e-8));
    }

    {
        // Mesh mesh = read_obj("./cnthd/test/data/sphere.obj");
        Mesh mesh = read_obj("./test/data/cell00038.obj20.obj");

        Vector center{find_center(mesh)};
        std::cout << "circumcenter: " << center << std::endl;

        // the center of mass
        Vector g{0,0,0};
        const std::size_t N{mesh.vertex.size()};
        for (const Vertex& v: mesh.vertex)
        {
            g.x += v.x() / N;
            g.y += v.y() / N;
            g.z += v.z() / N;
        }
        std::cout << "center of mass: " << g << std::endl;

    }
}

以下はmeshlabで作図してみたもの。赤は外接円の中心、青はメッシュの頂点の重心。入力がほぼ球面に近いことは分かっているので入力の大きさに合わせてループの回数を減らしても問題ないとは思われる。