Doesn’t Houdini already have a PolyBevel?

Houdini's PolyBevel has had a weird behaviour for the longest time. When beveling points (on polylines), and set to Round, the beveled arcs are not circular. If we put aside the various justifications for why that is, the question is then, can we make it circular?

Here's what I've managed to achieve, in a single Wrangle running over Primitives.

Wrangle Code

1// Parameters
2int usePscale = chi("use_pscale");
3string pscaleAttrib = chs("attribute");
4float offset = chf("offset");
5float res = chf("resolution");
6
7// Check if prim is open or closed
8int isClosed = primintrinsic(0, "closed", @primnum);
9
10int pts[] = primpoints(0, @primnum);
11int npts = len(pts);
12
13// Set the output prim type
14int rebuiltPrim;
15if (isClosed)
16    rebuiltPrim = addprim(0, "poly");
17else{
18    rebuiltPrim = addprim(0, "polyline");
19    addvertex(0, rebuiltPrim, pts[0]);
20    //setpointgroup(0, "transition", pts[0], 1);
21}
22
23// Loop over each point in prim and create fillets
24foreach (int index; int pt; pts){
25    // Skip endpoints if curve is open
26    if (!isClosed){
27        if (index == 0)
28            continue;
29        if (index == npts - 1)
30            continue;
31    }
32    int prevPt = pts[(index - 1) % npts];
33    int nextPt = pts[(index + 1) % npts];
34    
35    vector prevP = point(0, "P", prevPt);
36    vector nextP = point(0, "P", nextPt);
37    vector ptP = point(0, "P", pt);
38    
39    vector fromPrev = normalize(ptP - prevP);
40    vector toNext = normalize(nextP - ptP);    
41    vector up = normalize(cross(toNext, fromPrev));
42    vector right = normalize(cross(fromPrev, up));
43    
44    float signedAngleRad = atan2(dot(cross(fromPrev, toNext), up),
45                                dot(fromPrev, toNext));
46    //float signedAngleDeg = degrees(signedAngleRad);
47    //float arcAngle = signedAngleDeg;
48    
49    // Clamp offsets
50    float pscale = point(0, pscaleAttrib, pt);
51    if (usePscale)
52        offset *= pscale;
53    float distFromPrev = distance(ptP, prevP);
54    float distFromNext = distance(ptP, nextP);
55    float limit = min(distFromPrev, distFromNext) * 0.5;
56    offset = clamp(offset, 0, limit);
57    
58    // Start of fillet
59    vector firstP = ptP - fromPrev * offset;
60    
61    // Compute fillet radius
62    float interiorAngleRad = PI - signedAngleRad;
63    float hypo = offset / cos(interiorAngleRad * 0.5);
64    float radius = sqrt(hypo * hypo - offset * offset);
65    
66    // Find arc's center
67    vector arcCenterP = firstP + right * radius;
68    
69    // Compute number of points to generate
70    float arcLen = abs(signedAngleRad * radius);
71    int steps = floor(arcLen / res);
72    
73    float angleStep = signedAngleRad / steps;
74    
75    // Initialize relative offset from arc center
76    vector offsetFromArcCenter = firstP - arcCenterP;
77    
78    // Generate fillet points
79    for (int i = 0; i < steps + 1; i++){
80        matrix m = ident();
81        rotate(m, angleStep * i, up);
82        vector newOffset = offsetFromArcCenter * m;
83        int nextFilletPt = addpoint(0, arcCenterP + newOffset);
84        addvertex(0, rebuiltPrim, nextFilletPt);
85        if (i == 0)
86            setpointgroup(0, "transition", nextFilletPt, 1);
87        if (i == steps)
88            setpointgroup(0, "transition", nextFilletPt, 1);
89    }
90}
91// Add last point for open curves
92if (!isClosed){
93    addvertex(0, rebuiltPrim, pts[-1]);
94    //setpointgroup(0, "transition", pts[-1], 1);
95}
96// Remove original prim, use Attribute Transfer after to preserve attrib data
97removeprim(0, @primnum, 1);
98// Remove outgoing pscale
99removepointattrib(0, "pscale");

What a monster!

In the next section, I’ll try to break it down.

Breakdown

1. primintrinsic(0, "closed", @primnum) is used to handle open/closed curves appropriately:
   1. Endpoints of an open curve will be ignored
   2. All points on a closed curve will be beveled
2. primpoints() creates an array of points for each primitive, sorted by vertex order. This lets us retrieve “previous” and “next” point positions when computing directions.

int prevPt = pts[(index - 1) % npts];
int nextPt = pts[(index + 1) % npts];

vector prevP = point(0, "P", prevPt);
vector nextP = point(0, "P", nextPt);

‍

3. For each point that we want to bevel, I normalize its direction from the previous point (fromPrev), as well as the direction to the next point (toNext), then create a new vector from their cross product. This new vector will be an “axis of rotation” (up) when constructing the bevel.

vector ptP = point(0, "P", pt);
vector fromPrev = normalize(ptP - prevP);
vector toNext = normalize(nextP - ptP);
vector up = normalize(cross(toNext, fromPrev));

‍

4. Clamp the offset being applied to prevent overlaps, based on the distance to the closer neighbouring point.

float distFromPrev = distance(ptP, prevP);
float distFromNext = distance(ptP, nextP);
float limit = min(distFromPrev, distFromNext) * 0.5;
offset = clamp(offset, 0, limit);

‍

5. The radius of the arc is computed using some Pythagoras’, and the angle between fromPrev and toNext, using atan2(dot(cross(a, b), axis), dot(a, b)), which gives higher precision than acos(dot(a, b)), as well as the sign of the angle. Having the sign is not strictly necessary here, since up would just point in the opposite direction, but it does help in debugging.

float signedAngleRad = atan2(dot(cross(fromPrev, toNext), up),
                                dot(fromPrev, toNext));
                                
// Compute fillet radius
float interiorAngleRad = PI - signedAngleRad;
float hypo = offset / cos(interiorAngleRad * 0.5);
float radius = sqrt(hypo * hypo - offset * offset);

‍

At this point, I think it’s worth noting that this technique does not create fillets of constant radii, but rather, each fillet starts at the same distance from the original point. This means that tight angles will have smaller arcs, and vice versa.

‍

6. Using fromPrev and up, I can construct another orthogonal vector right, and use it to find the center of the arc.

// Find arc's center
vector arcCenterP = firstP + right * radius;

‍

7. Unlike PolyBevel, I prefer having a consistent density of points, rather than a fixed point count along the bevel. The number of points being created is thus the arc’s length divided by the Resolution parameter res.

// Compute number of points to generate
float arcLen = abs(signedAngleRad * radius);
int steps = floor(arcLen / res);

float angleStep = signedAngleRad / steps;

‍

8. To construct the actual bevel, I rotate the offsetFromArcCenter vector around the up axis at each step, and insert a point/vertex.

// Initialize relative offset from arc center
vector offsetFromArcCenter = firstP - arcCenterP;

// Generate fillet points
for (int i = 0; i < steps + 1; i++){
    matrix m = ident();
    rotate(m, angleStep * i, up);
    vector newOffset = offsetFromArcCenter * m;
    int nextFilletPt = addpoint(0, arcCenterP + newOffset);
    addvertex(0, rebuiltPrim, nextFilletPt);
}

‍

Extras

transition

What is this “transition” group I’m creating?

A typical use case for beveling points is to round sharp corners, then sweep along them to produce piping/ducts. This group marks the start/end of each fillet, so objects can be copied to those points to emulate joints.

pscale

Based on a suggestion from Probiner, I’ve also extended support to use a point attribute as a scale multiplier for the size of the bevel. You can specify which point attribute to use.

The existing limit of “half of shorter neighbouring edge” is still in effect. I’ll leave it as an exercise for the reader to implement a limit based on adjacent bevel sizes.

Attributes

Unfortunately, in the process of constructing the bevels, the original primitive is completely reconstructed. This means that groups and attributes are lost. While it is possible to use detailintrinsic() to retrieve and reassign them, remember that Group/Attribute Transfer/Copy exist as native SOPs, and don’t require as much maintenance (on my part).

‍