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.
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)