Ich habe einen 2D-Datensatz mit shape = (500, 500). Von einem bestimmten Ort (x_0, y_0) möchte ich die Entfernung für jedes Element / Pixel diesem bestimmten Ort zuordnen. Dazu bestimme ich alle eindeutigen Entfernungen von (x_0, y_0) und ordne sie mit ganzen Zahlen zu. Eine solche Karte für einen 6 x 6 Datensatz sieht folgendermaßen aus:

[9 8 7 6 7 8]
[8 5 4 3 4 5]
[7 4 2 1 2 4]
[6 3 1 0 1 3]
[7 4 2 1 2 4]
[8 5 4 3 4 5]

Dabei entsprechen die Ganzzahlen den eindeutigen Abständen, die im folgenden Array gespeichert sind:

[0.  1.  1.41421356  2.  2.23606798  2.82842712  3.  3.16227766  3.60555128  4.24264069]

Der Code, der diese Entfernungen bestimmt, lautet wie folgt:

def func(data, (x_0,y_0)):
  y, x = numpy.indices((data.shape))
  r = numpy.sqrt((x - x_0)**2 + (y - y_0)**2)

  float_values = numpy.unique(r.ravel())  # Unique already sorts the result 
  int_values = numpy.arange(float_values.shape[0]).astype(numpy.int) 

  for idx in range(float_values.shape[0])[::-1]:
    r[r == float_values[idx]] = int_values[idx] 

  return float_values, r

Die for Schleife ist ein Engpass. Die Bewerbung dauert zu lange. Gibt es eine Möglichkeit, die Leistung zu beschleunigen / zu steigern? Oder gibt es vielleicht eine völlig andere, aber schnellere Methode, um die von mir benötigte Ausgabe zu erhalten?

1
The Dude 17 Apr. 2018 im 22:59

3 Antworten

Beste Antwort

Hier ist ein vektorisierter Ansatz, der masking verwendet -

def func_mask_vectorized(data, (x_0, y_0)):
    # Leverage broadcasting with open meshes to create the squared distances/ids
    m,n = data.shape
    Y,X = np.ogrid[:m,:n]
    ids = (X-x_0)**2 + (Y-y_0)**2

    # Setup mask that will help us retrieve the unique "compressed" IDs
    # (similar to what return_inverse does).
    # This is done by setting 1s at ids places and then using that mask to 
    # assign range covered array, in effect setting up the unique compress. IDs.
    mask = np.zeros(ids.max()+1, dtype=bool)
    mask[ids] = 1    
    id_arr = mask.astype(int)
    id_arr[mask] = np.arange(mask.sum())
    r_out = id_arr[ids]

    # Finally extract out the unique ones among the IDs & get their sqrt values
    float_values_out = np.sqrt(np.flatnonzero(mask))
    return float_values_out, r_out

Benchmarking

Timing des vorgeschlagenen Setups mit Daten der Form (500,500) unter Verwendung des Zahlenbereichs 0-9, wie er auch im Beispiel aus der Frage verwendet wird, und Timing aller vollständigen Lösungen in diesem Abschnitt unten -

In [371]: np.random.seed(0)
     ...: data = np.random.randint(0,10,(500,500))
     ...: x_0 = 2
     ...: y_0 = 3

# Original soln
In [372]: %timeit func(data, (x_0,y_0))
1 loop, best of 3: 6.77 s per loop

# @Daniel's soln
In [373]: %timeit func_return_inverse(data, (x_0,y_0))
10 loops, best of 3: 23.9 ms per loop

# Soln from this post
In [374]: %timeit func_mask_vectorized(data, (x_0,y_0))
100 loops, best of 3: 5.02 ms per loop

Die Erweiterung für Fälle, in denen sich die Zahlen auf 100 oder sogar 1000 erstrecken können, ändert nicht viel daran, wie sich diese stapeln -

In [397]: np.random.seed(0)
     ...: data = np.random.randint(0,100,(500,500))
     ...: x_0 = 50
     ...: y_0 = 50

In [398]: %timeit func(data, (x_0,y_0))
     ...: %timeit func_return_inverse(data, (x_0,y_0))
     ...: %timeit func_mask_vectorized(data, (x_0,y_0))
1 loop, best of 3: 5.62 s per loop
10 loops, best of 3: 20.7 ms per loop
100 loops, best of 3: 4.28 ms per loop

In [399]: np.random.seed(0)
     ...: data = np.random.randint(0,1000,(500,500))
     ...: x_0 = 500
     ...: y_0 = 500

In [400]: %timeit func(data, (x_0,y_0))
     ...: %timeit func_return_inverse(data, (x_0,y_0))
     ...: %timeit func_mask_vectorized(data, (x_0,y_0))
1 loop, best of 3: 6.87 s per loop
10 loops, best of 3: 21.9 ms per loop
100 loops, best of 3: 5.05 ms per loop
2
Divakar 17 Apr. 2018 im 21:11

Verwenden Sie den Parameter return_inverse - von unique:

def func(data, (x_0,y_0)):
    y, x = numpy.indices(data.shape)
    r = (x - x_0)**2 + (y - y_0)**2
    float_values, r = numpy.unique(r, return_inverse=True)
    return float_values ** 0.5, r.reshape(data.shape)
1
Daniel 17 Apr. 2018 im 20:37

Es scheint, dass Ihr Indizierungsschema (Ganzzahlen in den Daten) in derselben Reihenfolge wie die Abstände liegt. Wenn dies immer der Fall ist, kann das Array von Entfernungen ohne den tatsächlichen Inhalt der Daten erzeugt werden.

Ich stütze diese Lösung auf eine Indexberechnung, bei der die x- und y-Pixel-Offsets jeder Position zur Ankerposition verwendet werden. Angenommen, "so" ist der kleinste Versatz und "ho" ist der größere Versatz, wobei "mo" der maximal mögliche Versatz in beide Richtungen ist:

Index = ho + (mo + 1) * lo - lo * (lo + 1) // 2

Um die Abstände im Array zu berechnen, müssen wir nur die Abmessungen der Matrix und die Position des Ankerpixels kennen.

import numpy as np
def distanceArray(x,y,cols,rows):
    maxDx  = max(x,cols-x)
    maxDy  = max(y,rows-y)
    maxD   = max(maxDx,maxDy)
    minD   = min(maxDx,maxDy)
    lo = np.arange(minD)[:,None]
    hi = np.arange(maxD)
    sqs = lo*lo + hi*hi
    unique = np.tri(*sqs.shape,maxD-minD, dtype=bool)[::-1,::-1]
    return np.sqrt(sqs[unique])

Wenn wir uns nur auf Pixelversätze von der Ankerposition konzentrieren, erhalten wir einen Bereich horizontaler und vertikaler Detlas, der durch die Grenzen der Datenform (maxDx und maxDy) bestimmt wird.

Für die Entfernungsberechnungen können wir die vertikalen / horizontalen Richtungen ignorieren und einen kleinen und einen großen (r) Bereich erstellen. (lo und hi von maxD und minD)

Um die gesamte Summe der Quadrate zu berechnen, können wir einen der beiden Bereiche in einen vertikalen Vektor (lo) transponieren und ihn dann nach Quadrieren ihrer Werte (hi * hi + lo * lo) zum anderen (hi) addieren. Dies ergibt eine 2D-Matrix mit allen Kombinationen von Quadratsummen (sqs).

Aus dieser Matrix ist das obere Dreieck eine Verdoppelung seines Gegenstücks. Also maskieren wir die doppelten Abstandspaare mit einer dreieckigen booleschen Matrix. (eindeutig) Durch Maskieren des oberen Dreiecks wird sichergestellt, dass die Reihenfolge der Quadratsummen, die aus der Maskierungsoperation hervorgehen, in der entsprechenden Reihenfolge liegt.

Schließlich enthalten die gefilterten sqs-Werte genau das, was wir benötigen, und zwar in der richtigen Reihenfolge. Wir können die kostspielige Quadratwurzelfunktion nur auf dieses Endergebnis anwenden.

Wenn Sie die Entfernungsberechnung nicht auf jedes Pixel anwenden, sollten Sie einige signifikante Leistungssteigerungen erzielen, da Sie die indizierten Entfernungen nur bei Bedarf verwenden können. Ich denke, es wäre unfair, die Leistung dieser distanceArray-Funktion mit den anderen Lösungen zu vergleichen (da sie nur einen Teil ihrer Arbeit leistet), aber da es auch Teil der Optimierung ist, etwas nicht tun zu müssen, wird das Endergebnis wahrscheinlich sein besser (um einen Faktor von ungefähr dem Fünffachen von Divakar in meinen nicht-wissenschaftlichen Tests).

Beachten Sie, dass Sie, wenn Sie die Abstände nur für eine kleine Teilmenge der Pixel verwenden, möglicherweise alle diese Berechnungen vermeiden und ein Wörterbuch als Cache verwenden möchten, um die Abstände "on demand" basierend auf den dX- und dY-Offsets (eingegeben auf und) zu berechnen bestelltes Tupel). Dies führt die absolute Mindestanzahl von Berechnungen durch und berechnet den Abstand nur einmal für ein bestimmtes Versatzpaar. Sie können diesen Cache sogar weiterhin für andere Ankerpositionen und Datenformen verwenden, da die Versatzpaare unabhängig von der Position des Ankers immer den gleichen Abstand erzeugen.

[BEARBEITEN] Um dieselbe Indizierung zu erhalten, die ich für distanceArray verwendet habe, können Sie Folgendes verwenden:

def offsets(x,y,cols,rows):
    mo   = max(x,cols-x-1,y,rows-y-1)+1

    dx   = abs(np.arange(cols)-x)
    dy   = abs(np.arange(rows)-y)[:,None]

    mo21 = 2 * mo - 1
    ly = dy*(mo21 - dy )//2  # mo*lo - lo*(lo+1)//2 when dy is lowest
    lx = dx*(mo21 - dx )//2  # mo*lo - lo*(lo+1)//2 when dx is lowest

    return np.maximum(dx,dy) + np.minimum(lx,ly)

offsets(3,3,6,6)

array([[9, 8, 6, 3, 6, 8],
       [8, 7, 5, 2, 5, 7],
       [6, 5, 4, 1, 4, 5],
       [3, 2, 1, 0, 1, 2],
       [6, 5, 4, 1, 4, 5],
       [8, 7, 5, 2, 5, 7]])
0
Alain T. 18 Apr. 2018 im 19:04