in vec2 texCoord;

out vec2 newPosition;

uniform int quads_x;
uniform int quads_y;

uniform GridMatrix{
   vec2 grid[400];
} g;

float distToLine(vec2 pt1, vec2 pt2, vec2 testPt){
    vec2 lineDir = pt2 - pt1;
    vec2 perpDir = vec2(lineDir.y, -lineDir.x);
    vec2 dirToPt1 = pt1 - testPt;
    return abs(dot(normalize(perpDir), dirToPt1));
}

bool pointInTriangle(vec2 A, vec2 B, vec2 C, vec2 X){
    float alpha, beta;

    // Calculate barycentric coordinates
    float det = ((B.y-C.y)*(A.x-C.x)+(C.x-B.x)*(A.y-C.y));
    alpha = ((B.y-C.y)*(X.x-C.x)+(C.x-B.x)*(X.y-C.y))/det;
    beta  = ((C.y-A.y)*(X.x-C.x)+(A.x-C.x)*(X.y-C.y))/det;

    if(alpha+beta >= 0.0 && alpha+beta <= 1.0 && alpha >= 0.0 && alpha <= 1.0 && beta >= 0.0 && beta <= 1.0)
        return true;
    return false;
}


vec2 bilinearInterpolation(vec2 point, vec2 local_bl, vec2 local_br, vec2 local_tl, vec2 local_tr, bool vert_parallel, bool hori_parallel){
    float s, t;

    bool vert_parallel = cross(vec3((local_tl-local_bl), 0), vec3(local_tr-local_br, 0)).z == 0.0;
    bool hori_parallel = cross(vec3((local_tl-local_tr), 0), vec3(local_bl-local_br, 0)).z == 0.0;

    vec3 A = vec3(local_bl.x, local_bl.y, 0);
    vec3 B = vec3(local_br.x, local_br.y, 0);
    vec3 C = vec3(local_tr.x, local_tr.y, 0);
    vec3 D = vec3(local_tl.x, local_tl.y, 0);
    vec3 M = vec3(point, 0);

    if(!vert_parallel && !hori_parallel){
        float a = (cross(C,D)+cross(D,B)+cross(A,C)+cross(B,A)).z;
        float b = (cross(M,C)+cross(D,M)+cross(B,M)+cross(M,A)+cross(C,A)+cross(B,D)+2.0*cross(A,B)).z;
        float c = (cross(M,B)+cross(A,M)+cross(B,A)).z;

        float t1 = (-b+sqrt(pow(b,2.0)-4.0*a*c))/(2.0*a);
        float t2 = (-b-sqrt(pow(b,2.0)-4.0*a*c))/(2.0*a);

        t = t1;
        if(t1 < 0.0 || t1 > 1.0)
            if(abs(t1) > abs(t2))
                t = t2;

        s = (M.x-A.x-t *(D.x-A.x))/(B.x-A.x+t*(C.x-D.x-B.x+A.x));
    }
    else if(vert_parallel && !hori_parallel){ // vertical edges parallel
            float dist1 = distToLine(local_tl, local_bl, point);
            float dist2 = distToLine(local_tr, local_br, point);
            s = dist1/(dist1+dist2);
            t = (M.y-A.y-s *(B.y-A.y))/(D.y-A.y+s*(C.y-D.y-B.y+A.y));
    }
    else if(!vert_parallel && hori_parallel){// horizontal edges parallel
            float dist1 = distToLine(local_bl, local_br, point);
            float dist2 = distToLine(local_tl, local_tr, point);
            t = dist1/(dist1+dist2);

            s = (M.x-A.x-t *(D.x-A.x))/(B.x-A.x+t*(C.x-D.x-B.x+A.x));
    }
    else{
            float dist1 = distToLine(local_bl, local_br, point);
            float dist2 = distToLine(local_tl, local_tr, point);
            t = dist1/(dist1+dist2);

            dist1 = distToLine(local_bl, local_tl, point);
            dist2 = distToLine(local_br, local_tr, point);
            s = dist1/(dist1+dist2);
    }
    return vec2(s, t);
}

vec2 staticWarp(vec2 point){
    int x = 0;
    int y = 0;
    while(y < quads_y){
        int curr_y = y*(quads_x+1);
        int curr_y_plus1 = (y+1)*(quads_x+1);
        vec2 bl = g.grid[curr_y_plus1+x]; // bottom left
        vec2 br = g.grid[curr_y_plus1+x+1]; // bottom right
        vec2 tl = g.grid[curr_y+x]; // top left
        vec2 tr = g.grid[curr_y+x+1]; // top right
        if(pointInTriangle(bl, br, tl, point) || pointInTriangle(tl, tr, br, point)){
             point = bilinearInterpolation(point, bl, br, tl, tr);
             point = vec2((x + point.x)/quads_x, (quads_y - y - 1 + point.y)/quads_y);
             return point;
        }
        x = x + 1;
        if ( x ==  quads_x){
            x = 0;
            y = y + 1;
        }
    }
    return vec2(-1.0, -1.0);
}


void main()
{
    vec2 point = texCoord;

    point = staticWarp(point);

    newPosition = point;
}