This page is about n-d array processing algorithms coded in ImgAlgos.PyAlgos. Algorithms are implemented on C++ and also have python callable interface using boost/python wrapper. All examples are shown for python level interface.

Common hood

Make object and set parameters

import psana
from ImgAlgos.PyAlgos import PyAlgos 

# create object:
alg = PyAlgos(windows=winds, mask=mask, pbits=0)
# where pbits - is a print info control bit-word:
# pbits = 0   - print nothing
#       + 1   - main results, list of peaks
#       + 2   - input parameters, index matrix of pixels for S/N algorithm
#       + 128 - tracking and all details in class PyAlgos.py
#       + 256 - tracking and all details in class AlgArrProc
#       + 512 - tracking and all details in class AlgImgProc

# set peak-selector parameters:
alg.set_peak_selection_pars(npix_min=5, npix_max=5000, amax_thr=0, atot_thr=0, son_min=10)

Define ROI using windows and/or mask

Region Of Interest (ROI)is defined by the set of rectangular windows on segments and mask, as shown in example below.

# List of windows
winds = None # entire size of all segments will be used for peak finding
winds = (( 0, 0, 185, 0, 388),              
         ( 1, 20,160, 30,300),
         ( 7, 0, 185, 0, 388))

# Mask
mask = None                   # (default) all pixels in windows will be used for peak finding
mask = det.mask()             # see class Detector.PyDetector
mask = np.loadtxt(fname_mask) # 
mask.shape = <should be the same as shape of data n-d array>

Evaluation of the background level, rms, and S/N ratio

When peak is found, its parameters can be precised for background level, noise rms, and signal over background ratio (S/N) can be estimated. All these values can be evaluated using pixels surrounding the peak on some distance. For all peak-finders we use the same algorithm. Surrounding pixels are defined by the ring with internal radial parameter r0 and ring width dr (both in pixels). The number of surrounding pixels depends on r0 and dr parameters as shown in matrices below. We use notation

+ central pixel with maximal intensity,
1 pixels counted in calculation of averaged background level and noise rms,
0 pixels not counted.

Matrices of pixels for r0=3 and 4 and different dr values

r0=3  dr=0.1 (4 pixels)   r0=3  dr=0.5 (12 pixels)   r0=3  dr=1 (24 pixels) 
0 0 0 0 0 0 0 0 0         0 0 0 0 0 0 0 0 0          0 0 0 0 1 0 0 0 0       
0 0 0 0 1 0 0 0 0         0 0 0 1 1 1 0 0 0          0 0 1 1 1 1 1 0 0       
0 0 0 0 0 0 0 0 0         0 0 0 0 0 0 0 0 0          0 1 0 0 0 0 0 1 0       
0 0 0 0 0 0 0 0 0         0 1 0 0 0 0 0 1 0          0 1 0 0 0 0 0 1 0       
0 1 0 0 + 0 0 1 0         0 1 0 0 + 0 0 1 0          1 1 0 0 + 0 0 1 1       
0 0 0 0 0 0 0 0 0         0 1 0 0 0 0 0 1 0          0 1 0 0 0 0 0 1 0       
0 0 0 0 0 0 0 0 0         0 0 0 0 0 0 0 0 0          0 1 0 0 0 0 0 1 0       
0 0 0 0 1 0 0 0 0         0 0 0 1 1 1 0 0 0          0 0 1 1 1 1 1 0 0       
0 0 0 0 0 0 0 0 0         0 0 0 0 0 0 0 0 0          0 0 0 0 1 0 0 0 0

r0=4  dr=0.2 (12 pixels)   r0=4  dr=0.3 (16 pixels)    r0=4  dr=0.5 (24 pixels)
0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 0 0 0 0      0 0 0 0 1 1 1 0 0 0 0       0 0 0 1 1 1 1 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0      0 0 1 0 0 0 0 0 1 0 0       0 0 1 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0       0 1 0 0 0 0 0 0 0 1 0
0 1 0 0 0 0 0 0 0 1 0      0 1 0 0 0 0 0 0 0 1 0       0 1 0 0 0 0 0 0 0 1 0
0 1 0 0 0 + 0 0 0 1 0      0 1 0 0 0 + 0 0 0 1 0       0 1 0 0 0 + 0 0 0 1 0
0 1 0 0 0 0 0 0 0 1 0      0 1 0 0 0 0 0 0 0 1 0       0 1 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0       0 1 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0      0 0 1 0 0 0 0 0 1 0 0       0 0 1 0 0 0 0 0 1 0 0
0 0 0 0 1 1 1 0 0 0 0      0 0 0 0 1 1 1 0 0 0 0       0 0 0 1 1 1 1 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0

Matrices of pixels for r0=5 and 6 and different dr values

r0=5  dr=0.05 (12 pixels)      r0=5  dr=0.5  (28 pixels)            
0 0 0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0      0 0 0 0 1 1 1 1 1 0 0 0 0
0 0 0 1 0 0 0 0 0 1 0 0 0      0 0 0 1 0 0 0 0 0 1 0 0 0
0 0 1 0 0 0 0 0 0 0 1 0 0      0 0 1 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0 0 0 1 0
0 1 0 0 0 0 + 0 0 0 0 1 0      0 1 0 0 0 0 + 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0      0 1 0 0 0 0 0 0 0 0 0 1 0
0 0 1 0 0 0 0 0 0 0 1 0 0      0 0 1 0 0 0 0 0 0 0 1 0 0
0 0 0 1 0 0 0 0 0 1 0 0 0      0 0 0 1 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0      0 0 0 0 1 1 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0 0 0 0

r0=6  dr=0.2 (12 pixels)        r0=6  dr=0.5 (28 pixels)      
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 0 0 0 0 0 0     0 0 0 0 0 1 1 1 1 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 1 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 1 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
0 1 0 0 0 0 0 0 0 0 0 0 0 1 0     0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
0 1 0 0 0 0 0 + 0 0 0 0 0 1 0     0 1 0 0 0 0 0 + 0 0 0 0 0 1 0
0 1 0 0 0 0 0 0 0 0 0 0 0 1 0     0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 1 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 1 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 1 1 0 0 0 0 0 0     0 0 0 0 0 1 1 1 1 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Hit finders

Hit finders return simple values for decision on event selection. Two algorithms are implemented in ImgAlgos.PyAlgos. They count number of pixels and intensity above threshold in the Region Of Interest (ROI) defined by windows and mask parameters in object constructor.

Number of pixels above threshold

number_of_pix_above_thr

npix = alg.number_of_pix_above_thr(data, thr=10)

Total intensity above threshold

intensity_of_pix_above_thr

intensity = alg.intensity_of_pix_above_thr(data, thr=12)

Peak finders

Two threshold "Droplet finder"

two-threshold peak-finding algorithm in restricted region around pixel with maximal intensity.

peak_finder_v1

peaks = alg.peak_finder_v1(nda, thr_low=10, thr_high=150, radius=5, dr=0.05)

Flood filling algorithm

define peaks for regoins of connected pixels above threshold

peak_finder_v2

peaks = alg.peak_finder_v2(nda, thr=10, r0=5, dr=0.05)

Local maximums search algorithm

define peaks in local maximums of specified rank (radius), for example rank=2 means 5x5 pixel region around central pixel.

peak_finder_v3

peaks = alg.peak_finder_v3(nda, rank=2, r0=5, dr=0.05)

makes a map of pixels with local maximums of requested rank for data ndarray and mask, pixel code in the map may have bits 0/1/2/4 standing for not-a-maximum / maximum-in-row / maximum-in-column / maximum-in-rectangular-region of radius=rank.
for each pixel with local maximal intensity in the region defined by the rank radius counts a number of pixels with intensity above zero, total positive intensity, center of gravity coordinates and rms,
using parameters r0(ex.=5.0), dr(ex.=0.05) evaluates background level, rms of noise, and S/N for the pixel with maximal intensity.
makes list of peaks which pass selector with parameters set in set_peak_selection_pars, for example

alg.set_peak_selection_pars(npix_min=5, npix_max=500, amax_thr=0, atot_thr=1000, son_min=6)

Demonstration for local maximum map

Test for 100x100 image with random normal distribution of intensities

Example of the map of local maximums found for rank from 1 to 5:

color coding of pixels:

blue=0 - not a local maximum
green=1 - local maximum in row
yellow=1+2 - local maximum in row and column
red=1+2+4 - local maximum in rectangular region of radius=rank.

Table for rank, associated 2-d region size, fraction of pixels recognized as local maximums for rank, and time consumption for this algorithm.

rank	2-d region	fraction	time, ms
1	3x3	0.1062	5.4
2	5x5	0.0372	5.2
3	7x7	0.0179	5.1
4	9x9	0.0104	5.2
5	11x11	0.0066	5.2

References

ImgAlgos.PyAlgos - code example in Sphinx documentation
Peak Finding - short announcement about peak finders
Hit and Peak Finders - examples in Chris' tutorial
Peak Finding Module - (depricated) psana module, it demonstaration examples and results
Psana Module Catalog - (depricated) peak finding psana modules
Psana Module Examples - (depricated) peak finding examples in psana modules
GUI for tuning peak finding - Chun's page in development