code.c

// choose initial point indices and radius
int O_ind = 3;
int P_ind = 2;
int Q_ind = 6;
float r = max(chf("r"), 0.05);

// get their positions
vector O = point(0, "P", O_ind);
vector P = point(0, "P", P_ind);
vector Q = point(0, "P", Q_ind);

// get M
float phi = acos(dot(normalize(O-P), normalize(Q-P)));
float t = r / sin(phi) / 2.;
vector M = P + t*(normalize(O-P) + normalize(Q-P));
int M_ind = addpoint(0, M);

// get projections from M to PO and PQ (M1 and M2)
vector T1 = O - P;
vector T2 = Q - P;
float t1 = dot(T1, M - P) / dot(T1, T1);
float t2 = dot(T2, M - P) / dot(T2, T2);
vector M1 = P + t1*T1;
vector M2 = P + t2*T2;

// initialize output polyline (point O, then bevel points, then point Q)
int out_poly = addprim(0, "polyline");

addvertex(0, out_poly, O_ind); // point O (you can ommit this if you only want the bevel)

// get vectors connecting M and its projections (called M1 and M2) to PO and PQ
vector MM1 = M1 - M;
vector MM2 = M2 - M;
// we need a vector perpendicular to those two vectors
vector perp = normalize(cross(normalize(MM1), normalize(MM2)));
// initialisation of rotation parameters
vector beta = acos(dot(normalize(MM1), normalize(MM2)));
matrix3 mRot = ident();
float inc = 0.01 * 1/r; // this is arbitrary
inc = beta/int(beta/inc);
// for each bevel segment (angle) we rotate MM1 around "perp" vector
// and add a point/vertex to the output polyline
for(float angle = 0; angle < beta; angle += inc){
    mRot = ident();
    rotate(mRot, angle, perp);
    vector MX = MM1 * mRot;
    vector MMX = M + MX;
    int MX_ind = addpoint(0, MMX);
    addvertex(0, out_poly, MX_ind); // bevel points
}
addvertex(0, out_poly, Q_ind); // point Q to finish the curve (you can ommit this if you only want the bevel)
removepoint(0, M_ind);

setdetailattrib(0, "phi", phi);
setdetailattrib(0, "beta", beta);