2009-12-15 18 views
5

Ich habe eine Liste von 300000 Listen (Faserbahnen), wobei jede Spur eine Liste von (x, y, z) Tupeln/Koordinaten:Effizienter Weg, Kreuzungen zu zählen?

tracks= 
[[(1,2,3),(3,2,4),...] 
[(4,2,1),(5,7,3),...] 
... 
] 

Ich habe auch eine Gruppe von Masken, wobei jede Maske definiert ist als eine Liste von (x, y, z) Tupeln/Koordinaten:

mask_coords_list= 
[[(1,2,3),(8,13,4),...] 
[(6,2,2),(5,7,3),...] 
... 
] 

ich versuche, für alle möglichen Paare von Masken zu finden, die:

  1. die Anzahl von Spuren, die jeweils überschneiden Maske-Maske-Paar (um eine Verbindung zu erstellen ectivity matrix)
  2. die Teilmenge der Spuren, die jede Maske schneiden, um 1 zu jedem (x, y, z hinzuzufügen) in der Teilmenge für jede Spur Koordinate (ein „Dichte“ Bild zu erzeugen)

ich mache zur Zeit Teil 1 wie folgt:

def mask_connectivity_matrix(tracks,masks,masks_coords_list): 
    connect_mat=zeros((len(masks),len(masks))) 
    for track in tracks: 
     cur=[] 
     for count,mask_coords in enumerate(masks_coords_list): 
      if any(set(track) & set(mask_coords)): 
       cur.append(count) 
      for x,y in list(itertools.combinations(cur,2)): 
       connect_mat[x,y] += 1 

und Teil 2 wie folgt:

def mask_tracks(tracks,masks,masks_coords_list): 
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks))) 
    for track in tracks: 
     for count,mask in enumerate(masks_coords_list): 
      if any(set(track) & set(mask)): 
       for x,y,z in track: 
        vox_tracks_img[x,y,z,count] += 1 

Sets Mit Kreuzungen finden hat diesen Prozess beschleunigt signifikant, aber beiden Teile stil Ich brauche eine Stunde, wenn ich eine Liste von 70 oder mehr Masken habe. Gibt es dafür einen effizienteren Weg als für jeden Track zu iterieren?

+0

Alle Antworten scheinen marginale Verbesserungen zu sein, aber ich denke, dass Sie mehr als das brauchen. – McPherrinM

+0

Wenn Sie einen Beispieldatensatz und die richtigen Antworten in einem Pastebin irgendwo veröffentlichen könnten, erhalten Sie möglicherweise mehr Hilfe. –

+0

Sehe ich das richtig, dass Schnittpunkte nur als zwei Koordinatentupel definiert sind, die gleich sind, und nicht als Linien zwischen den Schnittpunkten? – Svante

Antwort

3

Linearisieren Sie die Voxelkoordinaten und fügen Sie sie in zwei scipy.sparse.sparse.csc-Matrizen ein.

Sei v die Anzahl der Voxel, m die Anzahl der Masken und t die Anzahl der Spuren.
Sei M die Maske csc matrix, Größe (m x v), wobei eine 1 bei (i, j) bedeutet Maske i überlappt Voxel j.
Sei T die Spur csc-Matrix, Größe (t x v), wobei eine 1 bei (k, j) bedeutet, dass die Spur k das Voxel j überlappt.

Overlap = (M * T.transpose() > 0) # track T overlaps mask M 
Connected = (Overlap * Overlap.tranpose() > 0) # Connected masks 
Density[mask_idx] = numpy.take(T, nonzero(Overlap[mask_idx, :])[0], axis=0).sum(axis=0) 

Ich könnte auf dem letzten falsch sein, und ich bin nicht sicher, css_matrices kann durch Nicht-Null-& nehmen operiert werden. Sie müssen möglicherweise jede Spalte in einer Schleife herausziehen und in eine vollständige Matrix konvertieren.


Ich habe einige Experimente ausgeführt, um zu simulieren, was ich für eine vernünftige Menge an Daten hielt. Der Code unten dauert etwa 2 Minuten auf einem 2 Jahre alten MacBook. Wenn Sie csr_matrices verwenden, dauert es etwa 4 Minuten. Es gibt wahrscheinlich einen Kompromiss, abhängig davon, wie lang jede Spur ist.

from numpy import * 
from scipy.sparse import csc_matrix 

nvox = 1000000 
ntracks = 300000 
nmask = 100 

# create about 100 entries per track 
tcoords = random.uniform(0, ntracks, ntracks * 100).astype(int) 
vcoords = random.uniform(0, nvox, ntracks * 100).astype(int) 
d = ones(ntracks * 100) 
T = csc_matrix((d, vstack((tcoords, vcoords))), shape=(ntracks, nvox), dtype=bool) 

# create around 10000 entries per mask 
mcoords = random.uniform(0, nmask, nmask * 10000).astype(int) 
vcoords = random.uniform(0, nvox, nmask * 10000).astype(int) 
d = ones(nmask * 10000) 
M = csc_matrix((d, vstack((mcoords, vcoords))), shape=(nmask, nvox), dtype=bool) 

Overlap = (M * T.transpose()).astype(bool) # mask M overlaps track T 
Connected = (Overlap * Overlap.transpose()).astype(bool) # mask M1 and M2 are connected 
Density = Overlap * T.astype(float) # number of tracks overlapping mask M summed across voxels 
+0

Wenn der dtyp der Matrizen auf bool gesetzt ist, sind die Bits "> 0" nicht mehr nötig. –

+2

eigentlich nicht wahr. Zumindest für dünn besetzte Matrizen fördert Multiplikation sie zu einem Byte. (?) Ich hoffe, das bedeutet nicht, dass es auch Probleme gibt. –

+0

Danke dafür, beschleunigte mich auf unter einer Minute mit der durchschnittlichen Spurlänge um 10 und der durchschnittlichen Maskengröße um 500. – jbrown

0

Sie können wahrscheinlich beginnen, indem Sie die beiden Funktionen kombinieren, um beide Ergebnisse auf einmal zu erstellen. Es ist auch nicht notwendig, vor dem Looping eine Liste der Kombinationen zu erstellen, da es bereits ein Generator ist, und das könnte Ihnen Zeit sparen.

def mask_connectivity_matrix_and_tracks(tracks,masks,masks_coords_list): 
    connect_mat=zeros((len(masks),len(masks))) 
    vox_tracks_img=zeros((xdim,ydim,zdim,len(masks))) 
    for track in tracks: 
     cur=[] 
     for count,mask_coords in enumerate(masks_coords_list): 
      if any(set(track) & set(mask_coords)): 
       cur.append(count) 
       for x,y,z in track: 
        vox_tracks_img[x,y,z,count] += 1 
      for x,y in itertools.combinations(cur,2): 
       connect_mat[x,y] += 1 

Auch wird dies wahrscheinlich nie „schnell“ sein, wie in „fertig, bevor wir sterben“, so der beste Weg, schließlich ist es für Python mit Cython als c-Modul zu kompilieren.

0

Wenn Sie die einzelnen Masken Punkte gespeichert haben: (1,2,3), (1,2,4), (1,3,1) als ein Wörterbuch wie dieses: {1: [{2: set([3, 4])}, {3: set([1])}]}, könnten Sie am Ende in der Lage sein, schneller nach Übereinstimmungen zu suchen ... aber vielleicht nicht.

0

Eine geringfügige Optimierung (gleiches Big-O, sligthly kleinerer Multiplikator) kann durch Entfernen redundante Operationen gehabt werden:

  1. nicht set so viele Male auf jeder Spur und Maske nennen: nennt es einmal pro Spur und einmal pro Maske, Hilfs „parallel“ Listen von Sätzen einzurichten, dann Arbeit an den
  2. if any(someset): semantisch die gleichen wie if someset: aber etwas langsamer

Werde nicht einen dramatischen Unterschied machen, aber könnte minutiös helfen.

0

Lame noch eine weitere schrittweise Verbesserung vorzuschlagen, die gemacht werden könnte, ich weiß, aber:

Sets kleiner ganzer Zahlen können lange Ints als Bit-Vektoren mit Python modelliert werden. Angenommen, Sie ersetzen jedes Tupel durch eine kleine Integer-ID und konvertieren dann jede Spur und jede Gruppe von Maskenkoordinaten in eine Menge dieser kleinen IDs. Sie könnten diese Sätze als lange Intets darstellen, wodurch die Kreuzungsoperation etwas schneller wird (aber nicht asymptotisch schneller).

1

OK, ich denke, ich habe endlich etwas, das die Komplexität reduzieren wird. Dieser Code sollte wirklich fliegen im Vergleich zu dem, was Sie haben.

Es scheint zunächst so zu sein, dass Sie wissen müssen, welche Spuren mit welchen Masken übereinstimmen, die incidence matrix.

import numpy 
from collections import defaultdict 

def by_point(sets): 
    d = defaultdict(list) 
    for i, s in enumerate(sets): 
     for pt in s: 
      d[pt].append(i) 
    return d 

def calc(xdim, ydim, zdim, mask_coords_list, tracks): 
    masks_by_point = by_point(mask_coords_list) 
    tracks_by_point = by_point(tracks) 

    a = numpy.zeros((len(mask_coords_list), len(tracks)), dtype=int) 
    for pt, maskids in masks_by_point.iteritems(): 
     for trackid in tracks_by_point.get(pt,()): 
      a[maskids, trackid] = 1 
    m = numpy.matrix(a) 

Die adjacency matrix Sie suchen ist m * m.T.

Der Code, den Sie bisher haben, berechnet nur das obere Dreieck. Sie können triu verwenden, um gerade diese Hälfte zu greifen.

Die Voxelberechnung kann auch die Inzidenzmatrix verwenden.

vox_tracks_img = numpy.zeros((xdim, ydim, zdim, len(mask_coords_list)), dtype=int) 
    for trackid, track in enumerate(tracks): 
     for x, y, z in track: 
      vox_tracks_img[x, y, z, :] += a[:,trackid] 
    return am, vox_tracks_img 

Für mich läuft das in weniger als einer Sekunde für Datensätze mit Hunderten von Masken und Spuren.

Wenn Sie viele Punkte in Masken haben, aber keine Spuren haben, kann es sinnvoll sein, die Einträge für diese Punkte vor dem Eintritt in die Schleife aus masks_by_point zu löschen.

Verwandte Themen