geometrycomputational-geometryintersectionbeziercubic-bezier

Perpendicular line from bezier curve to point

Question

I need to get the the point Q of a cubic (2d) bezier curve B(t) where the line from point Q to another given point P intersects perpendicular with the bezier curve.

What I did until now (you can skip this)

Note that I think this ansatz is wrong. This is only included for completeness.

I tried to solve this with my (basic) knowledge of mathematics but I cannot finish it. This is what I have now (please don't too strict with the notation, I'm not very good in this):

The following formulas will be expressed as y(x) to get one result this has to be calculated for y(x) and x(y). The point P is the control point, Q is the point where the line g(x) from Q to P is perpendicular on the bezier curve B(t) = (x, y)T. The expression for the line g(x) can be retrieved by

where B(x) is the bezier curve in Cartesian coordinates, B'(x) the derivative (in Cartesian coordinates) and k is the intersection with the y axis. To get the slope of g(x) one has to solve

To calculate B(x) one has to solve B(t) for t and then plugging it back in B(t). So in each point on the bezier curve there is the relation

which also applies for the derivative B'(t).

The derivative of B(t) is (according to wikipedia)

Solving this for t (with wolfram alpha) results in

where a0 = (P1 - P0)x, a1 = (P2 - P1)x and a2 = (P3 - P2)x. Plugging the *ti*s back in B(t) results in (wolfram alpha for t1, wolfram alpha for t2, wolfram alpha for t3)

Now the next thing would be to use y = B'(x) and the second equation and eliminate x but I have no idea how and I don't even know if this is possible.

Solution

  • I found a implementation for a approximation for my problem made by mbostock on Github who implemented the idea from this online book which was acutally created by @Mike'Pomax'Kamermans about bezier curves. If you have to deal with bezier curves check this out. This explains most of the problems I had with bezier curves.

    The idea of the algorithm is the following one:

    1. Make a rough approximation:
      1. Calculate Qapprox, i by plugging in different tis in B(t) (mbostock uses t1 = 0, t2 = 1/8, t3 = 2/8, ...).
      2. Calculate the quadratic (save one square root calculation) distance between Qapprox, i and P.
      3. Save the Qapprox, i and the ti which is the closest (with the smallest distance).
    2. Make a more precise approximation:
      1. Choose a precision q.
      2. Calculate tbefore = ti - q.
      3. Check if the distance between between Qapprox, before = B(tbefore) and P is smaller than the distance between Qapprox, i and P, if yes set Qapprox, i = Qapprox, before and start at 2 (precise approximation), if no continue.
      4. Calculate tafter = ti + q.
      5. Check if the distance between between Qapprox, after = B(tafter) and P is smaller than the distance between Qapprox, i and P, if yes set Qapprox, i = Qafter and start at 2 (precise approximation), if no continue.
      6. Qapprox, i is the closest. If the precision is good enough stop here. If not decrease the precision q (mbostock divides by two) and start over again at 2 (precise approximation).

    As already mentioned there is the implementation of mbostock on Github. I paste the code here in case the link goes down. THIS IS NOT MY OWN CODE.

    var points = [[474,276],[586,393],[378,388],[338,323],[341,138],[547,252],[589,148],[346,227],[365,108],[562,62]];
var width = 960,
    height = 500;
var line = d3.svg.line()
    .interpolate("cardinal");
var svg = d3.select("body").append("svg")
    .attr("width", width)
    .attr("height", height);
var path = svg.append("path")
    .datum(points)
    .attr("d", line);
var line = svg.append("line");
var circle = svg.append("circle")
    .attr("cx", -10)
    .attr("cy", -10)
    .attr("r", 3.5);
svg.append("rect")
    .attr("width", width)
    .attr("height", height)
    .on("mousemove", mousemoved);
// adding coordinates display
var coords = svg.append("text");
function mousemoved() {
  var m = d3.mouse(this),
      p = closestPoint(path.node(), m);
  line.attr("x1", p[0]).attr("y1", p[1]).attr("x2", m[0]).attr("y2", m[1]);
  circle.attr("cx", p[0]).attr("cy", p[1]);
  coords.attr("x", (p[0] + m[0]) / 2).attr("y", (p[1] + m[1]) / 2).html("Q(" + Math.round(p[0]) + ", " + Math.round(p[1]) + ")");
}
function closestPoint(pathNode, point) {
  var pathLength = pathNode.getTotalLength(),
      precision = 8,
      best,
      bestLength,
      bestDistance = Infinity;
  // linear scan for coarse approximation
  for (var scan, scanLength = 0, scanDistance; scanLength <= pathLength; scanLength += precision) {
    if ((scanDistance = distance2(scan = pathNode.getPointAtLength(scanLength))) < bestDistance) {
      best = scan, bestLength = scanLength, bestDistance = scanDistance;
    }
  }
  // binary search for precise estimate
  precision /= 2;
  while (precision > 0.5) {
    var before,
        after,
        beforeLength,
        afterLength,
        beforeDistance,
        afterDistance;
    if ((beforeLength = bestLength - precision) >= 0 && (beforeDistance = distance2(before = pathNode.getPointAtLength(beforeLength))) < bestDistance) {
      best = before, bestLength = beforeLength, bestDistance = beforeDistance;
    } else if ((afterLength = bestLength + precision) <= pathLength && (afterDistance = distance2(after = pathNode.getPointAtLength(afterLength))) < bestDistance) {
      best = after, bestLength = afterLength, bestDistance = afterDistance;
    } else {
      precision /= 2;
    }
  }
  best = [best.x, best.y];
  best.distance = Math.sqrt(bestDistance);
  return best;
  function distance2(p) {
    var dx = p.x - point[0],
        dy = p.y - point[1];
    return dx * dx + dy * dy;
  }
}
    .disclaimer{
  padding: 10px;
  border-left: 3px solid #ffcc00;
  background: #fffddd;
}
path {
  fill: none;
  stroke: #000;
  stroke-width: 1.5px;
}
line {
  fill: none;
  stroke: red;
  stroke-width: 1.5px;
}
circle {
  fill: red;
}
rect {
  fill: none;
  cursor: crosshair;
  pointer-events: all;
}
    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.17/d3.min.js"></script>
<div class="disclaimer">
  This is not my own code. This is made by <a href="https://gist.github.com/mbostock/8027637">mbostock on Github</a> based on Mike Kamermans <a href="https://pomax.github.io/bezierinfo/#projections">Primer on Bézier Curves online book</a>.
</div>