LouConreux
commented
4 months ago 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:
-
Geometry conversion: transforming either a .geom or .data file in a pyFAI-friendly geometry instance.
-
The pyFAI optimizer that takes in input keV calibrated images and the detector geometry to optimize.
-
Transforming back the final pyFAI-ed geometry in the LCLS reference frame coordinates.
-
Deploying the calibrated geometry.
-
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:
- 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: _maxrings : defines the interval of ring index in which control points must be extracted _pts_perdeg : 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)-
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) -
The optimization relies on the scipy.optimize function called
fmin_slsqpsource 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
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
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:
Load Data
Load the calibrated images:
Here, we need to reshape the _calibavg (in (Np, Nx, Ny)) to the shape supported by PyFAI *(NpNx, Ny)**
Set Experiment Initial Conditions
Define the guessed dist, poni1, poni2 and the wavelength used.
Calling the Geometry fitting
Here, we voluntarily defined parameters to be way off the true center. We get:
Instantiating the PyFAI Calibrant
Here, this was Silver Behenate:
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.
Results
For n=5 loops, we obtained the following convergence:
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.
Next Steps
Finding why dist cannot be optimized ?
Finding best hyperparameters for _extractcp
In this optimization, the most crucial function is _extractcp which automatically extract control points on rings. We can play on the number of points we want on 3 hyperparameters: _maxrings selects the furthest rings where to extract control points _pts_perdeg 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).