Kugelkoordinaten mit gleichem Radius herausfinden

[g123]

Hallo Leute,

ich habe mit einem Numpy Array ein Voxel Gitter erstellt, in meinem array bedeutet Array[100][100][40] = 20: Der Voxel an der Stelle x=100 y=100 z=40 besitzt den Wert 20.

So nun möchte ich Messdaten in diesem array abspeichern, ich bekomme Messdate in dieser Form: x=10 y=10 z=10 r=3 s= 77, d.h. alle Voxel auf mit einem Radius 3 auf der Halbkugel um den Punkt (x,y,z) sollen den Wert 77 addiert bekommen. Jedoch kenn ich keine wirklich geeignete Funktion, mein Lösungsansatz wäre zurzeit:

Code: Alles auswählen

def Koordinaten(x, y, z, r, s):
    a = (100, 100, 100)
    koord = np.zeros(a)

    phi = 0
    eta = 0
    while phi <= 2pi:
        while eta <= pi:
            x = r * sin eta * cos phi
            y = r * sin eta * sin phi
            z = t * cos eta
            a[x][y][z] = s
            eta + = 0.1
    phi += 0.1

Finde es halt komisch dass ich per hand über eta und phi iteriere, gibt schnellere bzw. 'bessere' Wege ?

[ThomasL]

in welchem Wertebereich würde sich den die Radius Variable r bewegen?
Dein Code oben hat zwei Fehler, z = t * soll doch wohl z = r* sein
und die letzte Zeile phi*= muss noch 1x eingerückt werden

[ThomasL]

schau dir mal diesen Code an, mit
idx = adjust_center(index[radius], (x,y,z))
matrix[idx] += density
kannst du um Koordinate (x,y,z) im radius alle Werte um density verändern

Code: Alles auswählen

import numpy as np
import math


def create_index(max_radius):
    array = np.zeros((50, 50, 50))
    index_matrix = []

    for radius in range(max_radius):
        phi = 0.0
        while phi <= math.pi * 2:
            eta = 0.0
            while eta <= math.pi:
                x = round(25 + radius * math.sin(eta) * math.cos(phi))
                y = round(25 + radius * math.sin(eta) * math.sin(phi))
                z = round(25 + radius * math.cos(eta))
                array[x,y,z] = 1.
                eta += 0.01
            phi += 0.01
        where = np.where(array == 1.)

        for i in range(3):
            arr = where[i]
            arr -= 25

        index_matrix.append(where)
        
    return index_matrix


def adjust_center(matrix, center):
    x_coord = matrix[0] + center[0]
    y_coord = matrix[1] + center[1]
    z_coord = matrix[2] + center[2]
    return tuple([x_coord, y_coord, z_coord])

index = create_index(5)
matrix = np.zeros((9,9,9))

x = 4
y = 4
z = 4
radius = 3
density = 10

idx = adjust_center(index[radius], (x,y,z))
matrix[idx] += density

for i in range(9):
    for j in range(9):
        print(matrix[i,j])

[Sirius3]

Zuerst erzeugst Du Dir Arrays mit den Koordinaten.
Dann berechnest Du Dir eine Maske, welche Koordinaten innerhalb des Kugelradius liegen.
Zum Schluß änderst Du die Werte, die durch die Maske ausgewählt wurden:

Code: Alles auswählen

a = numpy.zeros((100, 100, 100))
sz, sy, sx = (numpy.arange(s) for s in a.shape)

x, y, z, r, s = 10, 10, 10, 3, 77
mask = (sx[None,None,:]-x)**2 + (sy[None,:,None]-y)**2 + (sz[:, None, None]-z)**2 <= r
a[mask] += s

[ThomasL]

Ich hatte mir schon gedacht, dass es eine reine Numpy Lösung gibt,
die Formel für die Maske hatte ich auch schon im Blick.
Müsste es aber nicht <= r**2 sein?

[Sirius3]

Stimmt, r**2

[g123]

Vielen Dank! @Sirius3 . Wie könnte ich denn hierbei den Polarwinkel einschränken, so dass er z.B. von 2/3pi - pi geht ?

[Sirius3]

@g123: Indem Du x und y in Winkel (atan2) umrechnest und einen ähnlichen Vergleich wie mit r machst.

[g123]

@Sirius3
atan2 ? Also meine Idee wäre:

mask = ((sz[None, None, :] - z) ** 2 + (sy[None, :, None] - y) ** 2 + (sx[:, None, None] - x) ** 2 == r ** 2) and (2/3 pi <= cos(sz[None, None, :] - z) <= pi)



Das Bild hat mir so ein bisschen als Anhaltspunkt gedient, ich hab versucht das eta von 2/3pi bis pi geht, so dass ich so einen 'Strahl' nach unten habe.
Und wie würde ich das machen, falls ein Indizee eine andere Auflösung hat.
Also ein Beispiel:

a=np.zeros((10, 10, 20))
Nehmen wir an jede Achse is 10 mm lang, dann wäre a[1][2][3], der Punkt bei x = 1mm y =2mm z = 1,5mm.

Müsste ich dann so etwas machen ?
mask = (((sz[None, None, :] - z * 2) * 2) ** 2 + (sy[None, :, None] - y) ** 2 + ((sx[:, None, None] - x)/voxelz) ** 2 == r ** 2) and (2/3 pi <= cos(sz[None, None, :] - z * 2) <= pi)

[g123]

Ich meinte arccos und dann noch durch r *

Und * 3 und nicht geteilt durch 3...

Diese mask Funktion rundet dann auch immer automatisch auf oder ab auf eine glatte Zahl?

[Sirius3]

Nein, es wird nicht gerunden, sondern nur getestet, ob ein Punkt innerhalb oder außerhalb Deines Kugelsegments liegt.

[g123]

Aber was ist denn wenn ich ie Punkte haben möchte, die am nächsten zur Kugeloberfläche liegen, und nicht alle die innerhalb liegen?

[ThomasL]

In der numpy Maske von Sirius weiter oben wird doch auf <= r**2 getest,
Grenze einfach diesen Bereich ein >=r**2-x and <=r**2
Für x musst du ermitteln welche Werte da Sinn machen.