I do believe that it is already invariant to linear transformations the way you want, i.e. we can evaluate the corners of an arbitrary parallelogram instead of a square and get a similar coverage estimate.

Similar maybe but it can't be the same surely? Just pick some function like f(x,y) = x-1 and start rotating it around your centre pixel, the average (s1+s2+s3+s4) will be the same (since it's a linear function) but there's no way those absolute values will remain constant.

You should be pretty close though. For a linear function you can just calculate the distance to the 0 line, which is invariant to any linear transformation that leaves that line where it is (which is what you want). This is just the function value divided by the norm of the gradient. Both of which you can estimate from those 4 points. This gives something like

    dx = (s2 - s1 + s4 - s3)
    dy = (s3 - s1 + s4 - s2)
    f  = (s1+s2+s3+s4)/4
    dist = f / sqrt(dx*dx + dy*dy)

My function approximate coverage of a square pixel, so indeed if you rotate a line around it at a certain distance that line will clip the corners at some angles and be clear of the pixel at other angles.