Deconvolution of the Thickness effect

This is the third part of this long journey in the thickness of the silicon sensor of the Pilatus detector: after characterizing the precise position of the Pilatus-1M on the goniometer of ID28 we noticed there are some discrepancies in the ring position and peak width as function of the position on the detector. In a second time the thickness of the detector was modeled and allowed us to define a sparse-matrix which represent the bluring of the signal with a point-spread function which actually depends on the position of the detector. This convolution can be revereted using techniques developped for inverse problems. Two pseudo-inversion algorithm have been tested, either limiting least-squares errors or with poissonian constrains.

We will now correct the images of the first notebook called “goniometer” with the sparse matrix calculated in the second one (called “raytracing”) and check if the pick-width is less chaotic and if peaks are actually thinner.

# %matplotlib widget
%matplotlib inline
# use `widget` instead of `inline` for better user-exeperience. `inline` allows to store plots into notebooks.

import os
import time
start_time = time.perf_counter()

import numpy
import fabio, pyFAI
from os.path import basename
from pyFAI.gui import jupyter
from pyFAI.calibrant import get_calibrant
from silx.resources import ExternalResources
from matplotlib.pyplot import subplots
from matplotlib.lines import Line2D
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

downloader = ExternalResources("thick", "http://www.silx.org/pub/pyFAI/testimages")
all_files = downloader.getdir("gonio_ID28.tar.bz2")
for afile in all_files:
    print(basename(afile))

gonio_ID28
det130_g45_0001p.npt
det130_g0_0001p.cbf
det130_g17_0001p.npt
det130_g0_0001p.npt
det130_g45_0001p.cbf
det130_g17_0001p.cbf

from scipy.sparse import save_npz, load_npz, csr_matrix, csc_matrix, linalg
#Saved in notebook called "raytracing"
csr = load_npz("csr.npz")

# Display all images as acquired
mask = numpy.load("mask.npy")
images = [fabio.open(i).data for i in sorted(all_files) if i.endswith("cbf")]

fig, ax = subplots(1, 3, figsize=(18, 6))
for i, img in enumerate(images):
    jupyter.display(img, ax=ax[i], label=str(i))

# Display all images once they have been corrected with least-squares minimization

fig, ax = subplots(1, 3, figsize=(18, 6))

msk = numpy.where(mask)
for i, img in enumerate(images):
    fl = img.astype("float32")
    fl[msk] = 0.0 # set masked values to 0, negative values could induce errors
    bl = fl.ravel()
    res = linalg.lsmr(csr.T, bl)
    print(res[1:])
    cor = res[0].reshape(img.shape)
    jupyter.display(cor, ax=ax[i], label=str(i))

(1, 22, 67.54663010219224, 11.13070997925794, 1.8245368369793549, 3.9484460058999806, 44886731.380958505)
(1, 22, 71.19787240560271, 12.08242996805402, 1.828034683996161, 3.9180960569347754, 38825917.68650676)
(1, 23, 30.375067157088367, 5.232227543184031, 1.8620364591591747, 3.9053502597275, 16554074.764745103)

It turns out the deconvolution is not that straight forwards and creates some negative wobbles near the rings. This phenomenon is well known in inverse methods which provide a damping factor to limit the effect. This damping factor needs to be adjusted manually to avoid this.

img = images[0]
fl = img.astype("float32")
fl[msk] = 0.0 # set masked values to 0
blured = fl.ravel()
fig, ax = subplots(1, 2, figsize=(8, 4))
im0 = ax[0].imshow(numpy.arcsinh(fl), cmap="inferno")
im1 = ax[1].imshow(numpy.arcsinh(fl), cmap="inferno")
ax[0].set_title("Original")
ax[1].set_title("Deconvolved")
def deconvol(damp):
    res = linalg.lsmr(csr.T, blured, damp=damp, x0=blured)
    restored = res[0].reshape(mask.shape)
    im1.set_data(numpy.arcsinh(restored))
    print("Number of negative pixels: %i"%(restored<0).sum())


interactive_plot = widgets.interactive(deconvol, damp=(0, 5.0))
display(interactive_plot)

#selection of the damping factor which provides no negative signal:
tot_neg = 100
damp = 0
while tot_neg>0:
    damp += 0.1
    tot_neg = 0
    deconvoluted = []
    for i, img in enumerate(images):
        fl = img.astype("float32")
        fl[msk] = 0.0 # set masked values to 0
        bl = fl.ravel()
        res = linalg.lsmr(csr.T, bl, damp=damp)
        neg=(res[0]<0).sum()
        print(i, damp, neg, res[1:])
        tot_neg += neg
        deconvoluted.append(res[0].reshape(img.shape))
print("damp:",damp)

0 0.1 12873 (2, 19, 4321777.074105922, 5.803905839281419, 1.6995848513433784, 3.3582577755005083, 41659950.06194745)
1 0.1 12250 (2, 19, 3757171.218986935, 6.13545179267901, 1.703297737531, 3.4429563190094186, 36395295.6587246)
2 0.1 17393 (2, 20, 1577644.8058258246, 2.224375524310791, 1.7392165507370456, 3.4856358181453064, 15078118.01837277)
0 0.2 11688 (2, 14, 7889537.446342854, 5.133911832121107, 1.4661895673156653, 2.702466733855008, 34991605.579435594)
1 0.2 10719 (2, 14, 6931107.204768687, 4.7184183946794445, 1.4712123116047053, 2.7256419459633983, 31156872.208067205)
2 0.2 14893 (2, 14, 2838311.780090173, 2.773752660454311, 1.4652037318093372, 2.7182470442201616, 12392352.91922903)
0 0.30000000000000004 9574 (2, 11, 10557036.27227694, 3.991240404131339, 1.3062473497865796, 2.262954193125589, 28202639.792946056)
1 0.30000000000000004 9043 (2, 10, 9364907.244862884, 11.482078593628, 1.2549593772240861, 2.2280694606631424, 25507026.746113963)
2 0.30000000000000004 12141 (2, 11, 3761519.519595857, 1.9463349070077876, 1.304574196287112, 2.229642941561226, 9872792.931475129)
0 0.4 8081 (2, 9, 12477540.081562467, 3.7332459027757765, 1.1904350950323244, 1.8837044446515787, 22441477.92458538)
1 0.4 7392 (2, 9, 11148104.812489938, 3.306798649652577, 1.1940614224199166, 1.8812512647701898, 20507005.251591522)
2 0.4 9534 (2, 9, 4419637.92908856, 1.8147912670290605, 1.183956480940489, 1.8820938989753015, 7805486.170852881)
0 0.5 6498 (2, 7, 13848158.643538255, 12.872335064169894, 1.0566656141446364, 1.632370808817917, 17878127.20875153)
1 0.5 5698 (2, 7, 12434455.716041615, 11.043407429189463, 1.0614658259363499, 1.6311538961850816, 16444500.750277305)
2 0.5 7015 (2, 8, 4886724.955397722, 0.979800043926933, 1.1181687534017324, 1.6274381364769144, 6193291.729488012)
0 0.6 4912 (2, 7, 14832502.623590384, 2.313604233809697, 1.0566656141446364, 1.4279098901784328, 14360842.637373138)
1 0.6 4047 (2, 7, 13364336.83169653, 1.9928757097329584, 1.0614658259363499, 1.4272622503605479, 13265387.8602033)
2 0.6 4768 (2, 7, 5220979.622287126, 1.0319321495647262, 1.0492603966918403, 1.4259634509194654, 4961402.176706391)
0 0.7 3451 (2, 6, 15549397.569375366, 4.576756768579681, 0.9839613211771707, 1.2647398976279356, 11668856.151128463)
1 0.7 2691 (2, 6, 14044355.899369746, 4.001389203582561, 0.9892162679551625, 1.267680033608625, 10809338.798894022)
2 0.7 3040 (2, 6, 5463815.332334352, 2.024441141325523, 0.9757176215849444, 1.2647003812458026, 4023716.0829573474)
0 0.7999999999999999 2260 (2, 6, 16080653.05498878, 1.3022158925671008, 0.9839613211771707, 1.1339599769976412, 9601236.440271325)
1 0.7999999999999999 1704 (2, 6, 14549644.60139676, 1.1412924894130865, 0.9892162679551625, 1.1360457291974302, 8911491.268508198)
2 0.7999999999999999 1833 (2, 6, 5643448.897663751, 0.5723827444096896, 0.9757176215849444, 1.1339565182550533, 3306188.2030229154)
0 0.8999999999999999 1351 (2, 5, 16481587.40573743, 5.802025301455812, 0.9054586638710889, 1.0389783960823007, 7999006.43706573)
1 0.8999999999999999 1045 (2, 5, 14931684.786443617, 5.197502843020406, 0.9126801430543429, 1.043353473858139, 7434800.689140501)
2 0.8999999999999999 1107 (2, 5, 5778838.587729241, 2.5971983123069617, 0.8966109242596816, 1.039114484165194, 2751636.3270156416)
0 0.9999999999999999 761 (2, 5, 16789606.418321703, 2.3708257037145213, 0.9054586638710889, 1.0318694569620497, 6743280.456540194)
1 0.9999999999999999 634 (2, 5, 15225571.475261556, 2.128207175861035, 0.9126801430543429, 1.0354908807978624, 6274138.06059209)
2 0.9999999999999999 668 (2, 5, 5882748.641113787, 1.0577594213654105, 0.8966109242596816, 1.0322711040251717, 2317861.1425666064)
0 1.0999999999999999 428 (2, 5, 17030240.41914222, 1.0318249422437455, 0.9054586638710889, 1.0265111477209021, 5747090.5271325065)
1 1.0999999999999999 349 (2, 5, 15455384.102332992, 0.9277503816463242, 0.9126801430543429, 1.0295553334569085, 5351440.56095096)
2 1.0999999999999999 379 (2, 5, 5963864.13808042, 0.45920313527835277, 0.8966109242596816, 1.0270355391194363, 1974247.0469713171)
0 1.2 246 (2, 4, 17221158.09127524, 11.655056407547853, 0.8174298454070358, 1.022561890341073, 4947158.483945605)
1 1.2 200 (2, 4, 15637847.046210708, 10.453365681655924, 0.8273308382836876, 1.024013715412026, 4609364.814528971)
2 1.2 233 (2, 5, 6028182.010611686, 0.21101549727639993, 0.8966109242596816, 1.0229529854181125, 1698642.5730693508)
0 1.3 149 (2, 4, 17374782.439631123, 6.530911111153145, 0.8174298454070358, 1.019310014781057, 4297295.369105349)
1 1.3 108 (2, 4, 15784749.312498542, 5.865688610012854, 0.8273308382836876, 1.0205667787432209, 4005781.914981552)
2 1.3 150 (2, 4, 6079911.413890551, 2.9289981207725035, 0.8083213310399351, 1.0197147638443347, 1474943.3366487934)
0 1.4000000000000001 89 (2, 4, 17499992.00231492, 3.7922736362752203, 0.8174298454070358, 1.0167071645200316, 3763527.0380031164)
1 1.4000000000000001 74 (2, 4, 15904532.089540195, 3.4098552362675725, 0.8273308382836876, 1.0178049205262003, 3509561.0590277207)
2 1.4000000000000001 99 (2, 4, 6122056.486261274, 1.698896604747032, 0.8083213310399351, 1.0171069152740568, 1291339.1415418575)
0 1.5000000000000002 56 (2, 4, 17603234.826595284, 2.273430545770306, 0.8174298454070358, 1.0145930979216815, 3320615.712793548)
1 1.5000000000000002 44 (2, 4, 16003334.465755593, 2.046078525383581, 0.8273308382836876, 1.0155596664987527, 3097497.9654209944)
2 1.5000000000000002 59 (2, 4, 6156796.634502143, 1.0175567462846833, 0.8083213310399351, 1.0149780826313985, 1139076.831542605)
0 1.6000000000000003 41 (2, 4, 17689266.368193828, 1.4025384807500054, 0.8174298454070358, 1.0128536575825917, 2949608.768257816)
1 1.6000000000000003 27 (2, 4, 16085688.994654555, 1.2632575975358207, 0.8273308382836876, 1.0137108236542425, 2752122.051642487)
2 1.6000000000000003 39 (2, 4, 6185737.780085941, 0.6272910785946967, 0.8083213310399351, 1.0132191259459935, 1011594.9346762811)
0 1.7000000000000004 32 (2, 4, 17761645.045293767, 0.887921383826259, 0.8174298454070358, 1.0114058923360862, 2636114.2972232825)
1 1.7000000000000004 17 (2, 4, 16154990.234904936, 0.8002662138462688, 0.8273308382836876, 1.0121874817702954, 2460140.1220960505)
2 1.7000000000000004 28 (2, 4, 6210080.748662998, 0.39688005355861533, 0.8083213310399351, 1.011749941354513, 903918.094338523)
0 1.8000000000000005 23 (2, 4, 17823069.947811823, 0.5754304408349723, 0.8174298454070358, 1.010188439412069, 2369082.9003932443)
1 1.8000000000000005 10 (2, 4, 16213814.779199643, 0.5189106914337613, 0.8273308382836876, 1.0109137012441618, 2211329.9924653866)
2 1.8000000000000005 13 (2, 4, 6230735.839631661, 0.2570692586754901, 0.8083213310399351, 1.010510772089116, 812230.9565479956)
0 1.9000000000000006 17 (2, 4, 17875614.40267762, 0.38092198424122997, 0.8174298454070358, 1.0091551656865132, 2139940.755542207)
1 1.9000000000000006 6 (2, 4, 16264142.981975865, 0.34366936164125106, 0.8273308382836876, 1.0098277754515839, 1997749.4619280791)
2 1.9000000000000006 7 (2, 4, 6248401.961468514, 0.17009821358881863, 0.8083213310399351, 1.0094563709486213, 733575.8504618482)
0 2.0000000000000004 16 (2, 4, 17920889.86534864, 0.25709016860313727, 0.8174298454070358, 1.0082708753590113, 1941967.0834511213)
1 2.0000000000000004 5 (2, 4, 16307514.761365455, 0.23204208791646103, 0.8273308382836876, 1.0088948373722588, 1813166.2508028413)
2 2.0000000000000004 2 (2, 4, 6263622.101374405, 0.11475802896945034, 0.8083213310399351, 1.0085520066962, 665636.2170065839)
0 2.1000000000000005 13 (2, 4, 17960162.535221867, 0.17661229041946044, 0.8174298454070358, 1.0075083564934988, 1769842.6976181953)
1 2.1000000000000005 5 (2, 3, 16345140.582826864, 11.292762837533477, 0.7304831518701634, 1.0081762999285853, 1652643.4412885439)
2 2.1000000000000005 1 (2, 4, 6276822.759159574, 0.07880884931028759, 0.8083213310399351, 1.0077706848111974, 606579.859818316)
0 2.2000000000000006 11 (2, 3, 17994437.441696446, 10.041258143748381, 0.7194137825276373, 1.0069199862565303, 1619319.2911782803)
1 2.2000000000000006 2 (2, 3, 16377981.557703085, 8.630970587121789, 0.7304831518701634, 1.0074587776446455, 1512235.1053357793)
2 2.2000000000000006 0 (2, 3, 6288342.360007805, 4.267213004501821, 0.7118559134639965, 1.007170666989057, 554944.300077617)
0 2.3000000000000007 10 (2, 3, 18024519.842516556, 7.758239862771203, 0.7194137825276373, 1.0063364404360613, 1486974.7093078333)
1 2.3000000000000007 1 (2, 3, 16406807.974364707, 6.670382049221476, 0.7304831518701634, 1.006831362071639, 1388760.6593225747)
2 2.3000000000000007 0 (2, 3, 6298451.987884327, 3.296405256684898, 0.7118559134639965, 1.0065634561247858, 509552.05210799567)
0 2.400000000000001 10 (2, 3, 18051060.59143047, 6.056381136939891, 0.7194137825276373, 1.0058235019735198, 1370030.0542053168)
1 2.400000000000001 1 (2, 3, 16432242.574793544, 5.20838364385342, 0.7304831518701634, 1.006279630679629, 1279635.9701282564)
2 2.400000000000001 0 (2, 3, 6307370.697269448, 2.5728897357227005, 0.7118559134639965, 1.0060300897384007, 469447.37468466524)
0 2.500000000000001 7 (2, 3, 18074590.037598077, 4.773149488950195, 0.7194137825276373, 1.005370255564005, 1266211.6782641541)
1 2.500000000000001 1 (2, 3, 16454792.908964928, 4.10568284854713, 0.7304831518701634, 1.005791924358184, 1182745.7549249923)
2 2.500000000000001 0 (2, 3, 6315276.948946852, 2.02745650919751, 0.7118559134639965, 1.005559085714414, 433848.5922897882)
0 2.600000000000001 6 (2, 3, 18095543.623578753, 3.795249163201464, 0.7194137825276373, 1.0049678102985675, 1173646.0889651652)
1 2.600000000000001 1 (2, 3, 16474875.777432902, 3.2651360363612842, 0.7304831518701634, 1.005358737935938, 1096346.304835549)
2 2.600000000000001 0 (2, 3, 6322317.2418242665, 1.6118778556085431, 0.7118559134639965, 1.0051411050715755, 402111.81565091526)
0 2.700000000000001 4 (2, 3, 18114281.40689112, 3.0426411510031306, 0.7194137825276373, 1.00460886661744, 1090779.2071267932)
1 2.700000000000001 0 (2, 3, 16492835.879176889, 2.618085221048346, 0.7304831518701634, 1.0049722613605143, 1018990.6816753438)
2 2.700000000000001 0 (2, 3, 6328612.693385167, 1.2920928586272646, 0.7118559134639965, 1.004768490374819, 373703.08908193855)
0 2.800000000000001 4 (2, 3, 18131103.08666679, 2.4580714919933397, 0.7194137825276373, 1.0042873872804914, 1016313.8121933455)
1 2.800000000000001 0 (2, 3, 16508960.170683466, 2.115398588092904, 0.7304831518701634, 1.004626030593109, 949470.7180944555)
2 2.800000000000001 0 (2, 3, 6334264.103373129, 1.043743278009079, 0.7118559134639965, 1.0044349158201462, 348176.8249268392)
0 2.9000000000000012 3 (2, 3, 18146259.672019884, 2.0001055033999595, 0.7194137825276373, 1.0039983448383105, 949160.6889750292)
1 2.9000000000000012 0 (2, 3, 16523489.019313192, 1.7215071617508866, 0.7304831518701634, 1.0043146593306773, 886771.6903944921)
2 2.9000000000000012 0 (2, 3, 6339355.8847040245, 0.8492052003811947, 0.7118559134639965, 1.0041351192288164, 325158.9706873846)
0 3.0000000000000013 2 (2, 3, 18159962.617622692, 1.6384287534816782, 0.7194137825276373, 1.0037375259880694, 888400.184191731)
1 3.0000000000000013 0 (2, 3, 16536624.937989092, 1.410379906834567, 0.7304831518701634, 1.0040336309220783, 830036.6280837732)
2 3.0000000000000013 0 (2, 3, 6343959.140351184, 0.6955873396943318, 0.7118559134639965, 1.0038646947993257, 304333.76927948784)
0 3.1000000000000014 0 (2, 3, 18172391.03201117, 1.3506352758902431, 0.7194137825276373, 1.0035013786646139, 833251.7353805491)
1 3.1000000000000014 0 (2, 3, 16548539.478812704, 1.1627714214660314, 0.7304831518701634, 1.0037791355678463, 778538.0099862082)
2 3.1000000000000014 0 (2, 3, 6348134.090582406, 0.5733633114904031, 0.7118559134639965, 1.0036199314030567, 285433.26938813186)
damp: 3.1000000000000014

#Display deconvoluted images:
fig, ax = subplots(1, 3, figsize=(18, 6))
for i,img in enumerate(deconvoluted):
    jupyter.display(img, ax=ax[i], label=str(i))

Validation using the goniometer refinement

At this point we are back to the initial state: can we fit the goniometer and check the width of the peaks to validate if it got better.

wavelength=0.6968e-10
calibrant = get_calibrant("LaB6")
calibrant.wavelength = wavelength
print(calibrant)

detector =  pyFAI.detector_factory("Pilatus1M")
detector.mask = mask

fimages = []
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[0], header={"angle":0}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[1], header={"angle":17}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[2], header={"angle":45}))

LaB6 Calibrant with 120 reflections at wavelength 6.968e-11

def get_angle(metadata):
    """Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
    return metadata["angle"]

from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
from pyFAI.gui import jupyter

gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)

for idx, fimg in enumerate(fimages):
    sg = gonioref2d.new_geometry(str(idx), image=fimg.data, metadata=fimg.header, calibrant=calibrant)
    print(sg.label, "Angle:", sg.get_position())

0 Angle: 0
1 Angle: 17
2 Angle: 45

# Display the control points freshly extracted on all 3 images:
fig,ax = subplots(1,3, figsize=(18, 6))
for lbl, sg in gonioref2d.single_geometries.items():
    idx = int(lbl)
    print(sg.extract_cp())
    jupyter.display(sg=sg, ax=ax[idx])

ControlPoints instance containing 6 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 6 groups of points:
# a ring 0: 348 points
# b ring 1: 341 points
# c ring 2: 316 points
# d ring 3: 129 points
# e ring 4: 48 points
# f ring 5: 2 points
ControlPoints instance containing 13 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 13 groups of points:
# g ring 0: 177 points
# h ring 1: 172 points
# i ring 2: 168 points
# j ring 3: 123 points
# k ring 4: 106 points
# l ring 5: 94 points
# m ring 6: 80 points
# n ring 7: 77 points
# o ring 8: 71 points
# p ring 9: 67 points
# q ring 10: 37 points
# r ring 11: 16 points
# s ring 12: 1 points
ControlPoints instance containing 26 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 26 groups of points:
# t ring 7: 37 points
# u ring 8: 54 points
# v ring 9: 66 points
# w ring 10: 64 points
# x ring 11: 62 points
# y ring 12: 61 points
# z ring 13: 57 points
#aa ring 14: 55 points
#ab ring 15: 55 points
#ac ring 16: 53 points
#ad ring 17: 53 points
#ae ring 18: 47 points
#af ring 19: 50 points
#ag ring 20: 49 points
#ah ring 21: 47 points
#ai ring 22: 47 points
#aj ring 23: 45 points
#ak ring 24: 45 points
#al ring 25: 43 points
#am ring 26: 43 points
#an ring 27: 42 points
#ao ring 28: 41 points
#ap ring 29: 41 points
#aq ring 30: 41 points
#ar ring 31: 19 points
#as ring 32: 3 points

#Refine the model:
gonioref2d.refine2()

Cost function before refinement: 3.3507684037890824e-08
[ 2.84547625e-01  8.86548910e-02  8.93078327e-02  5.48988066e-03
  5.54612874e-03  1.74619292e-02 -2.09754968e-05]
 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 2.007485066636828e-08
       x: [ 2.844e-01  8.866e-02  8.931e-02  5.490e-03  5.546e-03
            1.746e-02 -2.076e-05]
     nit: 8
     jac: [-9.441e-08  3.281e-08 -1.445e-07  2.570e-10  1.395e-07
           -5.724e-07 -4.626e-07]
    nfev: 72
    njev: 8
Cost function after refinement: 2.007485066636828e-08
GonioParam(dist=0.2843744779999335, poni1=0.08865741730339576, poni2=0.0893064450577343, rot1=0.00548956178434964, rot2=0.005546213297294622, scale1=0.017461900517488905, scale2=-2.0764946547014187e-05)
maxdelta on: dist (0) 0.2845476247783824 --> 0.2843744779999335

array([ 2.84374478e-01,  8.86574173e-02,  8.93064451e-02,  5.48956178e-03,
        5.54621330e-03,  1.74619005e-02, -2.07649465e-05])

# Create a MultiGeometry integrator from the refined geometry:


angles = []
dimages = []
for sg in gonioref2d.single_geometries.values():
    angles.append(sg.get_position())
    dimages.append(sg.image)

multigeo = gonioref2d.get_mg(angles)
multigeo.radial_range=(0, 63)
print(multigeo)
# Integrate the whole set of images in a single run:

res_mg = multigeo.integrate1d(dimages, 10000)
fig, ax = subplots(1, 2, figsize=(12,4))
ax0 = jupyter.plot1d(res_mg, label="multigeo", ax=ax[0])
ax1 = jupyter.plot1d(res_mg, label="multigeo", ax=ax[1])

# Let's focus on the inner most ring on the image taken at 45°:
for lbl, sg in gonioref2d.single_geometries.items():
    ai = gonioref2d.get_ai(sg.get_position())
    img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
    res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
    ax0.plot(*res, "--", label=lbl)
    ax1.plot(*res, "--", label=lbl)
ax1.set_xlim(29,29.3)
ax0.set_ylim(0, 4e10)
ax1.set_ylim(0, 1.5e10)
p8tth = numpy.rad2deg(calibrant.get_2th()[7])
ax1.set_title("Zoom on peak #8 at %.4f°"%p8tth)
l = Line2D([p8tth, p8tth], [0, 2e10])
ax1.add_line(l)
ax0.legend()
ax1.legend().remove()

MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)

#Peak FWHM

from scipy.interpolate import interp1d
from scipy.optimize import bisect

def calc_fwhm(integrate_result, calibrant):
    "calculate the tth position and FWHM for each peak"
    delta = integrate_result.intensity[1:] - integrate_result.intensity[:-1]
    maxima = numpy.where(numpy.logical_and(delta[:-1]>0, delta[1:]<0))[0]
    minima = numpy.where(numpy.logical_and(delta[:-1]<0, delta[1:]>0))[0]
    maxima += 1
    minima += 1
    tth = []
    FWHM = []
    for tth_rad in calibrant.get_2th():
        tth_deg = tth_rad*integrate_result.unit.scale
        if (tth_deg<=integrate_result.radial[0]) or (tth_deg>=integrate_result.radial[-1]):
            continue
        idx_theo = abs(integrate_result.radial-tth_deg).argmin()
        id0_max = abs(maxima-idx_theo).argmin()
        id0_min = abs(minima-idx_theo).argmin()
        I_max = integrate_result.intensity[maxima[id0_max]]
        I_min = integrate_result.intensity[minima[id0_min]]
        tth_maxi = integrate_result.radial[maxima[id0_max]]
        I_thres = (I_max + I_min)/2.0
        if minima[id0_min]>maxima[id0_max]:
            if id0_min == 0:
                min_lo = integrate_result.radial[0]
            else:
                min_lo = integrate_result.radial[minima[id0_min-1]]
            min_hi = integrate_result.radial[minima[id0_min]]
        else:
            if id0_min == len(minima) -1:
                min_hi = integrate_result.radial[-1]
            else:
                min_hi = integrate_result.radial[minima[id0_min+1]]
            min_lo = integrate_result.radial[minima[id0_min]]

        f = interp1d(integrate_result.radial, integrate_result.intensity-I_thres)
        tth_lo = bisect(f, min_lo, tth_maxi)
        tth_hi = bisect(f, tth_maxi, min_hi)
        FWHM.append(tth_hi-tth_lo)
        tth.append(tth_deg)
    return tth, FWHM

fig, ax = subplots()
ax.plot(*calc_fwhm(res_mg, calibrant), "o", label="multi")
for lbl, sg in gonioref2d.single_geometries.items():
    ai = gonioref2d.get_ai(sg.get_position())
    img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
    res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
    t,w = calc_fwhm(res, calibrant=calibrant)
    ax.plot(t, w,"-o", label=lbl)
ax.set_title("Peak shape as function of the angle")
ax.set_xlabel(res_mg.unit.label)
ax.legend()
pass

Conclusion on deconvolution using least-squares:

The results do not look enhanced compared to the initial fit: neither the peak position, nor the FWHM looks enhanced. Maybe there is an error somewhere.

Pseudo inverse with positivitiy constrain and poissonian noise (MLEM)

MLEM has provided similar results to least-sqares pseudo inversion with positivity constrains and the assumption of poisson noise.

# Function used for ONE MLEM iteration:

def iterMLEM_scipy(F, M, R):
    "Perform one iteration"
    #res = F * (R.T.dot(M))/R.dot(F)# / M.sum(axis=-1)
    norm = 1/R.T.dot(numpy.ones_like(F))
    cor = R.T.dot(M/R.dot(F))
    #print(norm.shape, F.shape, cor.shape)
    res = norm * F * cor
    res[numpy.isnan(res)] = 1.0
    return res

def deconv_MLEM(csr, data, thres=0.2, maxiter=1e6):
    R = csr.T
    msk = data<0
    img = data.astype("float32")
    img[msk] = 0.0 # set masked values to 0, negative values could induce errors
    M = img.ravel()
    #F0 = numpy.random.random(data.size)#M#
    F0 = R.T.dot(M)
    F1 = iterMLEM_scipy(F0, M, R)
    delta = abs(F1-F0).max()
    for i in range(int(maxiter)):
        if delta<thres:
            break
        F2 = iterMLEM_scipy(F1, M, R)
        delta = abs(F1-F2).max()
        if i%100==0:
            print(i, delta)
        F1 = F2
        i+=1
    print(i, delta)
    return F2.reshape(img.shape)

raw = images[0]
cor = deconv_MLEM(csr, images[0], maxiter=2000)

#Display one corrected image:
fig, ax = subplots(1, 2, figsize=(8,4))
im0 = ax[0].imshow(numpy.arcsinh(raw), origin="lower", cmap="inferno")
im1 = ax[1].imshow(numpy.arcsinh(cor), origin="lower", cmap="inferno")
ax[0].set_title("Original")
pass

/tmp/ipykernel_3816190/1273356881.py:6: RuntimeWarning: divide by zero encountered in divide
  norm = 1/R.T.dot(numpy.ones_like(F))
/tmp/ipykernel_3816190/1273356881.py:7: RuntimeWarning: invalid value encountered in divide
  cor = R.T.dot(M/R.dot(F))
/tmp/ipykernel_3816190/1273356881.py:9: RuntimeWarning: invalid value encountered in multiply
  res = norm * F * cor

0 261744.75
100 213.77441
200 33.487305
300 17.205078
400 7.6657715
500 4.576172
600 4.258789
700 3.6054688
800 2.5957031
900 1.807251
1000 1.651123
1100 1.4211426
1200 1.1898193
1300 0.9868164
1400 0.81848145
1500 0.68151855
1600 0.571228
1700 0.48199463
1800 0.40979004
1900 0.35058594
2000 0.3024292

deconvoluted = []
for i, img in enumerate(images):
    cor = deconv_MLEM(csr, img, maxiter=2000)
    deconvoluted.append(cor)

/tmp/ipykernel_3816190/1273356881.py:6: RuntimeWarning: divide by zero encountered in divide
  norm = 1/R.T.dot(numpy.ones_like(F))
/tmp/ipykernel_3816190/1273356881.py:7: RuntimeWarning: invalid value encountered in divide
  cor = R.T.dot(M/R.dot(F))
/tmp/ipykernel_3816190/1273356881.py:9: RuntimeWarning: invalid value encountered in multiply
  res = norm * F * cor

0 261744.75
100 213.77441
200 33.487305
300 17.205078
400 7.6657715
500 4.576172
600 4.258789
700 3.6054688
800 2.5957031
900 1.807251
1000 1.651123
1100 1.4211426
1200 1.1898193
1300 0.9868164
1400 0.81848145
1500 0.68151855
1600 0.571228
1700 0.48199463
1800 0.40979004
1900 0.35058594
2000 0.3024292
0 394157.62
100 159.6875
200 53.87256
300 19.878906
400 9.826172
500 4.1445312
600 2.4798584
700 1.6994629
800 1.2145691
900 0.89471436
1000 0.67407227
1100 0.5165558
1200 0.401062
1300 0.31455994
1400 0.25689697
1500 0.25
1600 0.25
1700 0.25
1800 0.25
1900 0.25
2000 0.25
0 144117.19
100 69.64014
200 61.936523
300 19.249756
400 16.606445
500 16.043457
600 11.894531
700 8.779297
800 6.486328
900 4.690918
1000 3.3764648
1100 2.5717773
1200 2.041504
1300 1.652832
1400 1.3605957
1500 1.1365967
1600 0.9613037
1700 0.82196045
1800 0.7095947
1900 0.6176758
2000 0.54241943

# # Re-normalization
# for img, cor in zip(images, deconvoluted):
#     msk = img>=0
#     ratio = img[msk].sum()/cor[msk].sum()
#     print(ratio)
#     cor *= ratio
#     cor[numpy.logical_not(msk)] = -1

fig, ax = subplots(1, 3, figsize=(18, 6))
for i,img in enumerate(deconvoluted):
    jupyter.display(img, ax=ax[i], label=str(i))

wavelength=0.6968e-10
calibrant = get_calibrant("LaB6")
calibrant.wavelength = wavelength
print(calibrant)

detector =  pyFAI.detector_factory("Pilatus1M")
detector.mask = mask

fimages = []
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[0], header={"angle":0}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[1], header={"angle":17}))
fimages.append(fabio.cbfimage.CbfImage(data=deconvoluted[2], header={"angle":45}))

LaB6 Calibrant with 120 reflections at wavelength 6.968e-11

def get_angle(metadata):
    """Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
    return metadata["angle"]

gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)

for idx, fimg in enumerate(fimages):
    sg = gonioref2d.new_geometry(str(idx), image=fimg.data, metadata=fimg.header, calibrant=calibrant)
    print(sg.label, "Angle:", sg.get_position())

0 Angle: 0
1 Angle: 17
2 Angle: 45

# Display the control points freshly extracted on all 3 images:

fig,ax = subplots(1,3, figsize=(18, 6))
for lbl, sg in gonioref2d.single_geometries.items():
    idx = int(lbl)
    print(sg.extract_cp())
    jupyter.display(sg=sg, ax=ax[idx])

ControlPoints instance containing 6 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 6 groups of points:
#at ring 0: 348 points
#au ring 1: 341 points
#av ring 2: 316 points
#aw ring 3: 129 points
#ax ring 4: 48 points
#ay ring 5: 2 points
ControlPoints instance containing 13 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 13 groups of points:
#az ring 0: 177 points
#ba ring 1: 172 points
#bb ring 2: 168 points
#bc ring 3: 123 points
#bd ring 4: 106 points
#be ring 5: 94 points
#bf ring 6: 80 points
#bg ring 7: 77 points
#bh ring 8: 71 points
#bi ring 9: 67 points
#bj ring 10: 37 points
#bk ring 11: 16 points
#bl ring 12: 1 points
ControlPoints instance containing 26 group of point:
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
Containing 26 groups of points:
#bm ring 7: 37 points
#bn ring 8: 54 points
#bo ring 9: 66 points
#bp ring 10: 64 points
#bq ring 11: 62 points
#br ring 12: 61 points
#bs ring 13: 57 points
#bt ring 14: 55 points
#bu ring 15: 55 points
#bv ring 16: 53 points
#bw ring 17: 53 points
#bx ring 18: 47 points
#by ring 19: 50 points
#bz ring 20: 49 points
#ca ring 21: 47 points
#cb ring 22: 47 points
#cc ring 23: 45 points
#cd ring 24: 45 points
#ce ring 25: 43 points
#cf ring 26: 43 points
#cg ring 27: 42 points
#ch ring 28: 41 points
#ci ring 29: 41 points
#cj ring 30: 41 points
#ck ring 31: 19 points
#cl ring 32: 3 points

gonioref2d.refine2()

Cost function before refinement: 3.062495458703491e-08
[ 2.84547625e-01  8.86548910e-02  8.93078327e-02  5.48988066e-03
  5.54612874e-03  1.74619292e-02 -2.09754968e-05]
 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 2.088574036367456e-08
       x: [ 2.844e-01  8.866e-02  8.931e-02  5.489e-03  5.547e-03
            1.746e-02 -2.068e-05]
     nit: 5
     jac: [ 1.727e-07 -1.658e-06  4.324e-06 -1.139e-06 -4.274e-07
           -4.582e-05  4.887e-07]
    nfev: 49
    njev: 5
Cost function after refinement: 2.088574036367456e-08
GonioParam(dist=0.2843995629872536, poni1=0.08865789000767253, poni2=0.08930799444232211, rot1=0.005488894123806088, rot2=0.005546761381912965, scale1=0.01746201634880945, scale2=-2.0683929290404354e-05)
maxdelta on: dist (0) 0.2845476247783824 --> 0.2843995629872536

array([ 2.84399563e-01,  8.86578900e-02,  8.93079944e-02,  5.48889412e-03,
        5.54676138e-03,  1.74620163e-02, -2.06839293e-05])

#Create a MultiGeometry integrator from the refined geometry:

angles = []
dimages = []
for sg in gonioref2d.single_geometries.values():
    angles.append(sg.get_position())
    dimages.append(sg.image)

multigeo = gonioref2d.get_mg(angles)
multigeo.radial_range=(0, 63)
print(multigeo)
# Integrate the whole set of images in a single run:

res_mg = multigeo.integrate1d(dimages, 10000)
fig, ax = subplots(1, 2, figsize=(12, 4))
ax0 = jupyter.plot1d(res_mg, label="multigeo", ax=ax[0])
ax1 = jupyter.plot1d(res_mg, label="multigeo", ax=ax[1])

# Let's focus on the inner most ring on the image taken at 45°:
for lbl, sg in gonioref2d.single_geometries.items():
    ai = gonioref2d.get_ai(sg.get_position())
    img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
    res = ai.integrate1d(img, 5000, unit="2th_deg", method=("full", "histogram","cython"))
    ax0.plot(*res, "--", label=lbl)
    ax1.plot(*res, "--", label=lbl)
ax1.set_xlim(29,29.3)
ax0.set_ylim(0, 1.5e12)
ax1.set_ylim(0, 7e11)
p8tth = numpy.rad2deg(calibrant.get_2th()[7])
ax1.set_title("Zoom on peak #8 at %.4f°"%p8tth)
l = Line2D([p8tth, p8tth], [0, 2e12])
ax1.add_line(l)
ax0.legend()
ax1.legend().remove()

MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)

# FWHM plot
fig, ax = subplots()
ax.plot(*calc_fwhm(res_mg, calibrant), "o", label="multi")
for lbl, sg in gonioref2d.single_geometries.items():
    ai = gonioref2d.get_ai(sg.get_position())
    img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
    res = ai.integrate1d(img, 5000, unit="2th_deg", method=("full", "histogram","cython"))
    t,w = calc_fwhm(res, calibrant=calibrant)
    ax.plot(t, w,"-o", label=lbl)
ax.set_title("Peak shape as function of the angle")
ax.set_xlabel(res_mg.unit.label)
ax.legend()
pass

Conclusion

The MLEM deconvolution provided thinner peak (~10% thinner), the fit for the geometry is in the same level of quality. Deconvoluted images are more noisy and there is a ring with zeros shadowing close to the highest intensities.

For the future (i.e. TODO): * Properly treat masked pixel as invalid ones * Provide some regulatization to induce less moise. * Implement MLEM in a more efficient way (cython, OpenCL, …)

print(f"Total execution time: {time.perf_counter()-start_time:.3f} s")

Total execution time: 185.878 s