Search code examples
c++computational-geometryraycastingconvex-hull

Can't understand the solution of HDU 2823


The following code snippet was taken from here. It is the solution of this problem HDU 2823.

#define eps 1e-9
double rc(point pp[],point qq[],int n,int m)    
{    
    int q=0;    
    int p=0;    
    for(int i=0;i<n;i++)     
        if(pp[i].y-pp[p].y<-eps)    
            p=i;    
    for(int i=0;i<m;i++)    
        if(qq[i].y-qq[q].y>eps)    
            q=i;    
    pp[n]=pp[0];    
    qq[m]=qq[0];    
    double tmp,ans=1e99;    
    for(int i=0;i<n;i++)    
    {    
        while((tmp=cross(pp[p+1],qq[q+1],pp[p])-cross(pp[p+1],qq[q],pp[p]))>eps)    
            q=(q+1)%m;    
        if(tmp<-eps)    
            ans=min(ans,dist_p_to_seg(qq[q],pp[p],pp[p+1]));    
        else    
            ans=min(ans,dist_seg_to_seg(pp[p],pp[p+1],qq[q],qq[q+1]));    
        p=(p+1)%n;    
    }    
    return ans;    

}    

pp[] and qq[] are two different convex hull. p is the highest point of pp convex hull and q is the lowest point of qq convex hull.

I can't seem to understand this line:

while((tmp=cross(pp[p+1],qq[q+1],pp[p])-cross(pp[p+1],qq[q],pp[p]))>eps) 
    q=(q+1)%m;

What is he trying to achieve?


Solution

  • The function cross(a, b, c) is finding the determinant of the following matrix,

    | a.x a.y 1 |
    | b.x b.x 1 |  = 2 * A
    | c.x c.y 1 | 
    

    Where A is the signed area of the triangle a, b, c. The sign of the determinant also tells us if 3 points are oriented clockwise or conter-clockwise. see here for an explanation

    Let's rewrite like this,

    triA ← cross(pp[p+1],qq[q+1],pp[p])
    triB ← cross(pp[p+1],qq[q],pp[p])
    
    // This is equivalent to,
    // just to make it a bit clearer
    triA ← cross(pp[p], pp[p+1],    qq[q+1])
    triB ← cross(pp[p], pp[p+1],    qq[q])
    

    So it checks if the triangle formed by one side of the hull pp with the lowest point on qq is smaller than the triangle formed by the same side and next highest point from qq.

    If yes, select the next point in qq as q and continue. —i.e. select a q so that the perpendicular distance of q from side <p, p+1> is minimized.

    enter image description here

    Once this is locally minimised for a given side <p, p+1>, repeat this for all sides of pp. At each step keep the minimum distance between the current sides.

    This is a somewhat shaky way of finding the minimum separation between two convexhulls. The intuition is correct —correct in the sense that it is simple to understand (the idea, not the code in question), quite general for convex polygons, and very useful a variety of problems (see the reference below); however, I have a feeling that this can be written in a much more efficient and easy to understand way.

    The intuition behind such ideas is nicely illustrated in this paper "Solving Geometric Problems with the Rotating Calipers" -- Toussaint G.