Scala: Orthogonal projection of a point onto a line

Question

I'm trying to make a function that takes 3 points as arguments. The first two of which represent two points on a line. The third one represents another point, outside of that line. Suppose a perpendicular through the third point on the line defined by the first two points. Now what I want to do, is calculate that intersection. I've come up with this procedure so far, but somehow it works only like 50% of the time. Could somebody figure out what I'm doing wrong here?

def calculateIntersection(p1: (Double, Double), p2: (Double, Double), c: (Double, Double)): (Double, Double) = {
    var intersection: (Double, Double) = null
    // CASE 1: line is vertical
    if(p1._1 == p2._1) {
      intersection = (p1._1, c._2)
    }
    // CASE 2: line is horizontal
    else if(p1._2 == p2._2) {
      intersection = (c._1, p1._2)
    }
    // CASE 3: line is neither vertical, nor horizontal
    else {
      val slope1: Double = (p2._2 - p1._2) / (p2._1 - p1._1) // slope of the line
      val slope2: Double = pow(slope1, -1) * -1 // slope of the perpendicular
      val intercept1: Double = p1._2 - (slope1 * p1._1) // y-intercept of the line
      val intercept2: Double = c._2 - (slope2 * c._1) // y-intercept of the perpendicular
      intersection = ((intercept2 - intercept1) / (slope1 - slope2), 
                     slope1 * ((intercept2 - intercept1) / (slope1 - slope2)) + intercept1)
    }
    intersection
}

There are at least three grave mistakes in this code: 1) You are treating horizontal lines separately. 2) You are treating vertical lines separately. 3) You are computing something with "slopes". No wonder it doesn't work. But that's hardly your fault: I've tried googling or finding a suitable duplicate, and I couldn't find anything. The amount of outrageous nonsense written on this topic is just staggering... Can anyone find a sane duplicate, please? — Andrey Tyukin, Jan 02 '19 at 17:05
What do you mean exactly? I thought working with the slope was the exact way to do it. My strategy was to calculate de equation of the line (given the two points), then extracting the slope out of that equation, because we know the slope of the perpendicular will be the opposite reciprocal of the slope of the line. Given the third point, and that slope, we can also find the equation of the perpendicular. So now I converted the two equations in order to find the intersection. — Frederik Vanclooster, Jan 02 '19 at 17:29
Possible duplicate of [Perpendicular on a line from a given point](https://stackoverflow.com/questions/1811549/perpendicular-on-a-line-from-a-given-point) — Ripi2, Jan 02 '19 at 17:51
I mean that if anyone tries to sell you "angles", "slopes" or "square roots" when dealing with objects that are completely covered by basic linear algebra, you should run as fast as possible. Remember, [Carmack wrote Quake 3 without even using square roots](https://en.wikipedia.org/wiki/Fast_inverse_square_root). If your sketch on the paper contains only points, lines, and planes, but the formula contains angles, slopes, and square roots, you should treat it with healthy suspicion. — Andrey Tyukin, Jan 02 '19 at 18:04

Andrey Tyukin · Answer 1 · 2019-01-02T19:03:07.067

Given the following definitions:

type Point = (Double, Double)

implicit class PointOps(p: Point) {
  def +(other: Point) = (p._1 + other._1, p._2 + other._2)
  def -(other: Point) = (p._1 - other._1, p._2 - other._2)
  def dot(other: Point) = p._1 * other._1 + p._2 * other._2
  def *(scalar: Double) = (p._1 * scalar, p._2 * scalar)
  def normSquare: Double = p._1 * p._1 + p._2 * p._2
}

meaning that

a + b        // is vector addition
a - b        // is vector subtraction
a dot b      // is the dot product (scalar product)
a * f        // is multiplication of a vector `a` with a scalar factor `f`
a.normSquare // is the squared length of a vector

you obtain the projection of a point p on line going through points line1 and line2 as follows:

/** Projects point `p` on line going through two points `line1` and `line2`. */
def projectPointOnLine(line1: Point, line2: Point, p: Point): Point = {
  val v = p - line1
  val d = line2 - line1
  line1 + d * ((v dot d) / d.normSquare)
}

Example:

println(projectPointOnLine((-1.0, 10.0), (7.0, 4.0), (6.0, 11.0)))

gives

(3.0, 7.0)

This works in 3D (or n-D) in exactly the same way.

Some math behind that (How to derive it from scratch)

(notation as above)

We have three points: l1 and l2 for the line and p for the target point. We want to project the point p orthogonally onto the line that goes through l1 and l2 (assuming l1 != l2).

Let d = l2 - l1 be the direction from l1 to l2. Then every point on the line can be represented as

l1 + d * t

with some scalar factor t. Now we want to find a t such that the vector connecting p and l1 + d * t is orthogonal to d, that is:

(p - (l1 + d * t)) dot d == 0

Recall that

(v1 + v2) dot v3 = (v1 dot v3) + (v2 dot v3)

for all vectors v1, v2, v3, and that

(v1 * s) dot v2 = (v1 dot v2) * s

for scalar factors s. Using this and the definition v.normSquared = v dot v, we obtain:

(p - l1 - d * t) dot d 
= (p - l1) dot d - (d dot d) * t
= (p - l1) dot d - d.normSquare * t

and this should become 0. Resolving for t gives:

t = ((p - l1) dot d) / d.normSquare

and this is exactly the formula that is used in the code.

(Thanks at SergGr for adding an initial sketch of the derivation)

Andrey I took the liberty of adding some description of the logic in words. Please take a look if you are OK with that. — SergGr, Jan 02 '19 at 18:18
@SergGr Thanks for trying. Adding mathematical / geometric explanations seems a bit difficult here, because there is no `TeX` support, and because I don't know where to get a halfway sane sketch to illustrate what's going on. I'll leave it as-is for now, and maybe add a sketch later. — Andrey Tyukin, Jan 02 '19 at 18:36

Scala: Orthogonal projection of a point onto a line

1 Answers1