Cross points of a line and a circle

Calendar Clock iconCalendar Clock icon

geometry

Table of contents

# Explanation

Let the center point of a circle CC, radius RR、2 points which a line ll go through p1p_1, p2p_2.
Let cross points of the line and the circle aa and bb.

# 1. Calculate cross point pp of a perpendicular from CC to ll and ll

Use projection formula

Vector2 project(Vector2 p1, Vector2 p2, Vector2 c) {
    Vector2 v = p2 - p1;
    return p1 + v * (v.dot(c - p1) / v.norm());
}

# 2. Calculate unit vector \vec

Devide a vector with its length.

Vector2 unitVec(Vector2 v) {
  return v / v.length();
}

# 3. Calculate the length of segment abab

pp devides abab into 2 equal parts.
∠Cpb\angle{Cpb} is a right angle.
Length of apap and pbpb is, by Pythegorean theorem,

ap=pb=R2−∣Cp⃗∣2ap = pb = \sqrt{R^2 - |\vec{Cp}|^2}

# 4. Calculate cross points

a,b=p⃗+e⃗⋅±apa, b = \vec{p} + \vec{e} \cdot \pm ap

# Entire Code

// C++ 14
#include <iostream>
#include <string>
#include <vector>
#include <list>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <unordered_map>
#include <math.h>

#define ll long long
#define Int int
#define loop(x, start, end) for(Int x = start; x < end; x++)
#define loopdown(x, start, end) for(int x = start; x > end; x--)
#define rep(n) for(int x = 0; x < n; x++)
#define span(a,x,y) a.begin()+x,a.begin()+y
#define span_all(a) a.begin(),a.end()
#define len(x) (x.size())
#define last(x) (*(x.end()-1))

using namespace std;

#define EPS 0.0000000001
#define fequals(a,b) (fabs((a) - (b)) < EPS)

class Vector2 {
public:
  double x, y;
  
  Vector2(double x = 0, double y = 0): x(x), y(y) {}
  
  Vector2 operator + (const Vector2 v) const { return Vector2(x + v.x, y + v.y); }
  Vector2 operator - (const Vector2 v) const { return Vector2(x - v.x, y - v.y); }
  Vector2 operator * (const double k) const { return Vector2(x * k, y * k); }
  Vector2 operator / (const double k) const { return Vector2(x / k, y / k); }
  
  double length() { return sqrt(norm()); }
  double norm() { return x * x + y * y; }
  double dot (Vector2 const v) { return x * v.x + y * v.y; }
  double cross (Vector2 const v) { return x * v.y - y * v.x; }
  
  bool parallel(Vector2 &other) {
    return fequals(0, cross(other));
  }
  
  bool orthogonal(Vector2 &other) {
    return fequals(0, dot(other));
  }
  
  bool operator < (const Vector2 &v) {
    return x != v.x ? x < v.x : y < v.y;
  }
  
  bool operator == (const Vector2 &v) {
    return fabs(x - v.x) < EPS && fabs(y - v.y) < EPS;
  }
};

ostream & operator << (ostream & out, Vector2 const & v) { 
  out<< "Vector2(" << v.x << ", " << v.y << ')';
  return out;
}

istream & operator >> (istream & in, Vector2 & v) { 
  double x, y;
  in >> x;
  in >> y;
  v.x = x;
  v.y = y;
  return in;
}
double R;
Vector2 C, p1, p2;

Vector2 project(Vector2 p1, Vector2 p2, Vector2 c) {
    Vector2 v = p2 - p1;
    return p1 + v * (v.dot(c - p1) / v.norm());
}

Vector2 unitVec(Vector2 v) {
  return v / v.length();
}

pair<Vector2, Vector2> crossPoints() {
  Vector2 p = project(p1, p2, C);
  Vector2 e = unitVec(p2 - p1);
  double len = sqrt(R*R - (p - C).norm());
  Vector2 a = p + e * (-len);
  Vector2 b = p + e * len;
  if (fabs(a.x) < EPS) a.x = 0.0;
  if (fabs(a.y) < EPS) a.y = 0.0;
  if (fabs(b.x) < EPS) b.x = 0.0;
  if (fabs(b.y) < EPS) b.y = 0.0;
  if (a < b) return make_pair(a, b);
  return make_pair(b, a);
}

Shun

Github Circle icon

Remote freelancer. A web and mobile application enginner.
Traveling around the world based on East Asia.
I'm looking forward to your job offers from all over the world!

Offer jobs or contact me!

Comments