Ausgangssituation

Jede Etappe wird durch zwei Geokoordinaten from = {lat: x.xxx, lng: x.xxxx} (Starpunkt) und to = {lat: x.xxx, lng: x.xxx} (Zielpunkt) beschrieben. Die Idee war es nun, die Verbindungslinie zwischen diesen beiden Punkten parabelformig zu krümmen. Die kürzeste Verbindungskurve zwischen zwei Punkten nennt man Geodäte oder im Spezialfall der sphärischen Geometrie auch Orthodrome. Da die Google Maps API Linien nur als Poligonzug (Array aus Geokoordinaten) darstellen kann, bestand die erste herausforderung darin, die Orthodrome in $N$ gleichlange Segmente zu unterteilen.

Segmentieren von Orthodromen

Berechnen der Länge auf der Einheitskugel und auf der Erdkugel:

    public static function distanceOnUnitSphere(
        LatLng $from, 
        LatLng $to
    )
    {
        $from = $from->toRad();
        $to   = $to->toRad();

        return acos(
            sin($from->lat) * sin($to->lat) +
            cos($from->lat) * cos($to->lat) * cos($to->lng - $from->lng)
        );
    }

    public static function distance(
        LatLng $from, 
        LatLng $to, 
        $radius = self::EARTH_RADIUS
     )
    {
        return self::distanceOnUnitSphere($from, $to) * $radius;
    }

Bruchteil zwischen zwei Geokoordinaten finden:

    public  static function fraction(
        LatLng $from, 
        LatLng $to, 
        $fraction
    ) 
    {
        if ($fraction == 0)
            return $from;

        if ($fraction == 1)
            return $to;

        // http://en.wikipedia.org/wiki/Slerp
        $from = $from->toRad();
        $to   = $to->toRad();
        $cosFromLat = cos($from->lat);
        $cosToLat   = cos($to->lat);

        // Computes spherical interpolation coefficients.
        $angle = self::distanceOnUnitSphere($from, $to);
        $sinAngle = sin($angle);

        $a = sin((1 - $fraction) * $angle) / $sinAngle;
        $b = sin($fraction * $angle) / $sinAngle;

        // Converts from polar to vector and interpolate.
        $x = $a * $cosFromLat * cos($from->lng) + $b * $cosToLat * cos($to->lng);
        $y = $a * $cosFromLat * sin($from->lng) + $b * $cosToLat * sin($to->lng);
        $z = $a * sin($from->lat) + $b * sin($to->lat);

        // Converts interpolated vector back to polar.
        $lat = atan2($z, sqrt($x ** 2 + $y ** 2));
        $lng = atan2($y, $x);

        return new LatLng($lat * Geo::RAD_TO_DEG, $lng * Geo::RAD_TO_DEG);
    }

Poligon aus Geokoordinaten unterteilen:

    public static function subdevidePoliLine(/*$steps, $point, $point, ...*/)
    {
        $points = func_get_args();
        $steps  = array_shift($points);

        $fraction = 1 / ($steps + 1);

        $result = [];

        for ($i = 0; $i < sizeof($points) - 1; $i++) {
            $from = $points[$i];
            $to   = $points[$i + 1];

            $result[] = $from;

            for ($j = 1; $j <= $steps; $j++) {
                $result[] = self::fraction($from, $to, $j * $fraction);
            }
        }

        $result[] = $to;
        return $result;
    }