28、几何算法-线段相交、凸包、球面弧长

1、判断线段是否相交
（1）标准方法：如果两直线不平行先求直线交点，再看这个交点是否分别在两个线段上。

（2）标准方法方法涉及到除法，计算机中一般方法：先判断以两个线段作为对角线的矩形是否相交，如果矩形相交再判断以一个线段为轴，另一个线段的两个端点是否分别位于顺时针和逆时针方向。

int lint(Point p1, Point p2, Point p3, Point p4)
{
    double z1, z2, z3, z4;
    int s1, s2, s3, s4;

    if(!( (MAX(p1.x, p2.x) >= MIN(p3.x, p4.x)) && (MAX(p3.x, p4.x) >= MIN(p1.x, p2.x))
        && (MAX(p1.y, p2.y) >= MIN(p3.y, p4.y))&& (MAX(p3.y, p4.y) >= MIN(p1.y, p2.y)) ))
        {
            return -1;
        }
    if( ( z1 = ((p3.x - p1.x)*(p2.y - p1.y)) - ((p3.y - p1.y)*(p2.x - p1.x)) ) < 0 )
        s1 = -1;
    else if(z1 > 0)
        s1 = 1;
    else
        s1 = 0;

    if( ( z2 = ((p4.x - p1.x)*(p2.y - p1.y)) - ((p4.y - p1.y)*(p2.x - p1.x)) ) < 0 )
        s2 = -1;
    else if(z1 > 0)
        s2 = 1;
    else
        s2 = 0;

    if( ( z3 = ((p1.x - p3.x)*(p4.y - p3.y)) - ((p1.y - p3.y)*(p4.x - p3.x)) ) < 0 )
        s3 = -1;
    else if(z1 > 0)
        s3 = 1;
    else
        s3 = 0;

    if( ( z4 = ((p2.x - p3.x)*(p4.y - p3.y)) - ((p2.y - p3.y)*(p4.x - p3.x)) ) < 0 )
        s4 = -1;
    else if(z1 > 0)
        s4 = 1;
    else
        s4 = 0;

    if((s1 * s2 <= 0) && (s3 * s4 <= 0))
        return 1;
    return -1; 
}

2、凸包
凸包：一个多边形内任意两点之间的连线都在多边形内部，即所有角都向外凸或平。

Jarvis’s March法：（1）选取y轴最低处点，为P0。（2）选取Pc，其他所有点和P0构成的直线顺时针到P0和Pc构成的直线。如下图z<0,则任何点Pi相比Pc符合这个条件。（3）如果z=0，则选离P0更远的点。（4）循环以Pc为P0，再次找Pc。

//P为原来所有节点构成的链表， polygon为最终所有外围点 
int cvxhull(const LIST_ATTRITIVE *P, LIST_ATTRITIVE *polygon)
{
    LIST_ELEMENT *element;
    Point *min, *low, *p0, *pi, *pc;
    double z, length1, length2;
    int count;//计数外围点的个数

    //找最低最左边的点为P0 
    min = list_data(list_head(p));
    for(element = list_head(p); element != NULL; element = list_next(element))
    {
        p0 = list_data(element);

        if(p0->y < min->y)
        {
            min = p0;
            low = list_data(element);
        }
        else
        {
            if( p0->y == min->y && p0->x < min->x)
            {
                    min = p0;
                    low = list_data(element);
            }
        }
    } 
    list_init(polygon, NULL);

    //先执行后判断，一直循环到再次找到P0这个点 
    p0 = low;
    do
    {
        list_ins_next(polygon, list_tail(polygon), p0);
        count = 0

        //找下一个点pc 
        for(element = list_head(p); element != NULL; element = list_next(element))
        {
            if((pi = list_data(element)) == p0)
                continue;
            count++;

            //每次大循环while，暂时把第二个点作为pc
            if(count == 1)
            {
                pc = list_data(element);
                continue;
            }

            if( ( z = ((pi.x - p0.x)*(pc.y - p0.y)) - ((pi.y - p0.y)*(pc.x - p0.x)) ) > 0 )
                pc = pi;

            //z==0时选取远的点
            else if(z == 0)
            {
                length1 = sqrt(pow(pi->x - p0->x, 2.0) + pow(pi->y - p0->y, 2.0));
                length2 = sqrt(pow(pc->x - p0->x, 2.0) + pow(pc->y - p0->y, 2.0));
                if(length1 > length2)
                    pc = pi
            } 
        } 

        p0 = pc;
    }while(p0 != low);

    return 0; 
}

3、球面弧长

void arclen(SPoint p1, SPoint p2, double *length)
{
    Point p1_rtc, p2_rtc;
    double alpha, dot;

    p1_rtc.x = p1.rho * sin(p1.phi) * cos(p1.theta);
    p1_rtc.y = p1.rho * sin(p1.phi) * sin(p1.theta);
    p1_rtc.z = p1.rho * cos(p1.phi);

    p2_rtc.x = p2.rho * sin(p2.phi) * cos(p2.theta);
    p2_rtc.y = p2.rho * sin(p2.phi) * sin(p2.theta);
    p2_rtc.z = p2.rho * cos(p2.phi);

    dot = (p1_rtc.x * p2_rtc.x) + (p1_rtc.y * p2_rtc.y) + (p1_rtc.z * p2_rtc.z);
    alpha = acos(dot / pow(p1.rho, 2.0));

    *length = alpha * p1.rho;
    return;
} 

4、球面弧长的应用-地球上两点的近视距离

地球半径约为3440065海里，北纬，西经为正。纬度（-90到+90需要转换成0到180，对应φ），经度（-180到+180需要转换成0到360，对应θ）