Zufällige 3D Vektoren innerhalb einer Kugel
-
bliblablid schrieb:
Ich habe jetzt nicht alle Post gelesen daher könnte es sein, dass mein Vorschlag schon gesagt wurde.
Ja, gleich als erste Antwort.
-
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.
Grüße
-
pumuckl schrieb:
KLar, nur wirfst du beim vorliegenden Fall, wo der innenradius = 1/2 Außenradius der Kugelschale ist, knapp 46% der Ergebnisse weg
Richtig. Aber wer sagt, dass das schlimm ist? Wenn am Anfang zur Initialisierung die Kugel in 2,17 ms statt in 1 ms mit zufälligen Punkten gefüllt wird, wen interessiert das?
Wenn das zu langsam ist, kann man sich immer noch was besseres ausdenken.
-
Muss es denn tatsächlich eine exakte Gleichverteilung in einer Kugel sein? Ich mein, die Sterne am Nachthimmel sind jetzt auch nicht wirklich gleichverteilt, sogar ziemlich weit davon entfernt...
-
Naja, es sollten zumindest nirgens Lücken zu sehen sein.
-
Wie wär's mit nem Würfel + Sichtschutz?
-
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.
-
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.
-
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.
-
@pumuckl:
Ich glaube das nennt man Skybox
-
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.
-
Alles klar. Alleine, dass du mir deinen 10000ten Post widmest, ist mir eine Ehre
-
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 normalisierenRadius 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.
-
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.
-
[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.
-
µ 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.
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
-
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
-
@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 r3)/(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 ).
