mean value coordinates

meanvaluecoordinates.cpp

#include <Eigen/Dense>
typedef Vector2d v2d;

int mean_value_coordinates( const vector<v2d> &cageCoords, const v2d &queryCoord, vector<double> &baryCoords)
{
	//////////////////////////////////////////////////////////////////////////////////
    // Input :
    //      1. polyCoords :  Coordinates of closed polygon in the Counter
    //                       clockwise direction. The input is not tested inside.
    //
    //       2. queryCoord:   the xyCoords of the query Point
    // Output:
    //       1:  baryCoords: baryCentric Coords of the query Point.
    //
    // Reference: Mean Value Coordinates for Arbitrary Planar Polygons:
    //            Kai Hormann and Michael Floater;
    // Written by:
    //            Chaman Singh Verma
    //            University of Wisconsin at Madison.
    //            18th March, 2011.
	// Modified by
	//            Jason
	//            The Chinese University of Hong Kong
	//            22nd, March, 2015
    /////////////////////////////////////////////////////////////////////////////////

    int nSize = cageCoords.size();
    assert( nSize );

    double dx, dy;

    vector<v2d>  s(nSize);
    for( int i = 0; i < nSize; i++)
    {
        dx  =   cageCoords[i][0] - queryCoord[0];
        dy  =   cageCoords[i][1] - queryCoord[1];
        s[i][0]  =   dx;
        s[i][1]  =   dy;
    }

    baryCoords.resize(nSize);
    for( int i = 0; i < nSize; i++)
        baryCoords[i] = 0.0;

    int ip, im;      // (i+1) and (i-1)
    double ri, rp, Ai, Di, dl, mu;  // Distance
    double eps = 10.0*std::numeric_limits<double>::min();

    // First check if any coordinates close to the cage point or
    // lie on the cage boundary. These are special cases.
    for( int i = 0; i < nSize; i++)
    {
        ip = (i+1)%nSize;
        ri = sqrt( s[i][0]*s[i][0] + s[i][1]*s[i][1] );
        Ai = 0.5*(s[i][0]*s[ip][1] - s[ip][0]*s[i][1]);
        Di = s[ip][0]*s[i][0] + s[ip][1]*s[i][1];
        if( ri <= eps)
        {
            baryCoords[i] = 1.0;
            return 0;
        }
        else if( fabs(Ai) <= 0 && Di < 0.0)
        {
            dx = cageCoords[ip][0] - cageCoords[i][0];
            dy = cageCoords[ip][1] - cageCoords[i][1];
            dl = sqrt(dx*dx + dy*dy);
            assert(dl > eps);
            dx = queryCoord[0] - cageCoords[i][0];
            dy = queryCoord[1] - cageCoords[i][1];
            mu = sqrt(dx*dx + dy*dy)/dl;
            assert( mu >= 0.0 && mu <= 1.0);
            baryCoords[i]  = 1.0-mu;
            baryCoords[ip] = mu;
            return 0;
        }
    }

    // Page #12, from the paper
    vector<double> tanalpha(nSize); // tan(alpha/2)
    for( int i = 0; i < nSize; i++)
    {
        ip = (i+1)%nSize;
        im = (nSize-1+i)%nSize;
        ri = sqrt( s[i][0]*s[i][0] + s[i][1]*s[i][1] );
        rp = sqrt( s[ip][0]*s[ip][0] + s[ip][1]*s[ip][1] );
        Ai = 0.5*(s[i][0]*s[ip][1] - s[ip][0]*s[i][1]);
        Di = s[ip][0]*s[i][0] + s[ip][1]*s[i][1];
        tanalpha[i] = (ri*rp - Di)/(2.0*Ai);
    }

    // Equation #11, from the paper
    double wi, wsum = 0.0;
    for( int i = 0; i < nSize; i++)
    {
        im = (nSize-1+i)%nSize;
        ri = sqrt( s[i][0]*s[i][0] + s[i][1]*s[i][1] );
        wi = 2.0*( tanalpha[i] + tanalpha[im] )/ri;
        wsum += wi;
        baryCoords[i] = wi;
    }

    if( fabs(wsum ) > 0.0)
    {
        for( int i = 0; i < nSize; i++)
            baryCoords[i] /= wsum;
    }

    return 0;
}