Kleine Sachen mit Python

'Ray Tracing' mit sektorierten Segmenten




Etwas schlichte Geometrie

Geraden lassen sich mathematisch auf viel­fältige Weise beschreiben. Meine Wahl zielt auf die Einfach­heit der Beschrei­bung: Ich beschreibe eine Gerade durch einen belie­bigen Punkt P0, Startpunkt genannt, und durch einen normierten Rich­tungs­vektor dg der Länge 1, der eine der beiden Rich­tun­gen der Geraden nach Be­lieben auszei­chnet. In der Skizze unten sind beide Ob­jekte blau einge­färbt.

Mit dieser Beschrei­bungs­weise gibt es keine Sonder­fälle wie die senk­rechte Gerade mit der Steigung ∞, es gibt keine häß­lichen Winkel und keine Winkel­funktionen - das Ganze ist einfach und un­mit­tel­bar anschau­lich, eben Geometrie.

Die Punkte Pg auf der Geraden g werden durch eine Para­meter­glei­chung definiert:

Der Parameter t hat die ein­fache geo­metrische Inter­preta­tion der Länge des Diffe­renz­vektors (Pg(t)-P0).

Gegeben sei nun ein belie­biger Punkt P1. Gesucht wird der zu P1 nächst­gele­gene Punkt P2, der auf der Gera­den g liegt.

Diesen Geradenpunkt P2 bestimme ich, indem ich den Diffe­renz­vek­tor (P1-P0) auf die Gerade g pro­jiziere, um den Vektor (P2-P0) zu erhalten und damit als End­punkt dieses Vek­tors den gesuch­ten Punkt P2, der als Geraden­punkt die Geraden­glei­chung mit einem Para­meter t2 erfüllt:

Den Parameter t2 berechne ich aus den Aus­gangs­größen mit­hilfe des Skalar­produk­tes zwischen Vek­toren zu:

Mein Wunsch war es nun, solch wunder­schön einfach aus­sehen­den Glei­chungen mög­lichst Eins zu Eins in Pro­gramm­code über­führen zu können. Wenn das gelingt, kann ich die ma­the­ma­ti­schen Glei­chungen leicht mit dem Pro­gramm­code verglei­chen und habe eine häß­liche Fehler­quelle weniger,

näm­lich die des fehler­behaf­tete Um­setzens von mathe­matischen Glei­chungen in Programm­code. - Und das gelingt mit Python, dem neben­stehen­den Bei­spiel zu ent­nehmen!

Vectors.dotOp bezeich­net hier das Ska­lar­pro­dukt zwi­schen zwei Vek­toren. In der Kürze liegt be­kannt­lich die Würze.

Und weiter: Das Auge erkennt sofort, dass der grüne Punkt P1 im Bild oben nicht auf der Gera­den g liegt. Wann liegt also ein Punkt P1 auf der Gera­den? Dazu berech­net man den zu P1 nächst­gele­genen Punkt P2 auf der Gera­den, und vergleicht diesen mit dem Aus­gangs­punkt P1. Wenn nun P1 = P2 erfüllt ist, genau dann liegt der Aus­gangs­punkt auf der Ge­raden.

Dazu noch ein Beispiel zur Imple­mentie­rung. Liegt ein Punkt etwa auf einer Gera­den? Die entsprechende Funktion steht nebenan.

'==' ist hier die implemen­tierte Ver­gleichs­ope­ra­tion für zwei Punkte.

Eine Software-Funktion baut auf der anderen auf. Effi­zienz spielt keine programm­gestal­tende Rolle, die Programm­sicher­heit ist hier das höhere Gut.

Und weiter: Strahlen werden durch den aus­ge­zeich­ne­ten Startpunkt und einen Richtung­vektor be­schrie­ben. Wie erhalte ich den Schnitt­punkt zweier Strahlen?

In der 'Hauptsache' wie folgt: Ich verlän­gere die Strahlen zu Gera­den und wenn deren Schnitt­punkt auf den beiden Aus­gangs­strah­len liegt, dann ist das der ge­wünschte Schnitt­punkt der beiden Strahlen.

Wie man am Code-Beispiel unten leicht erken­nen kann, wird aus der Haupt­sache schnell eine Ne­ben­sa­che, die Spezial­fälle NoPoint und In­fi­ni­ty­Points machen die Arbeit und den Test aus.

Aber dennoch: Man braucht keine Winkel­funktion atan2 und keinen Winkel­wirrwar, um das Problem mit einigen Skizzen und Tests zu durch­dringen.


 

Und weiter: Ein Segment, also eine Strecke wird durch zwei Punkte bestimmt. Wann liegt ein Punkt auf einer Strecke? Ich verlän­gere die Strecke zu­nächst zu zwei Strah­len, in beide Rich­tungen, sodass die Strecke auch auf beide Strah­len liegt. Wenn der Punkt nun auf den beiden so erzeug­ten Strah­len liegt, so liegt er auch auf der Strecke. Ganz einfach ist das.

Zuletzt: Wann schneidet ein Strahl eine Strecke? Ich verlän­gere die Strecke in beide Rich­tungen zu zwei Strahlen und unter­suche nun das Schnitt­verhal­ten der drei Strahlen.

Das Eine baut auf dem Ande­ren auf. Es entsteht eine Hie­rachie von ein­fachen Funk­tionen. Der Test wird so zum Kinder­spiel.

Mathematk macht Spaß. Nicht jeder hat in der Schule aufge­passt. Ich hatte einige gute Lehrer. Die moderne, koor­dinaten­freie Vek­tor­rech­nung war damals eine Offen­barung und wirk­lich etwas fürs Leben. Herrn Rum­berger sei gedankt.

Hier habe ich in einer Hand­reichung ein kleines mathe­mati­sches Grund­gerüst zusam­men getra­gen, erstellt mit Libre­Office Writer & Draw.

Vektoren, Punkte und mehr
(PDF-Format)
Formeln zur zwei­dimen­sio­nalen Geomerie


Wie hilft mir die Sprache Python?

Python hilft mir weitaus mehr als ich zunächst erwartet hatte, denn ich hatte Python bisher nur als klas­sen­lose Skript­sprache einge­setzt.

Und weiter: Geome­trische Objekte werden zu Klassen, einge­baute Funk­tionen lassen sich über­lagern, sodass sich etwa zum Ver­gleich von Punk­ten der Operator '==' (intern '__eq__') und zur Differenz­bildung von Punkten der Operator '-' (intern: '__sub__') verwenden lässt.

Ich hätte gerne für die Multi­plika­tion eines Vektors mit einem Skalar, etwa 2*v oder auch v*2, und für das Skalar­produkt zweier Vektorn, etwa v1*v2, eben auch das Operator-Symbol '*' (intern: '__rmul__') verwendet, Python hat sich aber bisher gestreubt. Deshalb gibt es halt die Funktion dotOp. Ich bin zufrieden damit.


Die Klasse Points ist ganz über­sicht­lich (siehe unten). Eine Instanz oder ein Objekt p mit den Koor­dinaten x=1, y=1 wird mit der Anwei­sung p=Points(1,1) erzeugt, wobei die interne Funk­tion '__init__' aus­geführt wird.


 

Ich schließe, ganz der konser­vativ Program­mierer, Zeilen mit einem Semi­kolon und längere Blöcke mit einem Kommen­tar ab, dann fällt der Code optisch nicht so ausein­ander. Auf Unter­striche in den Bezeich­nern ver­zichte ich und wo weißer Raum nicht nötig ist, lasse ich ihn weg. Das ist die dichte Packung des Codes, die ich mir bei Sci­Lab ange­wöhnt habe.

(Am Rande: Ich halte es immer noch eine recht abstruse Idee, die Syntax an weißen Raum zu binden. Ich habe mehrere Anläufe gebraucht, um mich daran zu gewöh­nen.)

Man verschiebt einen Punkt im Raum, hier in der Ebene, durch die Addi­tion eines Vektors, etwa P2=P1+v, dabei wird intern die Funk­tion '__add__' auf­ge­ru­fen.

Nützlich finde ich auch die Mög­lich­keit, bei (im Argu­ment) symme­trischen Funk­tionen einen klassi­schen Funktions­aufruf ver­wenden zu können; ein Bei­spiel für den Ab­stand zweier Punkte:

distance=Points.distanceOf(p1,p2);

sieht gut aus, denn die beiden Punkte tau­chen in der glei­chen Weise als Funk­tions­para­meter auf, aber

distance=p1.distanceOf(p2);

verletzt die Symme­trie der Funk­tion und entfernt sich von der mathe­matischen Schreib­weise zu sehr, der eine Punkt taucht als Objekt auf, der andere als Funk­tions­para­meter.

Sehr nützlich und hilf­reich ist die Mög­lich­keit, ein Modul sofort als Haupt­programm aus­füh­ren zu kön­nen. Daher braucht man bei über­schau­baren Modu­len kein ei­gen­stän­di­ges Test­programm, man kann den Test­code gleich im Modul un­ter­brin­gen.

So gehen Imple­mentie­rung und Test eine enge Ver­bindung ein. Und ich finde es eine sehr gute Idee, den Test­code beim Modul zu be­lassen. Bei den vielen Frei­heiten, die Python zulässt, ist das ein Schritt auf dem Weg zu robus­ter Soft­ware.

Ein wenig mehr der Typ­si­cher­heit wäre ei­gent­lich nicht schlecht. Ada nimmt da eine ex­tre­me Seite ein, Python bietet ein an­de­res Extrem.


 

Sind die Module umfang­reicher, ist es empfeh­lens­werter, für den Modul­test gleich ein Test­modul zu schrei­ben und sich dabei vom ‚Unit Testing Frame­work’ Pythons unter­stützen zu lassen. Man imple­men­tiert kleine Test­routinen, den Test­treiber stellt Python bereit. Das ist unkom­pli­ziert und herr­lich unauf­wändig.

Und der statischen Typ­sicher­heit kann man - wenn denn ge­wünscht - auch ein wenig auf die Sprünge hel­fen: Da gibt es ja die Funktions­annota­tionen und seit der Py­thon-Ver­sion 3.5 das Modul ‚typ­ing’. PyCharm meckert so­gar bei Un­stim­mig­keiten, ich bin begeis­tert.

 
class PolygonTrains(object): # Polygonzug

def __init__(self,ptL:List['Points'])->None:
   self.Points=[];
   self.Points+=ptL;
# end __init__

Ich habe mich aus Neugier ein paar Wochen lang mit Has­kell be­fasst; dessen star­kes Typ­system fand ich ele­gant, es hat mich beein­druckt. Aber ein paar Wochen rei­chen bei Has­kell leider nicht aus, um diese rein funk­tionale Pro­gram­mmier­sprache hinrei­chend zu beherr­schen und auch ele­gant ein­setzen zu kön­nen. Das war die Chance für Python, der syn­tak­tische wei­ßen Raum war nun nicht mehr so ganz befremd­lich!


Hilft mir das Werkzeug «PyCharm»?

Manchmal ist man einfach zu genüg­sam! Ich hatte mich zwar der­einst nach einer zeit­gemä­ßen Ent­wick­lungs­umge­bung für Python umge­schaut, wurde aber nicht so recht fün­dig und be­gnügte mich dann lange mit dem schlich­ten Edi­tor «Idle» aus dem Py­thon-Instal­lations­paket.

Eine Kollegin in Sachen Ada aus alte Tagen machte mich im letz­ten Jahr so neben­bei auf «PyCharm» auf­merk­sam. Was einen denn auch gleich erfreut: Eine ‚Com­munity Edi­tion’ des Werk­zeugs wird ins Netz gestellt, ein feiner Zug von «Jet­Brains».

Mein erster Eindruck von PyCharm war, dass da zuviel der Meckerei sei — und es gelang mir ein­fach nicht, die Schrift­größe des Edi­tors zu ändern, was mich doch sehr är­gerte.

Nach einigen Wochen habe ich mir das Werk­zeug noch ein­mal vor­genom­men: Und siehe da, bald waren die Mecke­reien weg­kon­figu­riert und die Schrift sperrte sich auch nicht mehr ge­gen ein Ver­größern; und mein häus­liches «Mercurial» wird sogar ohne Um­stände ins Werk­zeug einge­baut!

Es ist alles dabei, was das Ent­wick­ler­herz be­gehrt: Eine moder­ne, inte­grier­te Ent­wick­lungs­umge­bung haben wir da an der Hand, die nichts kostet und die wirk­lich hilft — nie, nie mehr das all­zu­schlichte Idle!

In der Tat, PyCharm macht mich zur Zeit mächtig an, ich habe mei­nen gesam­ten Python-Code einer stren­gen Revi­sion unter­zogen - ja, es ist schon er­schre­ckend, welch schrä­gen Code der Python-Inter­preter mit dem rich­tigen Ergeb­nis aus­führt. Aller­dings konnte ich mich nicht durch­ringen, end­lich den Stil-Stan­dard von Python (PEP 0008) anzu­wenden. Deshalb sieht mein Code immer noch so aus, wie ich ihn hier präsen­tiere, schlicht nicht stan­dard­kon­form - ich nenne es die hoch ver­dich­tete Form; Leer­raum und Unter­striche stö­ren in mei­nen Augen die Wohl­gestalt meines Python-Codes :-)


Die Anwendung – Ein Beispiel aus der Praxis

Ein Bewohner von Flachland möchte einen Teil eines Polygonzuges ge­kenn­zeich­net haben. Die Teilstücke des gegebenen Polygonzuges sollen sich nicht über­schnei­den dürfen. Der Bewohner spezifizierte den gewünschten Bereich auf mehr als zwei Seiten mit Hilfe von Pseudocode und Winkeln und Winkelfunktionen - ein großer, häßlicher Haufen undurchschaubarer Algorithmik!

Ich hielt dagegen: Ein geometrisches Problem müsse man mit Geometrie und nicht mit Algorithmen in Pseudocode angehen. Es brauchte einige Überzeugungsarbeit, für deren Vorbereitung ich vierzehn Tage Weihnachtsurlaub (2004 / 2005) gerne opferte. Dann war es geschafft!

Meine Spezifikation des Problems war es, eben den für den Bewohner sichtbaren Teil des Polygonzuges zu finden. Dieser eine Satz reicht vollständig aus, mehr braucht es nicht. Der geometrische Algorithmus als Anleitung für die Implementierung geht dann so: Man zerlege den Polygonzug in Seg­men­te (Teilstücke), auf dass alle Segmente Sektoren zugeordnet werden können. Dann entnehme man jedem Sektor das Segment, das dem Beobachter am nächsten ist – und schon hat man den sichtbaren Bereich in den Händen.

Mit etwas Geometrie, gegossen in Python, und den Listen von Python ist die Im­ple­men­tie­rung auch schnell erledigt.

Ich spendiere mir noch zwei Testfälle: Ein Polygonzug, der von unterschiedlichen Standorten betrachtet wird und ein zweiter, der für einen Sonderfall steht, hier liegen Segmente auf Sektorgrenzen.

Der rote Kringel mimt den Beobachter, der vom Beobachter aus der sichtbare Teil ist rot eingefärbt. Die Grafiken habe ich mit «GeoGebra» erstellt.


PolygonTrain points:

3 1 | 2 2 | 4 2 | 3 3 | 1 3 | 0 4 | 5 3 | 6 3 | 6 4 |


 


 

Observer position: 0 0

Visible part - points: 3

4 2 | 3 2 | 1 1


 


 

Observer position: 2 0

Visible part - points: 5

6 3 | 5 3 | 4 2 | 3 2 | 1 1


 


 

Observer position: 3 0

Visible part - points: 6

6 3 | 5 3 | 5 4 | 4 2 | 3 2 | 1 1


 


 

Observer position: 4, 0

Visible part - points: 7

6 3 | 5 3 | 5 4 | 4.0 3.5 | 4 2 | 3 2 | 1 1


 


 

Observer position: 5 0

Visible part - points: 7

6 3 | 5 3 | 5 4 | 3.4 3.2 | 4 2 | 3 2 | 1 1


 


 

Observer position: 6 0

Visible part - points: 7

6 3 | 5 3 | 4.71 3.86 | 3 3 | 4 2 | 3 2 | 1 1

PolygonTrain points:

3 1 | 2 2 | 4 2 | 1 3 | 0 4 | 5 3 | 6 3 | 6 4 |


 


 

Observer position: 4 0

Visible part - points: 5

6 3 | 5 3 | 4.0 3.2 | 4 2 | 2 2


Eine kleine Widmung

Ich widme die obige Bildkomposition aus vierzehn handgefertigten Einzelbildern
einem netten Arbeitskollegen, dem «unbeschwerten» Herrn Bernd M.

Es waren wirklich hitzige Diskussionen damals vor fast zehn Jahren.
Er hat das Ganze am Ende in Ada gegossen und mit hübschen Kommentaren und kunstvollen Grafiken verziert.


 
----------------------
-- local Nearest_Point
----------------------
-- this function returns the nearest point 
-- to point P on ray R
--  this point is 
--  a) the orthogonal projection of the 
--     point on the ray or
--  b) the start point of the ray 
--     (if no projection is possible 
--
-- a)               b)
--           R              R
--           ^              ^
--           |              |
--     P     |              |
--     *-----+ Pi           |
--           |              |
--           |              + Sp
--           * Sp       P   
--                      *
--     return Pi        return Sp
--
----------------------
                        
 


 

 

Ach ja, zu guter Letzt sei noch eine kleine Anmerkung nachgeschoben:
(Punkt+Punkt)=Punkt
und
(Punkt-Punkt)=Punkt
geht natürlich gar nicht!

Aber die komplexe Software läuft dennoch
— fehlerfrei und wie geschmiert!
Wohl auch in fünfzehn Jahren noch.

 



Der Quellcode für Python 3.5

geometry2d.html

geometry2d.py

test_geometry2d.html

test_geometry2d.py

visible_segments.html

visible_segments.py

 

Hand­reichung
          Vektoren, Punkte und mehr
         
 
(PDF-Format)
Formeln zur zweidimensionalen Geometrie



© 2015 Bernd Ragutt
Alle Rechte vorbehalten
 ... hier kann man hinschreiben letzte Änderung: 21. März 2018
Kruschtkiste