Zufällige 3D Vektoren innerhalb einer Kugel

314159265358979

Naja, es sollten zumindest nirgens Lücken zu sehen sein.

cooky451

Wie wär's mit nem Würfel + Sichtschutz?

Bashar

Würde das nicht alles ein bisschen komisch aussehen? Je nach Maßstab sollten Sterne entweder praktisch unendlich weit weg sein oder du mitten durch düsen. Eine relativ begrenzte Schale von Sternen, die das das Spielgebiet umschließt, würde dagegen doch mit Sicherheit auffallen.

314159265358979

Im Moment sind die Sterne mittendrin im Bereich, in dem sich das Raumschiff bewegen soll. Sieht eigentlich gut aus, bis auf die Tatsache, dass man zu den Kugeln hinfliegen kann und die dann halt sieht.

Bedenken habe ich da keine.

pumuckl

Wenn dein Rauschiff sich nur in einem begrenzen Bereich befindet und die nächsten Sterne nicht erreichen kann, dann ist die Näherung, dass alle Sterne unendlich weit entfernt sind, einiges realistischer. Die Lösung dafür ist einfach: Male das Bild der umgebenden Galaxis auf eine Kugelschale und projiziere die in entsprechend großem Abstand um dein Raumschiff. Und zwar so, dass es sich immer im Mittelpunkt der Kugel befindet.

hustbaer

@pumuckl:
Ich glaube das nennt man Skybox

SeppJ

314159265358979 schrieb:

Ich habe dein Verfahren, SeppJ, leider nicht verstanden, denn wir haben Integralrechnung in der Schule auch noch nicht gemacht. Es wäre schön, wenn mir das jemand programmieren könnte.

Ich habe im Moment gerade nicht die Mittel grossartig zu coden und an meinen schon fertigen Code fuer das Kugelproblem komme ich von hier gerade auch nicht dran. Wenn es dir morgen reicht, dann kann ich es dann machen, falls ich dran denke.

314159265358979

Alles klar. Alleine, dass du mir deinen 10000ten Post widmest, ist mir eine Ehre

krümelkacker

Das einfachste ist natürlich X/Y/Z aus einem die Kugel umschließenden Würfel gleichverteilt zu ziehen und das solange zu wiederholen, bis der Radius passt.

Sonst könnte man das vielleicht noch so machen:
Richtung und Radius unabhängig ziehen.

Gleichverteilte Richtung:
X/Y/Z/ Gauß-verteilt ziehen und normalisieren

Radius ziehen:
Die Wahrscheinlichkeitsdichte für den Radius muss mit dem Radius quadratisch zunehmen, also:
P(r) \propto r^2
da die Kugeloberfläche ja mit dem Radius auch quadratisch zunimmt. Sonst bekommt man hinterher keine Gleichverteilung.

*Rewind*

Das grundsätzliche Problem bei dieser Aufgabe ist, dass der Gültigkeitsbereich (und damit auch das Maximum) jeder Dimension sich in Abhängigkeit von den Werten der anderen zwei Dimensionen ergibt.
Ich habe zwar nichts von dem verstanden, was SeppJ erklärt hat, aber aus dem genannten Grund sehe ich keine andere Möglichkeit als den Würfel.

krümelkacker schrieb:

Radius ziehen:
Die Wahrscheinlichkeitsdichte für den Radius muss mit dem Radius quadratisch zunehmen, also:
P(r) \propto r^2
da die Kugeloberfläche ja mit dem Radius auch quadratisch zunimmt. Sonst bekommt man hinterher keine Gleichverteilung.

Wie bastelt man das dann im Endeffekt? Wenn Radius zwischen r1 und r2 liegen soll: Kann man einfach gleichverteilt zwischen r1*r2 und r2*r2 ziehen und daraus die Wurzel als Radius nehmen? Müsste doch eine quadratische Verteilung ergeben.

hustbaer

[Rewind] schrieb:

Das grundsätzliche Problem bei dieser Aufgabe ist, dass der Gültigkeitsbereich (und damit auch das Maximum) jeder Dimension sich in Abhängigkeit von den Werten der anderen zwei Dimensionen ergibt.

Nur wenn du mit X/Y/Z Koordinaten rechnest.
Wenn du dagegen mit Kugelkoordinaten (quasi Längengrad/Breitengrad/Radius) rechnest, dann sind alle drei Werte unabhängig.

Einer der Winkel hat immer den Bereich [0°...180°), der andere immer [0...360°), und der Radius eben immer [0...Radius].

Wenn du dann noch wie krümelkacker & SeppJ beschreiben haben die Verteilung anpasst, dann bekommst du auch Punkte raus, die im X/Y/Z Raum gleichverteilt sind.

fdfdg schrieb:

krümelkacker schrieb:

Radius ziehen:
Die Wahrscheinlichkeitsdichte für den Radius muss mit dem Radius quadratisch zunehmen, also:
P(r) \propto r^2
da die Kugeloberfläche ja mit dem Radius auch quadratisch zunimmt. Sonst bekommt man hinterher keine Gleichverteilung.

Wie bastelt man das dann im Endeffekt? Wenn Radius zwischen r1 und r2 liegen soll: Kann man einfach gleichverteilt zwischen r1*r2 und r2*r2 ziehen und daraus die Wurzel als Radius nehmen? Müsste doch eine quadratische Verteilung ergeben.

Stell eine Wahrscheinlichkeitsdichte Proportional zu r^2 auf und bestimme den Faktor mit der Bedingung, dass das Integral der Dichte von r1 bis r2 gleich 1 wird. Dann wendet man die Inversionsmethode an.

http://de.wikipedia.org/wiki/Inversionsmethode

µ schrieb:

fdfdg schrieb:

krümelkacker schrieb:

Radius ziehen:
Die Wahrscheinlichkeitsdichte für den Radius muss mit dem Radius quadratisch zunehmen, also:
P(r) \propto r^2
da die Kugeloberfläche ja mit dem Radius auch quadratisch zunimmt. Sonst bekommt man hinterher keine Gleichverteilung.

Wie bastelt man das dann im Endeffekt? Wenn Radius zwischen r1 und r2 liegen soll: Kann man einfach gleichverteilt zwischen r1*r2 und r2*r2 ziehen und daraus die Wurzel als Radius nehmen? Müsste doch eine quadratische Verteilung ergeben.

Stell eine Wahrscheinlichkeitsdichte Proportional zu r^2 auf und bestimme den Faktor mit der Bedingung, dass das Integral der Dichte von r1 bis r2 gleich 1 wird. Dann wendet man die Inversionsmethode an.

http://de.wikipedia.org/wiki/Inversionsmethode

Oder kurz: Ja, das stimmt so, fdfdg. Bis auf den Typo und die exakten Exponenten halt

Wen es interessiert, hab hier kurz was gebastelt: http://ideone.com/IsAj2

Werner Salomon

otze schrieb:

Für manche Dinge braucht es nicht C++11 sondern einfache Mathematik.

genau .. und es kommt öfter vor, als man denkt

Folgender Code sollte in etwa das sein, was PI braucht. Er sollte Positionen in einer Kugelschale zwischen zwei gegebenen Radien R1 und R2 erzeugen, die im Volumen annähernd gleich verteilt sind.

Im Gegensatz zu pumuckl und fdfdg bin ich aber der Meinung, dass man beim Radius die 3.Potenz nehmen muss, da der Radius im Volumen mit der 3.Potenz eingeht. Bei einer Kugelschale mit den Innenradius 1 und Außenradius 2 befinden sich 50% des Volumens jenseits eines Radius von ca. 1,651 bzw. = 3.Wurzel( 1 + (2^3 - 1^3)/2 ).

#include <boost/random/uniform_real.hpp>
#include <boost/random/linear_congruential.hpp> // minstd_rand

#include <iostream>
#include <cmath>

class Vector // nur für Demozwecke
{
public:
    typedef double value_type;
    Vector( double x, double y, double z )
        : m_x(x), m_y(y), m_z(z)
    {}
    friend Vector operator*( double factor, const Vector& v )
    {
        return Vector( v.m_x * factor, v.m_y * factor, v.m_z * factor );
    }
    friend double abs( const Vector& v )
    {
        return std::sqrt( v.m_x * v.m_x + v.m_y * v.m_y + v.m_z * v.m_z );
    }
    friend std::ostream& operator<<( std::ostream& out, const Vector& v )
    {
        return out << v.m_x << " " << v.m_y << " " << v.m_z;
    }
private:
    double m_x, m_y, m_z;
};

namespace
{
    const double PI = std::acos(-1.0);
} 

// -- Verteilung zur Generierung eines 'gleichverteilten' Einheitsvektors
template< typename V, typename T = typename V::value_type >
class Sphere
{
public:
    typedef T input_type;
    typedef V result_type;

    Sphere() 
        : m_sin_phi( -1., +1. )
        , m_beta( -PI, PI )
    {}
    template< typename G >
    result_type operator()( G& gen )
    {
        using std::cos;
        using std::sin;
        using std::sqrt;
        // siehe auch http://de.wikipedia.org/wiki/Kugelkalotte
        const T sin_phi = m_sin_phi( gen ); // sin( Breitengrad ) ist gleichverteilt, da H=2*PI*sin(phi)
        const T cos_phi = sqrt( 1 - sin_phi*sin_phi );
        const T beta = m_beta( gen ); // Längenwinkel
        return result_type( cos( beta )*cos_phi, sin( beta )*cos_phi, sin_phi );
    }
private:
    boost::uniform_real< T > m_sin_phi;
    boost::uniform_real< T > m_beta;
};

// --   Verteilung zur Generierung eines Radius im Intervall [x0, x1)
//      mit dem Ziel im Volumen einer Kugelschale eine 'Gleichverteilung' zu erreichen
template< typename T >
class Cubic
{
public:
    typedef T input_type;
    typedef T result_type;
    Cubic( result_type x0, result_type x1  )
        : m_distribution3( x0 * x0 *x0, x1 * x1 * x1 )
    {}
    template< typename G >
    result_type operator()( G& gen )
    {
        using std::pow;
        return pow( m_distribution3( gen ), result_type( 1/3. ) );
    }
private:
    boost::uniform_real< T > m_distribution3;
};

int main()
{
    using namespace std;
    const double R1 = 1.5;
    const double R2 = 4;
    boost::minstd_rand gen;
    Sphere< Vector > sphere;
    Cubic< double > cubic( R1, R2 );
    cout << fixed;
    cout.precision( 2 ); // 2 NK reichen hier
    for( int i=0; i<20; ++i )
    {
        Vector pos = cubic( gen ) * sphere( gen );
        cout << pos << "  " << abs(pos) << endl;
    }
    return 0;
}

Gruß
Werner

@Edit: in den Zeilen 42 und 58 den Typ 'Vector' gegen das Template-Argument V bzw. result_type ersetzt

SeppJ

@Werner: Bist du da sicher, dass die Verteilung deiner Längen stimmt? Dein Code ist mir zu verschlungen, um ihn jetzt noch zu entziffern. Du ziehst nur die dritte Wurzel, oder? Das würde nämlich nicht stimmen, da du dann nämlich entweder die Verteilung bei R1 als 0 hast oder die Verteilung verschoben (und somit die falsche Steigung) hast, je nachdem, was du da genau machst (sorry, wenn ich gerade keine 100 Zeilen Templates dekodieren möchte - habe heute Abend schon genug Denksport betrieben).

Ich habe jetzt stundenlang an der exakten Verteilung der Längen rumgetüftelt, ich kann nicht glauben, dass du so einen einfachen Code heraus bekommst, ich hingegen dieses Monstrum:

#include <iostream>
#include <vector>
#include <functional>
#include <numeric>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <complex>
using namespace std;

int main ()
{
  vector<double> bins(1000,0.);
  const double r1 = 2;
  const double r2 = 5;
  double r1_3 = r1*r1*r1;
  double len = r2-r1;
  double cut = r1*r1 / (pow(r2,3)/3. - r1_3/3.);  // Wert der Verteilung bei r1 -> Offset
   double fac1 = cut * len;   // Der Offset muss noch mit der Länge normiert werden
  for(unsigned i = 0; i < 1000000; ++i)
    {
      double a = rand() *1. / RAND_MAX;
      if (a < fac1)              // Gleichverteilungsanteil der Verteilung
        a = (a/fac1)*len + r1;   // Dies ist noch einfach...
      else                       // ...aber jetzt kommt die x^2 Funktion die da drüber gelegt ist: OMG!
        {
          a = (a-fac1) / (1-fac1); // a ist wieder Einheitsverteilung
          // Herleitung: Mathematica, Zettel und Stift und viel Geduld
          double fac2 = len*len*(2*r1+r2)*a;
          complex<double> fac3 = pow(-2*r1_3+fac2+sqrt(complex<double>(fac2*(-4*r1_3+fac2),0.)),1./3.);
          complex<double> result = 0.5 * (
                                          pow(2,2./3.)*fac3
                                          +2*r1*
                                          (-1. + pow(2,1./3.)*r1 / fac3 )
                                  );
          a = result.real();   // Zusatzwissen: Imaginärteile heben sich genau weg
          a += r1;
        }
      ++bins[a*100];
    }
  double factor = 1./(accumulate(bins.begin(), bins.end(), 0.)*0.01);
  transform (bins.begin(), bins.end(), bins.begin(), bind1st(multiplies<double>(),factor));
  for (unsigned i = 0; i < bins.size(); ++i)
    {
      cout << i / 100. << '\t' << bins[i] << '\n';
    }
}

Dieser Code ist erfolgreich geprüft (beziehungweise er ist mein Prüfcode, einfach mal die Werte plotten, wer es nachvollziehen möchte) gegen die theoretische Wahrscheinlichkeitsverteilung von

P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

Irgendwie muss man auch die komplexwertige Zwischenrechnung vereinfachen können, aber das würde ich selber gerne sehen und Mathematica schafft es nicht.

Die Längen kann man dann mit Punkten auf einer Einheitssphäre kombinieren, wie hier schon gezeigt wurde. Herleitung findet man hier:
http://mathworld.wolfram.com/SpherePointPicking.html
Das ist zum Glück ein ganzes Ende einfacher.

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat. Bei einer dünnen Schale sinkt der Anteil sehr schnell nahe 0. Und man hat natürlich auch milde Rechenkosten für die Abstandsberechnung.

Entschuldigung auch an alle Leser, dass mein Code ebenso unleserlich ist wie Werners, wenn auch aus weitaus weniger künstlerisch anspruchsvollen Gründen :p

P.S.: Herleitung ohne funktionierende LaTeX-Tags möchte ich niemandem antun. Im Prinzip ist es nichts kompliziertes, bloß eklig lang ohne CAS (und mit CAS auch ). Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, das einzige was glaube ich noch nicht genannt wurde ist, wie man das macht, wenn die Inverse der kumulierten Wahrscheinlichkeitsverteilung bei 0 nicht 0 ist. Das ist das was hier im ersten if-Zweig behandelt wird, indem eine Gleichverteilung passender Höhe "abgeschröpft" und einfach (nach ein bisschen Skalieren) durchgeleitet wird. Der zweite if-Zweig ist dann der Teil der übrig bleibt und mit herkömmlichen Mitteln behandelt werden kann (auch wenn die herkömmlichen Mittel hier ganz schön komplexe Rechnung erfordern). Die Tatsache, dass die Imaginärteile bei negativer Wurzel sich gegenseitig weg heben habe ich zugegebenermaßen nur durch Ausprobieren getestet, aber konstruktionsbedingt muss sich das auch so ergeben, weil die Methodik an sich bewiesen ist und in der Rechnung selbst bis dahin keine Fehler sind (weil der Computer gerechnet hat ).

krümelkacker

SeppJ: Probier's doch mal so zu trenne, wie ich es vorgeschlagen habe. Wenn Du X/Y/Z Gauß-verteilt ziehst und dann normalisierst, hast Du schon auf der Kugeloverfläche eine Gleichverteilung. Dann fehlt nur noch der Abstand zum Ursprung. Ich habe den letzten Schritt nicht zu Ende gerechnet, aber so kompliziert dürfte das nicht sein.

*Rewind*

@Werner: wie lange hast du für den Code gebraucht?

otze

Bashar schrieb:

otze schrieb:

Für manche Dinge braucht e snicht C++11 sondern einfache Mathematik.

"einfache"? Ich hab grad keine Ahnung, wie das mathematisch gehen könnte.

ich meinte in der Tat die ziehen und Wegwerfmethode. Bei den typischen schnellen Zufallszahlengeneratoren ist die analytische Lösung wahrscheinlich langsamer als wegwerfen und neu ziehen. In dem Setup hier ist das auch gut möglich. Ist auch nicht weniger mathematisch, läuft dann halt unter Akzeptanzkriterium.

SeppJ

krümelkacker schrieb:

Ich habe den letzten Schritt nicht zu Ende gerechnet, aber so kompliziert dürfte das nicht sein.

Ja, das hatte ich auch gedacht:

SeppJ schrieb:

hier ist das bei der parabelfoermigen Verteilung der r-Koordinate aber noch einfach

Dankste .

Die Gleichverteilung auf der Kugelfläche ist nämlich der einfache Teil, auch wenn es zuerst nicht so aus sieht.