Clever Geek Handbook

Construction de coque 3D convexe

Quelle? Pourquoi?



salut!



Je voudrais considérer un problème de géométrie de calcul, à savoir la construction d'une coque 3D convexe . Il me semble que ce n'est ni l'algorithme le plus compliqué ni le plus simple qui serait très intéressant et utile à analyser.



Si vous n'avez jamais fait face à une telle tâche, je pense qu'il sera intéressant pour vous d'en apprendre davantage, de voir ce que c'est.



Si vous venez d'entendre quelque chose sur les coques convexes, vous pouvez en savoir plus à leur sujet.



Si vous êtes un gourou des coques convexes, vous aimerez peut-être écouter à nouveau la solution à un problème intéressant.







Contenu



  • Quelle? Pourquoi?
  • Qu'est-ce qu'une coque convexe?
  • Mots communs
  • Description de l'algorithme
  • Asymptotiques
  • Implémentation d'algorithme
  • Mise en œuvre complète




J'exprime ma gratitude muji-4ok pour obtenir de l'aide pour rédiger et éditer l'article.



Qu'est-ce qu'une coque convexe?



La coque convexe d'un ensemble X est le plus petit ensemble convexe contenant l'ensemble X.



, . :

, , , , , . .





2D

,



, , , , , , "" , .





, , , 3D . , .



: , ?



— .



— . , , : .



, , .



, ( ) " ".



, 4 . 3 .





, , .



, , 3d , . , .



1.

, , , , .

, .



Z. , , "" Z . XY , . . , . P, , "" — f.



, . Q, f ( prQ). , () {P, Q} {Q, prQ} — — . , Q , . Q .





. R, , P, Q R (0, 0, 1) . , f — XY. , , . , , R .



, , , , ( - ).



, P, Q, R .



-! , .



2.

, , , . , .



: , . , , , , () ( ), , . , , , , .



: , . : .



1. ? , . . E. , , E — R. ( E) , E R. , , , — .



2. , , " ", E , , E .



3. , , . , , . , , , .



4. , . , . , "" , .



5. , "" , .



6. , . , .



, ( ), 3D .





, , . , . N, H , O(NH). H 2N, O(N^2).



, , , O(NlogN), , . 2D , O(NlogN), , , , .



, , : .





, , , , , . C++.



. . , , -. - . ( , ).



struct Point
{
    coordinate x;
    coordinate y;
    coordinate z;
    Point (coordinate x = 0, coordinate y = 0, coordinate z = 0) : x(x), y(y), z(z) {}
    Point operator- (const Point& other) const;
    Point operator+ (const Point& other) const;
    bool operator!= (const Point& other) const;
    bool operator== (const Point& other) const;
};


. 22 ( ).



coordinate Det2x2(coordinate a11, coordinate a12, coordinate a21, coordinate a22)
{
    return a11 * a22 - a12 * a21;
}


( A — — - - AB AC):



Point VectorProduct(const Point& A, const Point& B, const Point& C)
{
    Point a = B - A;
    Point b = C - A;
    return Point (Det2x2(a.y, a.z, b.y, b.z),
                  Det2x2(a.x, a.z, b.x, b.z),
                  Det2x2(a.x, a.y, b.x, b.y));
}


( ):



double GetLenghtVector(Point A, Point B = Point(0, 0, 0))
{
    Point vec = B - A;
    double lenght = std::sqrt(vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
    return lenght;
}


:



double GetAngle(const Point& n1, const Point& n2)
{
    double len_n1 = GetLenghtVector(n1);
    double len_n2 = GetLenghtVector(n2);
    double scalar_prod = n1.x * n2.x + n1.y * n2.y + n1.z * n2.z;
    if (scalar_prod == 0)
    {
        return 0;
    }
    return std::acos((scalar_prod) / (len_n1 * len_n2));
}


, .



. , . , , , , . , , .



struct Edge
{
    int first;
    int second;
    int flatness; // 
    bool is_close = false; //    ,     ?
    Edge(int first, int second, int flatness = -1, Point normal = Point(0, 0 , 0)) :
                    first(first), second(second), flatness(flatness) {}
};


Flatness — . 3 , ( ). — Another, . ( — if — .)



struct Flatness
{
    int first;
    int second;
    int third;
    Point normal; //   ,    
    Flatness(int first, int second, int third, Point normal) :
        first(first), second(second), third(third), normal(normal) {}
    int Another(int one, int two);
};


Class . 



class ConvexHull
{
    struct Flatness;
    struct Edge;

    std::vector<Point> points_; //   
    std::vector<Flatness> verge_; // 
    std::vector<Edge> edges_; // 

    int count_; //    

    int findMinZ() const;
    void findFirstFlatness();
    int returnIsEdgeInHull(int a, int b) const;
    void makeHull();
public:
    ConvexHull(const std::vector<Point>& points): points_(points), count_(points.size()) { makeHull(); }
};




: Z ( , )



int ConvexHull::findMinZ() const
{
    int min_id = 0;
    for (int i = 1; i < count_; ++i)
    {
        if (points_[i].z < points_[min_id].z ||
            (points_[i].z == points_[min_id].z && points_[i].y < points_[min_id].y) ||
            (points_[i].z == points_[min_id].z && points_[i].y == points_[min_id].y &&
                points_[i].x < points_[min_id].x))
        {
            min_id = i;
        }
    }
    return min_id;
}


, :



int ConvexHull::returnIsEdgeInHull(int a, int b) const
{
    for (int i = 0; i < edges_.size(); ++i)
    {
        if ((edges_[i].first == a && edges_[i].second == b) ||
            (edges_[i].first == b && edges_[i].second == a))
        {
            return i;
        }
    }

    return -1;
}


. (-)



. , : . - , , .



:



void ConvexHull::findFirstFlatness()
{
    int first_point, second_point, third_point; //    
    first_point = findMinZ();


, Z. "".



    double min_angle = 7; //         2pi,      7
    int min_id = -1;
    for (int i = 0; i < count_; ++i)
    {
        if (first_point == i) //        
        {
            continue;
        }
        Point start = points_[first_point];
        Point next = points_[i];
        double angle = GetAngle(start - next, next - Point(next.x, next.y, start.z));
        if (min_angle > angle)
        {
            min_angle = angle;
            min_id = i;
        }
    }
    second_point = min_id;


, .



    min_angle = 7;
    min_id = -1;
    for (int i = 0; i < count_; ++i)
    {
        if (first_point == i || second_point == i)
        {
            continue;
        }
        Point normal = VectorProduct(points_[first_point], points_[second_point], points_[i]);
        double angle = GetAngle(Point(0, 0, 1), normal);
        if (min_angle > angle)
        {
            min_angle = angle;
            min_id = i;
        }
    }
    third_point = min_id;


, , XY, (0, 0, 1). .



( ), .



    if (VectorProduct(points_[first_point], points_[second_point], points_[third_point]).z > 0)
    {
        std::swap (second_point, third_point);
    }
    Point new_normal = VectorProduct(points_[first_point], points_[second_point], points_[third_point]);
    verge_.push_back(Flatness(first_point, second_point, third_point, new_normal)); //  
    edges_.push_back(Edge(first_point, second_point, 0));
    edges_.push_back(Edge(second_point, third_point, 0));
    edges_.push_back(Edge(third_point, first_point, 0));
}


:

:



void ConvexHull::makeHull()
{
    findFirstFlatness();
    std::stack<int> stack;
    stack.push(0);
    stack.push(1);
    stack.push(2);


: , . , : .



    while (!stack.empty())
    {
        Point new_normal;
        Edge e = edges_[stack.top()]; //   ,  
        stack.pop();
        if (e.is_close) //       ,    
        {
            continue;
        }
        int min_id = -1;
        double min_angle = 7;

        for (int i = 0; i < count_; ++i)
        {
            int another = verge_[e.flatness].Another(e.first, e.second);
            if (i != another && e.first != i && e.second != i) // ,       
            {
                //    ,    i- 
                Point current_normal = VectorProduct(points_[e.second], points_[e.first], points_[i]);
                double angle = GetAngle(current_normal, verge_[e.flatness].normal);

                if (min_angle > angle)
                {
                    min_angle = angle;
                    min_id = i;
                    new_normal = current_normal;
                }
            }
        }


, , e , . , , is_closed = true, , , , — — , .



        if (min_id != -1) //    -     4 
        {
            e.is_close = true; //     ,        
            int count_flatness = verge_.size(); //   
            int first_edge_in_hull = returnIsEdgeInHull(e.first, min_id); //  -1,    
            int second_edge_in_hull = returnIsEdgeInHull(e.second, min_id);

            if (first_edge_in_hull == -1)
            {
                edges_.push_back(Edge(e.first, min_id, count_flatness));
                stack.push(edges_.size() - 1);
            }
            if (second_edge_in_hull == -1)
            {
                edges_.push_back(Edge(min_id, e.second, count_flatness));
                stack.push(edges_.size() - 1);
            }
            if (first_edge_in_hull != -1)
            {
                edges_[first_edge_in_hull].is_close = true;
            }
            if (second_edge_in_hull != -1)
            {
                edges_[second_edge_in_hull].is_close = true;
            }

            verge_.push_back(Flatness(e.first, e.second, min_id, new_normal));
        }

    } //  while
} //


: 
#include <iostream>
#include <vector>
#include <cmath>
#include <stack>
#include <iomanip>

using coordinate = int64_t;

struct Point;
coordinate Det2x2(coordinate a11, coordinate a12, coordinate a21, coordinate a22);
Point VectorProduct(const Point& A, const Point& B, const Point& C);
double GetLenghtVector(Point A, Point B);
double GetAngle(const Point& n1, const Point& n2);

struct Point
{
    coordinate x;
    coordinate y;
    coordinate z;
    Point(coordinate x = 0, coordinate y = 0, coordinate z = 0) : x(x), y(y), z(z) {}
    Point operator-(const Point& other) const
    {
        return Point(x - other.x, y - other.y, z - other.z);
    }
    Point operator+(const Point& other) const
    {
        return Point(x + other.x, y + other.y, z + other.z);
    }
    bool operator!= (const Point& other) const
    {
        return (x != other.x || y != other.y || z != other.z);
    }
    bool operator== (const Point& other) const
    {
        return (x == other.x && y == other.y && z == other.z);
    }
};

coordinate Det2x2(coordinate a11, coordinate a12, coordinate a21, coordinate a22)
{
    return a11 * a22 - a12 * a21;
}

//[AB, AC]
Point VectorProduct(const Point& A, const Point& B, const Point& C)
{
    Point a = B - A;
    Point b = C - A;
    return Point (Det2x2(a.y, a.z, b.y, b.z),
                  Det2x2(a.x, a.z, b.x, b.z),
                  Det2x2(a.x, a.y, b.x, b.y));
}

//vector AB
double GetLenghtVector(Point A, Point B = Point(0, 0, 0))
{
    Point vec = B - A;
    double lenght = std::sqrt(vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
    return lenght;
}

double GetAngle(const Point& n1, const Point& n2)
{
    double len_n1 = GetLenghtVector(n1);
    double len_n2 = GetLenghtVector(n2);
    double scalar_prod = n1.x * n2.x + n1.y * n2.y + n1.z * n2.z;
    if (scalar_prod == 0)
    {
        return 0;
    }
    return std::acos((scalar_prod) / (len_n1 * len_n2));
}

class ConvexHull
{
    struct Flatness
    {
        int first;
        int second;
        int third;
        Point normal; //   ,    
        Flatness(int first, int second, int third, Point normal) :
            first(first), second(second), third(third), normal(normal) {}
        int Another(int one, int two)
        {
            if ((one == first && two == second) || (one == second && two == first))
            {
                return third;
            }
            if ((one == first && two == third) || (one == third && two == first))
            {
                return second;
            }
            if ((one == third && two == second) || (one == second && two == third))
            {
                return first;
            }
            return -1; // error
        }
    };

    struct Edge
    {
        int first;
        int second;
        int flatness; // 
        bool is_close = false;
        Edge(int first, int second, int flatness = -1, Point normal = Point(0, 0 , 0)):
                    first(first), second(second), flatness(flatness) {}
    };

    std::vector<Point> points_;
    std::vector<Flatness> verge_;
    std::vector<Edge> edges_;

    int count_; //    

    int findMinZ() const;
    void findFirstFlatness();
    int returnIsEdgeInHull(int a, int b) const;
    void makeHull();
public:
    ConvexHull(const std::vector<Point>& points): points_(points), count_(points.size()) { makeHull();}
};

void ConvexHull::makeHull()
{
    findFirstFlatness();
    std::stack<int> stack;
    stack.push(0);
    stack.push(1);
    stack.push(2);

    while (!stack.empty())
    {
        Point new_normal;
        Edge e = edges_[stack.top()]; //   ,  
        stack.pop();
        if (e.is_close) //       ,    
        {
            continue;
        }
        int min_id = -1;
        double min_angle = 7;

        for (int i = 0; i < count_; ++i)
        {
            int another = verge_[e.flatness].Another(e.first, e.second);
            if (i != another && e.first != i && e.second != i) // ,       
            {
                //    ,    i- 
                Point current_normal = VectorProduct(points_[e.second], points_[e.first], points_[i]);
                double angle = GetAngle(current_normal, verge_[e.flatness].normal);

                if (min_angle > angle)
                {
                    min_angle = angle;
                    min_id = i;
                    new_normal = current_normal;
                }
            }
        }

        if (min_id != -1) //    -     4 
        {
            e.is_close = true; //     ,        
            int count_flatness = verge_.size(); //   
            int first_edge_in_hull = returnIsEdgeInHull(e.first, min_id);
            int second_edge_in_hull = returnIsEdgeInHull(e.second, min_id);

            if (first_edge_in_hull == -1)
            {
                edges_.push_back(Edge(e.first, min_id, count_flatness));
                stack.push(edges_.size() - 1);
            }
            if (second_edge_in_hull == -1)
            {
                edges_.push_back(Edge(min_id, e.second, count_flatness));
                stack.push(edges_.size() - 1);
            }
            if (first_edge_in_hull != -1)
            {
                edges_[first_edge_in_hull].is_close = true;
            }
            if (second_edge_in_hull != -1)
            {
                edges_[second_edge_in_hull].is_close = true;
            }

            verge_.push_back(Flatness(e.first, e.second, min_id, new_normal));
        }

    } //  while
} //  

int ConvexHull::findMinZ() const
{
    int min_id = 0;
    for (int i = 1; i < count_; ++i)
    {
        if (points_[i].z < points_[min_id].z ||
            (points_[i].z == points_[min_id].z && points_[i].y < points_[min_id].y) ||
            (points_[i].z == points_[min_id].z && points_[i].y == points_[min_id].y &&
                points_[i].x < points_[min_id].x))
        {
            min_id = i;
        }
    }
    return min_id;
}

void ConvexHull::findFirstFlatness()
{
    int first_point, second_point, third_point;
    first_point = findMinZ();
    double min_angle = 7;
    int min_id = -1;
    for (int i = 0; i < count_; ++i)
    {
        if (first_point == i)
        {
            continue;
        }
        Point start = points_[first_point];
        Point next = points_[i];
        double angle = GetAngle(start - next, next - Point(next.x, next.y, start.z));
        if (min_angle > angle)
        {
            min_angle = angle;
            min_id = i;
        }
    }
    second_point = min_id;

    min_angle = 7;
    min_id = -1;
    for (int i = 0; i < count_; ++i)
    {
        if (first_point == i || second_point == i)
        {
            continue;
        }
        Point normal = VectorProduct(points_[first_point], points_[second_point], points_[i]);
        double angle = GetAngle(Point(0, 0, 1), normal);
        if (min_angle > angle)
        {
            min_angle = angle;
            min_id = i;
        }
    }
    third_point = min_id;

    //  
    if (VectorProduct(points_[first_point], points_[second_point], points_[third_point]).z > 0)
    {
        std::swap (second_point, third_point);
    }
    Point new_normal = VectorProduct(points_[first_point], points_[second_point], points_[third_point]);
    verge_.push_back(Flatness(first_point, second_point, third_point, new_normal)); //  
    edges_.push_back(Edge(first_point, second_point, 0));
    edges_.push_back(Edge(second_point, third_point, 0));
    edges_.push_back(Edge(third_point, first_point, 0));
}

int ConvexHull::returnIsEdgeInHull(int a, int b) const
{
    for (int i = 0; i < edges_.size(); ++i)
    {
        if ((edges_[i].first == a && edges_[i].second == b) ||
            (edges_[i].first == b && edges_[i].second == a))
        {
            return i;
        }
    }

    return -1;
}


.



More articles:


All Articles