Analysis of PyFAI geometry optimization

[LouConreux]

Analysis of PyFAI geometry optimization

Context

Following the geometry calibration on experiment mfxx49820 run 8 (see #197), we would like to translate the custom-made geometry optimization with a more-globally used library pyFAI. Using pyFAI would allow more robustness and rapidity in finding the optimal geometry of an experiment.

Geometry format

LCLS is using the CrystFEL format, so first one need to parse the CrystFEL geometry in a format that PyFAI can digest. Then optimizing through PyFAI reduces to finding the optimal 6 parameters defining the coordinates of the detector reference frame compared to the origin being the sample position (3 translations and 3 rotations). See https://www.silx.org/doc/pyFAI/dev/geometry.html#geometry for more details. In the following, we will then try to optimize the geometry by playing on the 3 translation parameters: poni1, poni2 and dist (PONI standing for "Point of Normal Incidence" = the orthogonal projection of the sample position onto the detector plane) (We make the assumption that we are in a transmission setting meaning all rotations are 0).

Optimizing the Geometry

Ground Truth

To set the ground truth for making sure PyFAI finds the right geometry, we refer to the analysis of #197 where the center (in pixels) of the geometry was found to be

New center is (843.5532519182825, 831.0710902924419) \\
Detector distance inferred from powder rings: 280.47 mm

with all rotations being 0.

In PyFAI reference frame, this coordinate corresponds to approximately (x=0.5, y=5)

Optimization

PyFAI struggles to find the optimal geometry with only one shot. So, we will apply the optimization recursively until convergence.

Initialization

Geometry Format Conversion

First, need to convert the CrystFEL format into a PyFAI-friendly format:

geom_path = '/sdf/data/lcls/ds/mfx/mfxx49820/scratch/dorlhiac/geom/r0008.geom'
converter = CrystFEL_to_PyFAI(geom_path)
epix10k2M = converter.detector
epix10k2M.set_pixel_corners(converter.corner_array)

Load Data

Load the calibrated images:

exp = 'mfxx49820'
run = 8
detector = 'epix10k2M'
psi = PsanaInterface(exp=exp, run=run, det_type=detector)
print(f'Instantiated run {run} of exp {exp} looking at {detector}')
print("Raw Images")
psi.calibrate = True
raw_images = psi.get_images(100, assemble=False)
print(f'Images shape: {raw_images.shape}')
calib_avg = np.mean(raw_images, axis=0)
calib_avg_flat = np.reshape(calib_avg, (16*2*176,2*192))

Here, we need to reshape the calib_avg (in (Np, Nx, Ny)) to the shape supported by PyFAI (Np*Nx, Ny)

Set Experiment Initial Conditions

Define the guessed dist, poni1, poni2 and the wavelength used.

wavelength = psi.get_wavelength() * 1e-10
dist = 0.4
p1, p2, p3 = epix10k2M.calc_cartesian_positions()
poni1 = +p1.mean()*2
poni2 = +p2.mean()*4

Calling the Geometry fitting

from pyFAI.azimuthalIntegrator import AzimuthalIntegrator
from pyFAI.gui import jupyter
ai = AzimuthalIntegrator(dist=dist, detector=epix10k2M, wavelength=wavelength, poni1=poni1, poni2=poni2)
print(ai)
res = ai.integrate1d(calib_avg_flat, 1000)
ax = jupyter.plot1d(res)

Here, we voluntarily defined parameters to be way off the true center. We get:

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.000000e-01 m	PONI= -1.452228e-03, -5.401834e-04 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 400.000 mm	Center: x=-5.402, y=-14.522 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Instantiating the PyFAI Calibrant

Here, this was Silver Behenate:

from pyFAI.calibrant import CalibrantFactory, CALIBRANT_FACTORY
behenate = CALIBRANT_FACTORY("AgBh")
behenate.wavelength = wavelength

Optimization Loop

Below is what the optimization looks like. First, it instantiated the guessed geometry. Then it passes through a SingleGeometry object that extracts control points between raw data and calibrant data. Using those control points, a new fitting is finally set and we repeat until convergence.

from pyFAI.goniometer import SingleGeometry
geom_initial = pyFAI.geometry.Geometry(dist=dist, poni1=poni1, poni2=poni2, detector=epix10k2M, wavelength=wavelength)
# The SingleGeometry object (from goniometer) allows to extract automatically ring and calibrate
sg = SingleGeometry("test1", calib_avg_flat, calibrant=behenate, detector=epix10k2M, geometry=geom_initial)
sg.extract_cp(max_rings=5, pts_per_deg=0.2, Imin=3)
# We fix rot1 = rot2 = rot3 = 0 and wavelength should not be changed
sg.geometry_refinement.refine3(fix=["rot1", "rot2", "rot3", "wavelength"])
sg.geometry_refinement.param
geom_initial = pyFAI.geometry.Geometry(dist=sg.geometry_refinement.param[0], poni1=sg.geometry_refinement.param[1], poni2=sg.geometry_refinement.param[2], detector=epix10k2M, wavelength=wavelength)

Results

For n=5 loops, we obtained the following convergence:

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.169748e-01 m	PONI= -1.212360e-03, -5.071352e-04 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 416.975 mm	Center: x=-5.071, y=-12.124 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.210260e-01 m	PONI= -8.301785e-04, -3.989296e-04 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 421.026 mm	Center: x=-3.989, y=-8.302 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.217013e-01 m	PONI= 4.613455e-05, -3.942903e-05 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 421.701 mm	Center: x=-0.394, y=0.461 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.214463e-01 m	PONI= 3.905332e-04, 2.993149e-05 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 421.446 mm	Center: x=0.299, y=3.905 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.218728e-01 m	PONI= 4.825142e-04, 5.180493e-05 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 421.873 mm	Center: x=0.518, y=4.825 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

This pretty good for a start, not so far away from the ground truth center. However, PyFAI fails to find the best sample-detector distance which is dist=280mm. We can see that the more loops we perform, the more control points fit right the rings, which helps reaching convergence.

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 2.800000e-01 m	PONI= 5.000000e-04, 5.000000e-05 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 280.000 mm	Center: x=0.500, y=5.000 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Next Steps

Finding why dist cannot be optimized ?

Finding best hyperparameters for extract_cp

In this optimization, the most crucial function is extract_cp which automatically extract control points on rings.
We can play on the number of points we want on 3 hyperparameters:
max_rings selects the furthest rings where to extract control points
pts_per_deg selects the "spacing" or "resolution" between two extracted control points
Imin specifies a minimum intensity for which points below that intensity cannot be extracted

Testing robustness

Playing with different initialization conditions, setting the guessed geometry way way too far from the ground truth may result in divergence (maybe too few control points).

[LouConreux]

Update 04/24

Geometry Optimization Objective Workflow

The objective of this project is to eventually build a real-time Geometry Optimization that can be summarized in this big picture:

This workflow relies on the development of several processes:

1. Geometry conversion: transforming either a .geom or .data file in a pyFAI-friendly geometry instance.

2. The pyFAI optimizer that takes in input keV calibrated images and the detector geometry to optimize.

3. Transforming back the final pyFAI-ed geometry in the LCLS reference frame coordinates.

4. Deploying the calibrated geometry.

5. is already working, 2. is in development, 3. needs to be implemented, and 4. is already working.

pyFAI Optimization Description

pyFAI Implementation

An example tutorial of how pyFAI can be used to calibrate X-ray powder diffraction can be found [here].(https://pyfai.readthedocs.io/en/v2023.1/usage/tutorial/Recalib/Recalib_notebook.html)
This script relies on three major functions:

1. After defining the initial guessed geometry and the used calibrant, pyFAI automatically extracts control points on which the optimization will be performed. Those control points depends on the calibrant data and the initial guessed geometry. This extraction process can be tuned using 3 hyperparameters:
   max_rings : defines the interval of ring index in which control points must be extracted
   pts_per_deg : defines how many points can be extracted per unit degree. Higher values leads to higher resolution but more computation.
   Imin : defines the minimal intensity of points that can be extracted.

dist = 0.42
poni1 = 0
poni2 = 0
geom_initial = pyFAI.geometry.Geometry(dist=dist, poni1=poni1, poni2=poni2, detector=epix10k2M, wavelength=wavelength)
sg = SingleGeometry("initial", calib_max_flat, calibrant=behenate, detector=epix10k2M, geometry=geom_initial)
sg.extract_cp(max_rings=5, pts_per_deg=2, Imin=0*photon_energy)

2. From these control points, the score function to minimize is defined as:
   
   where theta_e depends on the 6 geometry parameters and the control points, and the theta_i are defined with the calibrant and wavelength data.
   From pyFAI:

def residu3(self, param, free, const, d1, d2, rings, weights=None):
        "Preform the calculation of $sum_(2\theta_e-2\theta_i)²$"
        param7 = self.calc_param7(param, free, const)
        delta_theta = self.tth(d1, d2, param7[:6]) - self.calc_2th(rings, param7[6])
        if weights is not None:
            delta_theta *= weights
        return numpy.dot(delta_theta, delta_theta)

3. The optimization relies on the scipy.optimize function called fmin_slsqp source
   which minimize the score function using Sequential Least Squares Programming

Score Function Mapping

Fixed Control Points

Starting from the guessing point dist=0.42, poni1=0pix, poni2=0pix,
we extract control points and maps the score function over a 2D subspace.

If we refine the mesh:

Thus, the optimization will not lead directly to global minimum which is around (poni1=5pix, poni2=0.5pix) as seen previously.
But at least the optimization is indicating the right direction. Hence, repeating the optimization until we reach convergence should work.

Variable Control Points

Here, we highlight the fact that in order to reach the global minimum, the control points have to be updated in between two iterations. We can see now the global minimum region if we map this new process:

Repeating the optimization loop a finite number of times should lead to convergence.

Hyperparameter Influence

Let's return to our first example with guessing point dist=0.28mm, poni1=0pix, poni2=0pix.
From this, we plot the calibrant data over the powder data.

We can see that the distance is way off the true value. In fact the optimal is more about dist=0.42. Hence, there is an offset of 0.14mm which turns out to be correponding to the camera offset of the experiment!

Max Rings

We can set it to 5 because after the 5th ring, we don't see that much the rings on the powder.

Points per degree

Hence for dist=0.42mm, poni1=0pix, poni2=0pix, if we zoom in one of the center panel:

We can see that the calibrant ring are not perfectly centered, so we extract points to perform the optimization. We can play on the pts_per_deg hyperparameter to refine the extraction (here 0.5, 1, 2 respectively).



Thus, it's preferable to have pts_per_deg large enough to have an accurate optimization, but not too large so the computation remains reasonably long. But what's problematic is that the extracted points are not fitting nicely the rings, the corresponding widths are too big. This could explain why only one iteration of the optimization is not immediately converging to the global minimum. So let's play now on Imin.

Imin

Plotting the same with Imin=0, 1, 2 and finally 3*photon_energy, we can filter out the background pixels.




It's still not improving that much the extraction, the noisy component is too high. Maybe it will need require some preprocessing on the powder data because it should fit nicely the rings like in the example:

Next Steps

Plotting the iteration succesively to understand how convergence is reached

Figure out the matrix relation between the LCLS pixel space and the pyFAI pixel space to change the reference frame.

A first idea solution can be investigated through repeating those operations: https://pyfai.readthedocs.io/en/v2023.1/geometry_conversion.html#geometry-definition-of-pyfai

[LouConreux]

Update 04/29

Extraction of Control Points & Convergence

A quick experiment was performed to measure the influence of the quality of control points on convergence. We start from a guessed geometry way off the true one:

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.409220e-01 m	PONI= 1.000000e-03, 1.000000e-03 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 440.922 mm	Center: x=10.000, y=10.000 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

We then extract the Control Points with:
sg.extract_cp(max_rings=5, pts_per_deg=1, Imin=0)
We get a lot of control points:

ControlPoints instance containing 5 group of point:
AgBh Calibrant with 49 reflections at wavelength 1.2909905595344899e-10
Containing 5 groups of points:
#at ring 0: 220 points
#au ring 1: 329 points
#av ring 2: 341 points
#aw ring 3: 347 points
#ax ring 4: 350 points

However, this results in poor quality control points which obviously lead to a suboptimal new geometry:

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.399303e-01 m	PONI= 9.736744e-04, 9.835861e-04 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 439.930 mm	Center: x=9.836, y=9.737 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

Since the noise component on the powder is very predominant, we need to refine the selection of points based on their intensity:
sg.extract_cp(max_rings=5, pts_per_deg=1, Imin=4*photon_energy)
Then, we get obviously a lot less of control points, but they are of definitely much better quality.

ControlPoints instance containing 5 group of point:
AgBh Calibrant with 49 reflections at wavelength 1.2909905595344899e-10
Containing 5 groups of points:
#ay ring 0: 83 points
#az ring 1: 5 points
#ba ring 2: 3 points
#bb ring 3: 0 points
#bc ring 4: 0 points

Giving this to the optimizer, we are able to make a much better significant improvement in a single step of optimization:

Detector Epix10k2M	 PixelSize= 1.000e-04, 1.000e-04 m	 BottomRight (3)
Wavelength= 1.290991e-10 m
SampleDetDist= 4.193749e-01 m	PONI= 4.912165e-04, 1.039317e-04 m	rot1=0.000000  rot2=0.000000  rot3=0.000000 rad
DirectBeamDist= 419.375 mm	Center: x=1.039, y=4.912 pix	Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.291Å

From this, we can deduce that the noise affects too much the extraction of high quality control points, thus leading to an ill-defined minimization problem.

Future Next Steps

It could be nice to see the what a preprocessing on the powder data can potentially add to the convergence!

Figure out the matrix relation between the LCLS pixel space and the pyFAI pixel space to change the reference frame.

A first idea solution can be investigated through repeating those operations: https://pyfai.readthedocs.io/en/v2023.1/geometry_conversion.html#geometry-definition-of-pyfai

[LouConreux]

Update 05/10

Reference Frame Definition

The geometry of an experiment can be defined using different frames. The first one involved is the psana frame used in metrology: the different translations x, y, z and rotation angles are measured and entered in .data file. Then those metrology data can be converted in Crystel format .geom file. At this point, those data can be chosen to be formatted in the CrystFEL format while staying in the psana frame or they could be transformed in the CrystFEL reference frame and format. Finally, pyFAI uses another system of coordinates. This update describes all three of them

CrystFEL or Laboratory Coordinate System

Following CrystFEL docs, CrystFEL's laboratory coordinate system defined as follows:
+z is the beam direction, and points along the beam (i.e. away from the source)
+y points towards the zenith (ceiling).
+x completes the right-handed coordinate system.
So, (z, y, x) is a direct reference frame, while traditional (x, y, z) is in fact indirect.

Psana Coordinate System

This is the coordinate system used in metrology measurements. According to the geometry_converter code, the psana coordinate system can be transformed into the CrystFEL coordinate system by doing those operations:

Hence, in this system:
+z is the opposite beam direction
+y points from right to left horizontally
+x points towards the ground

pyFAI Coordinate System

PyFAI's coordinate system is defined as such:

Hence, in this system:
+z is along the beam
+y points towards the zenith
+x points from left to right horizontally

Frame Transformations

To sum up, these are the three reference frames used and their respective translations:

Verification

From the psana .data file of exp mfxx49820, the beam center is shifted by X0=-600µm and Y0=-1900µm in psana
reference frame.

Thus, we would expect that the PONI being initiated at the center of the geometry would then be in pyFAI coordinates:
poni1 = 600µm = 6pix
poni2 = 1900µm = 19pix

This is indeed what we found after appropriate conversion!

Nota Bene

This frame definition took longer than expected because:

1. In reality, people at LCLS do not always define the reference frame they are using.
2. Moreover, the laboratory reference frame definition on the confluence doc is therefore wrong and needs to be updated. See Update 05/15
3. Finally, it seems that every scientist has its own preference when it comes to defining the detector orientation at the beginning of an experiment. Thus, the reference frame may remain the same, but in reality, the horizontal and vertical axis should be interchanged in consequence.

[LouConreux]

Update 05/15

Reference Frame Re-Definition

The geometry of an experiment can be defined using different frames. The first one involved is the psana frame used in metrology: the different translations x, y, z and rotation angles are measured and entered in .data file. Then those metrology data can be converted in Crystel format .geom file. At this point, those data can be chosen to be formatted in the CrystFEL format while staying in the psana frame or they could be transformed in the CrystFEL reference frame and format. Finally, pyFAI uses another system of coordinates. This update is correcting the previous one.

Psana Coordinate System

This is the coordinate system used in metrology measurements. According to the geometry_converter code, the psana coordinate system can be transformed into the CrystFEL coordinate system by doing those operations:

Hence, in this system:
+z is the opposite beam direction
+y points from right to left horizontally
+x points towards the ground

This coordinate system is convenient for calibrating the geometry since the pixel space in matplotlib definition has the same axis as psana. The top left corner of an image is the origin of the pixel space, the row-space being the vertical axis (i.e. the first dimension of an array) which in psana is the x-coordinate and the column-space being the horizontal axis (i.e. the second dimension of an array) which in psana is the y-coordinate.

CrystFEL or Laboratory Coordinate System

Following CrystFEL docs, CrystFEL's laboratory coordinate system defined as follows:
+z is the beam direction, and points along the beam (i.e. away from the source)
+y points towards the zenith (ceiling).
+x completes the right-handed coordinate system.
So, (z, y, x) is a direct reference frame, while traditional (x, y, z) is in fact indirect.

However, LCLS has chosen another reference for the writing of CrystFEL geometry file .geom.
In this system:
+z is the beam direction, and points along the beam
+y points towards the zenith (ceiling)
+x is defined going from right to left
So, (x, y, z) becomes a direct reference frame as well as (z, x, y)

pyFAI Coordinate System

PyFAI's coordinate system is defined as such:

Hence, in this system:
+z is along the beam
+y points towards the zenith
+x points from left to right horizontally

Frame Transformations

To sum up, these are the three reference frames used and their respective translations:

Verification

From the psana .data file of exp mfxx49820, the beam center is shifted by X0=-600µm and Y0=-1900µm in psana
reference frame.

Thus, we would expect that the PONI being initiated at the center of the geometry would then be in pyFAI coordinates:
poni1 = 600µm = 6pix
poni2 = -1900µm = -19pix

This is indeed what we found after appropriate conversion both from psana frame or LAB frame!

Available Format Converter Python Scripts

To go from one coordinate system to another, we have already at our disposal some conversion scripts.

geometry_convert .data to .geom and .geom to .data

This function implemented by Mikhail Dubrovin is able to convert a psana .data file into a CrystFEL format either in psana coordinate system or in LAB coordinate system. This latter conversion option is dictated by the cframe argument: cframe=0 for to psana coordinates or cframe=1 for to LAB coordinates.
This function can be called from the terminal as described here.
Besides, this function relies mainly on one function from the PSCalib repository called geometry_to_crystfel described here.
The conversion of coordinate systems is coherent with the definition of the different reference frames above.

pyFAI convert .geom to PyFAI pixel corner array

To be able to instantiate a pyFAI Detector object, one needs to provide it a flattened shape (N_mods * dim_{slow scan}, dim_{fast scan}) alongside with a pixel corner array corresponding of the (z, y, x) (respectively z along the beam, y the fast scan dimension, x the slow scan dimension) coordinates of each 4 pixel corners. This pixel corner array is then of shape (N_mods * dim_{slow scan}, dim_{fast scan}, 4, 3).

deploy_geometry

This script is particular taken from the geom_opt.py script on btx repository. From a given PSCalib.GeometryAccess object, it can generate both a .data file and a .geom file. For that, it relies on two functions:

1. save_pars_in_file
   This one simply takes the PSCalib.GeometryAccess object, reads its x0 y0 z0 rot_x rot_y rot_z tilt_x tilt_y tilt_z internal class parameters and exports them to .data file in the correct psana frame.
2. generate_geom_file
   This one is trickier. It is calling a different object called CompoundAreaCamera from the psgeom repository, in camera.py which does have an internal script from_crystfel_file or from_psana_file that allows to instantiate an ordered tree of PSCalib.GeometryObject objects (with parent/children relations). From that, there exists a to_crystfel_file function that writes a corresponding .geom file.
   However, the output of this function is not coherent with the reference frame definition (a geometry_convert .geom file and a generate_geom_file .geom file originated from the same .data file are not alike). Moreover, the axis definition in this .geom file makes no sense with respect to any of the 3 reference frames defined above.

Therefore, I would recommend not relying on this function and sticking to the well documented geometry_convert function.

GeomOpt vs PyFAI parameters

GeomOpt Optimization

Following analysis #197, it seems that the true center of the rings is:

New center is (843.5532519182825, 831.0710902924419)
Detector distance inferred from powder rings: 280.47 mm

Then, when calling deploy_geometry:

def deploy_geometry(self, outdir, pv_camera_length=None):
        """
        Write new geometry files (.geom and .data for CrystFEL and psana respectively) 
        with the optimized center and distance.
    
        Parameters
        ----------
        outdir : str
            path to output directory
        pv_camera_length : str
            PV associated with camera length
        """
        # retrieve original geometry
        run = self.diagnostics.psi.run
        geom = self.diagnostics.psi.det.geometry(run)
        top = geom.get_top_geo()
        children = top.get_list_of_children()[0]
        pixel_size = self.diagnostics.psi.get_pixel_size() * 1e3 # from mm to microns
    
        # determine and deploy shifts in x,y,z
        cy, cx = self.diagnostics.psi.det.point_indexes(run, pxy_um=(0,0), fract=True)
        dx = pixel_size * (self.center[0] - cx) # convert from pixels to microns
        dy = pixel_size * (self.center[1] - cy) # convert from pixels to microns
        dz = np.mean(-1*self.diagnostics.psi.det.coords_z(run)) - 1e3 * self.distance # convert from mm to microns
        geom.move_geo(children.oname, 0, dx=-dy, dy=-dx, dz=dz)

We can see that a PSCalib.GeometryAccess object is instantiated through the PsanaInterface object. This allows an easy manipulation of the geometry through low-level functions like geom.move_geo or geom.tilt_geo. The shift to be applied on the uncalibrated geometry is then defined as the difference in (x, y) coordinates between the true center in pixel and the current center of the detector in pixel (cy, cx) = psi.det.point_indexes(run, pxy_um=(0,0), fract=True) (this returns the pixel indexes of the point (0,0) in psana or LAB coordinates).

So starting from a 0-end.data file as such:

IP  0       CAMERA  0      -600    -1900  -300000      90      0      0    0.00000   0.00000   0.00000

We get an optimized r0008_end.data file as such:

IP  0       CAMERA  0       237    -1183  -282388       90      0      0    0.00000  0.00000  0.00000

PyFAI Optimization

Now let's draw a figure to understand the difference between looking for the optimal center of the rings and the optima PONI pyFAI coordinates:

Since the definition of PONI is relative to the geometry origin we are providing it, hence the optimal PONI found is actually the shift we need to perform on the uncalibrated geometry.
Hence, starting from a 0-end.data file as such:

IP  0       CAMERA  0      -600    -1900  -300000      90      0      0    0.00000   0.00000   0.00000

We get an optimized r0008_end.data file as such:

IP  0       CAMERA  0       188     -1234  -282388       90      0      0    0.00000  0.00000  0.00000

!!!

[LouConreux]

Update 06/09

Updated Geometry Optimization Workflow

The final endpoint of this analysis is to provide a central framework for dealing and handling geometry files for the calibration of the different LCLS expriments (mfx, mec, cxi...). This will help design a uniform geometry optimization workflow that aims to be applied for the calibration of all the different LCLS detectors (ePix10k2M, Jungfrau4M, ePix10ka...). This workflow will look like this:

Let's go into details:

Geometry Format Conversion

A geometry calibration starts with a psana .data file and is supposed to end with an optimal psana .data file. But pyFAI cannot comprehend this format, hence we need to have uniform scripts that are able to convert any geometry file format into another one, without affecting the geometry data. This will be handled by built-in functions from the PSCalib repository (.data to .geom and vice-versa) as well as some hand-coded functions (.geom to pyFAI format and vice-versa).

Optimization and Powder Preprocessing

PyFAI optimization will likely to need some preprocessing of the powder (as discussed in here: ) explaining why I am regrouping these two in one sub-category.

PyFAI provides a potential robust and fast geometry calibration workflow that was also described in here. PyFAI takes as an input a flattened pixel corner array (N_mods x dim(ss), dim(fs), 4, 3) where N_mods is the number of panels that the detector has, dim(ss) is the number of pixels in the slow-scan direction, dim(fs) the number of pixels in the fast-scan direction, 4 because a pixel has 4 corners, and 3 for the 3 coordinates. This pixel corner array is then used to instantiate a pyFAI Detector object where the optimization can be done given a powder data.

New Geometry Definition

Once the optimization is done, pyFAI will give 6 optimal parameters:
dist: distance IP-PONI where PONI is the Point of Normal Incidence.
poni1: slow-scan coordinate of PONI
poni2: fast-scan coordinate of PONI
rot1: rotation over slow-scan direction
rot2: rotation over fast-scan direction
rot3: rotation over beam direction

Nota Bene
These parameters in fact define the translations and rotations that need to be applied to correct the current detector frame.

According to pyFAI Geometry Definition, given those parameters, a pixel coordinate (x,y,z) in the detector frame can be translated into (X,Y,Z) coordinates in the chosen reference frame following this math:

So, by applying first the translations and then the rotations onto the current (x,y,z) coordinates, we can find the true (X,Y,Z) pixel coordinates. From that, we can write a script to write a CrystFEL format .geom file and finally convert that into an optimized psana .data file.

Analysis of mecl1016522 run 8

Following a small-science discussion, analysis of mecl1016522 run 8 was performed to determine the geometries of the 4 different ePix10ka detectors. LaB6 was the calibrant in this experiment.

Optimization of Quad 1

By analyzing the unassembled powder:

A first geometry guess was:

# Estimating initial geometry
poni1 = p1.mean()+525*0.0001
poni2 = p2.mean()+250*0.0001
dist = 0.14

Giving a first guess decent radial integration:

We then optimize without and with rotations and compare the resulting radial integrations:

We see little difference between those two. Indeed, pyFAI finds only small tilts in the rotation case, indicating that Quad 1 is relatively orthogonal to the beam direction.

Optimization of Quad 0

Looking at Quad 0, it seems to be slightly tilted, so we will focus on optimizing Quad 0 with rotations.

A first geometry guess was:

# Estimating initial geometry
poni1 = p1.mean()+510*0.0001
poni2 = p2.mean()+1180*0.0001
dist = 0.1329

Giving a first guess decent radial integration:

A first optimization shot was performed with the usual parameters for control points extraction:

# Extracting Control Points
sg0.extract_cp(max_rings=5, pts_per_deg=2, Imin=np.max(calib_max_flat)/100)

The resulting optimized geometry is really off as the radial integration indicates:

This can be explained by a poor quality control points extraction:

We can see that control points were only extracted on the lower left asic, leading to a bad fit:

Hence, we change the parameters of control point extraction to include more significant data in other asics:

# Extracting Control Points
sg0.extract_cp(max_rings=5, pts_per_deg=2, Imin=np.max(calib_max_flat)/400)

Optimization now leads to a way better radial integration:

Indeed, the control points were slightly more significant:

Thus, the fit is not perfect but still better than the previous one:

That's an example of why it is still unclear how to find an automated way of extracting control points. This could be done by some preprocessing of the powder beforehands (such as thresholding).

[LouConreux]

Update 07/18

Geometry Calibration Workflow Complete

Following last update, a new strategy has been developed for geometry calibration of experiments. This new strategy has been tested on mfxx49820 run 8 as well as mecl1016522 run 8. Since geometries in those experiments were rather flat, verification of the geometry calibration workflow will need some test on real tilted geometries to confirm the robustness of the workflow in the case where rotations truly matter.

Some points will be first discussed and then the overall workflow will be demonstrated in the case of the two previously mentionned experiments.

Discussions

Discussion 1: z-component, coffset, to be added or not to be added ?

When converting a psana .data file, the choice has been made to take into account the small tilt angles around x, y and z. This results in having a non-flat detector, i.e. the z-component is not a constant.

Then, when it comes to writing a corresponding crystFEL .geom file, the coffset parameter is therefore not a constant, the geometry_convert function is able to write that accordingly.

To ensure the whole loop makes sense, when writing the pyFAI corner arrays, choice has been made to center the z-component around z=0. This ensures keeping the fact that the detector is not flat, while making sure pyFAI is not troubled by an off-set in the z-component, meaning pyFAI will find the true distance between the PONI and the origin.

This far, zeroing-out the mean of the z-component has worked fine and has provided a robust way of treating different detector geometries.

Discussion 2: PONI and center of rings

While trying to solve a way of a writing new .geom files from the PONI parameters when rotations are not fixed to 0, it came to me that I needed to do some trigonometry in order to express the center of rings from the PONI parameters. Indeed, when calibrating, the geometry is sent to the center of the rings, not the PONI. Hence, when there are no rotations, center of rings = PONI, so no need to bother with that. But, when rotations come in, the PONI can be way-off the center of rings in terms of pixels (depending on the distance and rotations as we will see).

So, here is a scheme to explain the importance of differentiating PONI and center of the rings when in presence of rotations:

On this image, we can clearly see a shift between the PONI and the center of the rings, a little trigonometry then helps us to find the formula to express this center in terms of the PONI parameters:

This correction has been implemented in the function PONI_to_center inside the PyFAI_to_CrystFEL class.

Demonstrations

Experiment mfxx49820 run 8

This is what pyFAI found after optimization:

1 - Panels

2 - Radial Integration

After writing the new geom file:

1 - Panels

2 - Radial Integration

Experiment mecl1016522 run 8 Quad 2

This is what pyFAI found after optimization:

1 - Panels

2 - Radial Integration

After writing the new geom file:

1 - Panels

2 - Radial Integration

To go further

• It may need some refinement
• Discussions have to be closed in order to fix once and for all the strategy
• Optimization can now be investigated deeper

[LouConreux]

Big Update 09/12: Black-Box Optimization

Motivation

The previous project was about writing geometry conversion scripts to correct detector geometry after performing an optimization. Naturally, now we will be focusing on the optimization!

LCLS vs PyFAI

So far, LCLS tools to perform geometry conversion are able to find the beam center and the detector distance but separately. Furthermore, it is only implemented for the Silver Behenate AgBh calibrant. PyFAI offers a much more robust framework for calibrating detector geometry while being implemented for dozens of calibrant. Furthermore, pyFAI is able to not only find the beam center and the distance, pyFAI can also find the rotations around each axis to perfectly describe the detector geometry.

Description of pyFAI optimization

PyFAI optimization can be used in two cases. First, pyFAI has a calibration tool that consists in a GUI where the user can select manually the peak data that will provide ground for optimization. PyFAI is also able to recalibrate a geometry assuming the current one has been moved a bit. Hence, pyFAI is not a global optimizer without any user input. Project Control Point Extractor is actually focusing on extracting key peak data from powder to give to pyFAI in the same way of using the calibration tool. Project Black Box Optimization is focusing on using the recalibration functions coupled with exploration methods such as Bayesian Optimization. The recalibration tool mainly depends on two critical functions.

extract_cp

This function is using a Massif algorithm to extract peak data from the powder to provide ground for optimization. However, this function uses the geometry provided by the user to assign control points a ring index. Hence, this function makes the whole pyFAI optimization loop a geometry-dependent loop. One can help refining the extraction by increasing the maximum ring index until where to extract, increasing the resolution pts_per_deg or increasing the minimal intensity Imin for which a point can be considered as key data.

# Extracting control points
geom_initial = pyFAI.geometry.Geometry(dist=dist, poni1=poni1, poni2=poni2, detector=epix10k2M, wavelength=wavelength)
sg = SingleGeometry("optimization", powder, calibrant=behenate, detector=epix10k2M, geometry=geom_initial)
sg.extract_cp(max_rings=5, pts_per_deg=1, Imin=np.max(powder)/100)

refine3

Once data is collected, refine3 performs a least-squares method: pyFAI wants to minimize the difference between the diffraction angle of the extracted data and the expected angle given the calibrant data and the current geometry. Hence, refine3 is dependent on the control data points that extract_cp provides. If extract_cp outputs poor quality data, refine3 will not give a decent result.
This is where one tells pyFAI on which parameters to do the optimization. For example, here we are looking for all translations and rotations, we only fix the wavelength. Finally, refine3 outputs the residual (sum over the control points of squared differences between diffraction angles of data and calibrant) after finding the optimum.

# Continuing PyFAI optimization
score = sg.geometry_refinement.refine3(fix=["wavelength"])

Global Optimizer

Thanks to this previous discussion, we need to combine pyFAI robust optimization loop with some exploratory methods to hope finding a global optimizer. Since Paul-Emile previously worked on using Bayesian Optimization coupled with LCLS geometry tools to correct geometries, this project was all about translating the work of Paul-Emile to pyFAI. However, we still need to define a score function on which to apply the BO.

Score Function Investigation

Long Loop

Naturally, we started investigating the output of refine3 as a score. This makes one loop of Bayesian Optimization picking a geometry, extracting control points and calling refine3.

geom_initial = pyFAI.geometry.Geometry(dist=dist, poni1=poni1, poni2=poni2, detector=self.detector, wavelength=wavelength)
sg = SingleGeometry("extract_cp", powder_img, calibrant=calibrant, detector=self.detector, geometry=geom_initial)
sg.extract_cp(max_rings=5, pts_per_deg=1, Imin=Imin)
score = sg.geometry_refinement.refine3(fix=["wavelength"])

This loop has the advantage of being consistent with pyFAI recalibration tool. However, calling refine3 to find 6 parameters can take tens of seconds, hence becoming the bottleneck of time computation of this loop. Hence, we will refer this loop as being the long one.

When performing a grid search to capture if this score function will be suitable for a Bayesian Optimization, we identified pros and cons for this score.

For experiments where the geometry is more or less already centered (such as mfx), the long loop provides a nice landscape:

Here, we can clearly see that pyFAI will find the true geometry if the current one is laying within this zone of convergence. Note that this zone will shrink as the distance gets further away from the true one, and will always be centered around the beam center (PONI).

However, for experiments where no prior geometry is known (such as mec), we need to dezoom and look more global. This is where we started identify cons to this score.

It appears that we can identify the zone where the true geometry lies around (poni1=0.05, poni2=0.24). However, this zone does not longer correspond to the global minimum over all this space. Here, since we are looking at geometries really far away from the true one, pyFAI sometimes can not even find a single control point, which leads to an impossibility for computing a refine3 score. But that's fine as long as we set this to a default value such as infinity. However, remember that refine3 is dependent on the extract_cp output. To be more concise, when only a few control points are found, pyFAI computes a residual that will be extremely low due to the low number of control points. Hence, refine3 is stronlgy dependent on the number of control points that extract_cp provides as ground for optimization.

Reflecting on that, the long loop does not appear to be the ideal score function. Nonetheless, it helps drawing assumptions. If few control points are found, this is fairly because the current geometry is wrong. On the contrary, if a lot of control points are found, we hope to find this is because the geometry lies somehow close to the true one. Hopefully, we can extrapolate by assuming that the geometry where the number of CP is maximum is actually the true geometry!

Short Loop

Since applying this reasoning results in not calling refine3 anymore, this helps a lot reducing the computation time, hence the name of this subsection.
Luckily, this assumption has been verified for both mfx and mec!

Again, we may have concerns for the mec case on the ability to find the true geometry due to its very localized optimum. Since this landscape was now quite satisfactory, BO was tested and another interesting result was drawn. BO was run on mfxx49820 run 8 and run 15 to find the beam center using the default distance value stored in psana (0.3mm). True distance for run 8 was 0.282mm but for run 15 the detector was moved making the distance more close to 0.1mm. Here is the result of the BO for both case:

It turns out that if the distance is somehow close to the true one, the maximum for #CP is the true geometry. But if the distance is false, it is no longer true. Therefore, we can not hope to run BO over the 3 translations, we have to be smarter than that.

Short then Long Loop

This strategy combines the previous two by applying one BO on the short loop over different PONIs for different fixed distance. At the end of each BOs we compute refine3 of the best result, and take the result that has the lowest score. Hopefully, the BO results where the distance is wrong will give large numbers of false control points, hence large residuals. On the other hand, for distance close to the true one, both the final number of control points and the residual will be lowest.

This assumption was quickly verified: BO was run for mfxx49820 run 8 for distance 200mm, 300mm and 400mm.

WARNING:scripts.tasks:Final score: 7.701113797002203e-06
WARNING:scripts.tasks:Refined detector distance in m: 2.00e-01
WARNING:scripts.tasks:Refined detector center in m: -4.15e-03, 9.29e-04
WARNING:scripts.tasks:Total duration: 100.67898156796582 seconds

WARNING:scripts.tasks:Final score: 4.988275212597858e-08
WARNING:scripts.tasks:Refined detector distance in m: 2.82e-01
WARNING:scripts.tasks:Refined detector center in m: -7.86e-04, -6.76e-04
WARNING:scripts.tasks:Total duration: 95.27363136806525 seconds

WARNING:scripts.tasks:Final score: 1.0328984373059156e-05
WARNING:scripts.tasks:Refined detector distance in m: 3.88e-01
WARNING:scripts.tasks:Refined detector center in m: 4.71e-03, -3.44e-04
WARNING:scripts.tasks:Total duration: 104.26975293597206 seconds

Indeed, the lowest refine3 corresponds to the one having the distance closest to the true one!!