Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
// Copyright 2017 Todd Reed
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the “Software”), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
// BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
// DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// This code is part of the article at <http://toddreed.name/articles/line-fitting/>.
import Foundation
/// A line expressed in Hessian normal form: ρ = x cos θ + y sin θ.
///
/// `rho` is the length of the vector n from the origin to the line that is perpendicular to the
/// line. `theta` is the angle of between the x-axis and the vector n.
struct Line {
let rho: Double
let theta: Double
init(rho: Double, theta: Double) {
// Normalize rho and theta so that rho > 0 and theta in [0, 2π]
var r = rho
var t = theta
if (r < 0) {
r = -r
t += Double.pi
}
if (t < 0 || t > 2*Double.pi) {
t = NormalizeAngle(t)
}
if (r == 0 && t > Double.pi) {
t -= Double.pi
}
self.rho = r
self.theta = t
}
}
struct Point {
let x: Double
let y: Double
}
/// Normalizes an angle `theta` so that is lies in the 2π interval [center-π, center+π].
func NormalizeAngle(_ theta: Double, center: Double = Double.pi) -> Double
{
return theta-2*Double.pi*floor((theta+Double.pi-center)/(2*Double.pi))
}
/// Returns the best least-squares line fit to the provided data points using perpendicular
/// distances.
///
/// Returns nil if not enough points are provided or no unique solution exists.
func FitLine(points: [Point]) -> Line? {
if points.count < 2 {
return nil
}
var xSum = 0.0
var ySum = 0.0
for p in points {
xSum += p.x
ySum += p.y
}
let centroid = Point(x: xSum/Double(points.count), y: ySum/Double(points.count))
var xSquaredSum = 0.0
var ySquaredSum = 0.0
var xySum = 0.0
for p in points.map({ Point(x: $0.x-centroid.x, y: $0.y-centroid.y)}) {
xSquaredSum += p.x*p.x
ySquaredSum += p.y*p.y
xySum += p.x*p.y
}
let a = ySquaredSum-xSquaredSum
let b = xySum
if (a == 0.0 && b == 0.0) {
return nil
} else {
let beta = 2*b
let gamma = sqrt(a*a+beta*beta)
let theta = atan2(a-gamma, beta)
let rho = centroid.x*cos(theta) + centroid.y*sin(theta)
return Line(rho: rho, theta: theta)
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment