// Helper functions
// OSX seems to be using very strict compiler setting which demand the function
// prototype be defined before the function is implemented.  It seems like we
// only need the function prototype for helper methods, not kernels that are
// used directly.

float3 mtxMul( float16 matrix , float3 point );
float3 mtxMulAffine( float16 matrix , float3 point );
float3 lerp( const float3 a , const float3 b , const float factor );

//------------------------------------------------------------------------------
float3 calculateVertexNormal(
    unsigned int iAid,
    __global const float* iVerts,
    __global const unsigned int* iFrameIndexer,
    __global const unsigned int* iFramePairs);

//------------------------------------------------------------------------------
float16 calculateVertexFrame(
    unsigned int iAid,
    __global const float* iVerts,
    __global const unsigned int* iFrameIndexer,
    __global const unsigned int* iFramePairs);

//------------------------------------------------------------------------------
float3 mtxMul( float16 matrix , float3 point )
{
    const float4 homogeneousPoint = (float4)( point.xyz , 1.0f );
    const float4 transformedPoint = (float4)(
        dot( homogeneousPoint , matrix.s048c ) ,
        dot( homogeneousPoint , matrix.s159d ) ,
        dot( homogeneousPoint , matrix.s26ae ) ,
        dot( homogeneousPoint , matrix.s37bf )
        );
    return transformedPoint.xyz / transformedPoint.w;
}

//------------------------------------------------------------------------------
float3 mtxMulAffine( float16 matrix , float3 point )
{
    const float4 homogeneousPoint = (float4)( point.xyz , 1.0f );
    return (float3)(dot( homogeneousPoint , matrix.s048c ) ,
                    dot( homogeneousPoint , matrix.s159d ) ,
                    dot( homogeneousPoint , matrix.s26ae ));
}

//------------------------------------------------------------------------------
float3 lerp( const float3 a , const float3 b , const float factor )
{
    return a + factor * ( b - a );
}

//------------------------------------------------------------------------------
float3 calculateVertexNormal(
    unsigned int iAid,
    __global const float* iVerts,
    __global const unsigned int* iFrameIndexer,
    __global const unsigned int* iFramePairs)
{
    float3 normal = (float3)(0.0, 0.0, 0.0);

    const unsigned int s = iFrameIndexer[iAid];
    const unsigned int e = iFrameIndexer[iAid+1];

    if (s == e) {
        return normal;
    }

    float3 p = vload3(iAid, iVerts);

    for (unsigned int i = s; i < e; i+=2) {

        const float3 v0 = vload3(iFramePairs[i], iVerts) - p;
        const float3 v1 = vload3(iFramePairs[i+1], iVerts) - p;

        normal += cross(v0, v1); // area size influences normal (normalize to remove this)
    }

    return normalize(normal);
}

//------------------------------------------------------------------------------
float16
calculateVertexFrame(
    unsigned int iAid,
    __global const float* iVerts,
    __global const unsigned int* iFrameIndexer,
    __global const unsigned int* iFramePairs)
{
    float16 frame;

    /*
        const float3 vp = vload3(iAid, iVerts);
        frame.s0123 = (float4)(1.0f, 0.0f, 0.0f, 0.0f);
        frame.s4567 = (float4)(0.0f, 1.0f, 0.0f, 0.0f);
        frame.s89ab = (float4)(0.0f, 0.0f, 1.0f, 0.0f);
        frame.scdef = (float4)(vp, 1.0f);
    */

    const float3 normal = calculateVertexNormal(iAid, iVerts, iFrameIndexer, iFramePairs);

    const unsigned int s = iFrameIndexer[iAid];
    const float3 vPosition = vload3(iAid, iVerts);
    const float3 nPosition = vload3(iFramePairs[s], iVerts);

    float3 yvec = normalize(nPosition - vPosition);
    const float3 xvec = normalize(cross(normal, yvec));
    // Note: xvec and normal are orthogonal unit vectors, so their cross
    //       product will be a unit vector
    yvec = cross(xvec, normal);

    frame.s0123 = (float4)(xvec, 0.0f);
    frame.s4567 = (float4)(yvec, 0.0f);
    frame.s89ab = (float4)(normal, 0.0f);
    frame.scdef = (float4)(vPosition, 1.0f);

    return frame;
}

//------------------------------------------------------------------------------
__kernel void tangentConstraintKernel(
    __global float* iOutPos,                                // arg 0
    __global const float* iInPos,                           // arg 1
    __global const unsigned char* iPinned,                  // arg 2
    __global const unsigned int* iFrameIndexer,             // arg 3
    __global const unsigned int* iFramePairs,               // arg 4
    const float iInwardConstraint,                          // arg 5
    const float iOutwardConstraint,                         // arg 6
    const uint iAffectCount)                                // arg 7
{
    const unsigned int aid = get_global_id(0);
    if ( aid >= iAffectCount ) return;

    if (iPinned[aid])
        return;

    const float3 normal = calculateVertexNormal(aid, iInPos, iFrameIndexer, iFramePairs);

    float3 currentPos = vload3(aid, iOutPos);

    const float3 initialPos = vload3(aid, iInPos);
    const float3 displacement = currentPos-initialPos;

    const float d = dot(displacement, normal);

    if (d < 0.0) {
        const float3 offset = (iInwardConstraint*d)*normal;
        currentPos -= offset;
        vstore3(currentPos, aid, iOutPos);
    } else if (d > 0.0) {
        const float3 offset = (iOutwardConstraint*d)*normal;
        currentPos -= offset;
        vstore3(currentPos, aid, iOutPos);
    }

}

//------------------------------------------------------------------------------
__kernel void averageVertKernel(
    __global float* iOutPos ,                               // arg 0
    __global const float* iInPos,                           // arg 1
    __global const unsigned char* iPinned,                  // arg 2
    const float iLambda,                                    // arg 3
    __global const unsigned int* iIndexer,                  // arg 4
    __global const unsigned int* iConnections,              // arg 5
    __global const float* iWeights,                         // arg 6
     const uint iAffectCount)                               // arg 7
{
    const unsigned int aid = get_global_id(0);
    if ( aid >= iAffectCount ) return;

    const float3 initialPos = vload3(aid, iInPos);

    if (iPinned[aid]) {
        vstore3(initialPos, aid , iOutPos);
        return;
    }

    const unsigned int s = iIndexer[aid];
    const unsigned int e = iIndexer[aid+1];
    if (s == e) {
        vstore3(initialPos, aid , iOutPos);
        return;
    }

    float3 neighborSum = (float3)(0.0f, 0.0f, 0.0f);
    for (unsigned int c = s; c < e; ++c) {
        float3 np = vload3(iConnections[c], iInPos);
        neighborSum = neighborSum + ( iWeights[c] * np);
    }

    neighborSum /= (float)(e-s);
    const float3 pos = initialPos + iLambda * (neighborSum - initialPos);
    vstore3(pos, aid , iOutPos);
}

//------------------------------------------------------------------------------
__kernel void applyDisplacementKernel(
    __global float* iResultVerts,                           // arg 0
    __global const float* iSourceVerts,                     // arg 1
    __global float* iDisplacements,                         // arg 2
    const float3 iWeightedScale,                            // arg 3
    __global const unsigned int* iFrameIndexer,             // arg 4
    __global const unsigned int* iFramePairs,               // arg 5
    const uint iAffectCount)                                // arg 6
{
    const unsigned int aid = get_global_id(0);
    if ( aid >= iAffectCount ) return;

    const float3 displacement = vload3(aid, iDisplacements);

    const float3 dir = (float3)(displacement.x * iWeightedScale.x,
                                displacement.y * iWeightedScale.y,
                                displacement.z * iWeightedScale.z);

    const float16 frame = calculateVertexFrame(aid, iSourceVerts, iFrameIndexer, iFramePairs);
    vstore3(mtxMulAffine( frame, dir ), aid, iResultVerts);
}