Trigonometry: Compute a point location on an arc

solub · March 11, 2019, 12:37pm

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);

}

kll · March 11, 2019, 1:12pm

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);
}

solub · March 11, 2019, 1:33pm

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

jeremydouglass · March 12, 2019, 10:44pm

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

solub · March 13, 2019, 4:28pm

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.

jeremydouglass · March 13, 2019, 7:06pm

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

kll · March 14, 2019, 1:46am

or change from atan2 to PVector.heading

solub · March 14, 2019, 3:33pm

@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.

jeremydouglass · March 14, 2019, 8:30pm

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.

solub · March 14, 2019, 11:18pm

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])

jeremydouglass · March 16, 2019, 7:12am

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

github.com

processing/processing/blob/master/core/src/processing/core/PVector.java#L895-L898


      
            return angle;
          }

…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.

solub · March 16, 2019, 1:38pm

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).

jeremydouglass · March 16, 2019, 3:29pm

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()

github.com

jdf/processing.py/blob/master/runtime/src/jycessing/core.py#L325-L334


      
          for f in """
          acos append arrayCopy asin atan atan2 binary blendColor ceil
          concat cos createInput createOutput createReader day debug degrees dist
          exp expand floor hex hour join lerp  loadBytes log mag match matchAll
          minute month nf nfc nfp nfs norm pow radians reverse saveBytes saveStream
          saveStrings second shorten sin sort splice split splitTokens sq sqrt
          subset tan trim unbinary unhex year
          """.split():
              assert getattr(PApplet, f)
              setattr(__builtin__, f, getattr(PApplet, f))

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.

solub · March 16, 2019, 5:03pm

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.

jeremydouglass · March 16, 2019, 9:12pm

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)
    );
}

Topic		Replies	Views
Trigonometry help--find third point of the triangle Coding Questions	16	1191	March 15, 2019
Trig question, drawing an arrow head on an arc Coding Questions	14	1057	August 20, 2019
Programming advice (spiral from point to point) Coding Questions	16	499	November 27, 2022
Trig question re line/circle intersections Coding Questions	5	533	April 28, 2021
Arcs and angles with PI (wings) Coding Questions	17	281	September 23, 2023

Trigonometry: Compute a point location on an arc

Related Topics