Inspect a reconstructed point#

A reconstructed point is a stack of frames, one per sample depth. To examine how a feature varies with depth, choose a reference image, place square regions of interest (ROIs) on it, and plot the intensity of each region through depth. The inspection functions support both Jupyter notebooks and the Laue Portal.

Intensity traces use the stored pixels, after conversion to the output pixel type. In a signed output type, negative values are kept, so a trace can cross zero. Reconstruction scan format defines stored, computed, and raw values.

The examples inspect a point reconstructed with reconstruct_scan(); see Reconstruct a wire scan. Replace run/scan.h5 and the point ID with your own. The same functions accept a PointReader opened directly on a point file.

Choose a reference image#

Each point file contains three images with the shape of one frame:

Kind	Values	Meaning
`"sum_reconstructed"`	stored	Sum of the stored frames through depth; the default
`"first_raw"`	raw	The first scan frame the reconstruction read
`"sum_raw"`	raw	Sum of the scan frames, before filtering and normalization

reference_image() returns one of them as a ReferenceImage, containing the image array, its kind, and whether it contains raw or stored values.

from lauelab.reconstruct import ScanReader
from lauelab.reconstruct.inspection import reference_image

with ScanReader("run/scan.h5") as scan:
    point = scan.point("scan12_p1")
    background = reference_image(point, "sum_reconstructed")

assert background.values == "stored"
rows, columns = background.shape

Use the reference image to locate features and place ROIs. Traces are calculated from the reconstructed stack, independently of the reference image selected.

Place square ROIs#

An ROI is a square of N × N stored-image pixels. square_bounds() places it from a size and a click position and returns half-open bounds (y0, y1, x0, x1), which select frame[y0:y1, x0:x1]. Positions are zero-based stored-image pixel coordinates (x, y); an integer is a pixel centre.

from lauelab.reconstruct.inspection import bounds_center, square_bounds

spot = square_bounds(5, 64.2, 63.8, (rows, columns))
edge = square_bounds(4, 52.0, 70.0, (rows, columns))

assert spot == (62, 67, 62, 67)
assert bounds_center(spot) == (64.0, 64.0)
assert bounds_center(edge) == (51.5, 69.5)

Pixel (x, y) covers x - 0.5 to x + 0.5 on each axis, and NumPy reads it as frame[y, x]. ROI coordinates are stored-image pixels; unbinned detector coordinates are not used. The bounds select data[:, y0:y1, x0:x1], which is exactly N × N pixels for a positive integer size N.

For a click at (x, y), the centre of the square is the nearest legal centre. The legal centres of an odd size are the integers, and those of an even size are the half-integers. When a click is at equal distance from two legal centres, the lower coordinate is selected. On each axis this gives

\[x_0 = \lceil x - N/2 \rceil, \qquad x_1 = x_0 + N,\]

and bounds_center() returns the centre \(x_0 + (N - 1)/2\), which is derived from the bounds. square_bounds raises InputError when any bound lies outside the image or the click is not finite. Accepted bounds are used unchanged at every display scale. Each call creates new bounds, so existing ROIs retain their size.

ROIs are created during analysis and are not stored in the point file. The datasets /entry1/detector/startx, starty, binx, and biny of the point file, copied from the input, contain the acquisition readout region and binning that indexing uses; they are not ROIs. lauelab.reconstruct.inspection implements these rules, and the builders in lauelab.visualization display them. The fixtures under tests/data/reconstruction_contract/ are test cases for these rules, so the Laue Portal and lauelab use the same definition.

Trace intensity through depth#

depth_trace() sums the stored pixels of a region at every depth and returns a DepthTrace. With no bounds, it returns the full-frame trace from the reduction stored in the file, without reading pixels. roi_traces() does the same for several named regions.

from lauelab.reconstruct.inspection import depth_trace, roi_traces

with ScanReader("run/scan.h5") as scan:
    point = scan.point("scan12_p1")
    whole = depth_trace(point)
    traces = roi_traces(point, {"spot": spot, "edge": edge})

peak_depth_um = whole.depth_um[whole.values.argmax()]
spot_peak_um = traces["spot"].depth_um[traces["spot"].values.argmax()]

ROI reductions read one depth block at a time, using at most 64 MiB of pixel data by default. Set max_bytes on depth_trace or roi_traces to choose another budget; it must hold at least one plane of the ROI. The small returned traces and HDF5 caches are additional. Multiple ROIs are reduced in sequence.

values has shape (n_depths,): the sum of the stored pixels in the bounds at each depth. Integer pixels are summed exactly in numpy.int64; floating-point pixels are summed in numpy.float64. The full-frame trace is /entry1/reconstruction/depth_intensity of the point file. roi_traces with no ROIs returns an empty dictionary. For the horizontal axis, use depth_um for physical depth in µm or depth_index for the zero-based frame index.

A trace has two methods for display:

trace.normalized() divides the trace by its own maximum. Zero and the sign of each sample are kept. When the maximum is zero or negative, values is None and reason contains the cause.
trace.log_samples() returns the indices of positive samples for a logarithmic axis and counts the zero or negative samples omitted.

normalized = traces["spot"].normalized()
assert normalized.available and normalized.values.max() == 1.0

samples = traces["edge"].log_samples()
n_shown, n_omitted = len(samples.kept), samples.n_omitted

Plot#

The figure builders in lauelab.visualization take prepared data and display choices. A change of axis, colour, or normalization uses the prepared arrays and does not read the file again.

from lauelab.visualization import (
    RoiOverlay, plot_depth_trace, plot_reference_image, plot_roi_traces,
)

overlays = [RoiOverlay("spot", spot, "rgb(230,90,60)"), RoiOverlay("edge", edge, "rgb(60,140,230)")]
image_figure = plot_reference_image(background, overlays)
whole_figure = plot_depth_trace(whole, axis="depth", selected_index=int(whole.values.argmax()))
roi_figure = plot_roi_traces(
    traces, colors={"spot": "rgb(230,90,60)", "edge": "rgb(60,140,230)"}, normalized=True,
)

plot_reference_image() places the origin at the upper left, with x increasing to the right and y downward. Both axes use the same pixel scale, and ROI outlines follow pixel edges.

plot_depth_trace() plots intensity against physical depth in µm or depth index, with an optional marker at the selected depth. plot_roi_traces() plots one trace per ROI. If normalization is unavailable, the plot omits that trace and annotates the reason. On a logarithmic axis, it omits nonpositive samples and reports their count. The plotted line has a gap at each omitted sample. Supply matching names and colours to identify each ROI across the image and trace plots.

Every builder accepts layout_update and per-role trace_update; the roles are "image" and "roi" for the image, "trace" and "selected" for the single trace, and "roi" for the ROI traces. DEPTH_AXIS_OPTIONS, INTENSITY_OPTIONS, and REFERENCE_OPTIONS list the choices for controls.

Inspect an in-memory result#

ArrayPoint gives a (depth, y, x) array the same interface, so a result of Reconstructor.reconstruct_array() or of Reconstructor.reconstruct(..., return_images=True) can be inspected without writing a file. The array is used as the stored data, so a numpy.float64 stack remains numpy.float64. It provides the "sum_reconstructed" reference; requests for raw reference images raise InputError.

from lauelab.reconstruct.inspection import ArrayPoint

point = ArrayPoint(result.images, result.depth_um, "in-memory")
trace = depth_trace(point, spot)

Inspect existing per-depth files#

PerDepthReader supplies the same frame, region, trace, and reconstructed-reference operations for the individual HDF5 files written by Reconstructor, reconstruct, or export_per_depth. Give it the files for one point, excluding the summary text file. It orders them by embedded physical depth and checks that their shapes, dtypes, detector metadata and normalization agree. It reads stored values without rescaling them again.

from pathlib import Path
from lauelab.reconstruct import PerDepthReader
from lauelab.reconstruct.inspection import depth_trace, reference_image

paths = Path("run/export").glob("scan12_p1_*.h5")
with PerDepthReader(paths, point_id="scan12_p1") as point:
    whole = depth_trace(point)
    background = reference_image(point, "sum_reconstructed")
    trace = depth_trace(point, (60, 65, 50, 55))

These files have no embedded full-frame reduction, so depth_trace(point) reads one frame at a time. The calculation of the reconstructed reference also reads one frame at a time. Raw references are unavailable and raise InputError; retain the point file if you need those backgrounds. The reader opens at most one file at a time and closes it after each operation. Exiting the context prevents further reads; returned arrays remain usable.

Do your own analysis#

For custom analysis, use the NumPy arrays in the prepared objects or read a selection directly with PointReader. frame() reads one depth, and region() reads a rectangular region over the requested depth range.

import numpy as np

with ScanReader("run/scan.h5") as scan:
    point = scan.point("scan12_p1")
    y0, y1, x0, x1 = spot
    pixels = point.region(spot)                 # (n_depths, 5, 5) in the stored dtype
    brightest = point.frame(int(whole.values.argmax()))

profile = pixels.reshape(len(pixels), -1).max(axis=1)