ogrlinearring.cpp

mitab,读取MapInfo的地图文件

/* -------------------------------------------------------------------- */
/*      Copy in the raw data.                                           */
/* -------------------------------------------------------------------- */
    memcpy( pabyData, &nPointCount, 4 );

/* -------------------------------------------------------------------- */
/*      Copy in the raw data.                                           */
/* -------------------------------------------------------------------- */
    if( b3D )
    {
        nWords = 3 * nPointCount;
        for( i = 0; i < nPointCount; i++ )
        {
            memcpy( pabyData+4+i*24, &(paoPoints[i].x), 8 );
            memcpy( pabyData+4+i*24+8, &(paoPoints[i].y), 8 );
            if( padfZ == NULL )
                memset( pabyData+4+i*24+16, 0, 8 );
            else
                memcpy( pabyData+4+i*24+16, padfZ + i, 8 );
        }
    }
    else
    {
        nWords = 2 * nPointCount; 
        memcpy( pabyData+4, paoPoints, 16 * nPointCount );
    }

/* -------------------------------------------------------------------- */
/*      Swap if needed.                                                 */
/* -------------------------------------------------------------------- */
    if( OGR_SWAP( eByteOrder ) )
    {
        int     nCount;

        nCount = CPL_SWAP32( nPointCount );
        memcpy( pabyData, &nCount, 4 );

        for( i = 0; i < nWords; i++ )
        {
            CPL_SWAPDOUBLE( pabyData + 4 + 8 * i );
        }
    }
    
    return OGRERR_NONE;
}

/************************************************************************/
/*                              _WkbSize()                              */
/*                                                                      */
/*      Helper method for OGRPolygon.  NOT THE NORMAL WkbSize() METHOD! */
/************************************************************************/

int OGRLinearRing::_WkbSize( int b3D ) const

{
    if( b3D )
        return 4 + 24 * nPointCount;
    else
        return 4 + 16 * nPointCount;
}

/************************************************************************/
/*                               clone()                                */
/*                                                                      */
/*      We override the OGRCurve clone() to ensure that we get the      */
/*      correct virtual table.                                          */
/************************************************************************/

OGRGeometry *OGRLinearRing::clone() const

{
    OGRLinearRing       *poNewLinearRing;

    poNewLinearRing = new OGRLinearRing();
    poNewLinearRing->assignSpatialReference( getSpatialReference() );

    poNewLinearRing->setPoints( nPointCount, paoPoints, padfZ );

    return poNewLinearRing;
}

/************************************************************************/
/*                            isClockwise()                             */
/************************************************************************/

/**
 * Returns TRUE if the ring has clockwise winding.
 *
 * @return TRUE if clockwise otherwise FALSE.
 */

int OGRLinearRing::isClockwise() const

{
    double dfSum = 0.0;

    for( int iVert = 0; iVert < nPointCount-1; iVert++ )
    {
        dfSum += paoPoints[iVert].x * paoPoints[iVert+1].y
            - paoPoints[iVert].y * paoPoints[iVert+1].x;
    }

    dfSum += paoPoints[nPointCount-1].x * paoPoints[0].y
        - paoPoints[nPointCount-1].y * paoPoints[0].x;

    return dfSum < 0.0;
}

/************************************************************************/ 
/*                             reverseWindingOrder()                    */ 
/************************************************************************/ 

void OGRLinearRing::reverseWindingOrder() 

{ 
    int pos = 0; 
    OGRPoint tempPoint; 

    for( int i = 0; i < nPointCount / 2; i++ ) 
    { 
        getPoint( i, &tempPoint ); 
        pos = nPointCount - i - 1; 
        setPoint( i, getX(pos), getY(pos), getZ(pos) ); 
        setPoint( pos, tempPoint.getX(), tempPoint.getY(), tempPoint.getZ() ); 
    } 
} 

/************************************************************************/
/*                             closeRing()                              */
/************************************************************************/

void OGRLinearRing::closeRings()

{
    if( nPointCount < 2 )
        return;

    if( getX(0) != getX(nPointCount-1) 
        || getY(0) != getY(nPointCount-1)
        || getZ(0) != getZ(nPointCount-1) )
    {
        /* Avoid implicit change of coordinate dimensionality
         * if z=0.0 and dim=2
         */
        if( getCoordinateDimension() == 2 )
            addPoint( getX(0), getY(0) );
        else
            addPoint( getX(0), getY(0), getZ(0) );
            
    }
}

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

/**
 * Compute area of ring.
 *
 * The area is computed according to Green's Theorem:  
 *
 * Area is "Sum(x(i)*y(i+1) - x(i+1)*y(i))/2" for i = 0 to pointCount-1, 
 * assuming the last point is a duplicate of the first. 
 *
 * @return computed area.
 */

double OGRLinearRing::get_Area() const

{
    double dfAreaSum = 0.0;
    int i;

    for( i = 0; i < nPointCount-1; i++ )
    {
        dfAreaSum += 0.5 * ( paoPoints[i].x * paoPoints[i+1].y 
                             - paoPoints[i+1].x * paoPoints[i].y );
    }

    dfAreaSum += 0.5 * ( paoPoints[nPointCount-1].x * paoPoints[0].y 
                         - paoPoints[0].x * paoPoints[nPointCount-1].y );

    return fabs(dfAreaSum);
}

/************************************************************************/
/*                              isPointInRing()                         */
/************************************************************************/

OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint* poPoint) const
{
    if ( NULL == poPoint )
    {
        CPLDebug( "OGR", "OGRLinearRing::isPointInRing(const  OGRPoint* poPoint) - passed point is NULL!" );
        return 0;
    }

    CPLDebug( "OGR", "OGRLinearRing::isPointInRing(): passed point: (%.8f,%.8f)",
              poPoint->getX(), poPoint->getY() );

    const int iNumPoints = getNumPoints();

    // Simple validation
    if ( iNumPoints < 4 )
        return 0;

    const double dfTestX = poPoint->getX();
    const double dfTestY = poPoint->getY();

    // Fast test if point is inside extent of the ring
    OGREnvelope extent;
    getEnvelope(&extent);
    if ( !( dfTestX >= extent.MinX && dfTestX <= extent.MaxX
         && dfTestY >= extent.MinY && dfTestY <= extent.MaxY ) )
    {
        CPLDebug( "OGR", "OGRLinearRing::isPointInRing(): passed point is out of extent of ring" );
        return 0;
    }

	// For every point p in ring,
    // test if ray starting from given point crosses segment (p - 1, p)
    int iNumCrossings = 0;

    for ( int iPoint = 1; iPoint < iNumPoints; iPoint++ ) 
    {
        const int iPointPrev = iPoint - 1;

        const double x1 = getX(iPoint) - dfTestX;
        const double y1 = getY(iPoint) - dfTestY;

        const double x2 = getX(iPointPrev) - dfTestX;
        const double y2 = getY(iPointPrev) - dfTestY;

        if( ( ( y1 > 0 ) && ( y2 <= 0 ) ) || ( ( y2 > 0 ) && ( y1 <= 0 ) ) ) 
        {
            // Check if ray intersects with segment of the ring
            const double dfIntersection = ( x1 * y2 - x2 * y1 ) / (y2 - y1);
            if ( 0.0 < dfIntersection )
            {
                // Count intersections
                iNumCrossings++;
            }
        }
    }

    // If iNumCrossings number is even, given point is outside the ring,
    // when the crossings number is odd, the point is inside the ring.
    return ( ( iNumCrossings % 2 ) == 1 ? 1 : 0 );
}