dse_convex_hull.hpp

Gramham法求解凸包。从最基本数据结构定义开始实现

#ifndef __DSE_CONVEX_HULL_H__
#define __DSE_CONVEX_HULL_H__

#include "dse_classes.hpp"
#include <cmath>

namespace dse
{

    //用于凸包的平面点，包含相对极坐标角度量 r
    //##ModelId=443378C7032C
    class point2d
    {
    public:
        double x,y,r;
        //##ModelId=443378C70336
        point2d(const point2d &other)
        {
            x = other.x;
            y = other.y;
            r = other.r;
        }
        
        //##ModelId=443378C70338
        point2d()
        {
            x = y = r = 0;
        }
        //##ModelId=443378C70339
        point2d(int x, int y)
        {
            this->x = x;
            this->y = y;
            r = 0;
        }
        //##ModelId=443378C70341
        point2d(double x, double y)
        {
            this->x = x;
            this->y = y;
            r = 0;
        }

        //计算两个点之间的距离
        //##ModelId=443378C70344
        static double distance(point2d p, point2d q)
        {
            double x1 = q.x - p.x, y1 = q.y - p.y;
            return sqrt(x1*x1 + y1*y1);
        }

    };


    //凸包类
    //##ModelId=443378C70368
    class convex_hull
    {
        //##ModelId=443378C8000D
        vector<point2d> points_;
        //##ModelId=443378C80022
        vector<int> convex_points_index_;
        //##ModelId=443378C8002A
        size_t last_i_;
        //##ModelId=443378C80049
        bool finished_;

        //用于极坐标角度排序的比较的仿函数
        //##ModelId=443378C8008F
        class convex_points_less
        {
                //##ModelId=443378C80091
            const point2d base;
        public:
            //记录基准点
                //##ModelId=443378C8009B
            convex_points_less(const point2d &b)
                : base(b)
            {
            }

            //比较操作，规则如下：
            //如果极坐标角度相等，则让距离基准点远的点靠前
            //如果极坐标角度不相等，则让极坐标角度小的在前
                //##ModelId=443378C800A3
            bool operator()(const point2d &a, const point2d &b)
            {
                if(IS_ZERO(a.r - b.r))
                    return point2d::distance(a, base) > point2d::distance(b, base);

                return a.r < b.r;
            }
        };

        
    public:
        //##ModelId=443378C80053
        convex_hull()
        {
            finished_ = false;
            last_i_ = 0;
        }

        //清除所有数据
        //##ModelId=443378C80054
        void clear()
        {
            finished_ = false;
            last_i_ = 0;
            points_.clear();
            convex_points_index_.clear();
        }

        //取得点集中的第i个点
        //##ModelId=443378C80055
        point2d &get_point(size_t i)
        {
            return points_[i];
        }

        //取得整个点集的所有点的数目
        //##ModelId=443378C8005D
        size_t get_points_count()
        {
            return points_.size();
        }

        //取得凸包的第i个点
        //##ModelId=443378C8005E
        point2d &get_convex_point(size_t i)
        {
            return points_[convex_points_index_[i]];
        }

        //获取凸包顶点数
        //##ModelId=443378C80060
        size_t get_convex_points_count()
        {
            return convex_points_index_.size();
        }

        //向待求解的点集中添加一个点
        //##ModelId=443378C80061
        void add_point(double x, double y)
        {
            points_.push_back(point2d(x,y));

            //清除已经完成的标记
            finished_ = false;  
            convex_points_index_.clear();
        }

        //从待求解的点集中删除一个点
        //##ModelId=443378C80069
        void delete_point(size_t i)
        {
            points_.erase(i);

            //清除已经完成的标记
            finished_ = false;  
            convex_points_index_.clear();
        }

        //寻找基准点。原则：x 尽量小，或者 x 相等的时候 y 尽量小
        //##ModelId=443378C8006B
        point2d &find_base_point()
        {
            size_t idx = 0;
            for(size_t i = 0; i < points_.size(); i++)
                if ((points_[i].x < points_[idx].x) || (points_[i].x == points_[idx].x && points_[i].y < points_[idx].y))
                    idx = i;

            return points_[idx];
        }

        //对其他点相对基准点的极坐标角度排序
        //##ModelId=443378C80071
        void adjust_points_order(point2d &base)
        {
            //计算其他点相对基准点的极坐标角度（ r >= 0 ）
            for(size_t i = 0; i < points_.size(); i++)
            {
                double t = (base.y - points_[i].y) / point2d::distance(base, points_[i]);
                if(IS_ZERO(t)) points_[i].r = PI/2;
                else if(IS_ZERO(t-1)) points_[i].r = 0;
                else if(IS_ZERO(t+1)) points_[i].r = PI;
                else points_[i].r = acos(t);
                //DPRINTF((" r = %lf\n"), points_[i].r);
            }

            //保证基准点排序后处于第一个位置
            base.r = -1;

            //以基准点对所有点进行排序
            quick_sort(points_, 0, points_.size() - 1, convex_points_less(base));

            /*
            for(i = 0; i < points_.size(); i++)
                DPRINTF(" after sort: %lf,%lf\n", points_[i].x, points_[i].y);
            */

        }

        //进行准备工作
        //1、清理环境
        //2、寻找基本点
        //3、对其他点相对基本点的极坐标角度排序
        //4、把基准点放到凸包点集中
        //##ModelId=443378C80073
        void prepare()
        {
            convex_points_index_.clear();
            last_i_ = 0;
            finished_ = false;

            adjust_points_order(find_base_point());

            if(points_.size() != 0)
                convex_points_index_.push_back(0);
        }

        //求解凸包的工作是否完成
        //##ModelId=443378C8007B
        bool finished()
        {
            return finished_;
        }

        //使用向量积判断 p, q, r是否能构成凸包的边界。
        //##ModelId=443378C8007C
        bool is_valid_convex_point(point2d &p, point2d &q, point2d &r)
        {
            double x1 = q.x - p.x, y1 = q.y - p.y;
            double x2 = r.x - q.x, y2 = r.y - q.y;
            return (x1*y2 - x2*y1 >= 0) ;
        }

        //是否还有下一个点需要处理
        //##ModelId=443378C80080
        bool has_next_point()
        {
            return last_i_ + 1 < points_.size();
        }

        //返回待处理的下一个点的索引
        //##ModelId=443378C80085
        int get_next_point_index()
        {
            return last_i_ + 1;
        }

        //进行一步操作，寻找下一个合适的点
        //##ModelId=443378C80086
        bool search_next_point()
        {

            //如果没有点需要处理，则结束整个求解过程
            if(!has_next_point())
            {
                finished_ = true;
                return false;
            }


            //如果目前求得的凸包顶点小于2个，则直接把下一个点添加到凸包顶点集合中
            //TODO: 这一步在很多情况下可以优化掉。不过还需要增加其他判断
            //      有空的时候仔细证明一下，看看代码能否继续优化
            if(convex_points_index_.size() < 2)
            {
                convex_points_index_.push_back(++last_i_);
                return true;
            }

            //取出当前求得的凸包顶点集合中的最后两个点
            point2d &p1 = points_[convex_points_index_[convex_points_index_.size() - 2]];
            point2d &p2 = points_[convex_points_index_[convex_points_index_.size() - 1]];
            //取出要处理的下一个点
            point2d &p3 = points_[last_i_ + 1];
        

            //如果下一个点与凸包点集中最后一点的极坐标相同，
            //则可以忽略（因为排序的时候已经保证距离基准点最远的点排在最前面，
            //因此当前点必然在凸包内）
            //可以处理下一个点
            if(IS_ZERO(p2.r - p3.r))
            {
                last_i_++;
                return true;
            }

            //如果 p1 p2 p3 不能构成凸包的边界
            //则把已经求得的凸包点集合的最后一个点去掉（因为它必然不是凸包的顶点）
            if(!is_valid_convex_point(p1, p2, p3))
            {
               convex_points_index_.pop_back();
               return true;
            }
            else  //如果 p1 p2 p3 能构成凸包的边界则把p3添加到凸包点集合，并继续找下一个
            {
                convex_points_index_.push_back(++last_i_);
                return true;
            }

        }
        
        //求解凸包
        //##ModelId=443378C80087
        void solve()
        {
            //准备; 求解直道结束
            for(prepare();search_next_point();)
                ;
        }
    };
}
#endif//__DSE_CONVEX_HULL_H__