ugkpolygon.cpp

linux下一款GIS程序源码

int UGKPolygon::Centroid( UGKPoint *poPoint ) const

{
	// notdef ... not implemented yet.
    return UGKERR_FAILURE;

}

/************************************************************************/
/*                           PointOnSurface()                           */
/************************************************************************/

int UGKPolygon::PointOnSurface( UGKPoint * ) const

{
    // notdef ... not implemented yet.
    
    return UGKERR_FAILURE;
}



/************************************************************************/
/*                            getEnvelope()                             */
/************************************************************************/
 
void UGKPolygon::getEnvelope( UGKEnvelope * psEnvelope ) const

{
    UGKEnvelope         oRingEnv;
    
    if( nRingCount == 0 )
        return;

    papoRings[0]->getEnvelope( psEnvelope );

    for( int iRing = 1; iRing < nRingCount; iRing++ )
    {
        papoRings[iRing]->getEnvelope( &oRingEnv );

        if( psEnvelope->MinX > oRingEnv.MinX )
            psEnvelope->MinX = oRingEnv.MinX;
        if( psEnvelope->MinY > oRingEnv.MinY )
            psEnvelope->MinY = oRingEnv.MinY;
        if( psEnvelope->MaxX < oRingEnv.MaxX )
            psEnvelope->MaxX = oRingEnv.MaxX;
        if( psEnvelope->MaxY < oRingEnv.MaxY )
            psEnvelope->MaxY = oRingEnv.MaxY;
    }
}



/************************************************************************/
/*                               Equal()                                */
/************************************************************************/

UGKBool UGKPolygon::Equals( UGKGeometry * poOther ) const

{
    UGKPolygon *poOPoly = (UGKPolygon *) poOther;

    if( poOPoly == this )
        return TRUE;
    
    if( poOther->getGeometryType() != getGeometryType() )
        return FALSE;

    if( getNumInteriorRings() != poOPoly->getNumInteriorRings() )
        return FALSE;

    if( !getExteriorRing()->Equals( poOPoly->getExteriorRing() ) )
        return FALSE;
    

    for( int iRing = 0; iRing < getNumInteriorRings(); iRing++ )
    {
        if( !getInteriorRing(iRing)->Equals(poOPoly->getInteriorRing(iRing)) )
            return FALSE;
    }

    return TRUE;
}



/************************************************************************/
/*                             closeRings()                             */
/************************************************************************/

void UGKPolygon::closeRings()

{
    for( int iRing = 0; iRing < nRingCount; iRing++ )
        papoRings[iRing]->closeRings();
}


/************************************************************************/
/*                              get_Area()                              */
/************************************************************************/

/**
 * Compute area of polygon.
 *
 * The area is computed as the area of the outer ring less the area of all
 * internal rings. 
 *
 * @return computed area.
 */


double UGKPolygon::get_Area() const

{
    double dfArea = 0.0;

    if( getExteriorRing() != NULL )
    {
        int iRing;

        dfArea = getExteriorRing()->get_Area();

        for( iRing = 0; iRing < getNumInteriorRings(); iRing++ )
            dfArea -= getInteriorRing( iRing )->get_Area();
    }

    return dfArea;
}


int UGKPolygonLabelPoint(UGKPolygon *poPoly, UGKPoint *poLabelPoint)
{
    double        slope;
    UGKRawPoint   point1, point2;
    int           i, j, k, nfound;
    double        x, y, *xintersect, temp;
    double        hi_y, lo_y;
    int           wrong_order, n;
    double        len, max_len=0;
    double        skip;
    UGKEnvelope   oEnv;

    if (poPoly == NULL)
        return UGKERR_FAILURE;

    poPoly->getEnvelope(&oEnv);

    poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
    poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);

    //if(get_centroid(p, lp, &miny, &maxy) == -1) return(-1);

    if(UGKIntersectPointPolygon(poLabelPoint, poPoly) == TRUE) /* cool, done */
        return UGKERR_NONE;

    /* do it the hard way - scanline */

    skip = (oEnv.MaxY - oEnv.MinY)/NUM_SCANLINES;

    n=0;
    for(j=0; j<UGK_NUM_RINGS(poPoly); j++) 
    {
        /* count total number of points */
        n += UGK_GET_RING(poPoly, j)->getNumPoints();
    }

    xintersect = (double *)UGK_Calloc(n, sizeof(double));
    if (xintersect == NULL)
        return UGKERR_FAILURE;

    for(k=1; k<=NUM_SCANLINES; k++) 
    { 
        /* sample the shape in the y direction */
    
        y = oEnv.MaxY - k*skip; 

        /* need to find a y that won't intersect any vertices exactly */  
        hi_y = y - 1; /* first initializing lo_y, hi_y to be any 2 pnts on either side of lp->y */
        lo_y = y + 1;
        for(j=0; j<UGK_NUM_RINGS(poPoly); j++) 
        {
            UGKLinearRing *poRing = UGK_GET_RING(poPoly,j);

            if((lo_y < y) && (hi_y >= y)) 
                break; /* already initialized */
            for(i=0; i < poRing->getNumPoints(); i++) 
            {   
                if((lo_y < y) && (hi_y >= y)) 
                    break; /* already initialized */
                if(poRing->getY(i) < y)
                    lo_y = poRing->getY(i);
                if(poRing->getY(i) >= y)
                    hi_y = poRing->getY(i);
            }
        }

        n=0;
        for(j=0; j<UGK_NUM_RINGS(poPoly); j++) 
        {
            UGKLinearRing *poRing = UGK_GET_RING(poPoly,j);

            for(i=0; i < poRing->getNumPoints(); i++) 
            {
                if((poRing->getY(i) < y) && 
                   ((y - poRing->getY(i)) < (y - lo_y)))
                    lo_y = poRing->getY(i);
                if((poRing->getY(i) >= y) && 
                   ((poRing->getY(i) - y) < (hi_y - y)))
                    hi_y = poRing->getY(i);
            }      
        }

        if(lo_y == hi_y) 
            return UGKERR_FAILURE;
        else  
            y = (hi_y + lo_y)/2.0;    
    
        nfound = 0;
        for(j=0; j<UGK_NUM_RINGS(poPoly); j++)   /* for each line */
        {
            UGKLinearRing *poRing = UGK_GET_RING(poPoly,j);
            point1.x = poRing->getX(poRing->getNumPoints()-1);
            point1.y = poRing->getY(poRing->getNumPoints()-1);

            for(i=0; i < poRing->getNumPoints(); i++) 
            {
                point2.x = poRing->getX(i);
                point2.y = poRing->getY(i);
        
                if(EDGE_CHECK(point1.y, y, point2.y) == CLIP_MIDDLE) 
                {
                    if(point1.y == point2.y)
                        continue; /* ignore horizontal edges */
                    else
                        slope = (point2.x - point1.x) / (point2.y - point1.y);
          
                    x = point1.x + (y - point1.y)*slope;
                    xintersect[nfound++] = x;
                } /* End of checking this edge */
        
                point1 = point2;  /* Go on to next edge */
            }
        } /* Finished the scanline */
    
        /* First, sort the intersections */
        do 
        {
            wrong_order = 0;
            for(i=0; i < nfound-1; i++) 
            {
                if(xintersect[i] > xintersect[i+1]) 
                {
                    wrong_order = 1;
                    SWAP(xintersect[i], xintersect[i+1], temp);
                }
            }
        } while(wrong_order);
    
        /* Great, now find longest span */
        point1.y = point2.y = y;
        for(i=0; i < nfound; i += 2) 
        {
            point1.x = xintersect[i];
            point2.x = xintersect[i+1];
            /* len = length(point1, point2); */
            len = ABS((point2.x - point1.x));
            if(len > max_len) 
            {
                max_len = len;
                poLabelPoint->setX( (point1.x + point2.x)/2 );
                poLabelPoint->setY( y );
            }
        }
    }

    free(xintersect);

    /* __TODO__ Bug 673
     * There seem to be some polygons for which the label is returned
     * completely outside of the polygon's MBR and this messes the 
     * file bounds, etc.
     * Until we find the source of the problem, we'll at least validate
     * the label point to make sure that it overlaps the polygon MBR.
     */
    if( poLabelPoint->getX() < oEnv.MinX
        || poLabelPoint->getY() < oEnv.MinY
        || poLabelPoint->getX() > oEnv.MaxX
        || poLabelPoint->getY() > oEnv.MaxY )
    {
        // Reset label coordinates to center of MBR, just in case
        poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
        poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);

        // And return an error
        return UGKERR_FAILURE;
    }

    if(max_len > 0)
        return UGKERR_NONE;
    else
        return UGKERR_FAILURE;
}


/**********************************************************************
 *                   UGKIntersectPointPolygon()
 *
 * Instead of using ring orientation we count the number of parts the
 * point falls in. If odd the point is in the polygon, if 0 or even
 * then the point is in a hole or completely outside.
 *
 * Returns TRUE is point is inside polygon, FALSE otherwise
 *
 * Adapted version of msIntersectPointPolygon() from MapServer's mapsearch.c
 **********************************************************************/
UGKBool UGKIntersectPointPolygon(UGKPoint *poPoint, UGKPolygon *poPoly) 
{
    int i;
    UGKBool status = FALSE;

    for(i=0; i<UGK_NUM_RINGS(poPoly); i++) 
    {
        if (UGKPointInRing(poPoint, UGK_GET_RING(poPoly, i))) 
        {
            /* ok, the point is in a line */
            status = !status;
        }
    }

    return status;
}