构造射线算法，测试点是否在曲线内

admin1

[code]enum IncidenceType {
    kIncidenceToLeft = 0,
    kIncidenceToRight = 1,
    kIncidenceToFront =2,
    kIncidenceUnknown,
};

IncidenceType CurveIncidence(AcDbCurve* pCurve, double param, AcGeVector3d dir, AcGeVector3d normal )
{
    AcGeVector3d deriv1;
    pCurve->getFirstDeriv( param, deriv1 );

    if( deriv1.isParallelTo( dir ))
    { // need second degree analysis
        AcGeVector3d deriv2;
        pCurve->getSecondDeriv( param, deriv2 );

        if( deriv2.isZeroLength() || deriv2.isParallelTo( dir ) )
        {
            return kIncidenceToFront;
        }
        else if(deriv2.crossProduct( dir ).dotProduct( normal ) < 0 )
        {
            return kIncidenceToRight;
        }
        else
        {
            return kIncidenceToLeft;
        }
    }

    if( deriv1.crossProduct( dir ).dotProduct( normal ) < 0 )
    {
        return kIncidenceToLeft;
    }
    else
    {
        return kIncidenceToRight;
    }
}

Adesk::Boolean IsInsideCurve(AcDbCurve*pCurve, AcGePoint3d testPt )
{
    if(pCurve== NULL || !pCurve->isClosed())
    {    // cannot be inside
        return Adesk::kFalse;
    }

    AcDb2dPolyline* p2dPoly = AcDb2dPolyline::cast(pCurve);

    if( p2dPoly != NULL && p2dPoly->polyType() != AcDb::k2dSimplePoly )
    { // Not supported
        return Adesk::kFalse;
    }

    AcGePoint3d ptOnCurve;
    Acad::ErrorStatus es = pCurve->getClosestPointTo(testPt, ptOnCurve );

    if( testPt == ptOnCurve )
    {
        return Adesk::kTrue;
    }

    AcGePlane plane;
    AcDb::Planarity planarity; // check its planar
    es = pCurve->getPlane(plane, planarity );

    if( es != Acad::eOk || planarity != AcDb::kPlanar )
    { // 曲线不共面
        return Adesk::kFalse;
    }

    // make the test ray from the plane
    AcGeVector3d normal = plane.normal();
    AcGeVector3d testVector = normal.perpVector();
    AcDbRay ray;
    ray.setBasePoint(testPt);
    ray.setUnitDir(testVector);

    AcGePoint3dArray IntersectionPoints;
    // fire the ray at the curve
    es = pCurve->intersectWith(&ray, AcDb::kOnBothOperands, IntersectionPoints);

    if(es != Acad::eOk)
    {
        return Adesk::kFalse;
    }

    int numberOfInters = IntersectionPoints.length();

    if (numberOfInters == 0)
    { // must be outside
        return Adesk::kFalse;
    }

    int nGlancingHits = 0;
    double epsilon = 2e-6;

    for( int i=0;i < numberOfInters; i++ )
    { // get the first point, and get its parameter
        AcGePoint3d hitPt = IntersectionPoints;
        double hitParam;
        es = pCurve->getParamAtPoint(hitPt, hitParam );

        if(es != Acad::eOk)
        {
            return Adesk::kFalse;
        }

        double inParam = hitParam - epsilon;
        double outParam = hitParam + epsilon;
        IncidenceType inIncidence = CurveIncidence(pCurve, inParam, testVector, normal );
        IncidenceType outIncidence = CurveIncidence(pCurve, outParam, testVector, normal );

        if ((inIncidence == kIncidenceToRight && outIncidence == kIncidenceToLeft) ||
            ( inIncidence == kIncidenceToLeft && outIncidence == kIncidenceToRight ) )
        {
            nGlancingHits ++;
        }
    }

    return (( numberOfInters + nGlancingHits) % 2 == 1);
}

void utilsisinside()
{
    ads_name eName;
    ads_point pt;

    if( RTNORM != ads_entsel( "\nPlease pick a curve ", eName, pt) )
    {
        return;
    }

    AcDbObjectId curveId;
    acdbGetObjectId( curveId, eName );
    AcDbCurve* pCurve = NULL;
    acdbOpenObject(pCurve, curveId, AcDb::kForRead );

    if( pCurve == NULL )
    {
        return;
    }

    while( RTNORM == ads_getpoint( NULL, "\nPlease pick a point ", pt ) )
    {
        acdbUcs2Wcs( pt, pt, Adesk::kFalse );
        Adesk::Boolean isIn;
        isIn = IsInsideCurve( pCurve, asPnt3d( pt ) );
        ads_printf("\n%s", isIn ? "INSIDE OR ON" : "OUTSIDE" );
    }
    pCurve->close();
}[/code]