I render isosurfaces with marching cubes, (or perhaps marching squares as this is 2D) and I want to do set operations like set difference, intersection and union. I thought this was easy to implement, by simply choosing between two vertex scalars from two different implicit surfaces, but it is not.

For my initial testing, I tried with two spheres, and the set operation difference. i.e A - B. One sphere is moving and the other one is stationary. Here's the approach I tried when picking vertex scalars and when classifying corner vertices as inside or outside. The code is written in C++. OpenGL is used for rendering, but that's not important. Normal rendering without any CSG operations does give the expected result.

       void march(const vec2& cmin, //min x and y for the grid cell
                  const vec2& cmax, //max x and y for the grid cell
                  std::vector<vec2>& tri, 
                  float iso,
                  float (*cmp1)(const vec2&), //distance from stationary sphere
                  float (*cmp2)(const vec2&) //distance from moving sphere
)
{
  unsigned int squareindex = 0;
  float scalar[4];
  vec2 verts[8];
  /* initial setup of the grid cell */
  verts[0] = vec2(cmax.x, cmax.y);
  verts[2] = vec2(cmin.x, cmax.y);
  verts[4] = vec2(cmin.x, cmin.y);
  verts[6] = vec2(cmax.x, cmin.y);

  float s1,s2;
  /**********************************
   ********For-loop of interest******
   *******Set difference between ****
   *******two implicit surfaces******
   **********************************/
  for(int i=0,j=0; i<4; ++i, j+=2){
    s1 = cmp1(verts[j]);
    s2 = cmp2(verts[j]);
    if((s1 < iso)){ //if inside sphere1
      if((s2 < iso)){ //if inside sphere2
        scalar[i] = s2; //then set the scalar to the moving sphere
      } else {
        scalar[i] = s1; //only inside sphere1
        squareindex |= (1<<i); //mark as inside
      }
    }
    else {
      scalar[i] = s1; //inside neither sphere
    }
  }

  if(squareindex == 0)
    return;
  /* Usual interpolation between edge points to compute
     the new intersection points */
  verts[1] = mix(iso, verts[0], verts[2], scalar[0], scalar[1]);
  verts[3] = mix(iso, verts[2], verts[4], scalar[1], scalar[2]);
  verts[5] = mix(iso, verts[4], verts[6], scalar[2], scalar[3]);
  verts[7] = mix(iso, verts[6], verts[0], scalar[3], scalar[0]);

  for(int i=0; i<10; ++i){ //10 = maxmimum 3 triangles, + one end token
    int index = triTable[squareindex][i]; //look up our indices for triangulation
    if(index == -1)
      break;
    tri.push_back(verts[index]);
  }
}

This gives me weird jaggies:
It looks like the CSG operation is done without interpolation. It just "discards" the whole triangle. Do I need to interpolate in some other way, or combine the vertex scalar values? I'd love some help with this. A full testcase can be downloaded HERE