Calibration of the pixel position for a Pilatus detector

This tutorial summarizes the work done by Frederic Sulzman during his internship at ESRF during the summer 2015 entilted “Calibration for geometric distortion in multi- modules pixel detectors”.

The overall strategy is very similar to “CCD calibration” tutorial with some specificities due to the modular nature of the detector.

0. Image preprocessing

1. Peak picking

2. Grid assignment

3. Displacement fitting

4. Reconstruction of the pixel position

5. Saving into a detector definition file

Each module being made by lithographic processes, the error within a module will be assumeed to be constant. We will use the name “displacement of the module” to describe the rigide movement of the module.

This tutorial uses data from the Pilatus3 2M CdTe from the ID15 beam line of the ESRF. They provided not only the intersnip subject but also the couple of images uses to calibrate the detector.

This detector contains 48 half-modules, each bound to a single CdTe monocrystal sensor and is designed for high energy X-ray radiation detection. Due to the construction procedure, these half-modules could show a misalignment within the detector plane. While the manufacturer (Dectris) garanties a precision within a pixel (172µm), the miss-alignment of certain modules can be seen while calibrating Debye-Scherrer ring using refereance sample. So the aim of this work is to provide a detector description with a better precision better than the original detector.

This work will be performed on the image of a grid available: http://www.silx.org/pub/pyFAI/detector_calibration/Pilatus2MCdTe_ID15_grid_plus_sample_0004.cbf and the scattering of ceria (CeO2) at 72.1keV available here. http://www.silx.org/pub/pyFAI/detector_calibration/Pilatus2MCdTe_ID15_CeO2_72100eV_800mm_0000.cbf

It is a good exercise to calibrate all rings of the later image using the pyFAI-calib2 tool. A calibration close to perfection is needed to visualize the module miss-alignement we aim at correcting.

%matplotlib notebook

#many imports which will be used all along the notebook
import time
start_time = time.time()
import os
import pyFAI
import fabio
import glob
import numpy
from math import sin, cos, sqrt
from scipy.ndimage import convolve, binary_dilation
from scipy.spatial import distance_matrix
from scipy.optimize import minimize
from matplotlib.pyplot import subplots
from pyFAI.ext.bilinear import Bilinear
from pyFAI.ext.watershed import InverseWatershed
from silx.resources import ExternalResources

print("Using pyFAI verison: ", pyFAI.version)

Using pyFAI verison:  0.18.0

#Nota: Configure here your proxy if you are behind a firewall
#os.environ["http_proxy"] = "http://proxy.comany.com:3128"

downloader = ExternalResources("detector_calibration", "http://www.silx.org/pub/pyFAI/detector_calibration/")
ring_file = downloader.getfile("Pilatus2MCdTe_ID15_CeO2_72100eV_800mm_0000.cbf")
print(ring_file)
grid_file = downloader.getfile("Pilatus2MCdTe_ID15_grid_plus_sample_0004.cbf")
print(grid_file)

/tmp/detector_calibration_testdata_kieffer/Pilatus2MCdTe_ID15_CeO2_72100eV_800mm_0000.cbf
/tmp/detector_calibration_testdata_kieffer/Pilatus2MCdTe_ID15_grid_plus_sample_0004.cbf

rings = fabio.open(ring_file).data
img = fabio.open(grid_file).data
fig,ax = subplots(1,2, figsize=(10,5))
ax[0].imshow(img.clip(0,1000), interpolation="bilinear")
ax[0].set_title("grid")
ax[1].imshow(numpy.arcsinh(rings), interpolation="bilinear")
ax[1].set_title("rings")

Text(0.5, 1.0, 'rings')

Image processing

There are 3 pre-processing steps which are needed.

1. Define for each module a unique identifier which will be used later on during the fitting procedure

2. Define the proper mask: each module is the assembly of 4x2 sub-modules and there are (3) interpolated pixels between each sub-module, such “unreliable pixels should be masked out as well

3. Correct the grid image by the smoothed image to have a constant background.

4. Convolve the raw image with a typical hole shape to allow a precise spotting of the hole center.

# This is the default detector as definied in pyFAI according to the specification provided by Dectris:
pilatus = pyFAI.detector_factory("Pilatus_2m_CdTe")
print(pilatus)

mask1 = pilatus.mask
module_size = pilatus.MODULE_SIZE
module_gap = pilatus.MODULE_GAP
submodule_size = (96,60)

Detector Pilatus CdTe 2M         PixelSize= 1.720e-04, 1.720e-04 m

#1 + 2 Calculation of the module_id and the interpolated-mask:
mid = numpy.zeros(pilatus.shape, dtype=int)
mask2 = numpy.zeros(pilatus.shape, dtype=int)
idx = 1
for i in range(8):
    y_start = i*(module_gap[0] + module_size[0])
    y_stop = y_start + module_size[0]
    for j in range(3):
        x_start = j*(module_gap[1] + module_size[1])
        x_stop = x_start + module_size[1]
        mid[y_start:y_stop,x_start: x_start+module_size[1]//2] = idx
        idx+=1
        mid[y_start:y_stop,x_start+module_size[1]//2: x_stop] = idx
        idx+=1
        mask2[y_start+submodule_size[0]-1:y_start+submodule_size[0]+2,
              x_start:x_stop] = 1
        for k in range(1,8):
            mask2[y_start:y_stop,
              x_start+k*(submodule_size[1]+1)-1:x_start+k*(submodule_size[1]+1)+2] = 1

#Extra masking
mask0 = img<0
#Those pixel are miss-behaving... they are the hot pixels next to the beam-stop
mask0[915:922,793:800] = 1
mask0[817:820,747:750] = 1

fig,ax = subplots(1,3, figsize=(10,4))
ax[0].imshow(mid, interpolation="bilinear")
ax[0].set_title("Module Id")

ax[1].imshow(mask2+mask1+mask0, interpolation="bilinear")
ax[1].set_title("Combined mask")

nimg = img.astype(float)
nimg[numpy.where(mask0+mask1+mask2)] = numpy.nan


ax[2].imshow(nimg)#, interpolation="bilinear")
ax[2].set_title("Nan masked image")

Text(0.5, 1.0, 'Nan masked image')

# The Nan-masked image contains now only valid values (and Nan elsewhere). We will make a large median filter to
# build up a smooth image without gaps.
#
# This function is backported from future version of numpy ... it allows to expose a winbowed view
# to perform the nanmedian-filter
from numpy.lib.stride_tricks import as_strided
def sliding_window_view(x, shape, subok=False, readonly=True):
    """
    Creates sliding window views of the N dimensional array with the given window
    shape. Window slides across each dimension of `x` and extract subsets of `x`
    at any window position.
    Parameters
    ----------
    x : array_like
        Array to create sliding window views of.
    shape : sequence of int
        The shape of the window. Must have same length as the number of input array dimensions.
    subok : bool, optional
        If True, then sub-classes will be passed-through, otherwise the returned
        array will be forced to be a base-class array (default).
    readonly : bool, optional
        If set to True, the returned array will always be readonly view.
        Otherwise it will return writable copies(see Notes).
    Returns
    -------
    view : ndarray
        Sliding window views (or copies) of `x`. view.shape = x.shape - shape + 1
    See also
    --------
    as_strided: Create a view into the array with the given shape and strides.
    broadcast_to: broadcast an array to a given shape.
    Notes
    -----
    ``sliding_window_view`` create sliding window views of the N dimensions array
    with the given window shape and its implementation based on ``as_strided``.
    Please note that if readonly set to True, views are returned, not copies
    of array. In this case, write operations could be unpredictable, so the returned
    views are readonly. Bear in mind that returned copies (readonly=False) will
    take more memory than the original array, due to overlapping windows.
    For some cases there may be more efficient approaches to calculate transformations
    across multi-dimensional arrays, for instance `scipy.signal.fftconvolve`, where combining
    the iterating step with the calculation itself while storing partial results can result
    in significant speedups.
    Examples
    --------
    >>> i, j = np.ogrid[:3,:4]
    >>> x = 10*i + j
    >>> shape = (2,2)
    >>> np.lib.stride_tricks.sliding_window_view(x, shape)
    array([[[[ 0,  1],
             [10, 11]],
            [[ 1,  2],
             [11, 12]],
            [[ 2,  3],
             [12, 13]]],
           [[[10, 11],
             [20, 21]],
            [[11, 12],
             [21, 22]],
            [[12, 13],
             [22, 23]]]])
    """
    np = numpy
    # first convert input to array, possibly keeping subclass
    x = np.array(x, copy=False, subok=subok)

    try:
        shape = np.array(shape, np.int)
    except:
        raise TypeError('`shape` must be a sequence of integer')
    else:
        if shape.ndim > 1:
            raise ValueError('`shape` must be one-dimensional sequence of integer')
        if len(x.shape) != len(shape):
            raise ValueError("`shape` length doesn't match with input array dimensions")
        if np.any(shape <= 0):
            raise ValueError('`shape` cannot contain non-positive value')

    o = np.array(x.shape) - shape  + 1 # output shape
    if np.any(o <= 0):
        raise ValueError('window shape cannot larger than input array shape')

    if type(readonly) != bool:
        raise TypeError('readonly must be a boolean')

    strides = x.strides
    view_strides = strides

    view_shape = np.concatenate((o, shape), axis=0)
    view_strides = np.concatenate((view_strides, strides), axis=0)
    view = as_strided(x, view_shape, view_strides, subok=subok, writeable=not readonly)

    if not readonly:
        return view.copy()
    else:
        return view

%%time
#Calculate a background image using a large median filter ... takes a while
shape = (19,11)
print(nimg.shape)
padded = numpy.pad(nimg, tuple((i//2,) for i in shape), mode="edge")
print(padded.shape)
background = numpy.nanmedian(sliding_window_view(padded,shape), axis = (-2,-1))
print(background.shape)
fig,ax = subplots()
ax.imshow(background)
ax.set_title("Background image")

(1679, 1475)
(1697, 1485)
(1679, 1475)

CPU times: user 24.9 s, sys: 4.06 s, total: 29 s
Wall time: 29 s

fig,ax = subplots(1,2, figsize=(9,5))

normalized = (nimg/background)

low = numpy.nanmin(normalized)
high = numpy.nanmax(normalized)
print(low, high)
normalized[numpy.isnan(normalized)] = 0
normalized /= high

ax[0].imshow(normalized)
ax[0].set_title("Normalized image")

ax[1].hist(normalized.ravel(), 100, range=(0,1))
ax[1].set_title("Histogram of intensities in normalized image")

0.0 17.728813559322035

Text(0.5, 1.0, 'Histogram of intensities in normalized image')

For a precise measurement of the peak position, one trick is to convolve the image with a pattern which looks like a hole of the grid.

#print the profile of the normalized image: the center is difficult to measure due to the small size of the hole.
fig,ax = subplots(2)
ax[0].plot(normalized[:,545])
ax[1].plot(normalized[536,:])

[<matplotlib.lines.Line2D at 0x7ff81c79e710>]

#Definition of the convolution kernel
ksize = 5
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-2.5)))

kernel = fadeout(d)
mini=kernel.sum()
print(mini)

fig,ax = subplots(1,3)
ax[0].imshow(d)
ax[0].set_title("Distance array")

ax[1].plot(numpy.linspace(0,5,100),fadeout(numpy.linspace(0,5,100)))
ax[1].set_title("fade-out curve")

ax[2].imshow(kernel)
ax[2].set_title("Convolution kernel")

19.63857792789662

Text(0.5, 1.0, 'Convolution kernel')

my_smooth = convolve(normalized, kernel, mode="constant", cval=0)/mini
print(my_smooth.shape)
fig,ax = subplots(1,2)
ax[0].imshow(normalized.clip(0,1))
ax[0].set_ylim(1050,1100)
ax[0].set_xlim(300,350)
ax[1].imshow(my_smooth.clip(0,1))
ax[1].set_ylim(1050,1100)
ax[1].set_xlim(300,350)
numpy.where(my_smooth == my_smooth.max())

(1679, 1475)

(array([1065]), array([338]))

#mask out all pixels too close to any masked position

all_masks = numpy.logical_or(numpy.logical_or(mask0,mask1),mask2)
print(all_masks.sum())
big_mask = binary_dilation(all_masks, iterations=ksize//2+1+1)
print(big_mask.sum())
smooth2 = my_smooth.copy()
smooth2[big_mask] = 0
fig,ax = subplots()
ax.imshow(smooth2)

335009
782371

<matplotlib.image.AxesImage at 0x7ff81c647e80>

#Display the profile of the smoothed image: the center is easy to measure thanks to the smoothness of the signal
fig,ax = subplots(2)
ax[0].plot(my_smooth[:,545])
ax[1].plot(my_smooth[536,:])

[<matplotlib.lines.Line2D at 0x7ff81c62ca58>]

Peak picking

We use the watershed module from pyFAI to retrieve all peak positions. Those regions are sieved out respectively for:

• their size, it should be larger than the kernel itself

• the peaks too close to masked regions are removed

• the intensity of the peak

iw = InverseWatershed(my_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: 82126
Number of large enough regions : 41333

#Remove peaks on masked region
sieved_region = [i for i in regions if not big_mask[(i.index//nimg.shape[-1], i.index%nimg.shape[-1])]]
print("Number of peaks not on masked areea : %s"%len(sieved_region))

Number of peaks not on masked areea : 30001

# Histogram of peak height:
s = numpy.array([i.maxi for i in sieved_region])

fig, ax = subplots()
ax.hist(s, 100)

(array([1.0000e+00, 0.0000e+00, 1.0000e+01, 3.9940e+03, 2.2274e+04,
        1.4800e+03, 1.1700e+02, 3.4000e+01, 8.0000e+00, 2.0000e+00,
        2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00,
        1.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00, 0.0000e+00,
        1.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00,
        0.0000e+00, 1.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        1.0000e+00, 2.0000e+00, 2.0000e+00, 1.0000e+00, 0.0000e+00,
        1.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 1.0000e+00, 0.0000e+00,
        1.0000e+00, 3.0000e+00, 1.0000e+00, 2.0000e+00, 0.0000e+00,
        0.0000e+00, 2.0000e+00, 1.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 1.0000e+00, 0.0000e+00, 0.0000e+00, 2.0000e+00,
        0.0000e+00, 0.0000e+00, 2.0000e+00, 0.0000e+00, 0.0000e+00,
        2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 1.0000e+00, 1.0000e+00, 3.0000e+00, 3.0000e+00,
        7.0000e+00, 1.0000e+01, 1.0000e+01, 2.1000e+01, 1.9000e+01,
        2.9000e+01, 3.5000e+01, 5.7000e+01, 7.3000e+01, 9.9000e+01,
        1.1800e+02, 1.5900e+02, 1.9400e+02, 1.9500e+02, 1.9700e+02,
        1.9500e+02, 1.6300e+02, 1.6200e+02, 1.0200e+02, 8.1000e+01,
        6.6000e+01, 3.9000e+01, 7.0000e+00, 1.0000e+00, 1.0000e+00]),
 array([0.04771975, 0.05009548, 0.05247121, 0.05484694, 0.05722266,
        0.05959839, 0.06197412, 0.06434985, 0.06672558, 0.06910131,
        0.07147704, 0.07385276, 0.07622849, 0.07860422, 0.08097995,
        0.08335568, 0.08573141, 0.08810713, 0.09048286, 0.09285859,
        0.09523432, 0.09761005, 0.09998578, 0.1023615 , 0.10473723,
        0.10711296, 0.10948869, 0.11186442, 0.11424015, 0.11661588,
        0.1189916 , 0.12136733, 0.12374306, 0.12611879, 0.12849452,
        0.13087025, 0.13324597, 0.1356217 , 0.13799743, 0.14037316,
        0.14274889, 0.14512462, 0.14750035, 0.14987607, 0.1522518 ,
        0.15462753, 0.15700326, 0.15937899, 0.16175472, 0.16413044,
        0.16650617, 0.1688819 , 0.17125763, 0.17363336, 0.17600909,
        0.17838482, 0.18076054, 0.18313627, 0.185512  , 0.18788773,
        0.19026346, 0.19263919, 0.19501491, 0.19739064, 0.19976637,
        0.2021421 , 0.20451783, 0.20689356, 0.20926929, 0.21164501,
        0.21402074, 0.21639647, 0.2187722 , 0.22114793, 0.22352366,
        0.22589938, 0.22827511, 0.23065084, 0.23302657, 0.2354023 ,
        0.23777803, 0.24015376, 0.24252948, 0.24490521, 0.24728094,
        0.24965667, 0.2520324 , 0.25440813, 0.25678385, 0.25915958,
        0.26153531, 0.26391104, 0.26628677, 0.2686625 , 0.27103822,
        0.27341395, 0.27578968, 0.27816541, 0.28054114, 0.28291687,
        0.2852926 ]),
 <a list of 100 Patch objects>)

#sieve-out for peak intensity
int_mini = 0.1
peaks = [(i.index//nimg.shape[-1], i.index%nimg.shape[-1]) for i in sieved_region if (i.maxi)>int_mini]
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.1: 2075

# Finally the peak positions are interpolated using a second order taylor expansion
# in thevinicy of the maximum value of the signal:

#Create bilinear interpolator
bl = Bilinear(my_smooth)

#Overlay raw peak coordinate and refined peak positions

ref_peaks = [bl.local_maxi(p) for p in peaks]
fig, ax = subplots()
ax.imshow(img.clip(0,1000), interpolation="nearest")
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:")
print(ref_peaks[:10])

Refined peak coordinate:
[(918.3804241418839, 279.5866185426712), (918.55060338974, 338.37308943271637), (918.7212917506695, 397.05022401735187), (918.8901954367757, 455.8412373960018), (919.2437788248062, 515.3245315551758), (919.2116163223982, 544.6306998133659), (919.2832362055779, 574.1425551921129), (919.0927048921585, 603.6582366526127), (771.483876645565, 191.8469839990139), (771.4615589380264, 221.4316090941429)]

At this stage we have about 2000 peaks (with sub-pixel precision) which are visually distributed on all modules. Some modules have their peaks located along sub-module boundaries which are masked out, hence they have fewer ontrol point for the calculation. Let’s assign each peak to a module identifier. This allows to print out the number of peaks per module:

dt = numpy.dtype([('y', numpy.float64),
                  ('x', numpy.float64),
                  ('i', numpy.int64),
                 ])
yxi = numpy.array([i+(mid[round(i[0]),round(i[1])],)
                   for i in ref_peaks], dtype=dt)
print("Number of keypoint per module:")
for i in range(1,mid.max()+1):
    print("Module id:",i, "cp:", (yxi[:]["i"] == i).sum())

Number of keypoint per module:
Module id: 1 cp: 48
Module id: 2 cp: 30
Module id: 3 cp: 48
Module id: 4 cp: 46
Module id: 5 cp: 42
Module id: 6 cp: 47
Module id: 7 cp: 47
Module id: 8 cp: 30
Module id: 9 cp: 48
Module id: 10 cp: 48
Module id: 11 cp: 41
Module id: 12 cp: 39
Module id: 13 cp: 48
Module id: 14 cp: 30
Module id: 15 cp: 47
Module id: 16 cp: 48
Module id: 17 cp: 42
Module id: 18 cp: 48
Module id: 19 cp: 47
Module id: 20 cp: 30
Module id: 21 cp: 48
Module id: 22 cp: 47
Module id: 23 cp: 42
Module id: 24 cp: 47
Module id: 25 cp: 48
Module id: 26 cp: 30
Module id: 27 cp: 48
Module id: 28 cp: 47
Module id: 29 cp: 41
Module id: 30 cp: 50
Module id: 31 cp: 46
Module id: 32 cp: 30
Module id: 33 cp: 48
Module id: 34 cp: 42
Module id: 35 cp: 47
Module id: 36 cp: 48
Module id: 37 cp: 48
Module id: 38 cp: 28
Module id: 39 cp: 48
Module id: 40 cp: 42
Module id: 41 cp: 44
Module id: 42 cp: 48
Module id: 43 cp: 47
Module id: 44 cp: 26
Module id: 45 cp: 46
Module id: 46 cp: 42
Module id: 47 cp: 45
Module id: 48 cp: 48

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.

# pairwise distance calculation using scipy.spatial.distance_matrix

dist = distance_matrix(peaks_ref, peaks_ref)

fig, ax = subplots()
ax.hist(dist.ravel(), 100, range=(0,100))
ax.set_title("Pair-wise distribution function")

Text(0.5, 1.0, 'Pair-wise distribution function')

The histogram of the pair-distribution function has a first peak at 0 and the second peak between 29 and 30. Let’s start the fit with this value

Two other parameters correspond to the offset, in pixel for the grid index (X,Y) = (0,0). The easiest is to measure the smallest x and y for the first module.

#from pair-wise distribution histogram
step = 29
#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 peak: ", x_min, y_min)

offset for the first peak:  16.269793540239334 7.186804354190826

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.

The default rotation will be zero and will be fitted later on.

Once the indexes X,Y determined for eack peak, one can fit the parameter to properly align the grid with the first module. Those 4 parameters are step-size, x_min, y_min and angle

#Assign each peak to an index
dl = numpy.dtype([('y', numpy.float64),
                  ('x', numpy.float64),
                  ('i', numpy.int64),
                  ('Y', numpy.int64),
                  ('X', numpy.int64),
                 ])
indexed1 = numpy.zeros(len(first), dtype=dl)

for i,v in enumerate(first):
    Y = int(round((v["y"]-y_min)/step))
    X = int(round((v["x"]-x_min)/step))
    indexed1[i]["y"] = v["y"]
    indexed1[i]["x"] = v["x"]
    indexed1[i]["i"] = v["i"]
    indexed1[i]["Y"] = Y
    indexed1[i]["X"] = X
    print("peak id: %s %20s Y:%d (Δ=%.3f) X:%s (Δ=%.3f)"%
          (i,v, Y, (v["y"]-Y*step-y_min)/step, X, (v["x"]-X*step-x_min)/step))

peak id: 0 (65.94893547, 16.45467511, 1) Y:2 (Δ=0.026) X:0 (Δ=0.006)
peak id: 1 (65.94483702, 45.61336857, 1) Y:2 (Δ=0.026) X:1 (Δ=0.012)
peak id: 2 (66.03506869, 75.15414216, 1) Y:2 (Δ=0.029) X:2 (Δ=0.030)
peak id: 3 (65.89957739, 104.38019937, 1) Y:2 (Δ=0.025) X:3 (Δ=0.038)
peak id: 4 (66.12287057, 133.89461011, 1) Y:2 (Δ=0.032) X:4 (Δ=0.056)
peak id: 5 (66.1012452, 163.39174998, 1) Y:2 (Δ=0.032) X:5 (Δ=0.073)
peak id: 6 (66.10113759, 192.66102123, 1) Y:2 (Δ=0.032) X:6 (Δ=0.082)
peak id: 7 (66.11095764, 222.20659751, 1) Y:2 (Δ=0.032) X:7 (Δ=0.101)
peak id: 8 (36.49696946, 16.50531429, 1) Y:1 (Δ=0.011) X:0 (Δ=0.008)
peak id: 9 (36.54493174, 45.66864574, 1) Y:1 (Δ=0.012) X:1 (Δ=0.014)
peak id: 10 (36.57917181, 75.22501931, 1) Y:1 (Δ=0.014) X:2 (Δ=0.033)
peak id: 11 (36.48521471, 104.38776881, 1) Y:1 (Δ=0.010) X:3 (Δ=0.039)
peak id: 12 (36.5958997, 133.89308376, 1) Y:1 (Δ=0.014) X:4 (Δ=0.056)
peak id: 13 (36.56109262, 163.54238695, 1) Y:1 (Δ=0.013) X:5 (Δ=0.078)
peak id: 14 (36.62988734, 192.63348889, 1) Y:1 (Δ=0.015) X:6 (Δ=0.082)
peak id: 15 (36.66130301, 222.22433692, 1) Y:1 (Δ=0.016) X:7 (Δ=0.102)
peak id: 16 (7.18680435, 16.41273034, 1) Y:0 (Δ=0.000) X:0 (Δ=0.005)
peak id: 17 (7.23123039, 45.74501115, 1) Y:0 (Δ=0.002) X:1 (Δ=0.016)
peak id: 18 (7.29059145, 75.19019249, 1) Y:0 (Δ=0.004) X:2 (Δ=0.032)
peak id: 19 (7.19876912, 104.56766725, 1) Y:0 (Δ=0.000) X:3 (Δ=0.045)
peak id: 20 (7.28363913, 134.01557775, 1) Y:0 (Δ=0.003) X:4 (Δ=0.060)
peak id: 21 (7.25490889, 163.38282073, 1) Y:0 (Δ=0.002) X:5 (Δ=0.073)
peak id: 22 (7.36717248, 192.76562056, 1) Y:0 (Δ=0.006) X:6 (Δ=0.086)
peak id: 23 (7.37469381, 222.16684446, 1) Y:0 (Δ=0.006) X:7 (Δ=0.100)
peak id: 24 (183.49989897, 16.39247844, 1) Y:6 (Δ=0.080) X:0 (Δ=0.004)
peak id: 25 (183.50560439, 45.46551961, 1) Y:6 (Δ=0.080) X:1 (Δ=0.007)
peak id: 26 (183.60742423, 75.08415242, 1) Y:6 (Δ=0.083) X:2 (Δ=0.028)
peak id: 27 (183.51583004, 104.62201929, 1) Y:6 (Δ=0.080) X:3 (Δ=0.047)
peak id: 28 (183.66220078, 133.73585817, 1) Y:6 (Δ=0.085) X:4 (Δ=0.051)
peak id: 29 (183.65553248, 163.30831155, 1) Y:6 (Δ=0.085) X:5 (Δ=0.070)
peak id: 30 (183.60811904, 192.46339762, 1) Y:6 (Δ=0.083) X:6 (Δ=0.076)
peak id: 31 (183.75545964, 222.06266521, 1) Y:6 (Δ=0.089) X:7 (Δ=0.096)
peak id: 32 (154.25518942, 16.26979354, 1) Y:5 (Δ=0.071) X:0 (Δ=0.000)
peak id: 33 (154.21079713, 45.57556251, 1) Y:5 (Δ=0.070) X:1 (Δ=0.011)
peak id: 34 (154.28871083, 75.04928852, 1) Y:5 (Δ=0.072) X:2 (Δ=0.027)
peak id: 35 (154.21604921, 104.4366582, 1) Y:5 (Δ=0.070) X:3 (Δ=0.040)
peak id: 36 (154.37423909, 133.8660319, 1) Y:5 (Δ=0.075) X:4 (Δ=0.055)
peak id: 37 (154.35255641, 163.24378996, 1) Y:5 (Δ=0.075) X:5 (Δ=0.068)
peak id: 38 (154.32582822, 192.64246738, 1) Y:5 (Δ=0.074) X:6 (Δ=0.082)
peak id: 39 (154.46003968, 222.04596291, 1) Y:5 (Δ=0.078) X:7 (Δ=0.096)
peak id: 40 (124.69245848, 16.47663069, 1) Y:4 (Δ=0.052) X:0 (Δ=0.007)
peak id: 41 (124.74490786, 45.56257215, 1) Y:4 (Δ=0.054) X:1 (Δ=0.010)
peak id: 42 (124.81674041, 75.13362026, 1) Y:4 (Δ=0.056) X:2 (Δ=0.030)
peak id: 43 (124.75223885, 104.64793801, 1) Y:4 (Δ=0.054) X:3 (Δ=0.048)
peak id: 44 (124.84565204, 133.79748473, 1) Y:4 (Δ=0.057) X:4 (Δ=0.053)
peak id: 45 (124.85729225, 163.3060154, 1) Y:4 (Δ=0.058) X:5 (Δ=0.070)
peak id: 46 (124.89991033, 192.57832247, 1) Y:4 (Δ=0.059) X:6 (Δ=0.080)
peak id: 47 (124.96431521, 222.11254385, 1) Y:4 (Δ=0.061) X:7 (Δ=0.098)

The error in positionning each of the pixel is less than 0.1 pixel which is already excellent and will allow a straight forward fit.

The cost function for the first module is calculated as the sum of distances squared in pixel space. It uses 4 parameters which are step-size, x_min, y_min and angle

guess1 = [step, y_min, x_min, 0]

def cost1(param):
    """contains: step, y_min, x_min, angle for the first module
    returns the sum of distance squared in pixel space
    """
    step = param[0]
    y_min = param[1]
    x_min = param[2]
    angle = param[3]
    XY = numpy.vstack((indexed1["X"], indexed1["Y"]))
    rot = [[cos(angle),-sin(angle)],
           [sin(angle), cos(angle)]]
    xy_min = [[x_min], [y_min]]
    xy_guess = numpy.dot(rot, step*XY+xy_min)
    delta = xy_guess - numpy.vstack((indexed1["x"], indexed1["y"]))
    return (delta*delta).sum()

print("Before optimization", guess1, "cost=", cost1(guess1))
res1 = minimize(cost1, guess1, method = "slsqp")
print(res1)
print("After optimization", res1.x, "cost=", cost1(res1.x))
print("Average displacement (pixels): ",sqrt(cost1(res1.x)/len(indexed1)))

Before optimization [29, 7.186804354190826, 16.269793540239334, 0] cost= 250.36710826038234
     fun: 0.6337542147785734
     jac: array([ 2.15254724e-04, -1.18207186e-03,  1.05350465e-03, -3.25701334e-01])
 message: 'Optimization terminated successfully.'
    nfev: 51
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([2.93980817e+01, 7.13605736e+00, 1.63978976e+01, 8.76803451e-04])
After optimization [2.93980817e+01 7.13605736e+00 1.63978976e+01 8.76803451e-04] cost= 0.6337542147785734
Average displacement (pixels):  0.11490523403173132

At this step, the grid is perfectly aligned with the first module. This module is used as the reference one and all other are aligned along it, using this first fit:

#retrieve the result of the first module fit:
step, y_min, x_min, angle = res1.x
indexed = numpy.zeros(yxi.shape, dtype=dl)

rot =  [[cos(angle),-sin(angle)],
        [sin(angle), cos(angle)]]
irot =  [[cos(angle), sin(angle)],
         [-sin(angle), cos(angle)]]

print("cost1: ",cost1([step, y_min, x_min, angle]), "for:", step, y_min, x_min, angle)

xy_min = numpy.array([[x_min], [y_min]])
xy = numpy.vstack((yxi["x"], yxi["y"]))
indexed["y"] = yxi["y"]
indexed["x"] = yxi["x"]
indexed["i"] = yxi["i"]
XY_app = (numpy.dot(irot, xy)-xy_min) / step
XY_int = numpy.round((XY_app)).astype("int")
indexed["X"] = XY_int[0]
indexed["Y"] = XY_int[1]
xy_guess = numpy.dot(rot, step * XY_int + xy_min)

thres = 1.2
delta = abs(xy_guess - xy)
print((delta>thres).sum())
suspicious = indexed[numpy.where(abs(delta>thres))[1]]
print(suspicious)

cost1:  0.6337542147785734 for: 29.39808169136302 7.136057363507019 16.3978976414565 0.0008768034506612712
6
[(1654.4202804 , 954.44135761, 46, 56, 32)
 (1624.68617365, 807.47743243, 46, 55, 27)
 (1595.5792897 , 954.49445057, 46, 54, 32)
 (1448.30209509, 748.8784325 , 40, 49, 25)
 (1448.66381988, 895.86545633, 40, 49, 30)
 (1448.84506877, 954.56920397, 40, 49, 32)]

fig,ax = subplots()
ax.imshow(img.clip(0,1000))
ax.plot(indexed["x"], indexed["y"],".g")
ax.plot(suspicious["x"], suspicious["y"],".r")

[<matplotlib.lines.Line2D at 0x7ff811de84a8>]

Only 6 peaks have an initial displacement of more than 1.2 pixel, all located in modules 40 and 46. The visual inspection confirms their localization is valid.

There are 48 (half-)modules which have each of them 2 translations and one rotation. In addition to the step size, this represents 145 degrees of freedom for the fit. The first module is used to align the grid, all other modules are then aligned along this grid.

def submodule_cost(param, module=1):
    """contains: step, y_min_1, x_min_1, angle_1, y_min_2, x_min_2, angle_2, ...
    returns the sum of distance squared in pixel space
    """

    step = param[0]
    y_min1 = param[1]
    x_min1 = param[2]
    angle1 = param[3]

    mask = indexed["i"] == module
    substack = indexed[mask]

    XY = numpy.vstack((substack["X"], substack["Y"]))
    rot1 = [[cos(angle1), -sin(angle1)],
            [sin(angle1), cos(angle1)]]
    xy_min1 = numpy.array([[x_min1], [y_min1]])
    xy_guess1 = numpy.dot(rot1, step * XY + xy_min1)
    #This is guessed spot position for module #1
    if module == 1:
        "Not much to do for module 1"
        delta = xy_guess1 - numpy.vstack((substack["x"], substack["y"]))
    else:
        "perform the correction for given module"
        y_min = param[(module-1)*3+1]
        x_min = param[(module-1)*3+2]
        angle = param[(module-1)*3+3]

        rot = numpy.array([[cos(angle),-sin(angle)],
                           [sin(angle), cos(angle)]])
        xy_min = numpy.array([[x_min], [y_min]])
        xy_guess = numpy.dot(rot, xy_guess1+xy_min)
        delta = xy_guess - numpy.vstack((substack["x"], substack["y"]))

    return (delta*delta).sum()

guess145 = numpy.zeros(48*3+1)
guess145[:4] = res1.x
for i in range(1, 49):
    print("Cost for module #",i, submodule_cost(guess145, i))

Cost for module # 1 0.6337542147785734
Cost for module # 2 1.8137335142043005
Cost for module # 3 1.8517810791354326
Cost for module # 4 4.9127875221043515
Cost for module # 5 1.773711051583296
Cost for module # 6 7.48693792691493
Cost for module # 7 14.439344856965251
Cost for module # 8 5.219239827628602
Cost for module # 9 23.403682459536583
Cost for module # 10 5.689200096256752
Cost for module # 11 4.128382716158292
Cost for module # 12 8.41469167554883
Cost for module # 13 29.862717021161917
Cost for module # 14 12.629471494880102
Cost for module # 15 3.1375750005201866
Cost for module # 16 5.280069866042303
Cost for module # 17 8.497446247607941
Cost for module # 18 4.655804873784172
Cost for module # 19 6.008515385043476
Cost for module # 20 4.982313469667585
Cost for module # 21 5.13721994073221
Cost for module # 22 5.050522255483948
Cost for module # 23 1.437001287138772
Cost for module # 24 6.050459713222455
Cost for module # 25 38.26262874059206
Cost for module # 26 19.514969043053544
Cost for module # 27 2.4311336259686778
Cost for module # 28 5.491303543176383
Cost for module # 29 10.770599187417773
Cost for module # 30 14.051565356868426
Cost for module # 31 22.948589256445068
Cost for module # 32 18.498872034582465
Cost for module # 33 12.94939072065644
Cost for module # 34 15.380312194242064
Cost for module # 35 39.681054141224514
Cost for module # 36 37.298625739026946
Cost for module # 37 63.71626527040634
Cost for module # 38 23.55576001555159
Cost for module # 39 34.01669588160015
Cost for module # 40 43.81908020764373
Cost for module # 41 15.51614503136555
Cost for module # 42 8.057425004151995
Cost for module # 43 22.74618763379752
Cost for module # 44 10.802033729233333
Cost for module # 45 41.03620867139952
Cost for module # 46 46.276855734650006
Cost for module # 47 24.5621299905672
Cost for module # 48 34.78883325879872

On retrieves that the modules 40 and 46 have large errors. Module 37 as well.

The total cost funtion is hence the sum of all cost function for all modules:

def total_cost(param):
    """contains: step, y_min_1, x_min_1, angle_1, ...
    returns the sum of distance squared in pixel space
    """
    return sum(submodule_cost(param, module=i) for i in range(1,49))
total_cost(guess145)

778.6690275085203

%%time
print("Before optimization", guess145[:10], "cost=", total_cost(guess145))
res_all = minimize(total_cost, guess145, method = "slsqp")
print(res_all)
print("After optimization", res_all.x[:10], "cost=", total_cost(res_all.x))

Before optimization [2.93980817e+01 7.13605736e+00 1.63978976e+01 8.76803451e-04
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 0.00000000e+00 0.00000000e+00] cost= 778.6690275085203
     fun: 27.18365832628731
     jac: array([ 3.36764812e-01,  2.41279602e-04,  1.03421211e-02, -8.09069777e+00,
       -2.82526016e-04, -8.50439072e-04, -1.01323128e-02, -8.08238983e-05,
       -3.63349915e-04, -1.85797215e-02, -1.40905380e-04, -9.55820084e-04,
       -8.91995430e-03,  1.38282776e-05,  8.19206238e-04, -1.05245113e-01,
       -2.09569931e-04,  1.52707100e-03, -1.26422167e-01,  2.03394890e-03,
        8.81671906e-04, -3.47824097e-02,  1.43098831e-03,  1.91617012e-03,
       -2.55455971e-02, -3.88145447e-04, -1.01804733e-04, -1.31155491e-01,
       -3.17335129e-04,  2.88486481e-04, -1.11642122e-01,  1.31368637e-04,
       -3.75986099e-04,  8.54325294e-03, -1.75714493e-04, -1.26838684e-04,
       -1.08163357e-01,  3.51667404e-04,  3.56912613e-04, -9.72857475e-02,
       -7.34567642e-04, -7.43865967e-04, -3.60100269e-02, -1.92165375e-04,
       -1.22308731e-04, -6.77630901e-02,  2.69412994e-05,  3.32355499e-04,
       -1.52028084e-01, -7.98940659e-04, -1.35207176e-03, -3.32699060e-01,
        4.24385071e-04,  4.83512878e-04, -2.34975100e-01, -1.58643723e-03,
       -1.14440918e-05, -9.22925472e-02, -1.26242638e-03, -5.36918640e-04,
       -4.03177738e-02,  2.01225281e-04,  2.02655792e-04, -1.15622759e-01,
       -3.58581543e-04, -3.13043594e-04, -1.85351372e-01,  3.60727310e-04,
        2.01702118e-04, -1.06734037e-01,  4.86373901e-04,  7.61508942e-04,
       -8.16304684e-02,  8.49723816e-04,  2.74896622e-04, -1.41849279e-01,
        8.97884369e-04,  3.00168991e-04, -8.56087208e-02,  7.53402710e-04,
        3.68595123e-04,  1.32328033e-01, -1.43051147e-05, -1.49726868e-04,
        5.17210960e-02, -1.03306770e-03, -8.61167908e-04, -1.73690319e-01,
        9.89198685e-04,  4.65631485e-04, -9.13746357e-02,  4.30822372e-04,
        2.80380249e-04, -1.09273434e-01,  2.12192535e-04,  2.64644623e-04,
       -1.86563969e-01,  5.52654266e-04,  1.77383423e-04, -2.00356483e-01,
       -7.89165497e-05, -2.15530396e-04, -1.59254074e-01, -1.34491920e-03,
       -7.27415085e-04, -3.59580994e-01, -3.62873077e-04,  6.33001328e-04,
       -4.18133974e-01,  7.90119171e-04,  4.77075577e-04, -3.38363886e-01,
       -5.48362732e-06,  1.46865845e-04, -1.75536156e-01,  3.64065170e-04,
        8.61644745e-04, -2.99779892e-01,  6.23703003e-04,  2.30121613e-03,
       -2.43101120e-01,  7.79151917e-04,  9.93490219e-04, -3.10999632e-01,
        4.47273254e-04,  1.13582611e-03, -3.22041512e-01, -1.87778473e-03,
        2.86579132e-04, -3.00019503e-01,  2.96115875e-04, -5.55753708e-04,
       -3.44004393e-01, -1.32369995e-03, -2.36272812e-04, -3.29681873e-01,
       -1.74283981e-04,  6.97374344e-04, -3.31437826e-01, -5.36918640e-04,
        5.60283661e-05, -2.70724535e-01,  1.46508217e-03,  2.64167786e-04,
       -1.53188705e-01])
 message: 'Optimization terminated successfully.'
    nfev: 14738
     nit: 97
    njev: 97
  status: 0
 success: True
       x: array([ 2.94082848e+01,  7.10550236e+00,  1.63621431e+01,  8.76343987e-04,
       -2.61142002e-01, -1.72404666e-01,  2.75358939e-04, -4.65371706e-01,
       -1.49368876e-01,  9.23377934e-04, -1.40925206e+00, -3.92431812e-02,
        1.93786320e-03, -4.71619538e-01, -4.34506153e-01,  3.15955995e-04,
       -2.55783527e+00, -4.42186027e-01,  2.08129218e-03, -7.62295275e-01,
        5.00535589e-01,  1.34360338e-03, -8.36248315e-01,  2.10686060e-01,
        1.08362469e-03, -1.27272229e+00, -1.27829390e-01,  8.94580313e-04,
       -2.05296929e+00,  2.04524448e-01,  2.10557127e-03, -3.29262534e-01,
       -1.11556153e-01,  4.68865734e-04, -5.88218340e-01,  3.71081910e-02,
        6.46130657e-04, -8.24199310e-01,  4.59685359e-02,  1.14564520e-03,
       -9.99768822e-01,  1.20247034e-01,  1.31121860e-03, -9.06361732e-01,
        1.88553858e-01,  9.65779253e-04, -1.14461914e+00,  1.26388275e-01,
        1.27134440e-03, -8.21266979e-01, -5.18563664e-01,  4.32107560e-04,
       -2.16110379e+00,  1.33708683e-01,  1.53367349e-03, -5.67627925e-01,
        3.19370766e-01,  7.07441671e-04, -4.90631071e-01,  6.01221530e-02,
        7.15756603e-04, -2.58645439e-01, -6.04060774e-01, -3.01162879e-04,
       -1.66726458e+00,  6.64322082e-01,  1.62371556e-03, -1.58599237e+00,
        6.43415301e-01,  1.23645918e-03, -2.16623943e+00,  8.25942522e-01,
        1.65256754e-03, -8.55872469e-01,  1.36470630e-01,  9.44940792e-04,
       -1.04916637e+00,  8.15407272e-01,  1.75160691e-03,  2.14373439e-01,
       -6.90411962e-01, -5.25494569e-04, -4.87271205e-01,  3.35472535e-01,
        5.93598781e-04, -7.51957868e-01, -2.87512355e-01,  5.55801951e-04,
       -3.03293420e+00,  1.59247632e+00,  2.32811165e-03, -4.66935314e-01,
       -1.52659742e+00, -7.47172216e-04, -1.21117247e+00,  8.98551575e-01,
        1.43614378e-03,  2.05750982e-01, -1.38157673e+00, -6.38868427e-04,
       -1.40203507e+00,  1.03337865e+00,  1.55563831e-03, -9.47399558e-01,
       -5.13283043e-01,  6.36750009e-04, -8.52438513e-01, -6.12855414e-01,
        5.58099763e-04, -1.02281342e+00, -1.03557970e+00, -3.49555463e-05,
       -9.28000183e-01, -8.97165498e-01, -4.07212989e-05, -6.65844032e-01,
       -8.46731364e-01,  1.07502923e-04, -2.37944396e+00,  2.10370247e+00,
        2.43396212e-03, -5.84541323e-01, -4.18588167e-01,  3.37284557e-04,
       -1.08753981e+00,  6.65981023e-02,  6.11100460e-04, -9.52909811e-01,
        9.59526263e-01,  9.85626571e-04, -8.37807529e-01,  1.62657619e-01,
        5.46469118e-04, -3.42723669e-01, -1.67395622e+00, -3.58118168e-04,
       -1.50828980e+00,  7.38683145e-01,  1.28170883e-03, -1.40264511e+00,
        8.68409721e-01,  1.15128804e-03, -2.28992308e+00,  2.04723385e+00,
        1.83647125e-03])
After optimization [ 2.94082848e+01  7.10550236e+00  1.63621431e+01  8.76343987e-04
 -2.61142002e-01 -1.72404666e-01  2.75358939e-04 -4.65371706e-01
 -1.49368876e-01  9.23377934e-04] cost= 27.18365832628731
CPU times: user 25.7 s, sys: 44.7 ms, total: 25.8 s
Wall time: 25.8 s

for i in range(1,49):
    print("Module id: %d cost: %.3f Δx: %.3f, Δy: %.3f rot: %.3f°"%
          (i, submodule_cost(res_all.x, i), res_all.x[-2+i*3], res_all.x[-1+i*3], numpy.rad2deg(res_all.x[i*3])))

Module id: 1 cost: 0.683 Δx: 7.106, Δy: 16.362 rot: 0.050°
Module id: 2 cost: 0.473 Δx: -0.261, Δy: -0.172 rot: 0.016°
Module id: 3 cost: 0.650 Δx: -0.465, Δy: -0.149 rot: 0.053°
Module id: 4 cost: 0.613 Δx: -1.409, Δy: -0.039 rot: 0.111°
Module id: 5 cost: 0.566 Δx: -0.472, Δy: -0.435 rot: 0.018°
Module id: 6 cost: 0.666 Δx: -2.558, Δy: -0.442 rot: 0.119°
Module id: 7 cost: 0.696 Δx: -0.762, Δy: 0.501 rot: 0.077°
Module id: 8 cost: 0.341 Δx: -0.836, Δy: 0.211 rot: 0.062°
Module id: 9 cost: 0.566 Δx: -1.273, Δy: -0.128 rot: 0.051°
Module id: 10 cost: 0.661 Δx: -2.053, Δy: 0.205 rot: 0.121°
Module id: 11 cost: 0.616 Δx: -0.329, Δy: -0.112 rot: 0.027°
Module id: 12 cost: 0.495 Δx: -0.588, Δy: 0.037 rot: 0.037°
Module id: 13 cost: 0.763 Δx: -0.824, Δy: 0.046 rot: 0.066°
Module id: 14 cost: 0.353 Δx: -1.000, Δy: 0.120 rot: 0.075°
Module id: 15 cost: 0.645 Δx: -0.906, Δy: 0.189 rot: 0.055°
Module id: 16 cost: 0.492 Δx: -1.145, Δy: 0.126 rot: 0.073°
Module id: 17 cost: 0.440 Δx: -0.821, Δy: -0.519 rot: 0.025°
Module id: 18 cost: 0.779 Δx: -2.161, Δy: 0.134 rot: 0.088°
Module id: 19 cost: 0.748 Δx: -0.568, Δy: 0.319 rot: 0.041°
Module id: 20 cost: 0.423 Δx: -0.491, Δy: 0.060 rot: 0.041°
Module id: 21 cost: 0.520 Δx: -0.259, Δy: -0.604 rot: -0.017°
Module id: 22 cost: 0.676 Δx: -1.667, Δy: 0.664 rot: 0.093°
Module id: 23 cost: 0.503 Δx: -1.586, Δy: 0.643 rot: 0.071°
Module id: 24 cost: 0.822 Δx: -2.166, Δy: 0.826 rot: 0.095°
Module id: 25 cost: 0.575 Δx: -0.856, Δy: 0.136 rot: 0.054°
Module id: 26 cost: 0.320 Δx: -1.049, Δy: 0.815 rot: 0.100°
Module id: 27 cost: 0.642 Δx: 0.214, Δy: -0.690 rot: -0.030°
Module id: 28 cost: 0.599 Δx: -0.487, Δy: 0.335 rot: 0.034°
Module id: 29 cost: 0.498 Δx: -0.752, Δy: -0.288 rot: 0.032°
Module id: 30 cost: 0.830 Δx: -3.033, Δy: 1.592 rot: 0.133°
Module id: 31 cost: 0.687 Δx: -0.467, Δy: -1.527 rot: -0.043°
Module id: 32 cost: 0.422 Δx: -1.211, Δy: 0.899 rot: 0.082°
Module id: 33 cost: 0.616 Δx: 0.206, Δy: -1.382 rot: -0.037°
Module id: 34 cost: 0.393 Δx: -1.402, Δy: 1.033 rot: 0.089°
Module id: 35 cost: 0.447 Δx: -0.947, Δy: -0.513 rot: 0.036°
Module id: 36 cost: 0.587 Δx: -0.852, Δy: -0.613 rot: 0.032°
Module id: 37 cost: 0.603 Δx: -1.023, Δy: -1.036 rot: -0.002°
Module id: 38 cost: 0.295 Δx: -0.928, Δy: -0.897 rot: -0.002°
Module id: 39 cost: 0.558 Δx: -0.666, Δy: -0.847 rot: 0.006°
Module id: 40 cost: 0.498 Δx: -2.379, Δy: 2.104 rot: 0.139°
Module id: 41 cost: 0.468 Δx: -0.585, Δy: -0.419 rot: 0.019°
Module id: 42 cost: 0.569 Δx: -1.088, Δy: 0.067 rot: 0.035°
Module id: 43 cost: 0.633 Δx: -0.953, Δy: 0.960 rot: 0.056°
Module id: 44 cost: 0.413 Δx: -0.838, Δy: 0.163 rot: 0.031°
Module id: 45 cost: 0.671 Δx: -0.343, Δy: -1.674 rot: -0.021°
Module id: 46 cost: 0.542 Δx: -1.508, Δy: 0.739 rot: 0.073°
Module id: 47 cost: 0.537 Δx: -1.403, Δy: 0.868 rot: 0.066°
Module id: 48 cost: 0.593 Δx: -2.290, Δy: 2.047 rot: 0.105°

Analysis: Modules 40, 46 and 48 show large displacement but the fitting precedure allowed to reduce the residual cost to the same value as other modules.

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, except the first one which is left where it is.

pixel_coord = pyFAI.detector_factory("Pilatus2MCdTe").get_pixel_corners()
pixel_coord_raw = pixel_coord.copy()
for i in range(2, 49):
    # Extract the pixel corners for one module
    module_idx = numpy.where(mid == i)
    one_module = pixel_coord_raw[module_idx]
    #retrieve the fitted values
    dy = res_all.x[1+(i-1)*3]
    dx = res_all.x[2+(i-1)*3]
    angle = res_all.x[3+(i-1)*3]

    z = one_module[...,0]
    y = one_module[...,1]
    x = one_module[...,2]

    #apply the correction the other way around
    irot = [[cos(angle), sin(angle)],
           [-sin(angle), cos(angle)]]
    trans = [[dx * pilatus.pixel2],
              [dy * pilatus.pixel1]]
    xy_guess = numpy.vstack((x.ravel(), y.ravel()))
    xy_cor = numpy.dot(irot, xy_guess) - trans
    xy_cor.shape = ((2,)+x.shape)
    one_module[...,1] = xy_cor[1, ...] #y
    one_module[...,2] = xy_cor[0, ...] #x
    #Update the array
    pixel_coord_raw[module_idx] = one_module

Update the detector and save it in HDF5

pilatus.set_pixel_corners(pixel_coord_raw)
pilatus.mask = all_masks
pilatus.save("Pilatus_ID15_raw.h5")

/mntdirect/_scisoft/users/kieffer/.jupy37/lib/python3.7/site-packages/pyFAI/io/__init__.py:999: H5pyDeprecationWarning: The default file mode will change to 'r' (read-only) in h5py 3.0. To suppress this warning, pass the mode you need to h5py.File(), or set the global default h5.get_config().default_file_mode, or set the environment variable H5PY_DEFAULT_READONLY=1. Available modes are: 'r', 'r+', 'w', 'w-'/'x', 'a'. See the docs for details.
  self.h5 = h5py.File(self.filename)

displ = numpy.sqrt(((pixel_coord - pixel_coord_raw)**2).sum(axis=-1))
displ /= pilatus.pixel1 #convert in pixel units
fig, ax = subplots()
ax.hist(displ.ravel(), 100)
ax.set_title("Displacement of pixels versus the reference representation")
ax.set_xlabel("Error in pixel size (172µm)")

Text(0.5, 0, 'Error in pixel size (172µm)')

unmasked = numpy.logical_not(all_masks)
misaligned = numpy.vstack((pixel_coord_raw[unmasked][..., 2].ravel(), #x
                           pixel_coord_raw[unmasked][..., 1].ravel())) #y

reference = numpy.vstack((pixel_coord[unmasked][..., 2].ravel(), #x
                          pixel_coord[unmasked][..., 1].ravel())) #y

def global_displacement(param):
    """correct all pixel position with dy. dx, angle"""
    dy, dx, angle = param
    rot = [[cos(angle), -sin(angle)],
           [sin(angle), cos(angle)]]
    trans = [[dx],
              [dy]]
    xy_corr = numpy.dot(rot, misaligned) + trans
    delta = xy_corr - reference
    return (delta*delta).sum()
global_displacement([0,0,0])

0.16861754193929496

%%time
guess_global = [0,0,0]
print("Before optimization", guess_global, "cost=", global_displacement(guess_global))
res_global = minimize(global_displacement, guess_global, method = "slsqp")
print(res_global)
print("After optimization", res_global.x, "cost=", global_displacement(res_global.x))

Before optimization [0, 0, 0] cost= 0.16861754193929496
     fun: 0.02783531839357946
     jac: array([ 0.24096283, -0.02946085, -0.02225318])
 message: 'Optimization terminated successfully.'
    nfev: 45
     nit: 5
    njev: 5
  status: 0
 success: True
       x: array([-1.16044090e-04, -2.56376016e-05,  5.33220510e-04])
After optimization [-1.16044090e-04 -2.56376016e-05  5.33220510e-04] cost= 0.02783531839357946
CPU times: user 25.2 s, sys: 43.6 s, total: 1min 8s
Wall time: 10.7 s

dy, dx, angle  = res_global.x
rot = [[cos(angle), -sin(angle)],
       [sin(angle), cos(angle)]]
trans = [[dx],
          [dy]]

raw = numpy.vstack((pixel_coord_raw[..., 2].ravel(), #x
                    pixel_coord_raw[..., 1].ravel())) #y
aligned = numpy.dot(rot, raw) + trans

pixel_coord_ref = pixel_coord.copy()
pixel_coord_ref[...,1] = aligned[1,:].reshape(pixel_coord_ref.shape[:-1])
pixel_coord_ref[...,2] = aligned[0,:].reshape(pixel_coord_ref.shape[:-1])

displ = numpy.sqrt(((pixel_coord-pixel_coord_ref)**2).sum(axis=-1))
displ /= pilatus.pixel1 #convert in pixel units
fig, ax = subplots()
ax.hist(displ.ravel(), 100)
ax.set_title("Displacement of pixels versus the reference representation")
ax.set_xlabel("Pixel size (172µm)")

Text(0.5, 0, 'Pixel size (172µm)')

pilatus.set_pixel_corners(pixel_coord_ref)
pilatus.mask = all_masks
pilatus.save("Pilatus_ID15_ref.h5")

Conclusion

This tutorial presents the way to calibrate a module based detector using the Pilatus2M CdTe from ESRF-ID15. The HDF5 file generated is directly useable by any parts of pyFAI, the reader is invited in calibrating the rings images with the default definition and with this optimized definition and check the residual error is almost divided by a factor two.

To come back on the precision of the localization of the pixel: not all the pixel are within the specifications provided by Dectris which claims the misaliment of the modules is within one pixel.

Final comments:

Using Kabsch’s algorithm would have been faster then performing the fit …

print("Total execution time: ", time.time()-start_time)

Total execution time:  74.55711460113525