kernel MultigradKernel : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

  param:
    int numPoints; // Number of points that will be used
    float4 c0, c1, c2, c3, c4, c5, c6, c7, c8, c9;
    int2 p0, p1, p2, p3, p4, p5, p6, p7, p8, p9;


  void process(int2 pos) {
    float4 colors [10] = {c0,c1,c2,c3,c4,c5,c6,c7,c8,c9};
    int2 points [10] = {p0,p1,p2,p3,p4,p5,p6,p7,p8,p9};

    //--------------
    // 1. FIND THE 3 CLOSEST POINTS

    int vertexId [3] = {-1,-1,-1}; // Stores the ids of the closest vertices to avoid repetition
    float4 vertexColors [3]; // Will store the color of the 3 closest vertices
    vertexColors[0] = colors[0]; // Initialising...
    int2 vertexPos [3]; // Will store the position of the 3 closest vertices
    vertexPos[0] = points[0]; // Initialising...

    float vertexDistance [3] = {10000000000.0f,10000000000.0f,10000000000.0f}; // Distance to vertices, for iterating
    int vertices = min(numPoints, 3);
    float dist = 0.0f;
    //Point 1
    for(int i = 0; i<numPoints; i++){
      dist = length(float2(float(pos.x),float(pos.y))-float2(float(points[i][0]),float(points[i][1])));
      if(dist < vertexDistance[0]){
        vertexDistance[0] = dist;
        vertexColors[0] = colors[i];
        vertexPos[0] = points[i];
        vertexId[0] = i;
      }
    }
    //dst() = vertexColors[0];
    //dst() = vertexId[0] ;
    //dst() = float(vertexPos[0][0]);

    //Point 2
    if(numPoints > 1){
      for(int i = 0; i<numPoints; i++){
        dist = length(float2(float(pos.x),float(pos.y))-float2(float(points[i][0]),float(points[i][1])));
        if(dist < vertexDistance[1] && i != vertexId[0]){
          vertexDistance[1] = dist;
          vertexColors[1] = colors[i];
          vertexPos[1] = points[i];
          vertexId[1] = i;
        }
      }
    }

    //Point 3
    if(numPoints > 2){
      for(int i = 0; i<numPoints; i++){
        dist = length(float2(float(pos.x),float(pos.y))-float2(float(points[i][0]),float(points[i][1])));
        if(dist < vertexDistance[2] && i != vertexId[0] && i != vertexId[1]){
          vertexDistance[2] = dist;
          vertexColors[2] = colors[i];
          vertexPos[2] = points[i];
          vertexId[2] = i;
        }
      }
    }

    // 0,1 and 2 are ordered in their distance. Next step is to calculate the barycentric coordinates, in order to find out
    // if the pixel lays on the current triangle or not? Or otherwise we can simply take the ones with positive weights...

    //--------------
    // 2. GET THE VERTEX WEIGHTS

    float x1,x2,x3,y1,y2,y3, Px, Py;
    x1 = float(vertexPos[0][0]);
    x2 = float(vertexPos[1][0]);
    x3 = float(vertexPos[2][0]);
    y1 = float(vertexPos[0][1]);
    y2 = float(vertexPos[1][1]);
    y3 = float(vertexPos[2][1]);
    Px = float(pos.x);
    Py = float(pos.y);


    float w1, w2, w3; // Barycentric weights
    w1 = ( (y2-y3)*(Px-x3) + (x3-x2)*(Py-y3) ) / ( (y2-y3)*(x1-x3) + (x3-x2)*(y1-y3) );
    w2 = ( (y3-y1)*(Px-x3) + (x1-x3)*(Py-y3) ) / ( (y2-y3)*(x1-x3) + (x3-x2)*(y1-y3) );
    w3 = 1.0f - w1 - w2;

    float wc1, wc2, wc3; // Weights clamped and normalized
    wc1 = max(w1,0.0f);
    wc2 = max(w2,0.0f);
    wc3 = max(w3,0.0f);

    wc1 = wc1/(wc1+wc2+wc3);
    wc2 = wc2/(wc1+wc2+wc3);
    wc3 = wc3/(wc1+wc2+wc3);

    //--------------
    // 2. CALCULATE THE RESULTING COLOR

    float4 result = vertexColors[0];

    /*  Maybe clamp the weights to 0... and then normalize the resulting  */
    if(numPoints <= 1){
      result = vertexColors[0];
    }else if(numPoints == 2){
      result = vertexColors[0]; //TODO: DO THIS WITH A LINEAR GRADIENT!!!!!!!!
    }else{
      result = vertexColors[0] * wc1 + vertexColors[1] * wc2 + vertexColors[2] * wc3;
    }

    // Write the result to the output image
    dst() = result;
  }
};