I am trying to identify connected regions of pixels in an image stack. Since it is a stack, the input is quite large (on the order of 10 million pixels, although only about 1 million are bright), and looping through it pixel by pixel is extremely slow. Nonetheless, my first attempt works as follows:
1) Check if current pixel is bright, if so go to 2. If not go to 4.
2) Obtain slice of image containing all 26 neighbors. Set all neighbors that have not yet been scanned to be dark. Set all neighbors that fall outside of the image (e.g. indices of -1) to be dark. Set current pixel to dark.
3) If all pixels in slice are dark, assign current pixel a new label. If exactly one pixel in slice is bright, set current pixel label to be equal to this bright pixel's label. If more than one pixel is bright, set current pixel to lowest of bright pixel labels, and record the equivalence of these labels.
4) Advance one row. Once all rows in the column are scanned, advance one column. Once every pixel is in the slice is scanned, advance one slice in the stack. Once the entire image is scanned, advance to 5.
5) Loop through the list of equivalencies, reassigning labels in the output.
I don't think the processing in steps 2 and 3 are much slower than other, similar algorithms, but I could be wrong. I think the major bottleneck is simply looping through each pixel. I know the bottleneck is not #5 because I have never been patient enough to reach that step. Even just looping over the image indices with a print statement is quite slow. I know that with the same data I can do this analysis quickly, so I suspect there must be a smarter way to do this with a very different approach. I have read something vague about performing dilation and using intersections, but I can't imagine what the mechanism they were alluding to was.
Is there something fast and clever out there, or is this basically the only way to do it? How is the MATLAB algorithm so much faster?
for i in xrange(0,nPix):
for j in xrange(0,nPix):
for k in xrange(0,nFrames):
if img[i,j,k] != 0: #If the current pixel is bright
#Construct an array of already scanned
#neighbors. Pixels are scanned first in z, then y,
#then x.
#Construct indices allowing for periodic boundary
ir = np.array(range(i-1,i+2)).reshape(-1,1,1)%nPix
jr = np.array(range(j-1,j+2)).reshape(1,-1,1)%nPix
kr = np.array(range(k-1,k+2)).reshape(1,1,-1)%nFrames
pastNeighbors = img[ir,jr,kr] #Includes i-1:i+1
#Set all indices outside image and "future neighbors" to 0
pastNeighbors[2,:,:] = 0
pastNeighbors[1,2,:] = 0
pastNeighbors[1,1,2] = 0
pastNeighbors[1,1,1] = 0 #Pixel itself; not a neighbor
#Eliminate periodic boundary
if i-1 < 0: pastNeighbors[0,:,:] = 0
if i+1 == nPix: pastNeighbors[2,:,:] = 0
if j-1 < 0: pastNeighbors[:,0,:] = 0
if j+1 == nPix: pastNeighbors[:,2,:] = 0
if k-1 < 0: pastNeighbors[:,:,0] = 0
if k+1 == nFrames: pastNeighbors[:,:,2] = 0
if np.count_nonzero(pastNeighbors) == 0: #If all dark neighbors
label = label + 1 #Assign new label to pixel
out[i,j,k] = label
elif np.count_nonzero(pastNeighbors) == 1: #One bright neighbor
relX, relY, relZ = np.nonzero(pastNeighbors)
relX = relX - 1
relY = relY - 1
relZ = relZ - 1
out[i,j,k] = out[i + relX, j + relY, k + relZ]
else: #Mutliple bright neighbors
relX, relY, relZ = np.nonzero(pastNeighbors)
relX = relX - 1
relY = relY - 1
relZ = relZ - 1
equiv = out[i + relX, j + relY, k + relZ]
out[i, j, k] = np.min(equiv) #Assign one of the
#multiple labels
for i in xrange(0,equiv.size):
equivList = np.append(equivList, [[np.min(equiv),equiv[i]]], axis = 0)
