Distance to a bezier curve

Gallery
1

Hi everyone,

Following the discussion on this thread relative to pattern attractors, I wanted to implement a method to find the distance to a bezier curve.

You can download the project on this gitHub repository or below if you don’t mind the mess :sweat_smile:

There are plenty of stuff to improve but it is working and it may, one day, help someone :slight_smile:.

It someone is interested, I’ll explain in detail how this whole thing works.

Point[] P; 
Point M;
float deltaX, deltaY;
BezierDistanceFinder bdf;
ArrayList<Double> closest;

void setup() {
  size(500, 500);
  
  P = new Point[4];
  P[0] = new Point(100, 400, color(200, 200, 10), 10);
  P[1] = new Point(100, 200, color(200, 200, 10), 10);
  P[2] = new Point(400, 200, color(200, 200, 10), 10);
  P[3] = new Point(400, 400, color(200, 200, 10), 10);
  M = new Point(370, 270, color(255, 20, 20), 10);
  
  bdf = new BezierDistanceFinder(P[0].pos(), P[1].pos(), P[2].pos(), P[3].pos());
  closest = bdf.getTClosestTo(M.pos());
}


void draw() {
  background(20);

  // Draw the control lines
  stroke(50);
  strokeWeight(2);
  line(P[0].pos().x, P[0].pos().y, P[1].pos().x, P[1].pos().y);
  line(P[2].pos().x, P[2].pos().y, P[3].pos().x, P[3].pos().y);

  // Draw the bezier curves
  noFill();
  stroke(200);
  strokeWeight(2);
  bezier(P[0].pos().x, P[0].pos().y, P[1].pos().x, P[1].pos().y, P[2].pos().x, P[2].pos().y, P[3].pos().x, P[3].pos().y);


  // Draw the points of interest
  for (int i = 0; i < 4; i++) {
    P[i].show();
  }
  M.show();
  
  if (closest.size() == 0) {
    return;
  }
  for (int i = 0; i < closest.size(); i++) {
    fill(120, 170, 250);
    noStroke();
    double t = closest.get(i);
    
    ellipse(bezierPoint(P[0].pos().x, P[1].pos().x, P[2].pos().x, P[3].pos().x, (float)t), bezierPoint(P[0].pos().y, P[1].pos().y, P[2].pos().y, P[3].pos().y, (float)t), 10, 10);
  }
}


void mousePressed() {
  // Check if mouse is over a point
  for (int i = 0; i < 4; i++) {
    if (P[i].isHit(mouseX, mouseY)) {
      P[i].lock();
      deltaX = mouseX - P[i].pos().x;
      deltaY = mouseY - P[i].pos().y;
      return;
    }
  }

  if (M.isHit(mouseX, mouseY)) {
    M.lock();
    deltaX = mouseX - M.pos().x;
    deltaY = mouseY - M.pos().y;
  }
}


void mouseReleased() {
  for (int i = 0; i < 4; i++) {
    P[i].unlock();
  }
  M.unlock();
}


void mouseDragged() {
  for (int i = 0; i < 4; i++) {
    if (P[i].isLocked()) {
      P[i].move(mouseX - deltaX, mouseY - deltaY);
    }
  }
  
  if (M.isLocked()) {
    M.move(mouseX - deltaX, mouseY - deltaY);
  }
  
  bdf.changeParameters(P[0].pos(), P[1].pos(), P[2].pos(), P[3].pos());
  closest = bdf.getTClosestTo(M.pos());
}



class BezierDistanceFinder {

  private PVector P0, P1, P2, P3;


  BezierDistanceFinder(PVector p_P0, PVector p_P1, PVector p_P2, PVector p_P3) {
    P0 = p_P0;
    P1 = p_P1;
    P2 = p_P2;
    P3 = p_P3;
  }


  void changeParameters(PVector p_P0, PVector p_P1, PVector p_P2, PVector p_P3) {
    P0 = p_P0;
    P1 = p_P1;
    P2 = p_P2;
    P3 = p_P3;
  }


  ComplexNumber evaluateDotProduct(PVector M, ComplexNumber t) {
    ComplexNumber result, perX, perY, tanX, tanY;
    double divisor;

    perX = t.pow(3).mult(-P3.x).add(t.mult(-1.0).add(1.0).mult(t.pow(2)).mult(-3*P2.x)).add(t.mult(-1.0).add(1.0).pow(2).mult(t).mult(-3*P1.x)).add(t.mult(-1.0).add(1.0).pow(3).mult(-P0.x)).add(M.x);
    perY = t.pow(3).mult(-P3.y).add(t.mult(-1.0).add(1.0).mult(t.pow(2)).mult(-3*P2.y)).add(t.mult(-1.0).add(1.0).pow(2).mult(t).mult(-3*P1.y)).add(t.mult(-1.0).add(1.0).pow(3).mult(-P0.y)).add(M.y);
    tanX = t.pow(2).mult(3*P3.x).add(t.mult(-1.0).add(1.0).mult(t).mult(6*P2.x)).add(t.pow(2).mult(-3*P2.x)).add(t.mult(-1.0).add(1.0).mult(t).mult(-6*P1.x)).add(t.mult(-1.0).add(1.0).pow(2).mult(3*P1.x)).add(t.mult(-1.0).add(1.0).pow(2).mult(-3*P0.x));
    tanY = t.pow(2).mult(3*P3.y).add(t.mult(-1.0).add(1.0).mult(t).mult(6*P2.y)).add(t.pow(2).mult(-3*P2.y)).add(t.mult(-1.0).add(1.0).mult(t).mult(-6*P1.y)).add(t.mult(-1.0).add(1.0).pow(2).mult(3*P1.y)).add(t.mult(-1.0).add(1.0).pow(2).mult(-3*P0.y));
    divisor = -3*P0.x*P0.x+18*P0.x*P1.x-18*P0.x*P2.x+6*P0.x*P3.x-27*P1.x*P1.x+54*P1.x*P2.x-18*P1.x*P3.x-27*P2.x*P2.x+18*P2.x*P3.x-3*P3.x*P3.x-3*P0.y*P0.y+18*P0.y*P1.y-18*P0.y*P2.y+6*P0.y*P3.y-27*P1.y*P1.y+54*P1.y*P2.y-18*P1.y*P3.y-27*P2.y*P2.y+18*P2.y*P3.y-3*P3.y*P3.y;
    //divisor = 1;

    result = perX.mult(tanX).add(perY.mult(tanY)).divide(divisor);

    return result;
  }

  private ComplexNumber[] getNextGeneration(ComplexNumber[] old, PVector M) {
    ComplexNumber[] result = new ComplexNumber[5];

    result[0] = evaluateDotProduct(M, old[0]).divide(old[0].add(old[1].mult(-1)).mult(old[0].add(old[2].mult(-1))).mult(old[0].add(old[3].mult(-1))).mult(old[0].add(old[4].mult(-1)))).mult(-1).add(old[0]);
    result[1] = evaluateDotProduct(M, old[1]).divide(old[1].add(result[0].mult(-1)).mult(old[1].add(old[2].mult(-1))).mult(old[1].add(old[3].mult(-1))).mult(old[1].add(old[4].mult(-1)))).mult(-1).add(old[1]);
    result[2] = evaluateDotProduct(M, old[2]).divide(old[2].add(result[1].mult(-1)).mult(old[2].add(result[0].mult(-1))).mult(old[2].add(old[3].mult(-1))).mult(old[2].add(old[4].mult(-1)))).mult(-1).add(old[2]);
    result[3] = evaluateDotProduct(M, old[3]).divide(old[3].add(result[1].mult(-1)).mult(old[3].add(result[2].mult(-1))).mult(old[3].add(result[0].mult(-1))).mult(old[3].add(old[4].mult(-1)))).mult(-1).add(old[3]);
    result[4] = evaluateDotProduct(M, old[4]).divide(old[4].add(result[1].mult(-1)).mult(old[4].add(result[2].mult(-1))).mult(old[4].add(result[3].mult(-1))).mult(old[4].add(result[0].mult(-1)))).mult(-1).add(old[4]);

    return result;
  }


  private ComplexNumber[] getRootsOfDotProduct(PVector M) {
    ComplexNumber[] result = new ComplexNumber[5];

    result[0] = new ComplexNumber(1, 0);
    result[1] = new ComplexNumber(0.4, 0.9);
    result[2] = result[1].pow(2);
    result[3] = result[1].pow(3);
    result[4] = result[1].pow(4);

    for (int i = 0; i < 15; i++) {
      result = getNextGeneration(result, M);
    }

    return result;
  }


  ArrayList<Double> getTClosestTo(PVector M) {
    ComplexNumber[] roots = getRootsOfDotProduct(M);
    ArrayList<Double> tempResult = new ArrayList<Double>();
    ArrayList<Double> result = new ArrayList<Double>();
    double minDist = Double.POSITIVE_INFINITY;
    float threshold = 0.00001;

    for (int i = 0; i < 5; i++) {
      if (roots[i].Re() > 0 && roots[i].Re() < 1) {  
        tempResult.add(roots[i].Re());
      }
    }
    tempResult.add((double)0);
    tempResult.add((double)1);
    
    for (int i = 0; i < tempResult.size(); i++) {
      double t = tempResult.get(i);
      double dist = dist(M.x, M.y, bezierPoint(P[0].pos().x, P[1].pos().x, P[2].pos().x, P[3].pos().x, (float)t), bezierPoint(P[0].pos().y, P[1].pos().y, P[2].pos().y, P[3].pos().y, (float)t));
      if(dist < minDist) {
        minDist = dist;
      }
    }
    
    for (int i = 0; i < tempResult.size(); i++) {
      double t = tempResult.get(i);
      double dist = dist(M.x, M.y, bezierPoint(P[0].pos().x, P[1].pos().x, P[2].pos().x, P[3].pos().x, (float)t), bezierPoint(P[0].pos().y, P[1].pos().y, P[2].pos().y, P[3].pos().y, (float)t)) - minDist;
      if(abs((float)dist) < threshold) {
        result.add(t);
      }
    }
    
    return result;
  }
}



class Point {
  
  private PVector pos;
  private color col;
  private int size;
  private boolean locked;
 
  Point(float x, float y, color p_col, int p_size) {
    pos = new PVector(x, y);
    col = p_col;
    size = p_size;
    locked = false;
  }
  
  void move(float x, float y) {
   pos.x = x;
   pos.y = y;
  }
  
  void show() {
    fill(col);
    noStroke();
    ellipse(pos.x, pos.y, size, size);
  }
  
  PVector pos() {
    return pos;
  }
  
  void lock() {
    locked = true;
  }
  
  void unlock() {
    locked = false;
  }
  
  boolean isLocked() {
    return locked;
  }
  
  boolean isHit(float x, float y) {
    return (pos.x - x) * (pos.x - x) + (pos.y - y) * (pos.y - y) < (size * size) / 4.0;
  }
}



class ComplexNumber {

  private double rPart, iPart;
  
  
  ComplexNumber(double r, double i) {
    rPart = r;
    iPart = i;
  }
  
  
  ComplexNumber(double r) {
    rPart = r;
    iPart = 0;
  }
  
  
  double Re() {
    return rPart;
  }
  
  
  double Im() {
    return iPart;
  }
  
  
  ComplexNumber mult(ComplexNumber nb) {
    return new ComplexNumber(rPart * nb.Re() - iPart * nb.Im(), rPart * nb.Im() + iPart * nb.Re());
  }
  
  
  ComplexNumber mult(double lambda) {
    return new ComplexNumber(lambda * rPart, lambda * iPart);
  }
  
  
  ComplexNumber add(ComplexNumber nb) {
    return new ComplexNumber(rPart + nb.Re(), iPart + nb.Im());
  }
  
  ComplexNumber add(double nb) {
    return new ComplexNumber(rPart + nb, iPart);
  }
  
  
  ComplexNumber divide(ComplexNumber nb) {
    double divisor = nb.Re() * nb.Re() + nb.Im() * nb.Im();
    double dividendRe = rPart * nb.Re() + iPart * nb.Im();
    double dividendIm = iPart * nb.Re() - rPart * nb.Im();

    return new ComplexNumber(dividendRe / divisor, dividendIm / divisor);
  }
  
  
  ComplexNumber divide(double divisor) {
    return new ComplexNumber(rPart / divisor, iPart / divisor);
  }
  
  
  ComplexNumber pow(int n) {
    ComplexNumber result;
    
    if (n == 0) {
      return new ComplexNumber(1, 0);
    }
    
    if (n == 1) {
      return new ComplexNumber(rPart, iPart);
    }
    
    result = new ComplexNumber(rPart, iPart);   
    for (int i = 0; i < n - 1; i++) {
      result = result.mult(this);
    }
    return result;
  }
  
  
  void log() {
    println(rPart + " + " + iPart + "i");
  }
}
1 Like
2

Can you explain what you use evaluateDotProduct() for?

3

Hi,
I recently coded a progam to find the shortest distance between a point and an arbitrary line. I used a graphical calculation for that.
What I did was drawing circles around my point, getting bigger and bigger. At some point the circle will intersect with the line. Then the radius of the circle is the shortest distance wich will be at a right angle to the line.

4

Hi,
How were you testing for crossing between your line and your circle?
I first thought of doing something like this too but the goal here was to get really precise result quite fast.

5

Hey @tony,

Let’s consider the following figure:

BezierCurve

The goal is to find a time t0 ϵ [0, 1] for which the distance to the point M is the minimum.
At this time, the vector M->P(t0) is at a right angle with the bezier curve.

The problem can now be rewritten as finding t0 for wich M->P(t0) perpendicular to (dP/dt)(t0)

And that is the same as finding the roots of M->P(t).(dP/dt)(t) where . is the dot product. In this case, the function is a quintic function and the evaluateDotProduct() is evaluating this function for any given t.

Then, the Durand–Kerner method is used to solve the function.

Hope it is clear enough.

1 Like
6

I color-coded the line and the cirlce. Then I simply checked if the color in a pixel represents both and there you are.

7

So you had to check every pixels on the the screen?

8

I drew that on a hidden PGraphic but basically yes. It was not critical how long the calculation took…

9

The advantage here is that you can select the precision independently of the way it is draw.

If I get your method correctly the more zoomed out you are the less precise you are. So you need to be really zoomed in to get high precision.

After, it is just a matter of what you need in the end.