// Given a distance squared (dist2), say the distance to a line, and 
// the radius of influence squared (rad2) the function returns a Wyvill
// field value in the range 0 to 1.
// See: http://www.fundza.com/algorithmic/quadtree_blobby/index.html
float fieldValue(float dist2, float rad2) {
    if(dist2 >= rad2)
        return 0.0;
    float d4 = dist2 * dist2;
    float d6 = dist2 * d4;
    float r4 = rad2 * rad2;
    float r6 = rad2 * r4;
    return (-0.44444 * d6/r6 + 1.88889 * d4/r4 + -2.44444 * dist2/rad2 + 1);
    }
  
// Returns the distance from a point and a line.
// See: http://fundza.com/vectors/point2line/index.html
float distToLine(point p, point start, point end) {
    vector line_vec = start - end;
    vector pnt_vec = start - p;
    float  line_len = length(line_vec);
    vector line_unitvec = normalize(line_vec);
    vector pnt_vec_scaled = pnt_vec * (1.0/line_len);
    float  delta_t = dot(line_unitvec, pnt_vec_scaled);
    if(delta_t < 0)
        delta_t = 0.0;
    else if(delta_t > 1.0)
        delta_t = 1.0;
    point nearest = line_vec * delta_t;
    return distance(nearest, pnt_vec);
    }
    
#define MAX_VERTICES 2000
  
// Because arrays in OSL cannot be re-sized the shader declares
// a point array of a fixed arbitary length. Although a user attribute
// called "verts_data" references MAX_VERTICES of data another user
// attribure called "verts_count" specifies the number of valid
// vertices.
shader
WyvillLines(float mask_radius = 0.1,
            float  wyvill_radius = 0.2,
            int invertMask = 0,
            output float resultMask = 0,
            output float resultF = 0)
{
int len, accumulated_mask = 0;
float accumulated_fieldvalue = 0;
float wyvill_rad_rad = wyvill_radius * wyvill_radius;
point data[MAX_VERTICES];
  
int result = getattribute("user:verts_count", len);
  
if(result && len <= MAX_VERTICES) {
    result = getattribute("user:verts_data", data);
    if(result) {
        point p = transform("world", P);
        for(int n = 0; n < len; n += 2) {
            float dist = distToLine(p, data[n], data[n+1]);
            if(dist < mask_radius)
                accumulated_mask++;
            accumulated_fieldvalue += fieldValue(dist * dist, wyvill_rad_rad); 
            }
        }
    }
if(invertMask)
    resultMask = (accumulated_mask > 0) ? 1 : 0;
else
    resultMask = (accumulated_mask > 0) ? 0 : 1;
resultF = accumulated_fieldvalue;
}