Python Vektorisierung für Schleifen verschachtelt

Question

ich etwas Hilfe bei der Suche nach und das Verständnis eine pythonic Art und Weise schätzen würde die folgenden Array Manipulationen in verschachtelten for-Schleifen zu optimieren:

def _func(a, b, radius): 
    "Return 0 if a>b, otherwise return 1" 
    if distance.euclidean(a, b) < radius: 
     return 1 
    else: 
     return 0 

def _make_mask(volume, roi, radius): 
    mask = numpy.zeros(volume.shape) 
    for x in range(volume.shape[0]): 
     for y in range(volume.shape[1]): 
      for z in range(volume.shape[2]): 
       mask[x, y, z] = _func((x, y, z), roi, radius) 
    return mask 

Wo volume.shape (182, 218, 200) und roi.shape (3,) sind beide ndarray Arten; und radius ist ein int

Answer 1

Ansatz # 1

Hier ist ein vektorisiert Ansatz -

m,n,r = volume.shape 
x,y,z = np.mgrid[0:m,0:n,0:r] 
X = x - roi[0] 
Y = y - roi[1] 
Z = z - roi[2] 
mask = X**2 + Y**2 + Z**2 < radius**2 

Mögliche Verbesserung: Wir können wahrscheinlich den letzten Schritt mit numexpr Modul Speedup -

import numexpr as ne 

mask = ne.evaluate('X**2 + Y**2 + Z**2 < radius**2') 

Ansatz # 2

Wir können auch schrittweise die drei Bereiche entsprechend den Formparametern erstellen und die Subtraktion gegen die drei Elemente von roi im laufenden Betrieb durchführen, ohne tatsächlich die Gitter zu erzeugen, wie früher mit np.mgrid gemacht. Dies würde durch die Verwendung von broadcasting aus Effizienzgründen begünstigt werden. Die Umsetzung würde wie folgt aussehen -

m,n,r = volume.shape 
vals = ((np.arange(m)-roi[0])**2)[:,None,None] + \ 
     ((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2) 
mask = vals < radius**2 

Vereinfachte Version: Dank @Bi Rico für die Annahme, hier eine Verbesserung, wie wir np.ogrid können diese Operationen in etwas prägnanter Weise durchzuführen, wie so -

m,n,r = volume.shape  
x,y,z = np.ogrid[0:m,0:n,0:r]-roi 
mask = (x**2+y**2+z**2) < radius**2 

---

Runtime Test

Funktionsdefinitionen -

def vectorized_app1(volume, roi, radius): 
    m,n,r = volume.shape 
    x,y,z = np.mgrid[0:m,0:n,0:r] 
    X = x - roi[0] 
    Y = y - roi[1] 
    Z = z - roi[2] 
    return X**2 + Y**2 + Z**2 < radius**2 

def vectorized_app1_improved(volume, roi, radius): 
    m,n,r = volume.shape 
    x,y,z = np.mgrid[0:m,0:n,0:r] 
    X = x - roi[0] 
    Y = y - roi[1] 
    Z = z - roi[2] 
    return ne.evaluate('X**2 + Y**2 + Z**2 < radius**2') 

def vectorized_app2(volume, roi, radius): 
    m,n,r = volume.shape 
    vals = ((np.arange(m)-roi[0])**2)[:,None,None] + \ 
      ((np.arange(n)-roi[1])**2)[:,None] + ((np.arange(r)-roi[2])**2) 
    return vals < radius**2 

def vectorized_app2_simplified(volume, roi, radius): 
    m,n,r = volume.shape  
    x,y,z = np.ogrid[0:m,0:n,0:r]-roi 
    return (x**2+y**2+z**2) < radius**2 

Timings -

In [106]: # Setup input arrays 
    ...: volume = np.random.rand(90,110,100) # Half of original input sizes 
    ...: roi = np.random.rand(3) 
    ...: radius = 3.4 
    ...: 

In [107]: %timeit _make_mask(volume, roi, radius) 
1 loops, best of 3: 41.4 s per loop 

In [108]: %timeit vectorized_app1(volume, roi, radius) 
10 loops, best of 3: 62.3 ms per loop 

In [109]: %timeit vectorized_app1_improved(volume, roi, radius) 
10 loops, best of 3: 47 ms per loop 

In [110]: %timeit vectorized_app2(volume, roi, radius) 
100 loops, best of 3: 4.26 ms per loop 

In [139]: %timeit vectorized_app2_simplified(volume, roi, radius) 
100 loops, best of 3: 4.36 ms per loop 

So, wie immer broadcasting seine Magie für ein verrücktes zeigt fast 10,000x Speedup gegenüber dem ursprünglichen Code und mehr als 10x besser als Maschen zu schaffen, indem Sie on-the -fly Sendung gesendet!

Answer 2

Sagen Sie bitte zuerst eine xyzy Array bauen:

import itertools 

xyz = [np.array(p) for p in itertools.product(range(volume.shape[0]), range(volume.shape[1]), range(volume.shape[2]))] 

nun mit numpy.linalg.norm,

np.linalg.norm(xyz - roi, axis=1) < radius 

überprüft, ob der Abstand für jedes Tupel von roi kleiner als Radius.

Schließlich, nur reshape das Ergebnis zu den Dimensionen, die Sie benötigen.