Locate peak on 2-D area detector image

From APS Python Training for Bluesky Data Acquisition.

Objective

Measure a simulated diffraction spot with an EPICS area detector and find its centroid and width.

Note: The position, width, and peak intensity of the simulated is randomized and will be different each time this note book is run. In fact, the EPICS support applies a small amount of jitter to the position to better simulate a real signal.

Start the instrument package

Our instrument package is in the bluesky subdirectory here so we add that to the search path before importing it.

import pathlib, sys
sys.path.append(str(pathlib.Path.home() / "bluesky"))
from instrument.collection import *

/home/prjemian/bluesky/instrument/_iconfig.py
Activating auto-logging. Current session state plus future input saved.
Filename       : /home/prjemian/Documents/projects/BCDA-APS/bluesky_training/docs/source/howto/.logs/ipython_console.log
Mode           : rotate
Output logging : True
Raw input log  : False
Timestamping   : True
State          : active

I Fri-15:42:25 - ############################################################ startup
I Fri-15:42:25 - logging started
I Fri-15:42:25 - logging level = 10
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/session_logs.py
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/collection.py
I Fri-15:42:25 - CONDA_PREFIX = /home/prjemian/.conda/envs/bluesky_2023_2

Exception reporting mode: Minimal

I Fri-15:42:25 - xmode exception level: 'Minimal'
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/mpl/notebook.py
I Fri-15:42:25 - #### Bluesky Framework ####
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/framework/check_python.py
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/framework/check_bluesky.py
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/framework/initialize.py
I Fri-15:42:25 - RunEngine metadata saved in directory: /home/prjemian/Bluesky_RunEngine_md
I Fri-15:42:25 - using databroker catalog 'training'
I Fri-15:42:25 - using ophyd control layer: pyepics
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/framework/metadata.py
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/epics_signal_config.py
I Fri-15:42:25 - Using RunEngine metadata for scan_id
I Fri-15:42:25 - #### Devices ####
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/devices/area_detector.py
I Fri-15:42:25 - /home/prjemian/bluesky/instrument/devices/calculation_records.py
I Fri-15:42:27 - /home/prjemian/bluesky/instrument/devices/fourc_diffractometer.py
I Fri-15:42:28 - /home/prjemian/bluesky/instrument/devices/ioc_stats.py
I Fri-15:42:28 - /home/prjemian/bluesky/instrument/devices/kohzu_monochromator.py
I Fri-15:42:28 - /home/prjemian/bluesky/instrument/devices/motors.py
I Fri-15:42:28 - /home/prjemian/bluesky/instrument/devices/noisy_detector.py
I Fri-15:42:28 - /home/prjemian/bluesky/instrument/devices/scaler.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/shutter_simulator.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/simulated_fourc.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/simulated_kappa.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/slits.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/sixc_diffractometer.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/devices/temperature_signal.py
I Fri-15:42:29 - #### Callbacks ####
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/callbacks/spec_data_file_writer.py
I Fri-15:42:29 - #### Plans ####
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/plans/lup_plan.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/plans/peak_finder_example.py
I Fri-15:42:29 - /home/prjemian/bluesky/instrument/utils/image_analysis.py
I Fri-15:42:29 - #### Utilities ####
I Fri-15:42:29 - writing to SPEC file: /home/prjemian/Documents/projects/BCDA-APS/bluesky_training/docs/source/howto/20230414-154229.dat
I Fri-15:42:29 -    >>>>   Using default SPEC file name   <<<<
I Fri-15:42:29 -    file will be created when bluesky ends its next scan
I Fri-15:42:29 -    to change SPEC file, use command:   newSpecFile('title')
I Fri-15:42:29 - #### Startup is complete. ####

Count the detector : record an image.

RE(bp.count([adsimdet], md={"motive": "locate_image_peak"}))

Transient Scan ID: 970     Time: 2023-04-14 15:42:30
Persistent Unique Scan ID: '8923a036-6063-4115-8763-ad6ce324df07'
New stream: 'label_start_motor'
New stream: 'primary'
+-----------+------------+
|   seq_num |       time |
+-----------+------------+
|         1 | 15:42:31.6 |
+-----------+------------+
generator count ['8923a036'] (scan num: 970)

('8923a036-6063-4115-8763-ad6ce324df07',)

View the image

Since area detector images are typically quite large, they are never read into Python memory during data acquisition. To view the image here, we need to use the databroker which reads the run data. The run data includes the reference to the file that contains the image.

Get the run

First, get the most recent run from the catalog of runs.

run = cat[-1]
run

BlueskyRun
  uid='8923a036-6063-4115-8763-ad6ce324df07'
  exit_status='success'
  2023-04-14 15:42:30.957 -- 2023-04-14 15:42:31.658
  Streams:
    * primary
    * label_start_motor

Look at the run’s metadata. Confirm it has the {"motive": "locate_image_peak"} metadata we added. That was stored in the start document.

print(run.metadata["start"]["motive"])

locate_image_peak

Primary stream

The image was collected in the primary data stream. Get the primary stream data from the database.

dataset = run.primary.read()
dataset

<xarray.Dataset>
Dimensions:         (time: 1, dim_0: 100, dim_1: 1024, dim_2: 1024)
Coordinates:
  * time            (time) float64 1.682e+09
Dimensions without coordinates: dim_0, dim_1, dim_2
Data variables:
    adsimdet_image  (time, dim_0, dim_1, dim_2) uint8 17 8 19 21 ... 10 21 8 17

xarray.Dataset

• Dimensions:
  • time: 1
  • dim_0: 100
  • dim_1: 1024
  • dim_2: 1024
• Coordinates: (1)
  • time
    (time)
    float64
    1.682e+09

    array([1.681505e+09])

• Data variables: (1)
  • adsimdet_image
    (time, dim_0, dim_1, dim_2)
    uint8
    17 8 19 21 10 20 ... 20 10 21 8 17

    object :

    adsimdet

    array([[[[17,  8, 19, ..., 10, 19, 21],
             [19, 15, 23, ...,  8, 15, 14],
             [ 7, 26, 15, ..., 12, 10,  8],
             ...,
             [11, 21, 16, ..., 20, 16,  7],
             [ 9, 12, 19, ..., 10,  7, 14],
             [14, 20,  9, ..., 19, 22,  8]],
    
            [[25, 18, 25, ..., 17, 11, 15],
             [14, 13, 25, ..., 26, 16, 25],
             [ 8, 10, 20, ..., 19, 10, 18],
             ...,
             [18, 15,  9, ..., 18, 24, 12],
             [17, 18,  8, ..., 21, 21, 17],
             [24, 16, 20, ..., 18, 14, 22]],
    
            [[10, 13, 23, ...,  8, 12, 13],
             [18, 19, 22, ..., 10,  7, 23],
             [ 9, 11,  8, ..., 22, 19, 20],
             ...,
    ...
             ...,
             [ 9, 24, 14, ..., 16, 16, 24],
             [17, 13, 13, ..., 21, 21, 21],
             [ 7, 20, 11, ...,  9, 14, 15]],
    
            [[14, 24, 21, ..., 24, 23, 23],
             [19, 24, 22, ..., 25, 22, 24],
             [14, 15, 14, ..., 20, 25,  7],
             ...,
             [11, 15, 13, ..., 10, 13, 18],
             [23,  7, 12, ..., 11, 22, 15],
             [26, 15,  9, ..., 22, 25, 17]],
    
            [[13, 23, 22, ..., 16,  7, 10],
             [21, 11, 16, ..., 15, 22, 18],
             [23, 13,  7, ..., 11, 13, 10],
             ...,
             [ 9, 23, 21, ..., 20, 16, 21],
             [16,  7, 18, ..., 21, 25,  9],
             [20, 18, 24, ..., 21,  8, 17]]]], dtype=uint8)

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(Float64Index([1681504951.6499872], dtype='float64', name='time'))

• Attributes: (0)

Get the image data

Get the named image adsimdet_image from the primary stream.

image = dataset["adsimdet_image"]
image

<xarray.DataArray 'adsimdet_image' (time: 1, dim_0: 100, dim_1: 1024,
                                    dim_2: 1024)>
array([[[[17,  8, 19, ..., 10, 19, 21],
         [19, 15, 23, ...,  8, 15, 14],
         [ 7, 26, 15, ..., 12, 10,  8],
         ...,
         [11, 21, 16, ..., 20, 16,  7],
         [ 9, 12, 19, ..., 10,  7, 14],
         [14, 20,  9, ..., 19, 22,  8]],

        [[25, 18, 25, ..., 17, 11, 15],
         [14, 13, 25, ..., 26, 16, 25],
         [ 8, 10, 20, ..., 19, 10, 18],
         ...,
         [18, 15,  9, ..., 18, 24, 12],
         [17, 18,  8, ..., 21, 21, 17],
         [24, 16, 20, ..., 18, 14, 22]],

        [[10, 13, 23, ...,  8, 12, 13],
         [18, 19, 22, ..., 10,  7, 23],
         [ 9, 11,  8, ..., 22, 19, 20],
         ...,
...
         ...,
         [ 9, 24, 14, ..., 16, 16, 24],
         [17, 13, 13, ..., 21, 21, 21],
         [ 7, 20, 11, ...,  9, 14, 15]],

        [[14, 24, 21, ..., 24, 23, 23],
         [19, 24, 22, ..., 25, 22, 24],
         [14, 15, 14, ..., 20, 25,  7],
         ...,
         [11, 15, 13, ..., 10, 13, 18],
         [23,  7, 12, ..., 11, 22, 15],
         [26, 15,  9, ..., 22, 25, 17]],

        [[13, 23, 22, ..., 16,  7, 10],
         [21, 11, 16, ..., 15, 22, 18],
         [23, 13,  7, ..., 11, 13, 10],
         ...,
         [ 9, 23, 21, ..., 20, 16, 21],
         [16,  7, 18, ..., 21, 25,  9],
         [20, 18, 24, ..., 21,  8, 17]]]], dtype=uint8)
Coordinates:
  * time     (time) float64 1.682e+09
Dimensions without coordinates: dim_0, dim_1, dim_2
Attributes:
    object:   adsimdet

xarray.DataArray

'adsimdet_image'

• time: 1
• dim_0: 100
• dim_1: 1024
• dim_2: 1024

• 17 8 19 21 10 20 13 24 22 11 7 18 ... 21 22 23 10 24 20 20 10 21 8 17

  array([[[[17,  8, 19, ..., 10, 19, 21],
           [19, 15, 23, ...,  8, 15, 14],
           [ 7, 26, 15, ..., 12, 10,  8],
           ...,
           [11, 21, 16, ..., 20, 16,  7],
           [ 9, 12, 19, ..., 10,  7, 14],
           [14, 20,  9, ..., 19, 22,  8]],
  
          [[25, 18, 25, ..., 17, 11, 15],
           [14, 13, 25, ..., 26, 16, 25],
           [ 8, 10, 20, ..., 19, 10, 18],
           ...,
           [18, 15,  9, ..., 18, 24, 12],
           [17, 18,  8, ..., 21, 21, 17],
           [24, 16, 20, ..., 18, 14, 22]],
  
          [[10, 13, 23, ...,  8, 12, 13],
           [18, 19, 22, ..., 10,  7, 23],
           [ 9, 11,  8, ..., 22, 19, 20],
           ...,
  ...
           ...,
           [ 9, 24, 14, ..., 16, 16, 24],
           [17, 13, 13, ..., 21, 21, 21],
           [ 7, 20, 11, ...,  9, 14, 15]],
  
          [[14, 24, 21, ..., 24, 23, 23],
           [19, 24, 22, ..., 25, 22, 24],
           [14, 15, 14, ..., 20, 25,  7],
           ...,
           [11, 15, 13, ..., 10, 13, 18],
           [23,  7, 12, ..., 11, 22, 15],
           [26, 15,  9, ..., 22, 25, 17]],
  
          [[13, 23, 22, ..., 16,  7, 10],
           [21, 11, 16, ..., 15, 22, 18],
           [23, 13,  7, ..., 11, 13, 10],
           ...,
           [ 9, 23, 21, ..., 20, 16, 21],
           [16,  7, 18, ..., 21, 25,  9],
           [20, 18, 24, ..., 21,  8, 17]]]], dtype=uint8)

• Coordinates: (1)
  • time
    (time)
    float64
    1.682e+09

    array([1.681505e+09])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(Float64Index([1681504951.6499872], dtype='float64', name='time'))

• Attributes: (1)

  object :

  adsimdet

Image frame

This image object has rank of 4. We just want the image frame (the last two indices). Select the first item of each of the first two indices (time, frame number).

frame = image[0][0]
frame

<xarray.DataArray 'adsimdet_image' (dim_1: 1024, dim_2: 1024)>
array([[17,  8, 19, ..., 10, 19, 21],
       [19, 15, 23, ...,  8, 15, 14],
       [ 7, 26, 15, ..., 12, 10,  8],
       ...,
       [11, 21, 16, ..., 20, 16,  7],
       [ 9, 12, 19, ..., 10,  7, 14],
       [14, 20,  9, ..., 19, 22,  8]], dtype=uint8)
Coordinates:
    time     float64 1.682e+09
Dimensions without coordinates: dim_1, dim_2
Attributes:
    object:   adsimdet

xarray.DataArray

'adsimdet_image'

• dim_1: 1024
• dim_2: 1024

• 17 8 19 21 10 20 13 24 22 11 7 18 ... 13 20 21 9 11 24 23 9 16 19 22 8

  array([[17,  8, 19, ..., 10, 19, 21],
         [19, 15, 23, ...,  8, 15, 14],
         [ 7, 26, 15, ..., 12, 10,  8],
         ...,
         [11, 21, 16, ..., 20, 16,  7],
         [ 9, 12, 19, ..., 10,  7, 14],
         [14, 20,  9, ..., 19, 22,  8]], dtype=uint8)

• Coordinates: (1)
  • time
    ()
    float64
    1.682e+09

    array(1.68150495e+09)

• Indexes: (0)
• Attributes: (1)

  object :

  adsimdet

Show the image

frame.plot.pcolormesh()

<matplotlib.collections.QuadMesh at 0x7f50794f3880>

Center & size of spot

Determine the center and width of this spot in both directions. We’ll use our own function (from instrument/utils/image_analysis.py) since no package gives us all of this in one place. Our function is extracted and modified from the PeakStats support in bluesky’s BestEffortCallback and also uses a function from scipy.

Note that the centroid measure here is the channel number weighted by the intensity value (as computed by the scipy package). The fwhm measure here (from the BestEffortCallback) is the computed difference in channel number between the high and low side of the maximum at half the maximum-minimum value. The crossings measure shows the interpolated locations at which the intensity crossed this half value. The other measures should be self-explanatory.

analyze_image(frame)

================= =========================== =====================================================
measure           vertical (dim_1)            horizontal (dim_2)
================= =========================== =====================================================
centroid_position 476.14470403546034          265.72533623826075
fwhm              192.45267337704308          239.67524162973245
crossings         [379.91836735 572.37104072] [ 86.59934853 324.93690852 325.09049774 326.27459016]
maximum           (476, 52255)                (205, 46040)
center_position   501.04078251008247          421.7680222410236
minimum           (813, 15953)                (819, 15960)
================= =========================== =====================================================

Next Steps

Make these additions or improvements:

1. Repeat the same analysis for all area detector images in the catalog. (hint: for uid in catalog: run = catalog[uid] ...)

2. Computed centroid_position and center_position do not agree. Explain why.

3. Sort the catalog based on some metadata key such as scan_id (hint: run.metadata["start"][key])