# -*- coding: utf-8 -*-
"""
Created on Thu Jan 7 15:43:07 2016
Upgraded on Wed Sep 14 10:55:30 2016 - Version 2
- Improved algorithm: the first version can fail for some specific shapes of
clusters (worm-like shapes). For example with
Pos = np.array([
[0,0,0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,4,4,4,4,4,4,4,5,6,6,6,6,6,
[4,8,9,10,11,12,4,8,12,4,8,12,4,8,12, 0,1,2,3,4,5,6,7,8,12,12,8,9,10,11,12]
], dtype = np.int16).
- Bug fixed at the lines 78-81 of version 1 : the pixel identifier can go up
to dimx2*dimy2 (requiring more than 16 bits).
Demo
@author: lautrido
"""
import numpy as np
import time
###############################################################################
def findclu(dimx, dimy, Pos):
"""
FINDCLU finds the solitary and clustered pixels in a 2D array of pixels.
Recognition of these structures is commonly used in processing of ionizing
particles in pixelated detectors or of CCDs in astrophysics.
---PURPOSE---
The function provides the ordered lists of solitary and clustered pixels.
A Cluster is defined by the assembly of adjoining pixels (neighbors).
An isolated pixel is a pixel having no neighbors. Neighbors of the pixel
P[i,j] are defined as the 8 pixels of the sub-array P[i-1:i+1, j-1:j+1].
The function can rapidly process large 2D arrays having a large number of
solitaries and large sizes of clusters.
---INPUT---
dimx = number of rows of the array of pixels (int)
dimy = number of columns of the array of pixels (int)
Pos = numpy 2D array (of int16) of the indices of positions of the selected
pixels, ranked following [[list of row indexes], [list of column
indexes]]
(see demo script below)
---OUTPUT---
Sol = numpy 2D array (of int16) of the ordered indices of positions of the
solitary pixels, ranked following [[list of row indexes], [list of
column indexes]]
Clu = numpy 2D array (of int16) of the ordered indices of positions of the
clustered pixels, ranked following [[list of row indexes], [list of
column indexes]]. The cluster separator is a negative number on the
position x (first row of the output array Clu) and a no meaning
value at the corresponding position in y (second row of Clu).
(see demo script below)
---USAGE---
Sigl, Neig = findclu(dimx, dimy, Pos)
---AUTHOR---
P. Lautridou (2016) (Pascal.Lautridou@subatech.in2p3.fr)
Tested with python 3.5.1, numpy 1.11.1 - Matlab version also available
"""
import numpy as np
# print("Number of input pixels:",np.shape(Pos)[1])
dimx2, dimy2 = dimx+2, dimy+2
H = np.full((dimx2, dimy2), False, dtype=bool)
Pos = Pos + 1
H[Pos[0, :], Pos[1, :]] = True
# Search of the neighbors (cf. Scipy central: 100 numpy exercices)
V = (H[:-2, :-2] | H[:-2, 1:-1] | H[:-2, 2:] | H[1:-1, :-2] |
H[1:-1, 2:] | H[2:, :-2] | H[2:, 1:-1] | H[2:, 2:])
# Search of the solitaries
V = np.logical_and(np.logical_not(V), H[1:-1, 1:-1])
# 2D-array Sol = Coodinates [[raw],[col]] of the solitary pixels
Sol = np.int16((np.nonzero(V)))
print("Number of solitary pixels:", np.shape(Sol)[1])
# Search of the clusters
H[1:-1, 1:-1] = np.logical_xor(H[1:-1, 1:-1], V)
Hf = np.int32(np.flatnonzero(H))
nhit = np.shape(Hf)[0]
print("Number of clustered pixels:", nhit)
if nhit<2:
Clu = np.array([], np.int16)
return Sol, Clu
V = np.zeros((dimx2, dimy2), np.int32)
V[Hf // dimy2, Hf % dimy2] = Hf
for k in range(nhit):
i = (Hf[k] // dimy2) - 1
j = (Hf[k] % dimy2) - 1
V[i:i+3, j:j+3] = H[i:i+3, j:j+3] * V[i+1,j+1]
for k in range(nhit):
i = (Hf[k] // dimy2) - 1
j = (Hf[k] % dimy2) - 1
for m in range(i, i+3):
for n in range(j, j+3):
if (V[m,n] > 0) and (V[m,n] != V[i+1, j+1]):
B=np.int32(np.where(V==V[m,n]))
V[B[0],B[1]] = V[i+1, j+1]
Vf = V[V>0]
U, Nb = np.unique(Vf, return_counts=True)
Clu = np.array([], np.int32)
for i in U:
Clu=np.append(Clu, Hf[Vf==i])
Clu=np.append(Clu, [-1])
nclu=np.size(U)
print("Number of clusters:", nclu)
if nclu:
print("Largest cluster size:", np.amax(Nb), "pixels")
# 2D array Clu = Coodinates [[row],[col]] of the clustered pixels
Clu = np.int16([Clu // dimy2, Clu % dimy2])-1
return Sol, Clu
###############################################################################
# Demo with a 2D array of dimension (dimx, dimy),
# and the following coordinates (row,col) of the selected pixels:
# (0,0); (0,2); (1,4); (1,1); (2,0); (4,0); (4,2); (4,3); (4,4)
dimx, dimy = 5, 6
# coordinates array in format [[row],[col]]
Pos = np.array([[0,0,1,1,2,4,4,4,4], [0,2,4,1,0,0,2,3,4]], dtype = np.int16)
#Pos = np.array([[0,0,0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,4,4,4,4,4,4,4,5,6,6,6,6,6],
# [4,8,9,10,11,12,4,8,12,4,8,12,4,8,12,
# 0,1,2,3,4,5,6,7,8,12,12,8,9,10,11,12]], dtype = np.int16)
print(Pos)
#dimx, dimy = 8192, 8192
#Pos = np.int16([np.arange(dimx),np.arange(dimx)])
t1 = time.clock()
Sigl, Neig = findclu(dimx, dimy, Pos)
t2 = time.clock()
print("Runing time:", t2-t1)
print("Solitaries", Sigl)
print("Clusters", Neig)