Trigonometry: Compute a point location on an arc


#1

Hi all,

I would like to find a way to compute the location of p4 on the picture below.

All I have is:

  • the location of p1, p2 and p3
  • the radius (length of segment p2 - p3)
  • the angle between segment p1-p2 and segment p2-p3

Ideally, I would like to compute any location of p4 along the arc drawn on the picture.

My trigonometry is not on point and I would really appreciate some help !

void setup(){
    size(400, 400, P2D);
    strokeCap(SQUARE);
    textSize(12);
    fill(0);
    noFill();
    smooth(8);
}
    
void draw(){
    background(255);
    
    PVector p0 = new PVector(width/2 - 10, height/2);
    PVector p1 = new PVector(width/2 + 20, height/3);
    PVector p2 = new PVector(mouseX, mouseY);
    
    // Display points -> p1, p2 and p3
    strokeWeight(6);
    stroke(0);
    point(p0.x, p0.y);
    point(p1.x, p1.y);
    point(p2.x, p2.y);
    text("p0", p0.x + 10, p0.y);
    text("p1", p1.x - 23, p1.y);
    text("p2", p2.x + 10, p2.y);
    
    // Display lines between points
    strokeWeight(1);
    stroke(30, 90, 200);
    line(p0.x, p0.y, p1.x, p1.y);
    line(p1.x, p1.y, p2.x, p2.y);
    
    // Draw arc
    stroke(0, 205, 120);
    strokeWeight(10);
    arc(p1.x, p1.y, 40, 40, atan2(p1.y - p0.y, p1.x - p0.x), atan2(p2.y - p1.y, p2.x - p1.x));
    
    // Compute angle between (p0 - p1) and (p1 - p2)
    float angle = atan2(p2.y - p1.y, p2.x - p1.x) - atan2(p1.y - p0.y, p1.x - p0.x);
    
    // Display text
    text(degrees(angle), 15, 25);
        
    // PROBLEM: Compute p4 location
    float d = dist(p2.x, p2.y, p1.x, p1.y);
    float f = 0.5; // between 0 and 1
    float x = d * cos(angle*f);
    float y = d * sin(angle*f);
    
    stroke(255, 60, 60);
    strokeWeight(5);
    point(p1.x + x , p1.y + y);
    text("p4", p1.x + x + 10 , p1.y + y);

}

#2

try adding just half angle??
( i changed also the confusing numbering )

void setup() {
  size(400, 400, P2D);
  strokeCap(SQUARE);
  textSize(12);
  fill(0);
  noFill();
  smooth(8);
}

void draw() {
  background(255);

  PVector p0 = new PVector(width/2 - 10, height/2);
  PVector p1 = new PVector(width/2 + 20, height/3);
  PVector p2 = new PVector(mouseX, mouseY);

  // Display points -> p0, p1 and p2
  strokeWeight(6);
  stroke(0);
  point(p0.x, p0.y);
  point(p1.x, p1.y);
  point(p2.x, p2.y);
  text("p0", p0.x + 10, p0.y);
  text("p1", p1.x - 23, p1.y);
  text("p2", p2.x + 10, p2.y);

  // Display lines between points
  strokeWeight(1);
  stroke(30, 90, 200);
  line(p0.x, p0.y, p1.x, p1.y);
  line(p1.x, p1.y, p2.x, p2.y);

  // Draw arc
  stroke(0, 205, 120);
  strokeWeight(10);
  float ang_01 = atan2(p1.y - p0.y, p1.x - p0.x);
  float ang_12 = atan2(p2.y - p1.y, p2.x - p1.x);
  arc(p1.x, p1.y, 40, 40, ang_01, ang_12);
  float angle = ang_01 - ang_12;     // Compute angle between (p0 - p1) and (p1 - p2)
  text(degrees(angle), 15, 25);     // Display text

  float d = dist(p2.x, p2.y, p1.x, p1.y);
  float f = 0.5; // between 0 and 1
  float ang_13 = ang_01 - angle*f;
  float x = d * cos(ang_13);
  float y = d * sin(ang_13);
  PVector p3 = new PVector(p1.x+x, p1.y + y);

  strokeWeight(1);
  stroke(200, 0, 0);
  line(p1.x, p1.y, p3.x, p3.y);
  stroke(255, 60, 60);
  strokeWeight(5);
  point(p3.x, p3.y);
  text("p3", p3.x + 10, p3.y);
}


#3

Ah, it seems obvious now ! Thank you so much for the prompt reply @kll


#4

Here is a version with a small function that is all PVector built-ins – probably less efficient, but fairly clear to read.

/**
 * LerpObtuse - interpolate the rotation point at the end of an obtuse angle
 * out towards 180 degrees.
 * 2019-03-12 - Processing 3.4
 * discourse.processing.org/t/trigonometry-compute-a-point-location-on-an-arc/9191/2
 **/ 
 
PVector[] pts;
void setup() {
  size(300, 300);
  rectMode(CENTER);
  fill(0);
  pts = new PVector[]{
    new PVector(80, 280),   // p0 sets baseline
    new PVector(100, 200),  // p1 is midpoint
    new PVector(250, 100)   // p2 sets arm of arc
  };
}

void draw() {
  background(255);
  float amt = mouseX/(float)width;
  drawObtuse(pts);
  PVector solve = lerpObtuse(amt, pts);
  rect(solve.x, solve.y, 20, 20);
}

/**
 * Takes an obtuse angle of points p0-p1-p2 and
 * provides a lerp position between p2 rotated out towards the p0 baseline. 
 **/
PVector lerpObtuse(float amt, PVector[] pts) {
  PVector side1 = PVector.sub(pts[2], pts[1]);
  PVector side2 = PVector.sub(pts[1], pts[0]);
  PVector midpoint = side1.rotate(lerp(-PVector.angleBetween(side2, side1), 0, amt));
  return PVector.add(pts[1], midpoint); // convert relative to absolute position
}

/**
 * Visual illustration of the obtuse angle model.
 */
void drawObtuse(PVector[] pts) {
  for (PVector pt : pts) {
    rect(pt.x, pt.y, 5, 5);
  }
  line(pts[1].x, pts[1].y, pts[2].x, pts[2].y);
  line(pts[1].x, pts[1].y, pts[0].x, pts[0].y);
  pushMatrix();
  pushStyle();
  fill(255,0,0,64);
  stroke(255,0,0,64);
  translate(pts[1].x, pts[1].y);
  PVector side1 = PVector.sub(pts[2], pts[1]);
  PVector side2 = PVector.sub(pts[1], pts[0]);
  side2.setMag(height);
  line(0, 0, side1.x, side1.y);
  line(0, 0, side2.x, side2.y);
  arc(0, 0, 200, 200, side2.heading(), side1.heading());
  popStyle();
  popMatrix();
}

LerpObtuse--screenshot


#5

There’s still a problem with my approach.

When the angle between p2 and p1 ( ang_12 in @kll’s code) goes below -3.14 it instantly turns to 3.14 (from negative to positive PI). This interferes with the computation of the final angle ( ang_13 ) and causes p4 to be placed very far from where it should normally go.


#6

Can you use modulo (%) on the input to restrict its range?


#7

or change from atan2 to PVector.heading


#8

@jeremydouglass No, I’m afraid it’s not possible. The range of angle a2 is already restricted between 3.14 and -3.14. I would need something to prevent this angle to jump from below 3.14 to positive 3.14. Something like a sin() function but for PI (ex: …-3.12, -3.13, -3.14, -3.13, - 3.12, …up to 3.12, 3.13, 3.14, 3.13, 3.12, … down to -3.12…etc)

@kll Interesting. How would you compute ang_01 and ang_12 in you code with .heading() ? Also, wouldn’t that be much slower that to use atan2. Ideally I would like to compute a few million angles so speed is an issue.


#9

Perhaps use the range 0-TWO_PI instead?

Just to make sure that I understand what you actually want – do you want something like this, which “flips” when it crosses p1-p2?

PVector[] pts;
float amt = 0.5;

void setup() {
  size(300, 300);
  rectMode(CENTER);
  fill(0);
  pts = new PVector[]{
    new PVector(80, 280), // p0 sets baseline
    new PVector(100, 200), // p1 is midpoint
    new PVector(250, 100)   // p2 sets arm of arc
  };
}

void draw() {
  background(255);

  pts[2] = new PVector(mouseX, mouseY);
  drawObtuse(pts);
  PVector solve = lerpObtuse(amt, pts);
  rect(solve.x, solve.y, 20, 20);
}

/**
 * Takes an obtuse angle of points p0-p1-p2 and
 * provides a lerp position between p2 rotated out towards the p0 baseline. 
 **/
PVector lerpObtuse(float amt, PVector[] pts) {
  PVector side1 = PVector.sub(pts[2], pts[1]);
  PVector side2 = PVector.sub(pts[1], pts[0]);
  side2.rotate(-side1.heading());
  PVector midpoint = side1.rotate(lerp(side2.heading(), 0, amt));
  return PVector.add(pts[1], midpoint); // convert relative to absolute position
}

/**
 * Visual illustration of the obtuse angle model.
 */
void drawObtuse(PVector[] pts) {
  for (PVector pt : pts) {
    rect(pt.x, pt.y, 5, 5);
  }
  line(pts[1].x, pts[1].y, pts[2].x, pts[2].y);
  line(pts[1].x, pts[1].y, pts[0].x, pts[0].y);
  pushMatrix();
  pushStyle();
  fill(255, 0, 0, 64);
  stroke(255, 0, 0, 64);
  translate(pts[1].x, pts[1].y);
  PVector side1 = PVector.sub(pts[2], pts[1]);
  PVector side2 = PVector.sub(pts[1], pts[0]);
  side2.setMag(height);
  line(0, 0, side1.x, side1.y);
  line(0, 0, side2.x, side2.y);
  popStyle();
  popMatrix();
}

Or do you want to be able to track TWO_PI worth of winding in either direction, clockwise or counterclockwise? That would additionally require saving state so that you know which direction the mouse came from.


#10

I tried your previous example, it works great (thank you !) and yes, that’s what I would like to do.
It also answers my second question about heading() (thanks again).

The only concern I have with this solution is speed. I would like to compute millions of angles and converting the points locations to PVector() + calling sub(), rotate(), heading() for each one of them is very expensive.

side1 = PVector.sub(PVector(p1[0], p1[1]), PVector(p0[0], p0[1]))
side2 = PVector.sub(PVector(p2[0], p2[1]), PVector(p1[0], p1[1]))
side1.rotate(-side2.heading())
angle = side2.heading() - (side2.heading() - side1.heading())

That’s why I was wondering if a workaround involving atan2 or hard coded trigonometry would be possible. Something that would keep the calculation of the angle as simple as:

angle = atan2(p2[1]-p1[1], p2[0]-p1[0]) - atan2(p1[1]-p0[1], p1[0]-p0[0])


#11

Of course! .sub() is just subtraction, and .heading() is a simple one-liner, (float) Math.atan2(y, x):

…and .rotate() is:

float temp = x;
x = x*cos(theta) - y*sin(theta);
y = temp*sin(theta) + y*cos(theta);

So you could convert each step, starting with something like:

float side1x = pts[2].x - pts[1].x;
float side1y = pts[2].y - pts[1].y;
float side2x = pts[1].x - pts[0].x;
float side2y = pts[1].y - pts[0].y;
float side2xRot = side2x*cos(theta) - side2y*sin(theta)
float side2yRot = side2x*sin(theta) + side2y*cos(theta);

and see if that improves performance.

Are they all the same, with different orientations, or are they random… depending on what those angles area and how they are changing you could probably get away with reusing a lot of your intermediate calculation values – do them once outside the loop, then plug them in.


#12

Apologies, that was a silly question. Looking for the source code was the obvious way to go. Thanks for the kind reminder and the detailed answer, that helped a lot.

Recap:

##### with built-in functions #####

side1 = PVector.sub(PVector(p1.x, p1.y), PVector(p0.x, p0.y))
side2 = PVector.sub(PVector(p2.x, p2.y), PVector(p1.x, p1.y))
side1.rotate(-side2.heading())
angle = side2.heading() - (side2.heading() - side1.heading())

##############################


##### Hard-coded #####

#side1 & side2
s1x = p1.x - p0.x
s1y = p1.y - p0.y
s2x = p2.x - p1.x
s2y = p2.y - p1.y
                            
#side2.heading()
s2h = atan2(s2y, s2x)
                          
#side1.rotate(-side2.heading())
s1rotx = s1x * cos(-s2h) - s1y * sin(-s2h)
s1roty = s1x * sin(-s2h) + s1y * cos(-s2h)
                            
#side1.heading()
s1h = atan2(s1roty, s1rotx)
                            
#angle   
angle = s2h - (s2h - s1h)

###############################

I’m now wondering if calling the math module in Python mode would be faster than the Processing built-in trigonometric functions (math.cos() vs cos()).

Angles are random. I wish I could save a lot intermediate calculation but I’m afraid it’s not feasible in this case (computing angles between millions of points uniformly distributed).


#13

Absolutely – I think it is always worth trying and seeing what the timeit results are. It looks like Python mode maps cos against PApplet.cos()

and that uses the Java Math library. Keep in mind that Python mode runs on Jython runs on JVM, so the Java cos might actually be more “native” than Python math – but trying it out is the best thing.

If the millions of points are uniformly distributed you might (depending on the precision you need) try building a lookup table of points-angles and see if that is faster, e.g. a HashMap.


#14

I did some test and it turns out calling the functions from the math module in Python mode is indeed slower than using the bult-in functions from Processing.

Speed test
import timeit, math

p0 = PVector(18.1204, 140.3879)
p1 = PVector(94.1049, 31.9122)
p2 = PVector(108.1049, 231.9122)

def f():
    s1x = p1.x - p0.x
    s1y = p1.y - p0.y
    s2x = p2.x - p1.x
    s2y = p2.y - p1.y
                                
    s2h = atan2(s2y, s2x)
                            
    s1rotx = s1x * cos(-s2h) - s1y * sin(-s2h)
    s1roty = s1x * sin(-s2h) + s1y * cos(-s2h)
                                
    s1h = atan2(s1roty, s1rotx)
                                
    angle = s2h - (s2h - s1h)
    
    return angle
    

print timeit.timeit('f()', 'from __main__ import f', number = 1000000)

It takes around 5.6 seconds to run the test function a million times using Python math module against 3.6 seconds without. That’s about a 36% difference.

The example picture below took around 22 minutes to render with the built-in functions against 35 using Python math module.


#15

Thanks for sharing your test results – and the image!

If you need better performance and don’t want to try a lookup table, given your application, perhaps just solving the lerp point on a triangle, rather than an arc? That would simplify the math – no trig functions. You could also try it with no intermediate PVector objects – or no PVector input/output at all, if you like.

PVector lerpObtuse2(float amt, PVector[] pts) {
  float side1x = pts[2].x - pts[1].x;
  float side1y = pts[2].y - pts[1].y;
  float side2x = pts[1].x - pts[0].x;
  float side2y = pts[1].y - pts[0].y;
  float mx = lerp(side2x, side1x, amt);
  float my = lerp(side2y, side1y, amt);
  return new PVector(mx, my);
}

or, more compactly:

PVector lerpObtuse2(float amt, PVector[] pts) {
  return new PVector(
    lerp(pts[1].x - pts[0].x, pts[2].x - pts[1].x, amt), 
    lerp(pts[1].y - pts[0].y, pts[2].y - pts[1].y, amt)
    );
}