Code for comment: http://www.entagma.com/analytical-foam/#comment-40052

deform.vex

// deform code
vector4 planes[] = point(1, "planes", @id);

foreach(vector4 plane; planes) {
    vector n = set(plane.x, plane.y, plane.z);
    float d = dot(n, @P) + plane.w;
    // check if point is "above" intersection plane
    // since the plane normal always points away from each sphere's center
    // this will select all points intersecting a neighbor sphere
    if (d > 0) {
        // project back onto plane
        @P -= n * d;
    }
}

gather.vex

// gather code
float rp = @pscale;
int neighbors[] = nearpoints(0, @P, chf("../set_pscale/max") * 2);
vector4 planes[];

for(int i=1, n = len(neighbors); i<n; i++) {
    vector q = point(0, "P", neighbors[i]);
    float rq = point(0, "pscale", neighbors[i]);
    vector dir = q - @P;
    float dist = length(dir);
    if (dist < rp + rq) {
        dir = normalize(dir);
        float w = (pow(rp, 2) - pow(rq, 2) + pow(dist, 2)) / (2 * dist);
        // point on intersection plane
        vector isec = @P + dir * w;
        // store plane in normal form
        push(planes, set(dir.x, dir.y, dir.z, -dot(dir, isec)));
    }
}
p[]@planes = planes;
i@id = @ptnum;