# Bezier offset self intersections

An "offset", or parallel curve, is "a curve whose points are at a fixed normal distance from a given curve". It might happen that the offset curve intersect itself, as shown here (the most inner green curve). This kind of "loop" is called "local" and appears when the curvature of given curve is greater than $1/d$ (i.e. the offset distance).

I would like to compute the approximate cubic Bezier offset curve/curves to a given cubic Bezier curve ($B$) but the local loops are giving me troubles. I though I would split the original curve at parameters ($t$), where the curvature radius ($r$) falls below the given offset distance.

The signed curvature radius can be calculated as:

$$r = \frac{\|B'(t)\|^3}{B'(t) \times B''(t)}$$

The problem I have is I don't know how to find those parameter values. I thought I would somehow get some nice polynomial whose roots are the location along the original curve where the splits would occur. But I cannot figure out how to derive this polynomial from the curvature formula or whether it is even possible (there is a polynomial in denominator)?

Could someone, please, give me a helping hand here? My question is: where to trim the original curve (into possibly more parts) so that the curvature radius would be at least the offset distance everywhere?

Alternatively, there is a paper "Error Bounded Variable Distance Offset Operator for Free Form Curves and Surfaces" written by Elber and Cohen. In the chapter 4, it describes a way of trimming off the local loops by usage of tangents but unfortunately I don't really understand it. Perhaps someone could cast some light into this method instead?

• See stackoverflow.com/questions/4148831/… for a related discussion. In general this is a difficult problem and the only way to check the effectiveness of the algorithms in the literature is to try them out. Personally, I have never had to deal with cusps. – Biswajit Banerjee Oct 3 '16 at 20:54
• – Biswajit Banerjee Oct 3 '16 at 20:58
• @BiswajitBanerjee I'm not asking about offsetting Beziers, rather I would like to know where to find the parameter values at which the curvature radius drops below given offset distance. – Faaf Oct 3 '16 at 21:43
• One can try a Newton method to solve for $t$ but that can be problematic unless a good starting value is chosen. I suggest you try out the approach described in the appendix in Elber and Elaine Cohen's paper. But that approach is only for local self-intersections and not global ones. Will provide a more detailed answer if time permits. – Biswajit Banerjee Oct 4 '16 at 4:39
• @BiswajitBanerjee Yes, I'm only after local self-intersections. I'm trying to read the Appendix but it is bit over my head so I would be very thankful for any kind of explanation of it! – Faaf Oct 4 '16 at 13:39

Let the control points of the Bezier curve be $\mathbf{p}_1$, $\mathbf{p}_2$, $\mathbf{p}_3$, and $\mathbf{p}_4$. Then the parametrized Bezier curve is $$\mathbf{b}(t) = (1 - t)^3 \mathbf{p}_1 + 3 t(1 - t)^2 \mathbf{p}_2 + 3 t^2 (1 - t) \mathbf{p}_3 + t^3 \mathbf{p}_4 \,.$$ The slope (tangent) of $\mathbf{b}(t)$ is \begin{align} \mathbf{t}(t) & = \mathbf{b}^{'}(t) \\ &= - 3(1 - t)^2 \mathbf{p}_1 + 3 [(1 - t)^2 - 2 t(1 - t)] \mathbf{p}_2 + 3 [2t(1 - t) - t^2] \mathbf{p}_3 + 3 t^2 \mathbf{p}_4 \\ & = - 3(1 - t)^2 \mathbf{p}_1 + 3 (1 - t)(1 - 3t) \mathbf{p}_2 + 3 (2t - 3t^2) \mathbf{p}_3 + 3 t^2 \mathbf{p}_4 \,. \end{align} The normal to $\mathbf{b}(t)$ can be computed using the cross product of the tangent with the $\mathbf{e}_z$ vector where $\mathbf{e}_z = (0, 0, 1)$ and assuming $\mathbf{p}_i = (x_i, y_i, 0)$. Then \begin{align} \mathbf{n}(t) & = \mathbf{t}(t) \times \mathbf{e}_z = (t_2z_3-t_3z_2)\mathbf{e}_x-(t_1z_3-t_3z_1)\mathbf{e}_y+(t_1z_2-t_2z_1)\mathbf{e}_z \\ & = t_2\,\mathbf{e}_x - t_1\,\mathbf{e}_y + 0\,\mathbf{e}_z \\ \end{align} Note that $\lVert\mathbf{n}(t)\rVert_{} = \lVert\mathbf{t}(t)\rVert_{}$.

The convexity of $\mathbf{b}(t)$ is \begin{align} \mathbf{c}(t) = \mathbf{b}^{''}(t) & = 6 (1 - t) \mathbf{p}_1 + 3 [-(1 - 3t) - 3(1 - t)] \mathbf{p}_2 + 3 (2 - 6 t) \mathbf{p}_3 + 6 t \mathbf{p}_4 \\ & = 6 (1 - t) \mathbf{p}_1 - 6 (2 - 3t) \mathbf{p}_2 + 6 (1 - 3 t) \mathbf{p}_3 + 6 t \mathbf{p}_4 \\ \end{align} The offset curve is defined as $$\mathbf{b}_d(t) = \mathbf{b}(t) + \frac{\mathbf{n}(t)}{\lVert\mathbf{n}(t)\rVert_{}} \,d =: \mathbf{b}(t) + \widehat{\mathbf{n}}(t)\,d$$ where $d$ is the offset distance. The tangent to the offset curve is $$\mathbf{t}_d(t) = \mathbf{b}^{'}(t) + \cfrac{d \widehat{\mathbf{n}}}{d t}\,d = \mathbf{t}(t) + \cfrac{d \widehat{\mathbf{n}}}{d t}\,d\,.$$ The derivative of the unit normal is $$\cfrac{d \widehat{\mathbf{n}}}{d t} = \cfrac{d }{d t}\left(\frac{\mathbf{n}(t)}{\lVert\mathbf{n}(t)\rVert_{}}\right) = \frac{1}{\lVert\mathbf{n}(t)\rVert_{}} \cfrac{d \mathbf{n}}{d t} + \cfrac{d }{d t}\left(\frac{1}{\lVert\mathbf{n}(t)\rVert_{}}\right)\mathbf{n}(t)\,.$$ Note that $$\cfrac{d \mathbf{n}}{d t} = \cfrac{d }{d t}[\mathbf{t}(t) \times \mathbf{e}_z] = \cfrac{d \mathbf{t}}{d t} \times \mathbf{e}_z = \mathbf{b}^{''}(t) \times \mathbf{e}_z = \mathbf{c}(t) \times \mathbf{e}_z = c_2\,\mathbf{e}_x - c_1\,\mathbf{e}_y\,.$$ To find the derivative of the inverse of the norm of $\mathbf{n}(t)$, we note that $$\frac{1}{\lVert\mathbf{n}(t)\rVert_{}} = \frac{1}{\sqrt{n_j n_j}} = (n_j n_j)^{-1/2} \,.$$ Therefore, $$\cfrac{d }{d t}\left(\frac{1}{\lVert\mathbf{n}(t)\rVert_{}}\right) = -\frac{1}{2} (n_j n_j)^{-3/2} \left(2n_j \cfrac{d n_j}{d t}\right) = - \frac{1}{\lVert\mathbf{n}(t)\rVert_{}^3} \left(\mathbf{n}(t) \cdot \cfrac{d \mathbf{n}}{d t}\right)$$ and we have $$\cfrac{d \widehat{\mathbf{n}}}{d t} = \frac{1}{\lVert\mathbf{n}(t)\rVert_{}} \cfrac{d \mathbf{n}}{d t} - \frac{1}{\lVert\mathbf{n}(t)\rVert_{}^3} \left(\mathbf{n}(t) \cdot \cfrac{d \mathbf{n}}{d t}\right) \mathbf{n}(t) = \frac{1}{\lVert\mathbf{n}(t)\rVert_{}} \cfrac{d \mathbf{n}}{d t} - \frac{1}{\lVert\mathbf{n}(t)\rVert_{}^3} \left[\mathbf{n}(t) \otimes \mathbf{n}(t)\right] \cdot \cfrac{d \mathbf{n}}{d t}$$ or $$\cfrac{d \widehat{\mathbf{n}}}{d t} = \frac{1}{\lVert\mathbf{n}(t)\rVert_{}} \left[\boldsymbol{I} - \widehat{\mathbf{n}}(t) \otimes \widehat{\mathbf{n}}(t)\right] \cdot \cfrac{d \mathbf{n}}{d t}$$ where $\boldsymbol{I}$ is the identity matrix. We now have an algebraic expression for the offset tangent: $$\mathbf{t}_d(t) = \mathbf{t}(t) + \frac{d}{\lVert\mathbf{n}(t)\rVert_{}} \left[\boldsymbol{I} - \widehat{\mathbf{n}}(t) \otimes \widehat{\mathbf{n}}(t)\right] \cdot [\mathbf{c}(t) \times \mathbf{e}_z]$$ Define $\mathbf{m}(t) := \mathbf{c}(t) \times \mathbf{e}_z$ to get $$\mathbf{t}_d(t) = \mathbf{t}(t) + \frac{d}{\lVert\mathbf{n}(t)\rVert_{}} \left[\boldsymbol{I} - \widehat{\mathbf{n}}(t) \otimes \widehat{\mathbf{n}}(t)\right] \cdot \mathbf{m}(t) \,.$$

The Elber-Cohen approach uses the sign of $\mathbf{t}(t) \cdot \mathbf{t}_d(t)$ to find the locations of the cusps. This computation has to be done numerically.

Once the cusps have been found, the two segments are intersected numerically to find the crossing point.

## Curvature

Warning The algebra below needs to be checked for correctness.

The curvature of the Bezier curve is $$\frac{1}{d}\,\mathbf{e}_z = \kappa(t)\, \mathbf{e}_z = \frac{\mathbf{t}(t) \times \mathbf{c}(t)}{\lVert\mathbf{t}(t)\rVert_{}^3} \,.$$ Therefore, $$d = \frac{\lVert\mathbf{t}(t)\rVert_{}^3}{\left[\mathbf{t}(t) \times \mathbf{c}(t)\right]\cdot \mathbf{e}_z} = \frac{\lVert\mathbf{t}(t)\rVert_{}^3}{\mathbf{t}(t) \cdot \left[\mathbf{c}(t) \times \mathbf{e}_z\right]} = \frac{\lVert\mathbf{t}(t)\rVert_{}^3}{\mathbf{t}(t) \cdot \mathbf{m}(t)} \,.$$ Plugging the above relation for $d$ into the expression for $\mathbf{t}_d$, we have \begin{align} \mathbf{t}_d(t) & = \mathbf{t}(t) + \frac{\lVert\mathbf{t}(t)\rVert_{}^3}{\lVert\mathbf{n}(t)\rVert_{}\left[\mathbf{t}(t) \cdot \mathbf{m}(t)\right]} \left[\boldsymbol{I} - \widehat{\mathbf{n}}(t) \otimes \widehat{\mathbf{n}}(t)\right] \cdot \mathbf{m}(t) \\ & = \frac{\lVert\mathbf{n}\rVert_{} (\mathbf{t} \cdot \mathbf{m}) \,\mathbf{t} + \lVert\mathbf{t}\rVert_{}^3 \left[\boldsymbol{I} - \widehat{\mathbf{n}}\otimes\widehat{\mathbf{n}}\right] \cdot \mathbf{m}} {\lVert\mathbf{n}\rVert_{} (\mathbf{t} \cdot \mathbf{m})} \,. \end{align} Using $\lVert\mathbf{n}\rVert_{} = \lVert\mathbf{t}\rVert_{}$, we get $$\mathbf{t}_d = \frac{\lVert\mathbf{t}\rVert_{}^3 (\mathbf{t} \cdot \mathbf{m}) \,\mathbf{t} + \lVert\mathbf{t}\rVert_{}^3 \left[\lVert\mathbf{t}\rVert_{}^2\,\boldsymbol{I} - \mathbf{n} \otimes \mathbf{n}\right] \cdot \mathbf{m}} {\lVert\mathbf{t}\rVert_{}^3 (\mathbf{t} \cdot \mathbf{m})} = \frac{ (\mathbf{t} \cdot \mathbf{m}) \,\mathbf{t} + \left[\lVert\mathbf{t}\rVert_{}^2\,\boldsymbol{I} - \mathbf{n} \otimes \mathbf{n}\right] \cdot \mathbf{m}} { \mathbf{t} \cdot \mathbf{m}} \,.$$

The Elber-Cohen approach uses the equation $$\mathbf{t}(t) \cdot \mathbf{t}_d(t) = 0$$ to find the values of $t$ where trimming is needed. Thus, \begin{align} \mathbf{t}(t) \cdot \mathbf{t}_d(t) &= \frac{ (\mathbf{t} \cdot \mathbf{m}) \,(\mathbf{t} \cdot \mathbf{t}) + \left[\lVert\mathbf{t}\rVert_{}^2\,\boldsymbol{I} - \mathbf{n} \otimes \mathbf{n}\right] : (\mathbf{t} \otimes \mathbf{m})} {\mathbf{t} \cdot \mathbf{m}} \\ & =\frac{\lVert\mathbf{t}\rVert_{}^2 (\mathbf{t} \cdot \mathbf{m}) + \left[\lVert\mathbf{t}\rVert_{}^2\,\boldsymbol{I} - \mathbf{n} \otimes \mathbf{n}\right] : (\mathbf{t} \otimes \mathbf{m})} {\mathbf{t} \cdot \mathbf{m}} \\ & =\frac{2\lVert\mathbf{t}\rVert_{}^2 (\mathbf{t} \cdot \mathbf{m}) - (\mathbf{n} \cdot \mathbf{t}) (\mathbf{n} \cdot \mathbf{m})} {\mathbf{t} \cdot \mathbf{m}} \\ & = 2\lVert\mathbf{t}\rVert_{}^2 = 0 \end{align} where we have used $\mathbf{n} \cdot \mathbf{t} = 0$. The quantity $\lVert\mathbf{t}\rVert$ is never zero and hence the expression for $d$ cannot be used to solve the above quartic equation for $t$.

Instead, Elber and Cohen use the expression $$\mathbf{t} \cdot \mathbf{t}_d = \mathbf{t}\cdot \mathbf{t} + \frac{d}{\lVert\mathbf{t}\rVert}(\mathbf{t} \cdot \mathbf{m})$$ to find the sign of the vector product as $$\text{sign}[\mathbf{t} \cdot \mathbf{t}_d] = \text{sign}\left[1 + \frac{d}{\lVert\mathbf{t}\rVert^3}(\mathbf{t} \cdot \mathbf{m}]\right] = \text{sign}\left[1 + \kappa(t) d\right] \,.$$ This expression is used to justify the use of change in sign of the dot product of the two tangents in determining the cusps in the offset curve.

## R-Code

require("ggplot2")

# Install spatstat for polyline intersections
if (!require(spatstat)) {
install.packages("spatstat")
library(spatstat)
}

setwd(".")

Bezier <- function(t, x1, x2, x3, x4) {
B <- (1.0 - t)^3*x1 + 3.0*t*(1- t)^2*x2 + 3.0*t^2*(1-t)*x3 + t^3*x4
#print(paste("t = ", t, " x1 = ", x1, " B = ", B))
return(B)
}

BezierTangent <- function(t, x1, x2, x3, x4) {
Bp <- -3*(1 - t)^2*x1 + 3*(1 - t)*(1 - 3*t)*x2 + 3*(2*t - 3*t^2)*x3 + 3*t^2*x4
normBp <- sqrt(Bp[1]^2 + Bp[2]^2)
Bphat <- Bp/normBp
return(Bphat)
}

BezierNormal <- function(t, x1, x2, x3, x4) {
tangent <- BezierTangent(t, x1, x2, x3, x4)
normal <- c(tangent[2], -tangent[1])
return(normal)
}

BezierConvex <- function(t, x1, x2, x3, x4) {
Bpp <- 6*(1 - t)*x1 - 6*(2 - 3*t)*x2 + 6*(1 - 3*t)*x3 + 6*t*x4
normBpp <- sqrt(Bpp[1]^2 + Bpp[2]^2)
Bpphat <- Bpp/normBpp
return(Bpphat)
}

BezierOffset <- function(B, N, d) {
Bd = B + N*d
return(Bd)
}

BezierOffsetTangent <- function(t, x1, x2, x3, x4, d) {

# Compute tangent
tvec <- -3*(1 - t)^2*x1 + 3*(1 - t)*(1 - 3*t)*x2 + 3*(2*t - 3*t^2)*x3 + 3*t^2*x4

# Compute normal
nvec <- c(tvec[2], -tvec[1])

# Compute norm of normal
norm_nvec <- sqrt(nvec[1]^2 + nvec[2]^2)

# Compute nhat
nhat <- nvec/norm_nvec

# Compute convexity
cvec <- 6*(1 - t)*x1 - 6*(2 - 3*t)*x2 + 6*(1 - 3*t)*x3 + 6*t*x4

# Compute dn_dt
dn_dt <- c(cvec[2], -cvec[1])

# Compute nhat otimes nhat
nhat.nhat = matrix(c(nhat[1]*nhat[1], nhat[1]*nhat[2],
nhat[2]*nhat[1], nhat[2]*nhat[2]),
nrow = 2, byrow = TRUE)

# Compute (nhat o nhat).dn_dt
nhat.nhat.dn_dt = nhat.nhat %*% dn_dt

# Compute I.dn_dt
II = matrix(c(1, 0, 0, 1), nrow = 2, byrow = TRUE)
I.dn_dt = II %*% dn_dt

# Compute I.dn_dt/||n||
I.dn_dt_n = I.dn_dt/norm_nvec

# Compute dnhat_dt
dnhat_dt_old = (I.dn_dt - nhat.nhat.dn_dt)/norm_nvec

tdvec <- tvec + dnhat_dt_old*d

norm_tdvec <- sqrt(tdvec[1]^2 + tdvec[2]^2)
#print(tvec)

# Compute dn_dt
dn_dt = c(cvec[2], -cvec[1])

# Compute dn_dt/||n||
dn_dt_n = dn_dt/norm_nvec

# Compare dn_dt/||n||
#print(dn_dt_n - I.dn_dt_n)

# Compute nhat . dn_dt
nhat.dn_dt = nhat[1]*dn_dt[1] + nhat[2]*dn_dt[2]

# Compute (nhat . dn_dt) nhat
nhat.dn_dt_nhat = nhat*nhat.dn_dt

# Compare (nhat . dn_dt) nhat
#print(nhat.dn_dt_nhat - nhat.nhat.dn_dt)

# Compute (nhat . dn_dt) nhat/||n||
nhat.dn_dt_nhat_n = nhat.dn_dt_nhat/norm_nvec

# Compute dnhat/dt
dnhat_dt = dn_dt_n - nhat.dn_dt_nhat_n

# Compare dnhat_dt
#print(dnhat_dt - dnhat_dt_old)

# Compute td
td = tvec + dnhat_dt*d

# Compute tdhat
tdhat = td/sqrt(td[1]^2 + td[2]^2)

# Compare td and tdvec
#print(td - tdvec)

return(tdhat)
}

x1 = c(1,2)
x2 = c(1.45,3)
x3 = c(1.55,3)
x4 = c(1.7,2)

tlist = seq(from = 0, to = 1, length.out = 50)
B <- sapply(tlist, function(t) {Bezier(t, x1, x2, x3, x4)})
Bprime <- sapply(tlist, function(t) {BezierTangent(t, x1, x2, x3, x4)})
Bnorm <- sapply(tlist, function(t) {BezierNormal(t, x1, x2, x3, x4)})
Bpprime <- sapply(tlist, function(t) {BezierConvex(t, x1, x2, x3, x4)})

d = 0.1
Boffset <- BezierOffset(B, Bnorm, d)
Boffsetprime <- sapply(tlist, function(t) {BezierOffsetTangent(t, x1, x2, x3, x4, d)})

T.Td <- mapply(function(Tx, Ty, Tdx, Tdy) {
Tt = c(Tx, Ty)
Td = c(Tdx, Tdy)
Tx*Tdx + Ty*Tdy
}, Bprime[1,],
Bprime[2,],
Boffsetprime[1,],
Boffsetprime[2,])
print(T.Td)

# Find sign chnage indices
sign.T.Td = sign(T.Td)
sign.changes <- function(d) {
p <- cumsum(rle(d)$lengths) + 1 p[-length(p)] } indices = sign.changes(sign.T.Td) df = data.frame(x = B[1,], y = B[2,], label = "Curve") df = rbind(df, data.frame(x = c(x1[1], x2[1], x3[1], x4[1]), y = c(x1[2], x2[2], x3[2], x4[2]), label = "Control points")) df = rbind(df, data.frame(x = Boffset[1,], y = Boffset[2,], label = "Offset")) df_pts = rbind(data.frame(x = Boffset[1,indices[1]], y = Boffset[2,indices[1]]), data.frame(x = Boffset[1,indices[2]], y = Boffset[2,indices[2]])) df_pts$label = "Cusps"

segment_1 = data.frame(x0 = Boffset[1,1:(indices[1]-1)], x1 = Boffset[1,2:indices[1]],
y0 = Boffset[2,1:(indices[1]-1)], y1 = Boffset[2,2:indices[1]])
segment_2 = data.frame(x0 = Boffset[1,indices[2]:(ncol(Boffset)-1)],
x1 = Boffset[1,(indices[2]+1):ncol(Boffset)],
y0 = Boffset[2,indices[2]:(ncol(Boffset)-1)],
y1 = Boffset[2,(indices[2]+1):ncol(Boffset)])
segment_window = owin(c(min(Boffset[1,]), max(Boffset[1,])),
c(min(Boffset[2,]), max(Boffset[2,])))
segment_1.psp = as.psp(segment_1, window = segment_window)
segment_2.psp = as.psp(segment_2, window = segment_window)
crossings = as.data.frame(crossing.psp(segment_1.psp, segment_2.psp))
print(crossings)
df_cross = data.frame(x = crossings$x, y = crossings$y, label = "Crossing")

plt = ggplot() +
geom_path(data = df[which(df$label == "Curve" | df$label == "Offset"),],
aes(x = x, y = y, color = label), size = 1) +
geom_point(data = df[which(df$label == "Control points"),], aes(x = x, y = y, color = label), size = 5) + geom_point(data = df_pts, aes(x = x, y = y, color = label), size = 2) + geom_point(data = df_cross, aes(x = x, y = y, color = label), size = 2) print(plt) dev.copy(png, "Crossings.png") dev.off() df_tan = data.frame(x = B[1,], xend = B[1,]+Bprime[1,]/20.0, y = B[2,], yend = B[2,]+Bprime[2,]/20.0, label = "Tangent") plt = plt + geom_segment(data = df_tan, aes(x = x, xend = xend, y = y, yend = yend, color = label), size = 0.75, arrow = arrow(length = unit(0.2, "cm"))) df_tan_orig = data.frame(x = Boffset[1,], xend = Boffset[1,]+Bprime[1,]/20.0, y = Boffset[2,], yend = Boffset[2,]+Bprime[2,]/20.0, label = "OrigTangent") plt = plt + geom_segment(data = df_tan_orig, aes(x = x, xend = xend, y = y, yend = yend, color = label), size = 0.5, arrow = arrow(length = unit(0.2, "cm"))) df_tan_off = data.frame(x = Boffset[1,], xend = Boffset[1,]+Boffsetprime[1,]/20.0, y = Boffset[2,], yend = Boffset[2,]+Boffsetprime[2,]/20.0, label = "OffsetTangent") plt = plt + geom_segment(data = df_tan_off, aes(x = x, xend = xend, y = y, yend = yend, color = label), size = 0.5, arrow = arrow(length = unit(0.2, "cm"))) #dev.new() #print(plt) df_nor = data.frame(x = B[1,], xend = B[1,]+Bnorm[1,]/20.0, y = B[2,], yend = B[2,]+Bnorm[2,]/20.0, label = "Normal") plt = plt + geom_segment(data = df_nor, aes(x = x, xend = xend, y = y, yend = yend, color = label), size = 0.5, arrow = arrow(length = unit(0.2, "cm"))) #dev.new() #print(plt) df_con = data.frame(x = B[1,], xend = B[1,]+Bpprime[1,]/20.0, y = B[2,], yend = B[2,]+Bpprime[2,]/20.0, label = "Convexity") plt = plt + geom_segment(data = df_con, aes(x = x, xend = xend, y = y, yend = yend, color = label), size = 0.5, arrow = arrow(length = unit(0.2, "cm"))) + coord_fixed() dev.new() print(plt) dev.copy(png, "TangentsNormals.png") dev.off() df_T = data.frame(x = Bprime[1,], y = Bprime[2,], label = "Tangent") df_T = rbind(df_T, data.frame(x = Boffsetprime[1,], y = Boffsetprime[2,], label = "OffsetTangent")) plt1 = ggplot(data = df_T) + geom_path(aes(x = x, y = y, color = label), size = 1.0) coord_fixed() #dev.new() #print(plt1) df_TTd = data.frame(t = tlist, T.Td = T.Td) plt2 = ggplot(data = df_TTd) + geom_path(aes(x = t, y = T.Td), size = 1.0) dev.new() print(plt2) dev.copy(png, "T_Td.png") dev.off()  • Thanks a lot for extremely detailed answer but I got very confused at times. Please: 1. I think I followed how you derived that td(t) = t(t) + ... . m(t) (above the second image). However, I cannot see how solving it numerically is any easier than to solve curvature(t) = 1/distance, as both have the square root in denominator. What am I missing? 2. I didn't really understand the "Curvature" part. If I plug the curvature formula into t(t).td(t)=0 I get 2|t|^2=0, hence no real solutions, ok. So I should instead use t.td=t.t+d/|t|(t.m), ok. – Faaf Oct 10 '16 at 14:04 • But then it results in using the curvature in sign[1+k(t)d], which is what we trying to avoid, no? Ok, so it's there just to show that the two approaches (curvature=1/distance and t.td) lead to same result, right? 3. I guess the main problem I have is I cannot see how this helps finding the parameter values. I assume I cannot see which of the formulas are easier to solve. Could you please sketch how to practically solve the (which?) formula? I mean, should I use Newton method to solve t(t).td(t) where td(t) = t(t) + ... . m(t)? How many zeros I should expect to find? – Faaf Oct 10 '16 at 14:05 • There's a jump discontinuity in t.td. It is this jump that we need to find the cusps. In my example there are two cusps that can be found by just examining the sign of$t.td$. You don't need to be very accurate for that. Then you break up the curve into three segments, throw away the middle segment. and numerically intersect the two remaining segments to find$t$. – Biswajit Banerjee Oct 10 '16 at 19:28 • I'd added my R code to the bottom. That should help make clear why a sign based numerical approach is easier. – Biswajit Banerjee Oct 10 '16 at 19:34 • Just to check I understood: I take, say, 50 samples along the curve and as soon I find a sign chain in t.td I should bisect to find precise location of the cusp, correct? – Faaf Oct 10 '16 at 22:32 To find the$t$value where the curvature equals$1/d$, you can simply solve the following equation numerically$f(t) = \frac{B'(t)\times B^"(t)}{|B'(t)|^3} - \frac{1}{d} = 0$by using any root-finding algorithms (see link), such as bisection method or Newton's method. Having said this, it is recommended to go with algorithms that do not require computation of the function's derivative (i.e.,$f'(t)$), which will require the computation of curvature's derivative. It is also recommended to deal with curvature, instead of radius of curvature, in your computation as the latter will approach infinity when curvature approaches zero. • For bisection I would need to know two$t$values with where the function has opposite signs and for Newton I don't really know how to get that derivation or a good starting point. I was kind of hoping there would be a way to express the problem as a nice polynomial or perhaps there is some other technique on how to locate the$t$values where curvature equals$1/d\$... – Faaf Oct 4 '16 at 13:42
• In general, anything involving curvature cannot be represented as a polynomial as the |B'(t)|^3 term involves a square root. You can sample along the curve uniformly to compute f(t) values and check the sign change to obtain the two t values where the f(t) has opposite signs. It might not be efficient, but at least it is a plausible solution. – fang Oct 4 '16 at 17:23