se/SE/se-lib/Geometry.php

<?php

class GeometryPoint {
	public $x, $y;

	public function __construct($x, $y) {
		if (!(is_numeric($x) && is_numeric($y))) throw new Exception('Bad arguments');
		$this->x = (float)$x;
		$this->y = (float)$y;
	}

	public static function same(GeometryPoint $a, GeometryPoint $b) {
		return ($a->x == $b->x && $a->y == $b->y);
	}
}

class GeometryLine {
	public $a, $b;

	public function __construct(GeometryPoint $a, GeometryPoint $b) {
		$this->a = $a;
		$this->b = $b;
	}

	public function length() {
		return Geometry::distance($this->a, $this->b);
	}

	public function asText($reproject = null) {
		if ($reproject) list($a, $b) = GeometryObject::reproject([$this->a, $this->b], $reproject);
		else list($a, $b) = [$this->a, $this->b];
		return "LINESTRING({$a->x} {$a->y},{$b->x} {$b->y})";
	}
}

class GeometryPolygon {
	protected $points;

	public function __construct(GeometryPoint ...$points) {
		if (reset($points)->x != end($points)->x || reset($points)->y != end($points)->y) $points[] = reset($points);
		$this->points = $points;
	}

	public function asText($reproject = null) {
		if ($reproject) $points = GeometryObject::reproject($this->points, $reproject);
		else $points = $this->points;
		return 'POLYGON((' . implode(',', array_map(function ($point) {
			return "{$point->x} {$point->y}";
		}, $points)) . '))';
	}
}

class GeometryObject {
	protected $type, $points;
	public function type() { return $this->type; }
	public function numPoints() { return count($this->points); }
	public function points() { return $this->points; }

	public function __construct($type, $points) {
		$type = (string)$type;
		if (!($type && is_array($points))) throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument');
		foreach ($points as $point) if (get_class($point) != 'GeometryPoint') throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument');
		$this->type = $type;
		$this->points = $points;
	}

	public static function getFromText($text, $reproject = null) {
		if (!preg_match('/^([[:alpha:]]+)\((.*)\)$/', $text, $matches)) throw new Exception('Invalid WKT syntax');

		$type = (string)$matches[1];
		$points = array_map(function ($point) {
			list($x, $y) = explode(' ', $point, 2);
			return new GeometryPoint($x, $y);
		}, explode(',', $matches[2]));
		if ($reproject) $points = self::reproject($points, $reproject);

		return new GeometryObject($type, $points);
	}

	public static function reproject($points, $reproject) {
		if (!is_array($points)) throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument #1');
		Lib::loadClass('EpsgConversion');
		if (!method_exists('EpsgConversion', $reproject)) throw new Exception("Function EpsgConversion::{$reproject} not exists");

		return array_map(function ($point) use ($reproject) {
			if (get_class($point) != 'GeometryPoint') throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument #2');
			$_point = EpsgConversion::$reproject($point->x, $point->y);
			return new GeometryPoint($_point->x, $_point->y);
		}, $points);
	}
}

class GeometryLinearFunction {
	protected $a, $b, $x;

	public function __construct($a = null, $b = null, $x = null) {
		if (($a === null && $b !== null) || ($a !== null && $b === null) || ($a === null && $x === null) || ($a !== null && $x !== null)) throw new Exception('Bad arguments');
		$this->a = $a; $this->b = $b; $this->x = $x;
	}

	public function valueOfX($x) {
		if ($this->x !== null) return null;
		return $this->a * $x + $this->b;
	}

	public function valueOfY($y) {
		if ($this->x !== null) return $x;
		if ($this->a == 0) return null;
		return ($y - $this->b) / $this->a;
	}

	public static function getFromLine(GeometryLine $L) {
		$a = null; $b = null; $x = null;
		if ($L->a->x == $L->b->x) {
			$x = $L->a->x;
		} else {
			$a = ($L->a->y - $L->b->y) / ($L->a->x - $L->b->x);
			$b = $L->b->y - $a * $L->b->x;
		}
		return new GeometryLinearFunction($a, $b, $x);
	}

	public static function findCrossPoint(GeometryLinearFunction $f1, GeometryLinearFunction $f2) {
		if ($f1->x !== null) {
			if ($f2->x !== null) {
				if ($f1->x == $f2->x) return $f1;
				return null;
			}
			$x = $f1->x;
			$y = $f2->a * $x + $f2->b;
		} elseif ($f2->x !== null) {
			$x = $f2->x;
			$y = $f1->a * $x + $f1->b;
		} else {
			if ($f1->a == $f2->a) {
				if ($f1->b == $f2->b) return $f1;
				return null;
			}
			$x = ($f2->b - $f1->b) / ($f1->a - $f2->a);
			$y = $f1->a * $x + $f1->b;
		}
		return new GeometryPoint($x, $y);
	}

	public static function functionPerpendicularAtPoint(GeometryLinearFunction $f, GeometryPoint $p) {
		$a = null; $b = null; $x = null;
		if ($f->x !== null) {
			$a = 0;
			$b = $p->y;
		} elseif ($f->a == 0) {
			$x = $p->x;
		} else {
			$a = -1 * (1 / $f->a);
			$b = $p->y - $a * $p->x;
		}
		return new GeometryLinearFunction($a, $b, $x);
	}

	public static function pointsDistanced(GeometryLinearFunction $f, GeometryPoint $p, $distance) {
		if (!is_numeric($distance)) throw new Exception('Bad argument');
		$fpap = self::functionPerpendicularAtPoint($f, $p);
		if ($fpap->x !== null) {
			return [new GeometryPoint($p->x, $p->y - 1), new GeometryPoint($p->x, $p->y + 1)];
		} else {
			$d = $distance / sqrt(1 + pow($fpap->a, 2));
			$x1 = $p->x - $d;
			$x2 = $p->x + $d;
			$y1 = $fpap->valueOfX($x1);
			$y2 = $fpap->valueOfX($x2);
			return [new GeometryPoint($x1, $y1), new GeometryPoint($x2, $y2)];
		}
	}
}

class Geometry {
	public static function point($x, $y) { return new GeometryPoint($x, $y); }
	public static function line(GeometryPoint $a, GeometryPoint $b) { return new GeometryLine($a, $b); }
	public static function objectFromText($text, $reproject = null) { return GeometryObject::getFromText($text, $reproject); }
	public static function samePoint(GeometryPoint $a, GeometryPoint $b) { return GeometryPoint::same($a, $b); }

	public static function pointsToLines($points) {
		if (!is_array($points)) throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument');
		$lines = [];
		$lastPoint = null;
		foreach ($points as $point) {
			if (get_class($point) != 'GeometryPoint') throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument');
			if ($lastPoint) $lines[] = self::line($lastPoint, $point);
			$lastPoint = $point;
		}
		return $lines;
	}

	public static function closedPointOnLine(GeometryPoint $p, GeometryLine $L) {
		// Funkcja zwraca najbliży punkt należący do odcinka L względem punktu p
		$f1 = GeometryLinearFunction::getFromLine($L);
		$f2 = GeometryLinearFunction::functionPerpendicularAtPoint($f1, $p);
		$crossPoint = GeometryLinearFunction::findCrossPoint($f1, $f2);
		if (self::distance($crossPoint, $L->a) <= $L->length() && self::distance($crossPoint, $L->b) <= $L->length()) return $crossPoint;
		if (self::distance($crossPoint, $L->a) < self::distance($crossPoint, $L->b)) return $L->a;
		return $L->b;
	}

	public static function distance($a, $b) {
		if (!(is_object($a) && is_object($b))) throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - bad argument');

		if (get_class($a) == 'GeometryPoint' && get_class($b) == 'GeometryPoint') return round(sqrt(pow($a->x - $b->x, 2) + pow($a->y - $b->y, 2)), 7);
		if (get_class($a) == 'GeometryPoint' && get_class($b) == 'GeometryLine') return self::distance($a, self::closedPointOnLine($a, $b));
		if (get_class($a) == 'GeometryLine' && get_class($b) == 'GeometryPoint') return self::distance($b, $a);
		if (get_class($a) == 'GeometryLine' && get_class($b) == 'GeometryLine') {
			if (self::crossPoint($a, $b)) return 0;
			return min(self::distance($a, $b->a), self::distance($a, $b->b), self::distance($a->a, $b), self::distance($a->b, $b));
		}

		throw new Exception(__CLASS__ . "::" . __FUNCTION__ . ' - object not implemented');
	}

	public static function crossPoint(GeometryLine $A, GeometryLine $B) {
		if (!($p = GeometryLinearFunction::findCrossPoint(GeometryLinearFunction::getFromLine($A), GeometryLinearFunction::getFromLine($B)))) return null;
		switch (get_class($p)) {
			case 'GeometryPoint':
				if (self::distance($p, $A) == 0 && self::distance($p, $B) == 0) return [$p];
//				if ($p->x >= min($A->a->x, $A->b->x) && $p->x <= max($A->a->x, $A->b->x) && $p->y >= min($A->a->y, $A->b->y) && $p->y <= max($A->a->y, $A->b->y) &&
//					$p->x >= min($B->a->x, $B->b->x) && $p->x <= max($B->a->x, $B->b->x) && $p->y >= min($B->a->y, $B->b->y) && $p->y <= max($B->a->y, $B->b->y)) return [$p];
				break;
			case 'GeometryLinearFunction':
				if ($A->a->x > $A->b->x) {
					$tmp = $A->a; $A->a = $A->b; $A->b = $tmp;
					$reverse = true;
				} else $reverse = false;
				if ($B->a->x > $B->b->x) {
					$tmp = $B->a; $B->a = $B->b; $B->b = $tmp;
				}
				if ($A->a->x > $B->b->x || $A->b->x < $B->a->x) return null;
				if ($A->a->x > $B->a->x) $a = $A->a; else $a = $B->a;
				if ($A->b->x < $B->b->x) $b = $A->b; else $b = $B->b;
				if ($a->x == $b->x) return [$a]; elseif ($reverse) return [$b, $a]; else return [$a, $b];
			default: throw new Exception('Uknown error');
		}
		return null;
	}

	public static function lineToRectangle(GeometryLine $L, $size) {
		$f = GeometryLinearFunction::getFromLine($L);
		$p1 = GeometryLinearFunction::pointsDistanced($f, $L->a, $size);
		$p2 = GeometryLinearFunction::pointsDistanced($f, $L->b, $size);
		return new GeometryPolygon($p1[0], $p2[0], $p2[1], $p1[1]);
	}
}