CSG operations on implicit surfaces with marching cubes [SOLVED]

Question

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 circles, and the set operation difference. i.e A - B. One circle 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 circle
                  float (*cmp2)(const vec2&) //distance from moving circle
)
{
  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 circle1
      if((s2 < iso)){ //if inside circle2
        scalar[i] = s2; //then set the scalar to the moving circle
      } else {
        scalar[i] = s1; //only inside circle1
        squareindex |= (1<<i); //mark as inside
      }
    }
    else {
      scalar[i] = s1; //inside neither circle
    }
  }

  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

EDIT: Basically, my implementation of marching squares works fine. It is my scalar field which is broken, and I wonder what the correct way would look like. Preferably I'm looking for a general approach to implement the three set operations I discussed above, for the usual primitives (circle, rectangle/square, plane)

EDIT 2: Here are some new images after implementing the answerer's whitepaper:

1.Difference
2.Intersection
3.Union

EDIT 3: I implemented this in 3D too, with proper shading/lighting:

1.Difference between a greater sphere and a smaller sphere
2.Difference between a greater sphere and a smaller sphere in the center, clipped by two planes on both sides, and then union with a sphere in the center.
3.Union between two cylinders.

Answer 1

This is not how you mix the scalar fields. Your scalars say one thing, but your flags whether you are inside or not say another. First merge the fields, then render as if you were doing a single compound object:

for(int i=0,j=0; i<4; ++i, j+=2){
  s1 = cmp1(verts[j]);
  s2 = cmp2(verts[j]);
  s = max(s1, iso-s2); // This is the secret sauce
  if(s < iso) { // inside circle1, but not inside circle2
    squareindex |= (1<<i);
  }
  scalar[i] = s;
}

This article might be helpful: Combining CSG modeling with soft blending using Lipschitz-based implicit surfaces.