Triangle Voronoi Blender OSL

tri_voronoi.osl

/*
    OSL Triangle Voronoi Borders

    Based on a shadertoy by Thomas Hooper(tdhooper):
    https://www.shadertoy.com/view/ss3fW4
*/

vector fract(vector v) { return mod(v*.5,.5)*2.; }

vector vstep(vector a, vector b) {
    return vector(step(a.x,b.x),step(a.y,b.y),step(a.z,b.z));
}

vector sdTriEdges(vector p) {
    return vector(
        dot(p, vector(0.,-1.,0.)),
        dot(p, vector(.866025, .5,0.)),
        dot(p, vector(-.866025, .5,0.))
    );
}

float sdTri(vector p) {
    vector t = sdTriEdges(p);
    return max(t.x, max(t.y, t.z));
}

vector vstep(vector a, vector b) {
    return vector(step(a.x,b.x),step(a.y,b.y),step(a.z,b.z));
}

vector primaryAxis(vector p) {
    vector a = abs(p);
    return (1.-vstep(a, vector(a.y,a.z,a.x)))*vstep(vector(a.z,a.x,a.y), a)*sign(p);
}

vector sdgBorder(vector pt1, vector pt2) {
    vector tbRel = sdTriEdges(pt2 - pt1);    
    vector axis = primaryAxis(tbRel);
    
    vector gA = vector(0.,-1.,.0);
    vector gB = vector(.866025, .5, 0.);
    vector gC = vector(-.866025, .5, 0.);
    
    vector norA = gC * axis.x + gA * axis.y + gB * axis.z;
    vector norB = gB * -axis.x + gC * -axis.y + gA * -axis.z;
    
    vector dir = gA * axis.x + gB * axis.y + gC * axis.z;
    vector corner = dir * dot(dir, pt1 - pt2) * 2./3.;
            
    if (axis.x + axis.y + axis.z < 0.0) {
        corner = pt2 + corner;
        vector ca = corner + min(0., dot(corner, -norA)) * norA;
        vector cb = corner + max(0., dot(corner, -norB)) * norB;
        float side = step(dot(corner,vector( dot(dir,vector(0.,-1.,0.)), dot(dir,vector(1.,0.,0.)),0.0)), 0.);
        corner = mix(cb, ca, side);
    } else {
        corner = pt1 - corner;
        vector ca = corner + max(0., dot(corner, -norA)) * norA;
        vector cb = corner + min(0., dot(corner, -norB)) * norB;
        float side = step(dot(corner,vector( dot(dir,vector(0.,-1.,0.)), dot(dir,vector(1.,0.,0.)),0.0)), 0.);
        corner = mix(ca, cb, side);
    }
    
    vector nor = normalize(corner);
    float d = length(corner);

    return vector(abs(d), nor.x, nor.y);
}

// https://www.shadertoy.com/view/4djSRW
vector hash22(vector p)
{
	vector p3 = fract(vector(p.x,p.y,p.x) * vector(.1031, .1030, .0973));
    p3 += dot(p3, vector(p3.y,p3.z,p3.x)+33.33);
    return fract((vector(p3.x,p3.x,0.)+vector(p3.y,p3.z,0.))*vector(p3.z,p3.y,0.));
}

vector cellPoint(vector n, vector f, vector cell, float seed, float offset) {
    vector coord = n + cell;
    vector o = hash22( n + cell + vector(seed) );
    o = 0.5 + 0.5*sin(offset * 6.2831 + 6.2831 * o);
    vector pt = cell + vector(o.x,o.y,0.0) - f;
    return pt;
}

matrix voronoi(vector x, float seed, float offset)
{
    vector n = floor(x);
    vector f = fract(x);
    vector nor;

    //----------------------------------
    // first pass: regular voronoi
    //----------------------------------
	vector closestCell, closestPoint, id;

    int reach = 2;

    float closestDist = 8.0;
    for( int j = -reach; j <= reach; j++ )
    for( int i = -reach; i <= reach; i++ )
    {
        vector cell = vector(i, j, 0.0);
        vector pt = cellPoint(n, f, cell, seed,offset);
        float dist = sdTri(pt);

        if( dist < closestDist )
        {
            closestDist = dist;
            closestPoint = pt;
            closestCell = cell;
            id = n + cell;
        }
    }

    //----------------------------------
    // second pass: distance to borders
    //----------------------------------
    closestDist = 8.0;
    for( int j = -reach-1; j <= reach+1; j++ )
    for( int i = -reach-1; i <= reach+1; i++ )
    {
        vector cell = closestCell + vector(i, j, 0.);
        vector pt = cellPoint(n, f, cell, seed, offset);

        float dist = sdTri(closestPoint - pt);

        if( dist > 0.00001 ) {
            vector sdg = sdgBorder(closestPoint, pt);
            if (sdg.x < closestDist) {
                closestDist = sdg.x;
                nor = vector(sdg.y,sdg.z,0.);
            }
        }
    }
    return matrix(closestDist,closestPoint.x,closestPoint.y,nor.x,nor.y,id.x,id.y,0.,0.,0.,0.,0.,0.,0.,0.,0.);
}

shader triangle_voronoi(
    point Vector = P,
    float Seed = 0.,
    float Offset = 0.,
    output float Distance = 0.,
    output vector Coords = 0.,
    output vector Normal = 0.,
    output vector ID = 0.
){
    matrix vor = voronoi(Vector, Seed, Offset);
    Distance = vor[0][0];
    Coords = vector(vor[0][1],vor[0][2],0.);
    Normal = vector(vor[0][3],vor[1][0],0.);
    ID = vector(vor[1][1],vor[1][2],0.);
}