mathematische Frage von Anfänger

[nudelauflauf]

Hi Leute,

mein Problem ist folgendes:

Ich möchte in einem Feld beliebiger Größe (bspw: 2000x3000 Punkte) einen Kreis mit beliebigem Radius (bspw: r = 100) an einen bestimmten Ort setzen (bspw: x=1000, y=1000). Ich möchte den Kreis nicht angezeigt haben oder sonst was, ich brauch nun nur alle Koordinaten aller Punkte innerhalb dieses Kreises als Auflistung, meinetwegen als zusammenhängender String mit Zeilenumbruch dazwischen.

Für ein Quadrat oder Rechteck ist diese Berechnung ja nicht sonderlich schwer indem man einfach mit 2 Schleifen vorgeht welche einmal eine X-Koordinate und innerhalb dieser Schleife jeweils die Y-Koordinate hochzählt und immer die aktuellen Werte abspeichert. Wie funktioniert das ganze aber mit einem Kreis?

Wäre es evtl eine Möglichkeit einen Kreis zu berechnen, die Koordinaten des Umfangs zu speichern, jeweils nach Punkten mit gleichem Y-Wert zu suchen und dann alle X-Werte dazwischen hochzählen lassen?

Gibt es evtl eine einfachere Möglichkeit?


Vielen Dank schonmal!

[lunar]

@nudelauflauf Wofür benötigst Du eine solche Liste?!

Wie dem auch sei, Dein Ansatz, ausgehend von Koordinaten des Kreisumfangs entlang die Punkte entlang der Parallelen zur X-Achse zu berechnen, erscheint mir vernünftig. Um die natürlichen Koordinaten des Kreisrands zu berechne, kannst Du den Bresenham-Algorithmus verwenden.

[BlackJack]

@nudelauflauf: Du könntest für jeden Punkt bestimmen ob er innerhalb des Kreises liegt und dann in Dein Ergebnis aufnehmen. Das bekommt man sehr einfach heraus in dem man den Abstand vom aktuellen Punkt zum Kreismittelpunkt berechnet. Wenn der Abstand kleiner (oder gleich) dem Radius ist, dann ist der Punkt im Kreis (oder auf dem Rand).

[/me]

@nudelauflauf: Du könntest für jeden Punkt bestimmen ob er innerhalb des Kreises liegt und dann in Dein Ergebnis aufnehmen. Das bekommt man sehr einfach heraus in dem man den Abstand vom aktuellen Punkt zum Kreismittelpunkt berechnet. Wenn der Abstand kleiner (oder gleich) dem Radius ist, dann ist der Punkt im Kreis (oder auf dem Rand).

Das lässt sich dann noch optimieren indem du ein den Kreis umschließendes und ein vom Kreis genau umschlossenes Quadrat nimmst. Alles was außerhalb des äußeren Quadrats liegt gehört nicht dazu und alles was innerhalb des innenliegenden Quadrats liegt gehört zum Kreis. Du brauchst dann nur noch die Punkte im Bereich zwischen den Quadraten betrachten.

[nudelauflauf]

Ich hab jetzt mal n bisschen rumgespielt. der Code ist aber dreckslangsam.
Da ich den Mittelpunkt kenne und den Radius habe gehe ich einfach ähnlich der Lösung mit dem Quadrat vor.
Ich bau mir eine Linie in x-Richtung (l1) die immer 2 Punkte länger ist wie der Kreisdurchmesser, so stell ich sicher dass ich nie innerhalb bin. Diese Linie lass ich dann in einer Schleife nach unten laufen und bestimme immer den Schnittpunkt der Linie mit dem Kreis. Außer in 2 Fällen (ganz oben und unten am Kreis) bekomme ich dann immer 2 Punkte heraus und kann quasi alle Punkte dazwischen in meine Liste aufnehmen.
Das Problem was ich jetzt habe, ist, dass ich einen Punkt herausbekomme mit jeweils einer X- und Y-Koordinate (logisch).
Wie kann ich nun die X- Koordinate extrahieren??? Die Y-Koordinate kenn ich ja eh. Wenn ich aber die 2 X-Koordinaten zu den Schnittpunkten habe dann ist es ja ein leichtes das jeweils mit der zugehörigen Y-Koordinate zu verknüpfen und als String zu speichern.
Deshalb meine Fragen nun:
1. Wie bekomm ich die X-Koordinate/n aus den Punkten aus intersection?
2. Gibt es eine eventuell schnellere Lösung?

Code: Alles auswählen

from sympy.geometry import Point, Line, Circle, intersection

x, y, r = 2000, 2000, 200

p1= Point(x,y)
c = Circle(p1, r)


for t in range(-1, 2*r+1):
      p2, p3 = Point(x-r-1, y-r+t), Point(x+r+1, y-r+t)
      l1 = Line (p2, p3)
      p4 = intersection(l1, c)

[pillmuncher]

Gibt es eine eventuell schnellere Lösung?

Vermutlich diese hier (hab's nicht gemessen):

Code: Alles auswählen

from PIL import Image
from math import sqrt

def fill_circle(x0, y0, radius):
    f = 1 - radius
    ddF_x = 0
    ddF_y = -2 * radius
    x = 0
    y = radius
    while x < y:
        if f >= 0:
            y -= 1
            ddF_y += 2
            f += ddF_y
        x += 1
        ddF_x += 2
        f += ddF_x + 1
        x1 = x
        for y1 in xrange(x + 1, y):
            yield x0 + x1, y0 + y1
            yield x0 + x1, y0 - y1
            yield x0 - x1, y0 + y1
            yield x0 - x1, y0 - y1
            yield x0 + y1, y0 + x1
            yield x0 + y1, y0 - x1
            yield x0 - y1, y0 + x1
            yield x0 - y1, y0 - x1
    sq = sqrt(.5)
    for i in xrange(1, radius):
        j = i * sq
        yield x0, y0 + i
        yield x0, y0 - i
        yield x0 + i, y0
        yield x0 - i, y0
        yield x0 + j, y0 + j
        yield x0 - j, y0 - j
        yield x0 + j + 1, y0 - j
        yield x0 - j + 1, y0 + j
    yield x0, y0

L = list(fill_circle(250, 250, 200))
assert len(L) == len(set(L))

img = Image.new('RGB', (500, 500))
pixel = img.load()
for coor in L:
    pixel[coor] = 0xff, 0xff, 0xff
img.show()

Es basiert auf dem von lunar erwähnten Bresenham-Algorithmus. Nur, dass hier nicht nur Bogensegmente erzeugt werden, sondern auch die dazugehörigen "Kuchenstücke". Jeder Punkt wird dabei nur genau einmal generiert.

Und um lunars Frage zu wiederholen: Wofür benötigst Du eine solche Liste?!

[mutetella]

Ok, als Ahnungsloser sollte man eigentlich seine Klappe halten, aber zufällig rutschte gerade dieser Blogeintrag an mir vorbei.
Sollte das darin erwähnte PIL dasselbe sein, das pillmuncher in seinem Beispiel verwendet, ist die Sache ja vielleicht für Euch interessant.

mutetella

[nudelauflauf]

Danke schonmal. Ich hab mir noch eine andere Lösung ausgedacht welche ich ersteinmal ausprobieren möchte. Werde diese dann mit der obigen auf Speed vergleichen.

Zur Frage für was:
Ich möchte ein Streamingfile generieren, welches ich einem FIB übergebe. Dieser soll dann an jedem dieser Punkte einmal "feuern".

[nudelauflauf]

Ich denke das ist eine gute und schnelle Lösung geworden. Ich nehme mir den ersten Schnittpunkt und verschiebe eine kurze Linie nur in eine Richtung um einen Schnittpunkt zu finden. Da der Kreis ja symmetrisch ist kann ich aus einem Schnittpunkt schnell den 2. bestimmen und die Punkte zwischen Schnittpunkt und X-Koordinate des Mittelpunkts aufnehmen. Das ganze spiegel ich dann einfach 2 Mal und schon hab ich alle Punkte erschlagen.

Code: Alles auswählen

  
    from sympy.geometry import Point, Line, Circle, intersection  
    p1= Point(x,y)
    c = Circle(p1, r)
    q=x-r
    k=1
    out_str=''

    for t in range(0, r):
        p2, p3 = Point(q-2, y+t), Point(q+5, y+t)
        l1 = Line (p2, p3)
        pl = intersection(l1, c)
        p4 = pl.pop(0)
        q = int(round(p4.x))
        q2 = 2*x-q
        z=0
        while q+z < x:
            out_str = out_str + ' ' + str(q+z) +' ' + str(y-t) + '\n'
            out_str = out_str + ' ' + str(q+z) +' ' + str(y+t) + '\n'
            out_str = out_str + ' ' + str(q2-z) +' ' + str(y-t) + '\n'
            out_str = out_str + ' ' + str(q2-z) +' ' + str(y+t) + '\n'
            k+=4
            z+=1

[BlackJack]

@nudelauflauf: Das wiederholte zusammensetzen einer immer länger werdenden Zeichenkette mit ``+`` ist ineffizient, weil potentiell jedes mal eine neue Zeichenkette angelegt wird in welche die Daten der beiden alten Zeichenketten kopiert werden. Man sammelt deshalb besser die ganzen Teilzeichenketten in einer Liste und fügt sie am Ende mit der `join()`-Methode auf Zeichenketten zusammen.

Zusammensetzen von Werten und Zeichenketten mittels ``+`` und `str()` ist auch eher „unpythonisch” und erinnert mehr an BASIC. Python bietet die `format()`-Methode oder den ``%``-Operator auf Zeichenketten.

Der `pop()`-Aufruf auf `pl` ist unnötig. Hier könnte man einfach auf das erste Element zugreifen statt es zu entfernen und damit alle folgenden kopieren zu müssen. An dieser Stelle könnte man den zweiten Punkt aber eigentlich auch benutzen, statt ihn durch Spiegelung selbst noch einmal zu berechnen. Andererseits könnte man die Linie auch kürzer machen, so das grundsätzlich nur ein Schnittpunkt mit dem Kreis besteht, wenn man sowieso nur einen haben möchte.

Die ``while``-Schleife könnte man durch eine ``for``-Schleife ersetzen, und zwar über die gesamte Breite der Linie und sich damit eine der Spiegelungsberechnungen, also zwei Zeilen, *in* der Schleife sparen.

Wozu ist `k` gut? Anzahl der Punkte plus eins? Warum plus eins? Und das manuelle mitzählen könnte man sich sparen, wenn man die Teilzeichenketten in einer Liste sammelt und dann deren Länge ermittelt.

Es werden für meinen Geschmack zu viele Objekte an durchnummerierte einbuchstaben Namen gebunden.

[nudelauflauf]

Erneut vielen Dank. Wie gesagt ich bin gerade am Anfang meiner Python(bzw. Programmier)-Karriere und hab leider wenig Ahnung über effizientes Arbeiten. Warum ich direkt einen String daraus machen ist da vor die Koordinaten immer noch die Belichtungszeit muss. Da im obenstehenden Versuch die mittlere Zeile in X-Richtung Doppelt und die mittlere Zeile in Y-Richtung nicht vorhanden waren musste ich ein paar Fallunterscheidungen einführen.

Den Pop Aufruf habe ich mal wie empfohlen rausgenommen.

Die Linie ist ja schon sehr kurz (nur 7 Pixel lang). Ich denke diese wird aber als Gerade interpretiert, denn es kommen wenn vorhanden immer beide Schnittpunkte heraus. Wenn sich das aber bewahrheitet kann man die Linie tatsächlich auch nur 2 Pixel lang machen.
Den 2. Punkt benutzen will ich nicht da ich diesen wieder auslesen, runden und in Integer umwandeln müsste. Da ist eine kurze Berechnung (q2 = 2*x-q) denke ich deutlich schneller.

Mit k zähle ich tatsächlich mit wieviele Punkte ich ingesamt habe. +1 weil bevor out_str in eine Datei geschrieben wird, eine weitere Zeile eingefügt wird.

Habe nun auch versucht das ganze hier mit einer Liste abzuarbeiten. Die Performance ist auch deutlich besser.

Der Code sieht nun folgendermaßen aus:

Code: Alles auswählen

def circlecreate(n, r, x, y):                                          

    from sympy.geometry import Point, Line, Circle, intersection

    fobj = open('stream.txt', 'w')
    print >>fobj, 's16'
    print >>fobj, 1
    out_lis=[]

    p1= Point(x,y)
    c = Circle(p1, r)
    q=x-r
    dur = duration (9.6) #Berechnet Belichtungszeit, noch statisch, kann aber variabel werden
    for t in range(0, r+1):
        p2, p3 = Point(q-2, y+t), Point(q+5, y+t)
        l1 = Line (p2, p3)
        pl = intersection(l1, c)
        if bool(pl)==True:
            p4 = pl[0]
            q = int(round(p4.x))
            q2 = 2*x-q
            z=0
            while q+z <= x:
                if q+z==x:
                    if t==0:
                        out_lis.extend(['%d %d %d' %(dur, q+z, y-t)])               
                        z+=1
                    elif t!=0:
                        out_lis.extend(['%d %d %d' %(dur, q+z, y-t)])
                        out_lis.extend(['%d %d %d' %(dur, q+z, y+t)])
                        z+=1
                elif q+z != x:
                    if t==0:
                        out_lis.extend(['%d %d %d' %(dur, q+z, y-t)])
                        out_lis.extend(['%d %d %d' %(dur, q2-z, y-t)])
                        z+=1

                    elif t!=0:
                        out_lis.extend(['%d %d %d' %(dur, q+z, y-t)])
                        out_lis.extend(['%d %d %d' %(dur, q+z, y+t)])
                        out_lis.extend(['%d %d %d' %(dur, q2-z, y-t)])
                        out_lis.extend(['%d %d %d' %(dur, q2-z, y+t)])
                        z+=1
    k = len(out_lis) +1
    print >>fobj, k
    print >>fobj, duration(n), x, y                               
    for daten in out_lis:
        print >>fobj, daten                                                
    fobj.close()

[Sirius3]

Hallo Nudelauflauf,

das Berechnen eines Kreises ist nun ja nicht gerade so schwierig:

Code: Alles auswählen

center = (70, 30)
radius = 50
result = []
for y in xrange(0,200):
    x = (radius**2-(y-center[1])**2)
    if x>=0:
        x = int(round(x**0.5))
        x1 = max(0, center[0]-x)
        x2 = min(200, center[0]+x)
        for x in xrange(x1,x2+1):
            result.append('%d %d'%(x,y))
print '\n'.join(result)

So sparst Du Dir auch das 4fache (9fache) Schreiben in Deine Liste.

Auf True muß man nicht explizit prüfen, es würde ein

Code: Alles auswählen

if pl:

reichen.

Grüße
Sirius

[BlackJack]

@nudelauflauf: Das ``print`` mit der ``>>``-Syntax in einer Datei umzulenken ist relativ ungewöhnlich. Diese Syntax ist AFAIK nur noch in der Sprache (in den 2er-Versionen) weil man die Rückwärtskompatibilität nicht verletzen wollte.

Importe schreibt man üblicherweise alle oben in den Quelltext und nicht in die Funktionen. Dann sieht man oben gleich wovon ein Modul abhängt und es gibt auch keine unnötigen Wiederholungen wenn Module von mehr als einer Funktion benötigt werden.

Das erhöhen von `z` passiert in *jedem* Zweig am Ende, also kann man das auch *einmal* hinter die Abfragen schreiben.

Wenn in einem ``elif`` genau das Gegenteil von der Bedingung im ``if`` geprüft wird, dann braucht man keine zweite Prüfung sondern schlicht ein ``else``.

Und wie gesagt, könnte man an der Stelle der ``while``-Schleife auch gleich eine ``for``-Schleife über die X-Koordinaten *aller* Punkte der horizontalen Linie schreiben und damit den Schleifenkörper vereinfachen.

Dateien kann man mit der ``with``-Anweisung öffnen und damit sicherstellen, dass sie in jedem Fall auch wieder geschlossen werden, wenn der Block verlassen wird.

`list.extend()` mit einer Liste mit einem Element aufzurufen ist eine umständliche Art `list.append()` auszudrücken. Du hättest da eher die Fälle mit mehreren Koordinaten in einer Liste mit *einem* `extend()`-Aufruf erledigen können.

[kbr]

Gibt es eine eventuell schnellere Lösung?

Aus Neugier habe ich mal ausprobiert, wie schnell sich die Koordinaten der Punkte auf einem Kreisoktanten mittels verschiedener Verfahren berechnen lassen:

Code: Alles auswählen

import time

def get_circleoctant_points(radius):
    # midpoint circle algorithm
    x = 0
    y = radius
    d = 1 - radius
    points = []
    while x < y:
        points.append((x, y))
        if d < 0:
            d += 2 * x + 3
            x += 1
        else:
            d += 2 * (x - y) + 5
            x += 1
            y -= 1
    return points

def get_circleoctant_points_so(radius):
    # second order differences
    x = 0
    y = radius
    d = 1 - radius
    dE = 3
    dSE = -2 * radius + 5
    points = []
    while y > x:
        points.append((x, y))
        if d < 0:
            d += dE
            dSE += 2
        else:
            d += dSE
            dSE += 4
            y -= 1
        dE += 2
        x += 1
    return points

def get_circleoctant_points_float(radius):
    # real arithmetic computation
    x = 0
    y = radius
    r = radius ** 2
    points = []
    while x < y:
        points.append((x, y))
        y2 = r - x ** 2
        y = int(round(y2 ** 0.5))
        x += 1
    return points

cl1 = time.clock()
p = get_circleoctant_points(40)
cl2 = time.clock()
p2 = get_circleoctant_points_so(40)
cl3 = time.clock()
q = get_circleoctant_points_float(40)
cl4 = time.clock()

# print p
# print p2
# print q

print "time p:  %2.6f" % (cl2 - cl1,)
print "time p2: %2.6f" % (cl3 - cl2,)
print "time q:  %2.6f" % (cl4 - cl3,)

Second order differences gewinnt. Aufgrund der Kreis-Symmetrie ist es nun leicht, die Koordinaten der Punkte innhalb des Kreises zu berechnen:

Code: Alles auswählen

def get_quadrant_row_points((x, y)):
    return [(x, y) for x in xrange(-x, x+1)]

def get_row_points((x, y)):
    pts = get_quadrant_row_points((x, y))
    pts.extend(get_quadrant_row_points((x, -y)))
    pts.extend(get_quadrant_row_points((y, x)))
    pts.extend(get_quadrant_row_points((y, -x)))
    return pts
    
def get_circle_points(octant_points):
    pts = []
    for point in octant_points:
        pts.extend(get_row_points(point))
    return pts

octant_points = get_circleoctant_points_so(40)
circle_points = get_circle_points(octant_points)
# nun ggf. den Ursprung verschieben

Dies läuft etwa so schnell wie die Lösung von pillmuncher (gleichfalls Bresenham, wie auch von lunar erwähnt). Ob sympy schneller ist, habe ich nicht ausprobiert.

[nudelauflauf]

Vielen Dank an alle. Die Lösungen waren echt gut, und auch super schnell. Da ich hinter der Lösung stehen muss und diese auch erklären können muss hab ich mich für eine Lösung ähnlich der von Sirius und gegen Bresenham entschieden =)

Bin heute morgen auf die Idee mit der Kreisgleichung gestoßen und hab dann im nachhinein festgestellt, dass diese hier ja schon als Vorschlag drin war

Mein finaler Code sieht dann so aus, falls es jemanden noch interessiert:

Code: Alles auswählen

def circlecreate(r, x, y):                                          
    out_lis3=[]
    
    dur = duration (9.6)                                              
    for y1 in range(0, r):
        x1 = int(round((r*r - y1*y1)**(0.5)))
        for q in range (-x1, x1):
            if y1==0:
                out_lis3.append('%d %d %d' %(dur, q+x, y+y1))
            else:
                out_lis3.append('%d %d %d' %(dur, q+x, y+y1))
                out_lis3.append('%d %d %d' %(dur, q+x, y-y1))
    return out_lis3

[EyDu]

Du solltest noch ein wenig zusammenfassen:

Code: Alles auswählen

FORMAT = '%d %d %d'

def circlecreate(r, x, y):                                          
    out_lis3=[]
   
    dur = duration (9.6)                                              
    for y1 in range(0, r):
        x1 = int(round((r*r - y1*y1)**(0.5)))
        for q in range (-x1, x1):
            out_lis3.append(FORMAT % (dur, q+x, y+y1))

            if y1:
                out_lis3.append(FORMAT % (dur, q+x, y-y1))

    return out_lis3

Auch sind deine Namen alle sehr seltsam. Du kürzt zum Beispiel zu viel ab (warum out_lis und nicht out_list) oder bentutz nichtssagende Namen. Es wäre interessant, was in out_lis3 drin ist und nicht, dass es sich um eine Liste handelt, die irgend ein ergebnis liefert. Und warum die 3 am Ende von out_lis3 oder ab die Einsen am Ende von x1 und y1?

Ich habe mir das Problem nicht genauer angeschaut, aber ist dir klar, dass range(m, n) von m bis n-1 reicht?

[nudelauflauf]

Die 3 ist da weil es mein dritter Versuch ist, ist aber inzwischen weg. Die 1en hinter x und y sind da, weil ich hab x und y übergebe und den Kreis damit später wenn ich die Liste erstelle ja dorthin verschiebe wo ich den Kreis gern hätte.
Deine Zusammenfassung macht nicht so viel Sinn wenn ich das richtig verstehe. Ich hab die Unterscheidung drin, damit ich wenn ich die Punkte mit -x1, x1, = r nicht doppelt aufnehme. Wenn x,x1 kleiner als r sind muss ich immer +y1 und -y1 in meine Liste aufnehmen um den ganzen Kreis zu bekommen.

[BlackJack]

(Etwas zu spät, aber ich schick's doch noch ab. )

@nudelauflauf: Die Namen sind immer noch verbesserungswürdig. `out_lis3`? Das `t` wäre zuviel gewesen, aber eine unötige `3` ging dann doch noch? `dur` würde ich auch als `duration` ausschreiben und der Funktion einen passenderen Namen geben. Funktionen werden üblicherweise mit Tätigkeiten benannt und `duration` ist Keine.

Code: Alles auswählen

from math import sqrt


def create_circle(radius, center_x, center_y):
    result = list()
    duration = calclulate_duration(9.6)
    for y in xrange(0, radius):
        width = int(round(sqrt(radius**2 - y**2)))
        for x in xrange(-width + center_x, width + center_x):
            result.append('%d %d %d' % (duration, x, center_y + y))
            if y:
                result.append('%d %d %d' % (duration, x, center_y - y))
    return result

PS: Die Zusammenfassung macht Sinn und sie macht letztlich genau das selbe wie Dein Code. 0 ist in Python als Wahheitswert „Falsch” und jede andere Zahl ist „Wahr”.

[nudelauflauf]

Ahhhh, jetzt verstehe ich. Das mit dem besser benennen macht natürlich Sinn. Zwar weiß man wahrscheinlich wenn man es selbst geschrieben hat bestens Bescheid, aber für Aussenstehende ist das mit Sicherheit besser zu durchschauen später Mal.


Thx

[JonasR]

Ahhhh, jetzt verstehe ich. Das mit dem besser benennen macht natürlich Sinn. Zwar weiß man wahrscheinlich wenn man es selbst geschrieben hat bestens Bescheid, aber für Aussenstehende ist das mit Sicherheit besser zu durchschauen später Mal.


Thx

Schau dir deinen Code in einem Jahr nochmal an