4. [Geometrie] Effektiv Senkrechte eines Punktes auf einer Linie

[Geometrie] Effektiv Senkrechte eines Punktes auf einer Linie

  • L

    Hallo Leute,

    ich bekomme 2 Punkte rein (x1, y1) und (x2, y2), welche eine (endliche) Linie ergeben sollen.

    Dann bekomme ich einen dritten Punkt (x3, y3), welcher sich irgendwo befindet. Ich würde gerne diesen dritten Punkt senkrecht auf die Linie projezieren und diesen projezierten neuen Punkt (x4, y4) zurückgeben. (BZW -1 wenn die senkrechte die Linie nicht trifft).

    Ich würde gerne ohne Trigonometrie auskommen, also nur Vektorrechnung etc. Gibt es da eine effiziente Möglichkeit?

    Ich danke schonmal!

    VG
    Lespaul

    0
  • Mod

    Skalarprodukt?

    0
  • W

    Hallo Lespaul,

    so könnte es gehen:

    #include <cassert>
#include <complex>
#include <iostream>
#include <tuple>

typedef std::complex< double > Vektor; // spart Tipparbeit
double X( const Vektor& v ) { return v.real(); }
double Y( const Vektor& v ) { return v.imag(); } // .. syntaktischer Zucker

// --   das Skalarprodukt
double operator,( const Vektor& a, const Vektor& b )
{
    return X(a)*X(b) + Y(a)*Y(b);
}

// ... effizient; zumindest was die Anzahl der Code-Zeilen angeht ;-)
std::tuple< Vektor, bool > fusspunkt( const Vektor& p1, const Vektor& p2, const Vektor& p3 )
{
    auto r = p2 - p1;
    assert( ("Die Punkte p1 und p2 duerfen nicht identisch sein", (r,r)) );
    auto t_f = ((p3 - p1),r)/(r,r);
    return std::make_tuple( p1 + t_f * r, 0 <= t_f && t_f <= 1 );
}
enum Idx { Value, IsValid };

int main() 
{
    using namespace std;
    for( Vektor p1, p2, p3; cin >> p1 >> p2 >> p3; )
    {
        auto res = fusspunkt( p1, p2, p3 );
        if( get< IsValid >( res ) )
            cout << "Der gesuchte Fusspunkt ist: " << get< Value >( res ) << endl;
        else
            cout << "Fusspunkt " << get< Value >( res ) << " liegt ausserhalb der Strecke" << endl;
    }
    return 0;
}

    Ein Dialog mit diesem Programm sieht z.B. so aus:

    (-1,3) (5,6) (5,1)
    Der gesuchte Fusspunkt ist: (3,5)
    (-1,3) (5,6) (10,1)
    Fusspunkt (7,7) liegt ausserhalb der Strecke

    Gruß
    Werner

    0
  • L

    Hallo Werner,

    Werner Salomon schrieb:

    so könnte es gehen:

    Danke erstmal. Ich habe das mal bei mir schnell hingedeichselt:

    double fussPunkt(double px, double py, double x1, double y1, double x2, double y2, double &x3, double &y3)
{
  double rx = x1 - x2;
  double ry = y1 - y2;;
  double dotn = crossProduct(rx, ry, rx, ry);

  double x31 = px - x1;
  double y31 = py - y1;

  //auto t_f = ((p3 - p1),r)/(r,r);
  double dot = crossProduct( x31, y31, rx, ry);
  double t_f = dot / dotn;

  //return std::make_tuple( p1 + t_f * r, 0 <= t_f && t_f <= 1 );
  x3 = x1 + t_f * rx;
  y3 = y1 + t_f * ry;

  if(t_f < 0)
  {
    return -1.0;
  }
  else if(1.0 > t_f)
  {
    return -2.0;
  }
  else
  {
//    return getDistance(px, py, x3, y3);
  }
  return -1.0;
}

    mit dem Kreuzprodukt:

    inline static double crossProduct(double x1, double y1, double x2, double y2) { return (x1*x2 + y1*y2); }

    Dabei habe ich zufällig folgende Koordinaten benutzt:
    p: (68.3595504,-133.728197) // zu projezierender Punkt
    l1: (68.3599538,-133.728215) // Linienpunkt 1
    l2: (68.3593208,-133.729684) // Linienpunkt 2

    fussPunkt: 68.35989705, -133.7283463 // Fusspunkt

    Google Maps sagt dann bei diesen Werten, wie man hier sehen kann:
    http://www.bilder-upload.eu/show.php?file=c89b54-1430295376.png

    dass der Fusspunkt nicht senrecht auf p zugeht. Wenigstens liegt der Punkt auf der Linie. Ich habe bestimmt ein Vorzeichen verdreht oder so? Finde jedoch nichts??

    Kann es vlt. an dem spärischen Koordinatensystem liegen? Das mit Kreuzprodukt etc. gilt ja für ebene Vektorrechnung.
    Aber die Abstände im Beispiel sind ja unbedeutend klein.

    😕

    Ich danke schonmal für alle Hinweise!

    0
  • L

    Erwartet hätte ich es irgendwo an der gelben Linie wie hier mit paint eingezeichnet:
    http://www.bilder-upload.eu/show.php?file=0696c6-1430295968.png

    0
  • L

    Wenn du die Punkte bzw Linien in einem orthonormalen Koordinatensystem aufträgst, ist deine Lösung auch senkrecht.

    Längen- und Breitengrad sind nunmal sowas wie Kugelkoordinaten (oder Ellipsoidenkoordinaten, falls es das gibt 😉 )...

    Um etwas umrechnen wirst du da nicht herumkommen. (Falls lon/lat die Anwendung ist... einen Punkt als Ursprung nehmen und dann in north/east den Fusspunkt berechnen...)

    0
  • W

    Hallo Lespaul,

    lagalopex hat völlig Recht. Bei 68Grad nördlicher Breite liegen die Längengrade weniger als halb so dicht beieinander wie die Breitengrade. Wenn Du das als kartesisches Koordinatensystem interpretierst, so ist dies auf keinen Fall mehr winkeltreu; daher dieses von Dir beschriebene Verhalten. Deine Rechnung ist völlig in Ordnung.

    Um das Problem zu lösen gibt es zwei Möglichkeiten. Die erste Möglichkeit besteht darin, diese 'Verzerrung' wieder rückgängig zu machen. Einfach durch einen Stretch-Faktor - das sieht dann so aus:

    Vektor strechtLon( const Vektor& v, double factor )
{
    return Vektor( X(v), Y(v)*factor );
}
std::tuple< Vektor, bool > fusspunkt_kugel( const Vektor& p1, const Vektor& p2, const Vektor& p3 )
{
    static double GRAD2RAD = std::acos(-1.0) / 180.;    // PI/180; erfordert #include <cmath>
    auto const factor = std::cos( 0.5 * (X(p1) + X(p2)) * GRAD2RAD ); // X(p_i) ist der Breitengrad in [Grad]
    auto r = strechtLon(p2,factor) - strechtLon(p1,factor);
    assert( ("Die Punkte p1 und p2 duerfen nicht identisch sein", (r,r)) );
    auto t_f = ((strechtLon(p3,factor) - strechtLon(p1,factor)),r) / (r,r);
   //  std::cout << "t_f = " << t_f << std::endl;
    return std::make_tuple( p1 + t_f * (p2 -p1), 0 <= t_f && t_f <= 1 );
}

    Aus der Funktion ' fusspunkt ' wird ' fusspunkt_kugel ' und ich habe auch die precision der Ausgabe auf 8 erhöht. Dann erhält man mit Deinen Ausgangsdaten als Ergebnis:

    (68.359724,-133.72875)

    .. und das ist schon genau da, wo Du es auch erwartet hast.

    Diese Lösung hat aber einige Nachteile. Sie liefert nur 'gute' Ergebnisse, wenn Die drei Punkten nicht allzu weit von einander entfernt sind. Dann solltest Du unbedingt den Polarregionen fern bleiben und auf keinen Fall darf eine der Verbindungslinien der vier Punkte die Datumsgrenze kreuzen.

    Um eine korrekte Lösung Lösung zu erhalten, so muss man einen anderen Algorithmus wählen. Dann lässt sich aber auch ein bisschen Trigonometrie nicht vermeiden. Schließlich besteht der Input nur aus Winkeln und die Koordinaten des Ergebnisses sollen ja auch wieder Winkel sein. Außerdem brauchen wir eine 'richtige' Vektor-Klasse, die ich aus eine alten Thread übernehme.
    BTW.:

    lespaul schrieb:

    mit dem Kreuzprodukt:

    inline static double crossProduct(double x1, double y1, double x2, double y2) { return (x1*x2 + y1*y2); }

    das ist nicht das Kreuzprodukt, sondern das Skalarprodukt. Das Kreuzprodukt verwende ich für die 2.Lösung (s. operator%)

    #include "vector.h" // s.: https://www.c-plusplus.net/forum/p1108802#1108802
#include <cassert>    
#include <cmath>    
#include <complex>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <tuple>

namespace
{
    double const PI = std::acos(-1.0);
    double const GRAD2RAD = PI/180.;
    double const RAD2GRAD = 180./PI;
}

typedef std::complex< double > Koordinate; // spart Tipparbeit
double Lat( const Koordinate& k ) { return k.real(); }  // geo. Breite in [Grad]
double Lon( const Koordinate& k ) { return k.imag(); }  // geo. Länge in [Grad]

typedef basic_Vector< double > Vector;

Vector make_vector( const Koordinate& k )
{
    using std::sin;
    using std::cos;
    auto r = cos( Lat(k) * GRAD2RAD );
    return Vector( r * cos( Lon(k) * GRAD2RAD ), r * sin( Lon(k) * GRAD2RAD ), sin( Lat(k) * GRAD2RAD ) );
}
Koordinate make_geo_coordinate( const Vector& v )
{
    using std::asin;
    using std::atan2;
    using std::abs;
    assert( ("Der Vektor zur Umrechnung in Koordinaten  muss die Laenge 1 haben!", abs(abs(v) - 1.0) <= 1.E-8 ) );
    return Koordinate( asin( Z(v) ) * RAD2GRAD, atan2( Y(v), X(v) ) * RAD2GRAD );
}

std::tuple< Koordinate, bool > fusspunkt_kugel( const Koordinate& k1, const Koordinate& k2, const Koordinate& k3 )
{
    static auto const EPS = std::sqrt( std::numeric_limits< double >::min() );
    auto p1 = make_vector( k1 );
    auto p2 = make_vector( k2 );
    auto p3 = make_vector( k3 );
    auto p12 = p1 % p2;
    auto f = p12 % p3 % p12;
    auto f_len = abs(f);
    if( f_len < EPS )
        throw std::invalid_argument( "Fusspunkt nicht eindeutig bei diesen Eingabeparametern" );
    f /= f_len; // |f| == 1
    return std::make_tuple( make_geo_coordinate( f ), p1*p2 <= f*p2 && p1*p2 <= p1*f );
}
enum Idx { Value, IsValid };

int main()
{
    using namespace std;
    cout.precision( 8 );
    for( Koordinate p1, p2, p3; cin >> p1 >> p2 >> p3; )
    {
        auto res = fusspunkt_kugel( p1, p2, p3 );
        if( get< IsValid >( res ) )
            cout << "Der gesuchte Fusspunkt ist: " << get< Value >( res ) << endl;
        else
            cout << "Fusspunkt " << get< Value >( res ) << " liegt ausserhalb der Strecke" << endl;
    }
    return 0;
}

    .. falls es nicht kompiliert; es fehlen noch die Funktionen X, Y und Z mit Zugriff auf die einzelnen Koordinaten.

    Dieser Algorithmus liefert auch dann noch korrekte Werte, wenn die Linie über einen der Pole verläuft. Das Ergebnis für Deine Eingabewerte ist hier übrigens identisch zu oben

    Gruß
    Werner

    0
  • B

    @lespaul
    In welchem Koordinatensystem sind deine XY Koordinaten? Irgentwie fehlt mir da die Höhe. Vieleicht kann man da eine Schätzrechnung machen.

    Eigentlich müsstest du die ellipsoiden Koordinaten mittels Gauss Krüger Transformation umrechnen und danach den Fusspunkt bestimmen. Den Fusspunkt müsstest du dann mittels der inversen Gauss Krüger Transformation wieder zurückrechnen. Aber das ist keine einfache Rechnung...

    0
  • L

    Ich werde die Effizienz abwägen. Da die Koordinaten sowieso in Mercator verfügbar waren, ich diese für mich in WGS84 umgerechnet hatte, werde ich mir diesen Zwischenschritt sparen und die Mercator Koordinaten ohne Konvertierung auf den ersten "nicht sphärischen" Algorithmus von Werner anwenden.

    Vorteile:
    + keine Zwischenkonvertierung
    + keine Trigonometrie

    Nachteil:
    - Mercator ist ein bissl ungewohnt für mich

    Vielen Dank nochmal allerseits

    0
  • L

    Bitte ein Bit schrieb:

    @lespaul
    In welchem Koordinatensystem sind deine XY Koordinaten?

    Mercator und WGS84.

    Bitte ein Bit schrieb:

    @lespaul
    Irgentwie fehlt mir da die Höhe.

    Ja die Höhe habe ich noch nicht. Weil mir eig. auch (noch) die Höhen-Daten fehlen.

    Ist es richtig, dass die "Kugel-Variante" später in der 3.Dimension einfacher nachgerüstet werden kann, als "nicht sphärische" Systeme?

    0
  • W

    lespaul schrieb:

    Ist es richtig, dass die "Kugel-Variante" später in der 3.Dimension einfacher nachgerüstet werden kann, als "nicht sphärische" Systeme?

    IMHO ist es eher umgekehrt! In der 'kartesischen' Variante, musst Du nur aus dem 2D-Vektor einen 3D-Vektor machen. Den Code selbst musst Du nicht ändern, um eine eine evt. Höhenlage der gegebenen Punkte zu berücksichtigen. Nur aufpassen bei der 'Einheit' der Höhe - die muss zu den anderen Angaben passen!

    Auf der Kugeloberfläche mit Höhen zu rechnen ist wesentlich schwieriger. Da musst Du erst mal definieren (rein mathematisch) was eine 'Linie' zwischen zwei Punkten mit unterschiedlichne Höhen überhaupt sein soll. Eine Gerade ist das ja nicht mehr. Und was genau bedeutet dann 'steht senkrecht' drauf? .. für die Verbindung von P3 zum gesuchten Fußpunkt.

    Also ich könnte Dir das nicht so einfach hinschreiben. Und ich denke, dass dies sehr davon abhängt, was für ein Problem Du mit diesem Vorgehen lösen möchtest. Das hast Du uns ja noch nicht verraten 😉

    Gruß
    Werner

    0
Anmelden zum Antworten