# Welcome to ImageD11 : Peaksearch

Peaksearching in two dimensional images can be carried out in a relatively automated way using the `peaksearch.py` script.
Giving the command "peaksearch.ph -h" should print the following help message:

## options:

`-h`, `--help` show this help message and exit
`-n STEM`,`--namestem=STEM` Name of the files up the digits part, eg mydata in mydata0000.edf
`-F FORMAT`, `--format=FORMAT` Image File format, eg edf or bruker
`-S OMEGASTEP`, `--step=OMEGASTEP` Step size in Omega when you have no header info
`-T OMEGA`, `--start=OMEGA` Start position in Omega when you have no header info
`-f FIRST`, `--first=FIRST` Number of first file to process, default=0
`-l LAST`, `--last=LAST` Number of last file to process
`-o OUTFILE`, `--outfile=OUTFILE` Output filename, gets extension .spt
`-d DARK`, `--darkfile=DARK` Dark current filename, to be subtracted, default=None
`-D DARKOFFSET`, `--darkfileoffset=DARKOFFSET` Constant to subtract from dark to avoid overflows, default=100
`-s SPLINE`, `--splinefile=SPLINE` Spline file for spatial distortion, default=/data/opid11/inhouse/Frelon2K/spatial2k.spline
`-p PERFECT`, `--perfect_images=PERFECT` Ignore spline Y|N, default=N
`-O FLOOD`, `--flood=FLOOD` Flood file, default=None
`-t THRESHOLDS`, `--threshold=THRESHOLDS` Threshold level, you can have several

Usually data images are collected as a numerical series while the sample rotates and are named using the convention stem0000.edf, stem0001.edf etc.
The first three options `"-n"`, `"-f"` and `"-l"` specify the series of images in terms of stem name, start and end number. Currently the script will automatically assume ".edf" or ".edf.gz" extensions and corresponding format. If you want for example ".tif.gz" say `"-F .tif.gz"`. It should accept the following formats:

".edf",".edf.gz",".edf.bz2",
".tif",".tif.gz",".tif.bz2",
".mccd",".mccd.gz",".mccd.bz2"

The `"-d"` option is optional, for a dark current image that can be subtracted. If you use the Frelon2K CCD camera at ESRF the please allow for the offset of 1000 ADU if you do not decide to subtract a dark current.
The `"-D"` option alters the offset to try to avoid overflow when subtraction gives a negative number (it is though to be marginally faster to remain with unsigned short data than going over to int)
The `"-s"` option specifies a `fit2d` spline file characterising the spatial distortion of the camera. The default is probably OK for recent data from ID11, but more generally you will want to change this option. Be sure to use the pixel size data from this file when computing the detector calibration parameters. In the case where no spatial corrections are needed then please be sure to specify the `"-p Y"` option to indicate there is no need for spatial distortion correction (eg a nice image plate). Corrections are applied to peaks positions, NOT to the raw images.
Threshold levels are used to determine which pixels contribute to diffraction spots. A segmentation of the image will be made where all pixels above the threshold are the collected into connected objects. Try to choose levels which are significantly above the background level ; you can always have some help using `Imageviewer`. You may use more than one level if you like, eg "-t 1050 -t 5000" for levels 1050 ADU and 5000 ADU.
The output file contains all of the information contained in the header as well as the obtained peak positions. You can specify a more logical name for this file if you do not like "peaks.out".
Some examples:

`peaksearch.py -n plat_a -f 50 -l 90 -o peaks_50_90.out -s spatial.spline -t 5000`
This will peaksearch images:

- `plat_a0050.edf`
- `plat_a0051.edf`
- ...
- `plat_a0090.edf`

A spatial distortion computed from the file spatial.spline (in the current directory) will be applied. A threshold level of 5000 is used to decide which pixels contribute to peaks. Example output:

``

# File plat_a0050.edf
# Processed on Wed Nov 30 17:09:34 2005
# Spatial correction from /data/opid11/inhouse/Frelon2K/spatial2k.spline
# SPLINE X-PIXEL-SIZE 46.776480
# SPLINE Y-PIXEL-SIZE 48.081500
# HeaderID = EH:000001:000000:000000
# y1 = -3.88889e-05
# Instron strain = 0.002
[... lots of header information snipped ...]
# Omega = -39
#Threshold level 1200.000000
# Number_of_pixels Average_counts x y xc yc sig_x sig_y cov_xy
119 2371.050420 377.755365 1467.049632 372.961306 1468.529932 3.122391 2.406393 -0.363226
44 1600.363636 424.880254 1857.962651 420.424252 1852.106077 1.570551 2.092762 -0.350194
41 1495.170732 1118.681952 1952.197416 1104.657513 1933.297804 1.19973 6 2.575913 -0.063616
# File plat_a0050.edf
[ ...etc... ]

The numbers in the output file have the following meanings:

- `Number_of_pixels` is the number of pixels in the connected object.
- `Average_counts` is the average counts in each pixel in that object. If the camera has a constant offset (in this case about 1000) then you need to subtract that number before computing intensities. In this case intensity is approximately number_of_pixels*(Average_counts-1000.0)
- `x` and `y` are the position of the centre of mass of the peak in the uncorrected image (first moments of the intensity distributions in the blobs is the mean)
- `xc` and `yc` are the x and y positions corrected for spatial distortion.
- `sig_x` and `sig_y` are the second moments of the intensity distribution in the blob. Something like the width, but off by some factor of twopi perhaps (FIXME)
- `cov_xy` is the covariance (the third one of the second moments). It ranges between -1 and 1 with a value of 0 for a circular peak and +1 and -1 refering to elliptical shapes rotated by 90 degrees from each other. (FIXME a picture would help)

The definition of `x` and`y` are in terms of the `fast` and `slow` array indices as the image comes into memory, so it depends on how the image was stored in the file and the routine which read it in as to what you might finally get.
This should later be better defined by the use of a better image file format, like `ImageCIF`, which defines such things.
For most experiments at ID11 the position in the rotation scan should be specified by a value "Omega" in the image headers. If this is not the case there is a script called recoveromega.py which will read a text file containing lines with omega values followed by filenames to recover the appropriate information.

## How does it work?

For now, we are talking about the most up to date SVN version (0.8.1), a new release will be made soon.
The algorithm used is based on a disjoint set, which is described in "Introduction to Algorithms" by Cormen, Leiserson., Rivest and Stein. An image is scanned row by row and each pixel is compared to a threshold value. If the pixel is above the threshold, it will be labelled as a peak. To determine the labels the pixel is compared to the previous pixel, and the pixels on the previous row. If one of these pixels is already labelled, the current pixel takes the same label. When there is a disagreement about labelling, the two labels are made the same using the "disjoint set"

oooooooooooo o = pixel below threshold
oo1ooooo2ooo 1,2 = labels
o111ooo222oo X pixel where label 1 and label 2 must be made equivalent
oo1111X

Within the sourcecode the routine "connectedpixels" is a compiled extension in C which takes a data image and computes an image of peak assignment labels. It optionally applies the dark and flood corrections when thresholding. The connectivity in the 2D image is therefore:

0 | 0 | 0
0 | X | 0
0 | 0 | 0

To compute the properties of each connected object the routine "blobproperties" is then called with the data image and label image (this could be more efficient). It scans through the image forming properties on a pixel by pixel basis.
In order to merge peaks on adjacent frames, the label images are compared. When exactly overlapping pixels are both labelled, then these labels should be made equivalent. The code which does it is called "bloboverlaps".

Files which are produced for each threshold are like:

- name_t500.spt add ".spt" to users request : 2D peaks, header info
- name_t500.flt (add "_t500.flt" to users request) : 3D peaks for threshold at 500

### Using peaksearch for programmers

For the full pleasures of 3D peaksearching in a program designed for fast 2D, one should read the ImageD11.labelimage code. The quick way to get an idea what thresholding will do is this:

- data = fabio.openimage("my_lovely_data.edf")
- blobs = numpy.zeros(data.shape, numpy.int)
- threshold = 1500.0
- npks = connectedpixels.connectedpixels(data, blobs, threshold, verbose=0)

On exit from connectedpixels the npks is the number of peaks found and blobs is an array of integer peak assignments (the argument is supposed to be modified during the function call).
You can then have a look at blobs to see where the peaks are. In practice the peaksearch code will be doing a dark and flood correction too, which can add to the confusion.