Distortion of Eiger2 CdTe detector from ID11¶
The detector’s distortion is characterized by a Tungsten plate with a grid pattern of holes provided by Gavin Vaughan. The calibration experiment has been performed by Marie Ruat with a conventional source. Those files are not available for share, please get in contact with Marie for the experimental details.
Preprocessing¶
Data were acquired with and without the grid on from of the detector. 100 frames were recorded in each case withthe same exposure time. The first step is to filter out outliers using a pixel-wise median filter on the stack.
Note: a faily recent version of FabIO is needed to read those HDF5 files produced by LImA.
[1]:
%matplotlib nbagg
import time
start_time = time.perf_counter()
#load many libraries ...
from matplotlib.pyplot import subplots
import matplotlib.colors as colors
import numpy
import pyFAI
from collections import namedtuple
from scipy.ndimage import convolve, binary_dilation, label, distance_transform_edt
from scipy.spatial import distance_matrix
from pyFAI.ext.watershed import InverseWatershed
from pyFAI.ext.bilinear import Bilinear
[2]:
!pyFAI-average -m median W200um_40kVp_5mA_T4_E8_FF_noretrigger_0000.h5 -F numpy -o flat.npy
!pyFAI-average -m median W200um_40kVp_5mA_T4_E8_GRID_noretrigger_0000.h5 -F numpy -o grid.npy
Median reduction finished
Median reduction finished
[3]:
# A couple of constants and compound dtypes ...
dt = numpy.dtype([('y', numpy.float64),
('x', numpy.float64),
('i', numpy.int64),
])
Alignment = namedtuple("Alignment", "points RMSD rotation center_ref center_set matrix")
mpl = {"origin":"lower",
"cmap":"viridis",
"interpolation":"nearest"}
print("Using pyFAI verison: ", pyFAI.version)
Using pyFAI verison: 0.20.0-beta1
[4]:
flat = numpy.load("flat.npy")
grid = numpy.load("grid.npy")
fig, ax = subplots(figsize=(8,8))
ax.imshow(grid/flat, **mpl)
ax.set_title("Raw grid (normalized)")
pass
Define the right mask¶
As we want to measure the position of the grid with a sub-pixel precision, it is of crutial important to discard all pixels wich have been interpolated.
The Eiger2 4M is built of 8 modules 500 kpixels, each of them consists of the assemble of 8 chips of 256x256 pixel each. Some pixels are systematically masked out as they are known to be noisy. Some pixels are also missing at the junction of the sub-modules. Finally the CdTe sensors is made of single crystals of the semi-conductor which cover 512x512 pixels. Hence one sensor covers half a module or 4 chips.
This is best demonstrated by the pixel-wise standard deviation along a stack of images like the one acquired for the flatfield.
A Poissonian detector should have a variance equal to the average signal. Thus plotting the standard deviation squared over the median highlights: * Noisy pixels which should be discarded for quantitative analysis std²>>median * Interpolated pixels which have only half/quater of the expected noise (std²<<median).
The detector has also an internal map of invalid pixel which are set to the maximum value of the range.
[5]:
!pyFAI-average -m std W200um_40kVp_5mA_T4_E8_FF_noretrigger_0000.h5 -F numpy -o flat_std.npy
Std reduction finished
[6]:
std = numpy.load("flat_std.npy")
fig, ax = subplots(figsize=(8,8))
ax.imshow((std**2/flat).clip(0,2), **mpl)
ax.set_title("Poissonianness map of the detctor (should be 1)")
pass
[7]:
fig, ax = subplots()
ax.hist((std**2/flat).ravel(), 100, range=(0,2))
ax.set_title("Histogram of std²/median values")
pass
This test of Poissonian-ness is not enough to discriminate between interpolated pixels and non interpolated ones. One needs to build the mask by other methods.
[8]:
# This is the default detector as definied in pyFAI according to the specification provided by Dectris:
eiger2 = pyFAI.detector_factory("Eiger2CdTe_4M")
width = eiger2.shape[1]
module_size = eiger2.MODULE_SIZE
module_gap = eiger2.MODULE_GAP
submodule_size = (256,256)
[9]:
#Calculate the default mask
mask = eiger2.calc_mask()
# Mask out the interpolated pixels along X
for j in [256, 772]:
for i in range(j, eiger2.max_shape[1],
eiger2.module_size[1] + eiger2.MODULE_GAP[1]):
mask[:,i-2:i + 2] = 1
# Mask out the interpolated pixels along Y
for i in range(256, eiger2.max_shape[0],
eiger2.module_size[0] + eiger2.MODULE_GAP[0]):
mask[i-2:i + 2, :] = 1
# Finally mask out invalid/miss-behaving pixels known from the detector
mask[ flat>(flat.max()*0.99) ] = 1
[10]:
pois = (std**2/flat).clip(0,2)
pois[numpy.where(mask)] = numpy.NAN
fig, ax = subplots(figsize=(8,8))
ax.imshow(pois, **mpl)
ax.set_title("Histogram of std²/median values (masked)")
pass
[11]:
normalized = grid/flat
normalized[numpy.where(mask)] = numpy.NaN
fig, ax = subplots(figsize=(8,8))
ax.imshow(normalized, **mpl)
ax.set_ylim(1100,1400)
ax.set_xlim(1000,1300)
ax.set_title("Zoom on one chip")
pass
We can see we have between 8 and 16 grid spots per sub-module. This is enough as 3 are needed to localize precisely the position of the chip.
Grid spot position measurement.¶
Let’s measure the position of the grid spots precisely. For this we perform a convolution with a kernel which looks like the spot itself. By zooming onto one spot, it is roughly 10 pixels wide, the kernel needs to be of odd size. The fade-out function is tuned to set precisely the spot diameter.
The masked values have been set to NaN, this ensures any spot close to a masked region get discarded automatically.
[12]:
#Definition of the convolution kernel
ksize = 11
y,x = numpy.ogrid[-(ksize-1)//2:ksize//2+1,-(ksize-1)//2:ksize//2+1]
d = numpy.sqrt(y*y+x*x)
#Fade out curve definition
fadeout = lambda x: 1/(1+numpy.exp(5*(x-5)))
kernel = fadeout(d)
mini=kernel.sum()
print("Integral of the kernel: ", mini)
fig,ax = subplots(1,3, figsize=(8,4))
ax[0].imshow(d)
ax[0].set_title("Distance array")
ax[1].plot(numpy.linspace(0,8,100),fadeout(numpy.linspace(0,8,100)))
ax[1].set_title("fade-out curve")
ax[2].imshow(kernel)
ax[2].set_title("Convolution kernel")
pass
Integral of the kernel: 78.56685275622786
[13]:
smooth = convolve(normalized, kernel, mode="constant", cval=0)/mini
fig,ax = subplots(1,2, figsize=(8,4))
ax[0].imshow(normalized)
ax[0].set_ylim(1100,1400)
ax[0].set_xlim(1000,1300)
ax[1].imshow(smooth)
ax[1].set_ylim(1100,1400)
ax[1].set_xlim(1000,1300)
ax[0].set_title("Grid (zoomed)")
ax[1].set_title("Smoothed grid")
pass
[14]:
# Calculate a mask with all pixels close to any gap is discared
#big_mask = numpy.isnan(my_smooth)
#big_mask = binary_dilation(mask, iterations=ksize//2+2) #even a bit larger
big_mask = binary_dilation(mask, iterations=ksize) #Extremely conservative !
Peak position measurement¶
The center of spot is now easily measured by segmenting out the image
[15]:
iw = InverseWatershed(smooth)
iw.init()
iw.merge_singleton()
all_regions = set(iw.regions.values())
regions = [i for i in all_regions if i.size>mini]
print("Number of region segmented: %s"%len(all_regions))
print("Number of large enough regions : %s"%len(regions))
Number of region segmented: 661404
Number of large enough regions : 4043
[16]:
# Histogram of peak height:
s = numpy.array([i.maxi for i in regions])
fig, ax = subplots()
ax.hist(s, 100)
[16]:
(array([3.047e+03, 5.000e+00, 0.000e+00, 2.000e+00, 0.000e+00, 1.000e+00,
0.000e+00, 1.000e+00, 0.000e+00, 1.000e+00, 0.000e+00, 2.000e+00,
0.000e+00, 1.000e+00, 0.000e+00, 0.000e+00, 1.000e+00, 0.000e+00,
1.000e+00, 1.000e+00, 1.000e+00, 1.000e+00, 1.000e+00, 3.000e+00,
1.000e+00, 1.000e+00, 3.000e+00, 1.000e+00, 3.000e+00, 3.000e+00,
1.000e+00, 2.000e+00, 0.000e+00, 2.000e+00, 1.000e+00, 4.000e+00,
2.000e+00, 2.000e+00, 0.000e+00, 2.000e+00, 1.000e+00, 0.000e+00,
1.000e+00, 1.000e+00, 1.000e+00, 2.000e+00, 1.000e+00, 1.000e+00,
1.000e+00, 2.000e+00, 0.000e+00, 2.000e+00, 0.000e+00, 2.000e+00,
0.000e+00, 2.000e+00, 0.000e+00, 2.000e+00, 0.000e+00, 2.000e+00,
1.000e+00, 1.000e+00, 3.000e+00, 0.000e+00, 1.000e+00, 2.000e+00,
0.000e+00, 2.000e+00, 2.000e+00, 0.000e+00, 1.000e+00, 2.000e+00,
2.000e+00, 1.000e+00, 2.000e+00, 4.000e+00, 4.000e+00, 1.500e+01,
1.300e+01, 1.900e+01, 1.800e+01, 2.300e+01, 1.900e+01, 1.900e+01,
1.800e+01, 3.400e+01, 4.800e+01, 5.200e+01, 5.300e+01, 5.600e+01,
6.700e+01, 7.700e+01, 5.600e+01, 8.700e+01, 6.800e+01, 4.900e+01,
4.600e+01, 4.200e+01, 1.100e+01, 7.000e+00]),
array([0.00389302, 0.01098051, 0.01806801, 0.0251555 , 0.03224299,
0.03933049, 0.04641798, 0.05350547, 0.06059297, 0.06768046,
0.07476795, 0.08185545, 0.08894294, 0.09603043, 0.10311793,
0.11020542, 0.11729291, 0.12438041, 0.1314679 , 0.1385554 ,
0.14564289, 0.15273038, 0.15981788, 0.16690537, 0.17399286,
0.18108036, 0.18816785, 0.19525534, 0.20234284, 0.20943033,
0.21651782, 0.22360532, 0.23069281, 0.2377803 , 0.2448678 ,
0.25195529, 0.25904279, 0.26613028, 0.27321777, 0.28030527,
0.28739276, 0.29448025, 0.30156775, 0.30865524, 0.31574273,
0.32283023, 0.32991772, 0.33700521, 0.34409271, 0.3511802 ,
0.35826769, 0.36535519, 0.37244268, 0.37953018, 0.38661767,
0.39370516, 0.40079266, 0.40788015, 0.41496764, 0.42205514,
0.42914263, 0.43623012, 0.44331762, 0.45040511, 0.4574926 ,
0.4645801 , 0.47166759, 0.47875508, 0.48584258, 0.49293007,
0.50001757, 0.50710506, 0.51419255, 0.52128005, 0.52836754,
0.53545503, 0.54254253, 0.54963002, 0.55671751, 0.56380501,
0.5708925 , 0.57797999, 0.58506749, 0.59215498, 0.59924248,
0.60632997, 0.61341746, 0.62050496, 0.62759245, 0.63467994,
0.64176744, 0.64885493, 0.65594242, 0.66302992, 0.67011741,
0.6772049 , 0.6842924 , 0.69137989, 0.69846738, 0.70555488,
0.71264237]),
<BarContainer object of 100 artists>)
[17]:
#Sieve-out for peak intensity
int_mini = 0.5
peaks = [(i.index//width, i.index%width) for i in regions if (i.maxi)>int_mini and
not big_mask[(i.index//width, i.index%width)]]
print("Number of remaining peaks with I>%s: %s"%(int_mini, len(peaks)))
peaks_raw = numpy.array(peaks)
Number of remaining peaks with I>0.5: 804
About 800 spots were found as valid out of a maximum of 1024 (64 chips with 16 spots per chip)
Those spot positions are interpolated using a second order taylor expansion in the region around the maximum value of the smoothed image.
[18]:
#Use a bilinear interpolator to localize/refine the maxima
bl = Bilinear(smooth)
ref_peaks = [bl.local_maxi(p) for p in peaks]
[19]:
#Overlay raw peak coordinate and refined peak positions
fig, ax = subplots(figsize=(8,8))
ax.imshow(normalized, **mpl)
peaks_ref = numpy.array(ref_peaks)
ax.plot(peaks_raw[:,1], peaks_raw[:, 0], ".r")
ax.plot(peaks_ref[:,1],peaks_ref[:, 0], ".b")
ax.set_title("Extracted peak position (red: raw, blue: refined)")
print("Refined peak coordinate:", ref_peaks[:10], "...")
Refined peak coordinate: [(1849.8662119656801, 988.3285617828369), (1780.0930171236396, 114.93192630261183), (1781.5487900078297, 652.5473097264767), (1781.8086816519499, 719.7614358514547), (1782.0510528832674, 786.9847892848775), (1782.2430633455515, 854.1485713273287), (1782.4817058444023, 921.467636436224), (1782.6218681633472, 988.386855006218), (1783.066966727376, 1055.5694026947021), (1783.2517661750317, 1122.8249076157808)] ...
At this stage we have about 800 peaks (with sub-pixel precision) which are visually distributed on all modules and on all chips. We could have expected 16*64=1024 hence most of the spots were properly located.
Let’s assign each peak to a module identifier. This allows to print out the number of peaks per module:
[20]:
#Module identification
mid, cnt = label(numpy.isfinite(normalized), structure=numpy.ones((3,3), dtype=int))
print(cnt, "chips have been labeled")
64 chips have been labeled
[21]:
# Fill the gaps in module identification array
#From http://stackoverflow.com/questions/3662361/fill-in-missing-values-with-nearest-neighbour-in-python-numpy-masked-arrays
def fill(data, invalid=None):
"""
Replace the value of invalid 'data' cells (indicated by 'invalid')
by the value of the nearest valid data cell
Input:
data: numpy array of any dimension
invalid: a binary array of same shape as 'data'. True cells set where data
value should be replaced.
If None (default), use: invalid = np.isnan(data)
Output:
Return a filled array.
"""
if invalid is None:
invalid = numpy.isnan(data)
ind = distance_transform_edt(invalid, return_distances=False, return_indices=True)
return data[tuple(ind)]
filled_mid = fill(mid, invalid=mid==0)
[22]:
fig,ax = subplots(1, 2, figsize=(10,5))
ax[0].imshow(mid, **mpl)
ax[0].set_title("Chip identification number")
ax[1].imshow(filled_mid, **mpl)
ax[1].set_title("Filled-gaps version")
pass
[23]:
yxi = numpy.array([i+(mid[round(i[0]),round(i[1])],)
for i in ref_peaks], dtype=dt)
print("Number of keypoint per chip:")
for i in range(1, cnt+1):
print(f"Chip id: {i:2d} \t Number of spots: {(yxi[:]['i'] == i).sum():2d}")
Number of keypoint per chip:
Chip id: 1 Number of spots: 12
Chip id: 2 Number of spots: 12
Chip id: 3 Number of spots: 12
Chip id: 4 Number of spots: 16
Chip id: 5 Number of spots: 16
Chip id: 6 Number of spots: 16
Chip id: 7 Number of spots: 12
Chip id: 8 Number of spots: 12
Chip id: 9 Number of spots: 9
Chip id: 10 Number of spots: 9
Chip id: 11 Number of spots: 9
Chip id: 12 Number of spots: 12
Chip id: 13 Number of spots: 12
Chip id: 14 Number of spots: 11
Chip id: 15 Number of spots: 9
Chip id: 16 Number of spots: 9
Chip id: 17 Number of spots: 12
Chip id: 18 Number of spots: 12
Chip id: 19 Number of spots: 12
Chip id: 20 Number of spots: 16
Chip id: 21 Number of spots: 15
Chip id: 22 Number of spots: 16
Chip id: 23 Number of spots: 15
Chip id: 24 Number of spots: 12
Chip id: 25 Number of spots: 12
Chip id: 26 Number of spots: 12
Chip id: 27 Number of spots: 12
Chip id: 28 Number of spots: 16
Chip id: 29 Number of spots: 16
Chip id: 30 Number of spots: 16
Chip id: 31 Number of spots: 16
Chip id: 32 Number of spots: 12
Chip id: 33 Number of spots: 8
Chip id: 34 Number of spots: 9
Chip id: 35 Number of spots: 8
Chip id: 36 Number of spots: 14
Chip id: 37 Number of spots: 16
Chip id: 38 Number of spots: 16
Chip id: 39 Number of spots: 16
Chip id: 40 Number of spots: 15
Chip id: 41 Number of spots: 10
Chip id: 42 Number of spots: 11
Chip id: 43 Number of spots: 12
Chip id: 44 Number of spots: 16
Chip id: 45 Number of spots: 16
Chip id: 46 Number of spots: 16
Chip id: 47 Number of spots: 15
Chip id: 48 Number of spots: 16
Chip id: 49 Number of spots: 9
Chip id: 50 Number of spots: 9
Chip id: 51 Number of spots: 9
Chip id: 52 Number of spots: 12
Chip id: 53 Number of spots: 12
Chip id: 54 Number of spots: 11
Chip id: 55 Number of spots: 12
Chip id: 56 Number of spots: 12
Chip id: 57 Number of spots: 9
Chip id: 58 Number of spots: 9
Chip id: 59 Number of spots: 9
Chip id: 60 Number of spots: 12
Chip id: 61 Number of spots: 12
Chip id: 62 Number of spots: 12
Chip id: 63 Number of spots: 15
Chip id: 64 Number of spots: 16
Grid assignment¶
The calibration is performed using a regular grid, the idea is to assign to each peak of coordinates (x,y) the integer value (X, Y) which correspond to the grid corrdinate system.
The first step is to measure the grid pitch which correspond to the distance (in pixels) from one peak to the next. This is easily obtained from a pair-wise distribution function.
[24]:
# pairwise distance calculation using scipy.spatial.distance_matrix
dist = distance_matrix(peaks_ref, peaks_ref)
fig, ax = subplots()
ax.hist(dist.ravel(), 200, range=(0,200))
ax.set_title("Pair-wise distribution function")
ax.set_xlabel("Distance (pixel)")
ax.set_ylabel("Occurances")
pass
The histogram of the pair-distribution function has a first peak at 0 and the second peak between 66 and 67. The value of the step size is taken as the average of this second peak in the histogram as it correspond to the first neighbour distance.
Two other parameters correspond to the offset, in pixels, for the grid index (X,Y) = (0,0). The easiest is to measure the smallest x and y for the first chip.
The grid looks pretty well aligned with the detector, so we will not take into account the rotation of the grid with the detector.
[25]:
#from pair-wise distribution histogram
valid_pairs = dist[numpy.logical_and(60<dist, dist<70)]
step = valid_pairs.mean()
print("Average step size", step, "±", valid_pairs.std(), " pixels, condidering ", valid_pairs.size, "paires")
Average step size 67.20312873609967 ± 0.12603174907696846 pixels, condidering 2770 paires
[26]:
#Refinement of the step size when considering only intra-chip distances.
distances = []
for i in range(1, cnt+1):
locale = yxi[yxi[:]["i"] == i]
xy = tmp = numpy.vstack((locale["x"], locale["y"])).T
ldist = distance_matrix(tmp, tmp)
distances.append(ldist[numpy.logical_and(60<ldist, ldist<70)])
valid_pairs = numpy.concatenate(distances)
step = valid_pairs.mean()
print("Average step size", step, "±", valid_pairs.std(), " pixels, condidering ", valid_pairs.size, "paires")
Average step size 67.20212525166619 ± 0.11948927150151667 pixels, condidering 2290 paires
It is interesting to note the standard deviation (0.12) corresponds to the precision of our measurement.
Any tweaking when defining the kernel prior to the convolution should be checked against this error.
Finally, it will not be possible to calibrate with a precision better than 0.12 pixel.
[27]:
#work with the first module and fit the peak positions
first = yxi[yxi[:]["i"] == 1]
y_min = first[:]["y"].min()
x_min = first[:]["x"].min()
print("offset for the first spot: ", x_min, y_min)
offset for the first spot: 51.99926513293758 32.647715866565704
The grid looks very well aligned with the axes which makes this step easier but nothing garanties it is perfect, so the rotation of the grid has to be measured as well. We will use Kabsch’s algoritm for this: https://en.wikipedia.org/wiki/Kabsch_algorithm
[28]:
def kabsch(P, R):
"Align P on R"
P = numpy.ascontiguousarray(P)
R = numpy.ascontiguousarray(R)
size = P.shape[0]
assert R.shape[0] == size
centroid_P = P.mean(axis=0)
centroid_R = R.mean(axis=0)
centered_P = P - centroid_P
centered_R = R - centroid_R
C = numpy.dot(centered_P.T, centered_R)
V, S, W = numpy.linalg.svd(C)
d = (numpy.linalg.det(V) * numpy.linalg.det(W)) < 0.0
if d:
S[-1] = -S[-1]
V[:, -1] = -V[:, -1]
# Create Rotation matrix U
U = numpy.dot(V, W)
rotated_P = numpy.dot(centered_P, U)
aligned_P = rotated_P + centroid_R
rmsd = numpy.sqrt(((centered_R - rotated_P)**2).sum()/size)
angle = 180*numpy.arctan2(U[1,0],U[0,0])/numpy.pi
return Alignment(aligned_P, rmsd, angle, centroid_R, centroid_P, U)
reference_1 = numpy.empty((first.size, 2))
measured_1 = numpy.empty((first.size, 2))
measured_1[:, 0] = first[:]["x"]
measured_1[:, 1] = first[:]["y"]
reference_1[:, 0] = numpy.round((first[:]["x"]-x_min)/step)*step+x_min
reference_1[:, 1] = numpy.round((first[:]["y"]-y_min)/step)*step+y_min
%time kabsch(measured_1, reference_1)
CPU times: user 443 µs, sys: 30 µs, total: 473 µs
Wall time: 447 µs
[28]:
Alignment(points=array([[ 52.21972564, 99.72317216],
[119.16605393, 99.82767768],
[186.46198903, 99.9937887 ],
[ 51.78194341, 32.47353632],
[119.17009389, 32.68479293],
[186.48768901, 32.82819913],
[119.2490204 , 234.25738509],
[ 51.92672117, 234.31193451],
[186.78477529, 234.17332045],
[ 51.96370067, 167.07745319],
[119.12344307, 167.02715893],
[186.0815291 , 167.03242584]]), RMSD=0.2011512982256994, rotation=0.1268046482418138, center_ref=array([119.20139038, 133.45090374]), center_set=array([119.49699086, 133.7740458 ]), matrix=array([[ 0.99999755, -0.00221316],
[ 0.00221316, 0.99999755]]))
[29]:
# Print alignment info for all chips:
kabsch_results = {}
raw_distances = []
distances = []
for i in range(1, cnt+1):
local = yxi[yxi[:]["i"] == i]
reference = numpy.empty((local.size, 2))
measured = numpy.empty((local.size, 2))
measured[:, 0] = local[:]["x"]
measured[:, 1] = local[:]["y"]
reference[:, 0] = numpy.round((local[:]["x"]-x_min)/step)*step+x_min
reference[:, 1] = numpy.round((local[:]["y"]-y_min)/step)*step+y_min
raw_distances.append(numpy.sqrt(((reference-measured)**2).sum(axis=-1)))
res = kabsch_results[i] = kabsch(measured, reference)
print(f"Chip: {i:02d} \t RMSD: {res.RMSD:.4f} \t Angle: {res.rotation:.4f}°\t Displacement: {res.center_ref-res.center_set}")
distances.append(numpy.sqrt(((reference-res.points)**2).sum(axis=-1)))
Chip: 01 RMSD: 0.2012 Angle: 0.1268° Displacement: [-0.29560048 -0.32314206]
Chip: 02 RMSD: 0.1510 Angle: 0.1656° Displacement: [-0.38391344 -1.22423316]
Chip: 03 RMSD: 0.1083 Angle: 0.2022° Displacement: [-0.23829583 -2.22024286]
Chip: 04 RMSD: 0.1306 Angle: 0.1836° Displacement: [-0.21810019 -2.90728647]
Chip: 05 RMSD: 0.1365 Angle: 0.1232° Displacement: [ 0.03213908 -3.53275582]
Chip: 06 RMSD: 0.1334 Angle: 0.1232° Displacement: [-0.06032228 -4.13097843]
Chip: 07 RMSD: 0.1475 Angle: 0.1536° Displacement: [-0.17562643 -4.64086857]
Chip: 08 RMSD: 0.1410 Angle: 0.1508° Displacement: [-0.34166363 -5.26099384]
Chip: 09 RMSD: 0.1927 Angle: 0.2098° Displacement: [ 0.49720431 -0.35965719]
Chip: 10 RMSD: 0.2274 Angle: 0.1516° Displacement: [ 0.43206139 -1.2323675 ]
Chip: 11 RMSD: 0.1829 Angle: 0.1413° Displacement: [ 0.48729789 -2.22309062]
Chip: 12 RMSD: 0.1436 Angle: 0.1747° Displacement: [ 0.50151207 -2.97743933]
Chip: 13 RMSD: 0.1093 Angle: 0.1599° Displacement: [ 0.68731817 -3.60915833]
Chip: 14 RMSD: 0.1450 Angle: 0.1404° Displacement: [ 0.51759231 -4.22358055]
Chip: 15 RMSD: 0.1604 Angle: 0.1154° Displacement: [ 0.50287889 -4.72060995]
Chip: 16 RMSD: 0.1761 Angle: 0.1194° Displacement: [ 0.39368985 -5.35868294]
Chip: 17 RMSD: 0.1366 Angle: 0.1373° Displacement: [ 1.00461638 -0.86454067]
Chip: 18 RMSD: 0.1178 Angle: 0.1796° Displacement: [ 1.05850903 -1.54616142]
Chip: 19 RMSD: 0.1333 Angle: 0.1285° Displacement: [ 1.57215237 -2.12885442]
Chip: 20 RMSD: 0.0857 Angle: 0.1542° Displacement: [ 1.66369684 -2.71096796]
Chip: 21 RMSD: 0.1137 Angle: 0.0904° Displacement: [ 1.47368641 -3.4015823 ]
Chip: 22 RMSD: 0.1381 Angle: 0.0841° Displacement: [ 1.41603634 -3.75261809]
Chip: 23 RMSD: 0.1437 Angle: 0.1674° Displacement: [ 1.10109899 -4.73906987]
Chip: 24 RMSD: 0.1398 Angle: 0.1837° Displacement: [ 0.97095115 -5.37045418]
Chip: 25 RMSD: 0.1166 Angle: 0.1324° Displacement: [ 1.68764075 -0.85490263]
Chip: 26 RMSD: 0.1177 Angle: 0.1185° Displacement: [ 1.69238379 -1.52111382]
Chip: 27 RMSD: 0.1055 Angle: 0.1307° Displacement: [ 2.26789772 -2.11232468]
Chip: 28 RMSD: 0.0909 Angle: 0.1354° Displacement: [ 2.29945194 -2.70842143]
Chip: 29 RMSD: 0.1050 Angle: 0.0733° Displacement: [ 1.88217578 -3.4212234 ]
Chip: 30 RMSD: 0.0839 Angle: 0.0652° Displacement: [ 1.80781409 -3.80468263]
Chip: 31 RMSD: 0.1050 Angle: 0.1635° Displacement: [ 1.84852818 -4.82808066]
Chip: 32 RMSD: 0.1269 Angle: 0.1861° Displacement: [ 1.73064256 -5.46682588]
Chip: 33 RMSD: 0.0862 Angle: 0.1965° Displacement: [ 2.37680974 -0.34603374]
Chip: 34 RMSD: 0.1233 Angle: 0.1742° Displacement: [ 2.38897066 -1.13349797]
Chip: 35 RMSD: 0.0483 Angle: 0.1682° Displacement: [ 2.80260071 -1.87052911]
Chip: 36 RMSD: 0.0919 Angle: 0.1821° Displacement: [ 2.83955264 -2.65322456]
Chip: 37 RMSD: 0.0837 Angle: 0.1706° Displacement: [ 2.82376951 -3.21370844]
Chip: 38 RMSD: 0.0725 Angle: 0.1691° Displacement: [ 2.82901261 -4.04532675]
Chip: 39 RMSD: 0.0955 Angle: 0.1839° Displacement: [ 2.81983799 -4.89687442]
Chip: 40 RMSD: 0.2603 Angle: 0.2233° Displacement: [ 2.80277945 -5.75942135]
Chip: 41 RMSD: 0.0663 Angle: 0.1500° Displacement: [ 2.94866733 -0.24353623]
Chip: 42 RMSD: 0.0657 Angle: 0.1755° Displacement: [ 3.08040408 -1.08261991]
Chip: 43 RMSD: 0.0795 Angle: 0.1534° Displacement: [ 3.54815692 -1.83588844]
Chip: 44 RMSD: 0.0957 Angle: 0.1752° Displacement: [ 3.63888685 -2.56507256]
Chip: 45 RMSD: 0.0908 Angle: 0.1744° Displacement: [ 3.68243121 -3.18435017]
Chip: 46 RMSD: 0.0824 Angle: 0.1811° Displacement: [ 3.67889941 -4.0305737 ]
Chip: 47 RMSD: 0.0923 Angle: 0.2053° Displacement: [ 3.7640197 -4.89743891]
Chip: 48 RMSD: 0.1995 Angle: 0.1995° Displacement: [ 3.7140228 -5.81118489]
Chip: 49 RMSD: 0.0507 Angle: 0.1030° Displacement: [ 4.24836704 -0.15775402]
Chip: 50 RMSD: 0.1115 Angle: 0.1149° Displacement: [ 4.34989515 -0.55983273]
Chip: 51 RMSD: 0.1024 Angle: 0.1895° Displacement: [ 4.1699786 -1.6241195]
Chip: 52 RMSD: 0.0841 Angle: 0.1909° Displacement: [ 4.2606929 -2.44430756]
Chip: 53 RMSD: 0.0557 Angle: 0.1228° Displacement: [ 4.43171405 -3.41100148]
Chip: 54 RMSD: 0.0763 Angle: 0.1238° Displacement: [ 4.43648478 -3.94109092]
Chip: 55 RMSD: 0.0687 Angle: 0.2055° Displacement: [ 4.31061586 -4.46464599]
Chip: 56 RMSD: 0.1348 Angle: 0.1873° Displacement: [ 4.24341779 -5.29907885]
Chip: 57 RMSD: 0.1909 Angle: 0.0670° Displacement: [4.59114254 0.05261615]
Chip: 58 RMSD: 0.0919 Angle: 0.0797° Displacement: [ 4.70254294 -0.42917319]
Chip: 59 RMSD: 0.0816 Angle: 0.2176° Displacement: [ 5.07864071 -1.50543293]
Chip: 60 RMSD: 0.0880 Angle: 0.1947° Displacement: [ 5.12083489 -2.29531844]
Chip: 61 RMSD: 0.0921 Angle: 0.1047° Displacement: [ 4.89500604 -3.27207393]
Chip: 62 RMSD: 0.0813 Angle: 0.0870° Displacement: [ 4.96782182 -3.80671605]
Chip: 63 RMSD: 0.0884 Angle: 0.1909° Displacement: [ 5.07659518 -4.33725322]
Chip: 64 RMSD: 0.0915 Angle: 0.1956° Displacement: [ 5.02915659 -5.17905991]
[30]:
ndistances = numpy.concatenate(distances)
nraw_distances = numpy.concatenate(raw_distances)
fig, ax = subplots()
ax.hist(nraw_distances, bins=100, label="raw positions")
ax.hist(ndistances, bins=100, label="refined positions")
ax.set_title("Histogram of displacement")
ax.set_xlabel("Distance in pixels")
ax.set_ylabel("Number of spots")
ax.legend()
pass
Reconstruction of the pixel position¶
The pixel position can be obtained from the standard Pilatus detector. Each module is then displaced according to the fitted values.
[31]:
%%time
pixel_coord = pyFAI.detector_factory("Eiger2CdTe_4M").get_pixel_corners().astype(numpy.float64)
pixel_coord_raw = pixel_coord.copy()
for module in range(1, cnt+1):
# Extract the pixel corners for one module
module_idx = numpy.where(filled_mid == module)
one_module = pixel_coord_raw[module_idx]
#retrieve the fitted values
res = kabsch_results[module]
#z = one_module[...,0]
y = one_module[...,1].ravel()/eiger2.pixel1
x = one_module[...,2].ravel()/eiger2.pixel2
xy_initial = numpy.vstack((x, y)).T
xy_centered = xy_initial - res.center_set
xy_rotated = numpy.dot(xy_centered, res.matrix)
xy_aligned = xy_rotated + res.center_ref
one_module[...,1] = (xy_aligned[:,1] * eiger2.pixel1).reshape(one_module.shape[:-1])
one_module[...,2] = (xy_aligned[:,0] * eiger2.pixel2).reshape(one_module.shape[:-1])
#Update the array
pixel_coord[module_idx] = one_module
CPU times: user 1.47 s, sys: 116 ms, total: 1.58 s
Wall time: 1.58 s
[32]:
# displacement for every pixel corner (before/after global displacement):
displ_refined = numpy.sqrt(((pixel_coord - pixel_coord_raw)**2).sum(axis=-1))/eiger2.pixel1
res = kabsch(pixel_coord.reshape((-1, 3)), pixel_coord_raw.reshape((-1, 3)))
new_pixel_coord = res.points.reshape(pixel_coord_raw.shape)
displ_aligned = numpy.sqrt(((new_pixel_coord - pixel_coord_raw)**2).sum(axis=-1))/eiger2.pixel1
[33]:
fig, ax = subplots()
ax.hist(displ_refined.ravel(), 100, label="refined before alignment")
ax.hist(displ_aligned.ravel(), 100, label="refined after alignement")
ax.set_title("Displacement of pixel corners versus a regular pixel layout")
ax.set_xlabel("Distance in pixels")
ax.set_ylabel("Number of corners of pixels")
ax.legend()
pass
Validation of the distortion¶
To validate the new pixel layout, we can use the new grid to calculate the spot postition in space and look how well aligned they are.
First we build a function which performes the bilinear interpolation of any detector coordinate (return a 3D position). This function is then used to calculate the position for the original grid and for the corrected grid.
As previouly, all spot distances are calculated and histogrammed. The standard deviation is used to evaluate how much was gained.
[34]:
def intepolate_3d(yx, coord=pixel_coord_raw):
y,x = yx
X = int(x)
Y = int(y)
pixel = coord[Y,X]
#print(pixel)
dx = x - X
dy = y - Y
res = pixel[0]*(1.0-dx)*(1.0-dy)+\
pixel[3]*dx*(1.0-dy)+\
pixel[1]*(1.0-dx)*dy+\
pixel[2]*dx*dy
return res
intepolate_3d((0.99,0.01))
[34]:
array([0.00000000e+00, 7.42500035e-05, 7.50000036e-07])
[35]:
raw_pixel_64 = pixel_coord_raw.astype("float64")/eiger2.pixel1
new_pixel_64 = new_pixel_coord.astype("float64")/eiger2.pixel1
spot3d_raw = numpy.array([intepolate_3d(i, coord=raw_pixel_64) for i in ref_peaks])
spot3d_ref = numpy.array([intepolate_3d(i, coord=new_pixel_64) for i in ref_peaks])
dist_raw = distance_matrix(spot3d_raw, spot3d_raw)
valid_raw = dist_raw[numpy.logical_and(65<dist_raw, dist_raw<70)]
dist_ref = distance_matrix(spot3d_ref, spot3d_ref)
valid_ref = dist_ref[numpy.logical_and(65<dist_ref, dist_ref<70)]
[36]:
fig,ax = subplots()
h_raw = numpy.histogram(valid_raw, 100, range=(66, 68))
h_ref = numpy.histogram(valid_ref, 100, range=(66, 68))
x = 0.5*(h_raw[1][1:]+h_raw[1][:-1])
ax.plot(x, h_raw[0], label="raw")
ax.plot(x, h_ref[0], label="ref", alpha=0.8)
ax.legend()
ax.set_title("Spot distance")
ax.set_xlabel("Distance in pixels")
ax.set_ylabel("occurance")
pass
[37]:
print(f"Distance before: {valid_raw.mean():.4f} after {valid_ref.mean():.4f},\nDeviation before {valid_raw.std():.4f}, after {valid_ref.std():.4f}")
Distance before: 67.2031 after 67.2000,
Deviation before 0.1260, after 0.1193
[38]:
#Saving of the result as a distortion file usable in pyFAI
eiger2.set_pixel_corners(new_pixel_coord.astype(numpy.float32))
eiger2.mask = eiger2.calc_mask() + flat>(flat.max()*0.9)
eiger2.save("Eiger2CdTe_4M_ID11.h5")
Conclusion¶
The distortion measured on the Eiger2 CdTe 4M detector for ID11 looks so small that it hardly deserves correction for module misalignment.
[39]:
print(f"Total execution time: {time.perf_counter()-start_time:.3f}s")
Total execution time: 43.484s