【c++】ObjectARX 最近点对及凸包算法代码源码

admin1

#include <cmath>
#include <cfloat>
#include <algorithm>
#include <vector>

/* 二维坐标定义 */
struct Point
{
        double x;
        double y;
};

/* 坐标数组 */
typedef std::vector<Point> PointArray;

#define MAXN 10000
Point points[MAXN];     // 顶点数组
int S[MAXN];            // graham算法使用
int D[MAXN*2];          // melkman算法使用
int n;                  // 顶点的个数

/* 计算两点之间距离 */
double dist(Point& a,Point& b)
{
        double x = a.x - b.x;
        double y = a.y - b.y;
        return sqrt(x*x + y*y);
}

/* 计算两点之间距离的平方和 */
double dist_2(Point& a,Point& b)
{
        double x = a.x - b.x;
        double y = a.y - b.y;
        return (x*x + y*y);
}

/* 即a-c是否在a-b的左边 */
/* a-b是堆栈S中的两个点，c则是需要判断的点 */
int isleft(Point& a,Point& b,Point& c)
{
        return (b.x-a.x)*(c.y-a.y) - (c.x-a.x)*(b.y-a.y);       
}

// 浮点数比较
int cmp_real(double x,double y)
{
        double d = x - y;
        if(d < 0) return 1; // d < 0
        //if(d > DBL_EPSILON) return 1; // d > 0
        return 0;   // d == 0
}

int cmp_y(int i,int j)
{
        return cmp_real(points[i].y, points[j].y);
}

// 按照极角的大小递增排序
// 如果角度相同，那么取距离p0最远的点
int cmp_angle(Point& a,Point& b)
{
        // 如果大于0，表示0-b在0-a的右边
        // 如果小于0，表示0-b在0-a的左边
        if(isleft(points[0], a, b) > 0) return 1;
        if(isleft(points[0], a, b) < 0) return 0;
        if(dist(points[0], a) < dist(points[0], b)) return 1;
        return 0;
}

// 按照先x后y排序
int cmp_xy(Point& a,Point& b)
{
        double dx = a.x - b.x;
        double dy = a.y - b.y;
        if(dx < 0) return 1;
        if(dx < DBL_EPSILON  && dy < 0) return 1; // dx == 0
        return 0;
}

//找第一个顶点,做为算法的起始顶点
int getStartPoint()
{
        int k = 0;
        for(int i = 0; i < n; ++i)
        {
                double xi = points[i].x;
                double yi = points[i].y;

                double xk = points[k].x;
                double yk = points[k].y;

                if((yi < yk) || ((yi == yk) && (xi < xk)))
                {
                        k=i;
                }
        }
        return k;
}

/* 采用Graham算法计算凸包，并返回在栈S中的位置:[0, top] */
void graham(int& top)
{
        //找第一个顶点,做为算法的起始顶点
        int k = getStartPoint();
        Point tmp = points[0];
        points[0] = points[k];
        points[k] = tmp;

        acutPrintf(_T("\n第1点坐标(%.3f, %.3f)"), points[0].x, points[0].y);

        // 将1-n-1的点按照极角[0, PI]进行排序
        std::sort(points + 1, points + n, cmp_angle);

        for(int i = 0; i <= 2; ++i)
        {
                S[i] = i;   //将前3个顶点压栈
        }

        top = 2;
        for(int i=3; i<n; i++)
        {
                /*如果新的点,与最近入栈中的2点构成了一个"凹"角, 则将栈顶元素出栈. 直到把栈                  检查完*/
                while(top > 0 && isleft(points[S[top-1]], points[S[top]], points[i])<0)               
                {
                        top--;
                }

                //acutPrintf(_T("top: %d"), top);
                top++;        // 将新点压栈
                S[top] = i;
        }
}

/* 采用melkman算法计算凸包,并返回在数组D中的区间位置(bot, top) */
void melkman(int& bot, int& top)
{
        // 将1-n-1的点按照xy进行排序
        std::sort(points, points + n, cmp_xy);

        int i,t;

        bot=n-1; top=n;
        D[top++]=0; D[top++]=1;

        for(i=2;i<n;i++)
        {
                //寻找第三个点 要保证3个点不共线！！
                if(isleft(points[D[top-2]],points[D[top-1]],points[i])!=0) break;
                D[top-1]=i;                        //共线就更换顶点
        }

        D[bot--]=i; D[top++]=i;                //i是第三个点 不共线！！

        //此时队列中有3个点，要保证3个点a,b,c是成逆时针的，不是就调换ab
        if(isleft(points[D[n]],points[D[n+1]],points[D[n+2]])<0)
        {
                t=D[n];
                D[n]=D[n+1];
                D[n+1]=t;
        }

        for(i++;i<n;i++)
        {
                //如果成立就是i在凸包内，跳过
                if(isleft(points[D[top-2]],points[D[top-1]],points[i])>0 &&
                                isleft(points[D[bot+1]],points[D[bot+2]],points[i])>0) continue;

                while(isleft(points[D[top-2]],points[D[top-1]],points[i])<=0) top--;
                D[top++]=i;

                while(isleft(points[D[bot+1]],points[D[bot+2]],points[i])<=0) bot++;
                D[bot--]=i;
        }

        // 凸包构造完成，D数组里bot+1至top-1内就是凸包的序列(头尾是同一点)
}

// 最近点对算法,采用分治算法
int CPY[MAXN], CPlen; // 数组Y'以及长度
double closestPoint_recursion(int L, int R, int& p1, int& p2)
{
        //acutPrintf(_T("\nL=%d, R=%d"), L, R);

        if(L >= R) return DBL_MAX;
        if(R-L == 1)
        {
                p1 = L;
                p2 = R;
                return dist_2(points[L], points[R]);
        }

        int mid = (L+R)>>1;
        int i,j;
        //for(i=mid; i>=L /*&& !cmp_real(points[i].x, points[mid].x)*/; i--);

        int lp1, lp2;
        double d1 = closestPoint_recursion(L, mid, lp1, lp2);

        //for(i=mid; i<=R /*&& !cmp_real(points[i].x, points[mid].x)*/; i++);
        int rp1, rp2;
        double d2 = closestPoint_recursion(mid, R, rp1, rp2);

        double ret;
        if((d1 - d2) < 0) // d1 < d2
        {
                p1 = lp1;
                p2 = lp2;
                ret = d1;
        }
        else             // d1 >= d2
        {
                p1 = rp1;
                p2 = rp2;
                ret = d2;
        }

        CPlen = 0;
        for(i=L; i<=R; i++)
        {
                double dx = (points[i].x - points[mid].x);
                if((dx*dx - ret) < 0)
                {
                        CPY[CPlen] = i;
                        CPlen++;
                }
        }

        std::sort(CPY, CPY+CPlen, cmp_y);

        for(i=0; i<CPlen; i++)
        {
                int cnt = 1;
                // 至多有8个点cnt = [0, 7]
                for(j=i+1; j<CPlen && cnt<=7; j++)
                {
                        // 不需要开方
                        double d = dist_2(points[CPY[i]], points[CPY[j]]);
                        if((d - ret) < 0) // d < ret
                        {
                                ret = d;
                                p1 = CPY[i];
                                p2 = CPY[j];
                        }
                        cnt++;
                }
        }
        return ret;
}

//分治算法求最近点对
void closestPoint(Point& pt1, Point& pt2)
{
        // 预排序(先x后y)
        std::sort(points, points+n, cmp_xy);
        //分治算法计算最近点对
        int p1, p2;
        double d = closestPoint_recursion(0, n-1, p1, p2);
        //得到最近点对
        pt1.x = points[p1].x; pt1.y= points[p1].y;
        pt2.x = points[p2].x; pt2.y = points[p2].y;
}

// graham算法求凸包
bool graham_ConvexHull(const PointArray& V, PointArray& H)
{
        // 至少有3个点数据
        if(V.size() < 3)
        {
                return false;
        }

        // 从V中读取数据到全局变量points中
        for(int i=0;i<V.size();i++)
        {
                points[i].x = V[i].x;
                points[i].y = V[i].y;
        }

        // 使用GrahamScan算法计算凸包，保存在stack数组中
        // 返回数组的最大索引位置top
        int top;
        graham(top);

        // 从S数组中获取凸包
        for(int i=0; i<=top; ++i)
        {
                Point pt;
                pt.x = points[S[i]].x;
                pt.y = points[S[i]].y;
                H.push_back(pt);
        }

        return true;
}

// melkman算法求凸包
bool melkman_ConvexHull( const PointArray& V, PointArray& H )
{
        / 至少有3个点数据
                if(V.size() < 3)
                {
                        return false;
                }

        // 从V中读取数据到全局变量points中
        for(int i=0;i<V.size();i++)
        {
                points[i].x = V[i].x;
                points[i].y = V[i].y;
        }

        // 使用melkman算法计算凸包，保存在D数组中
        // 返回数组的底部bot和顶部top
        int bot, top;
        melkman(bot, top);

        // 从stack数组中获取凸包
        for(int i=bot+1;i<top-1;i++)
        {
                Point pt;
                pt.x = points[D[i]].x;
                pt.y = points[D[i]].y;
                H.push_back(pt);
        }

        return true;
}