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.
In [1]:
%matplotlib nbagg
In [2]:
import os
import time
start_time = time.time()
In [3]:
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
In [4]:
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
In [5]:
from scipy.sparse import save_npz, load_npz, csr_matrix, csc_matrix, linalg
#Saved in notebook called "raytracing"
csr = load_npz("csr.npz")
In [6]:
# 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=(9,3))
for i, img in enumerate(images):
jupyter.display(img, ax=ax[i], label=str(i))
In [7]:
# Display all images once they have been corrected with least-squares minimization
fig, ax = subplots(1, 3, figsize=(9,3))
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.53766751120409, 11.129642745894946, 1.8245371636630272, 3.948378116039215, 44886758.46784883)
(1, 22, 71.19252340495744, 12.081752507142443, 1.8280343579375222, 3.918061149282325, 38825906.04977599)
(1, 23, 30.369625768691204, 5.231363709428408, 1.8620379956606512, 3.905307157273793, 16554013.80570858)
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.
In [8]:
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)
interactive(children=(FloatSlider(value=2.5, description='damp', max=5.0), Output()), _dom_classes=('widget-in…
In [9]:
#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 12872 (2, 19, 4321779.35150566, 5.803641113519436, 1.6995850617640857, 3.3582417538199376, 41659969.98554161)
1 0.1 12250 (2, 19, 3757169.8995306906, 6.135248689092134, 1.7032968276952476, 3.442943267728284, 36395286.6620507)
2 0.1 17393 (2, 20, 1577640.5341739412, 2.2241199080898126, 1.739217853033561, 3.485600761960406, 15078086.2824897)
0 0.2 11689 (2, 14, 7889540.65802813, 5.133807239271567, 1.466189729926675, 2.7024543620078028, 34991613.76254249)
1 0.2 10719 (2, 14, 6931105.172705148, 4.7182678428707705, 1.4712118254370397, 2.7256366365716107, 31156865.95377047)
2 0.2 14893 (2, 14, 2838307.0872005536, 2.7735637679209546, 1.4652046267714205, 2.7182425846299467, 12392344.981435211)
0 0.30000000000000004 9574 (2, 11, 10557039.34545076, 3.9912227744921145, 1.3062471672647964, 2.262955649019415, 28202640.945378628)
1 0.30000000000000004 9043 (2, 10, 9364905.024348436, 11.481832262884733, 1.2549590447570944, 2.2280699454584427, 25507023.638021868)
2 0.30000000000000004 12141 (2, 11, 3761515.5620772634, 1.9462379637685716, 1.3045749273103469, 2.2296428063638745, 9872793.073419824)
0 0.4 8081 (2, 9, 12477542.683871275, 3.7332217465162136, 1.1904350950323244, 1.8837040069820077, 22441476.459564008)
1 0.4 7392 (2, 9, 11148102.743700482, 3.3067429654305527, 1.1940612227496195, 1.8812516614353345, 20507004.367669377)
2 0.4 9534 (2, 9, 4419634.770566149, 1.8147208128711052, 1.1839571857508362, 1.8820936392361118, 7805488.867632971)
0 0.5 6498 (2, 7, 13848160.766432375, 12.872262447930522, 1.056665726961105, 1.6323711374276255, 17878124.950661585)
1 0.5 5699 (2, 7, 12434453.925741194, 11.043285919766197, 1.061465769783203, 1.6311527529962528, 16444501.00257793)
2 0.5 7015 (2, 8, 4886722.454829522, 0.979766475272105, 1.11816965959639, 1.6274367495221003, 6193295.071873871)
0 0.6 4912 (2, 7, 14832504.344398007, 2.3135919536183502, 1.056665726961105, 1.427910148735896, 14360840.296283662)
1 0.6 4047 (2, 7, 13364335.327901691, 1.9928542874095296, 1.061465769783203, 1.4272615124530288, 13265388.592819687)
2 0.6 4768 (2, 7, 5220977.629033409, 1.0319042322071525, 1.049261078367721, 1.4259628536520232, 4961405.468468706)
0 0.7 3451 (2, 6, 15549398.97210384, 4.576723349921045, 0.9839615029057753, 1.2647394559310077, 11668853.971670033)
1 0.7 2691 (2, 6, 14044354.646574955, 4.001336977544828, 0.9892160269374843, 1.2676799039844624, 10809339.68532789)
2 0.7 3040 (2, 6, 5463813.724084949, 2.0243967600633885, 0.9757181713768514, 1.2647002172393345, 4023719.085641064)
0 0.7999999999999999 2260 (2, 6, 16080654.210058544, 1.302206612960886, 0.9839615029057753, 1.1339596776754033, 9601234.499255516)
1 0.7999999999999999 1704 (2, 6, 14549643.555596681, 1.1412774956890175, 0.9892160269374843, 1.1360456401101637, 8911492.15813172)
2 0.7999999999999999 1833 (2, 6, 5643447.58223499, 0.5723704440227377, 0.9757181713768514, 1.1339564121991543, 3306190.8648090656)
0 0.8999999999999999 1351 (2, 5, 16481588.367830684, 5.801990783829077, 0.90545869678515, 1.038978183742405, 7999004.730084728)
1 0.8999999999999999 1045 (2, 5, 14931683.907623399, 5.19744533510426, 0.912679816517955, 1.0433532186310783, 7434801.520768994)
2 0.8999999999999999 1107 (2, 5, 5778837.497164339, 2.5971514652494854, 0.8966114893200743, 1.039114381836669, 2751638.6587414625)
0 0.9999999999999999 761 (2, 5, 16789607.22888084, 2.3708116466893925, 0.90545869678515, 1.0318692834712748, 6743278.965397153)
1 0.9999999999999999 634 (2, 5, 15225570.730472477, 2.1281834623003464, 0.912679816517955, 1.0354906769480396, 6274138.816129405)
2 0.9999999999999999 668 (2, 5, 5882747.725411044, 1.0577405584624433, 0.8966114893200743, 1.0322710220962337, 2317863.179796034)
0 1.0999999999999999 428 (2, 5, 17030241.109451465, 1.0318188409636586, 0.90545869678515, 1.0265110034106488, 5747089.222152081)
1 1.0999999999999999 349 (2, 5, 15455383.465454841, 0.9277399854902225, 0.912679816517955, 1.02955516705276, 5351441.238055776)
2 1.0999999999999999 379 (2, 5, 5963863.3601903515, 0.4591950175373597, 0.8966114893200743, 1.0270354721669281, 1974248.8293181416)
0 1.2 246 (2, 4, 17221158.68506743, 11.655045380956649, 0.8174296995727469, 1.0225618622515276, 4947157.3378323335)
1 1.2 200 (2, 4, 15637846.496804776, 10.453251889498636, 0.8273302259051267, 1.0240137734462216, 4609365.418140597)
2 1.2 233 (2, 5, 6028181.3427790245, 0.21101179228079514, 0.8966114893200743, 1.0229529297366875, 1698644.138311225)
0 1.3 149 (2, 4, 17374782.95507173, 6.530904821138379, 0.8174296995727469, 1.0193099905429288, 4297294.358381241)
1 1.3 108 (2, 4, 15784748.83459533, 5.8656242641779235, 0.8273302259051267, 1.0205668300020216, 4005782.4531853963)
2 1.3 150 (2, 4, 6079910.8350604735, 2.928954729219734, 0.808321736603232, 1.0197147168348824, 1474944.7177167924)
0 1.4000000000000001 89 (2, 4, 17499992.453455716, 3.7922699319513766, 0.8174296995727469, 1.0167071434076267, 3763526.142293161)
1 1.4000000000000001 74 (2, 4, 15904531.670609849, 3.409817596831565, 0.8273302259051267, 1.0178049660305932, 3509561.539649996)
2 1.4000000000000001 99 (2, 4, 6122055.980254809, 1.698871515444936, 0.808321736603232, 1.0171068750683203, 1291340.366034497)
0 1.5000000000000002 56 (2, 4, 17603235.22442668, 2.2734282996870334, 0.8174296995727469, 1.0145930793774571, 3320614.9148289496)
1 1.5000000000000002 44 (2, 4, 16003334.095898438, 2.046055824293337, 0.8273302259051267, 1.0155597070975697, 3097498.395944294)
2 1.5000000000000002 59 (2, 4, 6156796.188729298, 1.0175417576053376, 0.808321736603232, 1.0149780478562425, 1139077.9227852263)
0 1.6000000000000003 41 (2, 4, 17689266.72140641, 1.4025370821602559, 0.8174296995727469, 1.0128536411720097, 2949608.0539421244)
1 1.6000000000000003 27 (2, 4, 16085688.665984478, 1.2632435223356266, 0.8273302259051267, 1.0137108600521891, 2752122.4385273485)
2 1.6000000000000003 39 (2, 4, 6185737.384630948, 0.6272818582146866, 0.808321736603232, 1.0132190955717184, 1011595.9118964572)
0 1.7000000000000004 32 (2, 4, 17761645.36083233, 0.8879204915768787, 0.8174296995727469, 1.0114058777159978, 2636113.654651129)
1 1.7000000000000004 17 (2, 4, 16154989.941085158, 0.8002572656161739, 0.8273302259051267, 1.0121875812380794, 2460140.471015919)
2 1.7000000000000004 28 (2, 4, 6210080.395627306, 0.3968742303173911, 0.808321736603232, 1.0117499145949043, 903918.9735781325)
0 1.8000000000000005 23 (2, 4, 17823070.231281802, 0.5754298588698396, 0.8174296995727469, 1.0101884263082217, 2369082.3198262774)
1 1.8000000000000005 10 (2, 4, 16213814.515093967, 0.5189048717740923, 0.8273302259051267, 1.010913790272343, 2211330.308483282)
2 1.8000000000000005 13 (2, 4, 6230735.522656641, 0.2570654924916575, 0.808321736603232, 1.010510748334062, 812231.751225542)
0 1.9000000000000006 17 (2, 4, 17875614.65864606, 0.3809215968934904, 0.8174296995727469, 1.0091551538772998, 2139940.228755556)
1 1.9000000000000006 6 (2, 4, 16264142.7433859, 0.3436654974690589, 0.8273302259051267, 1.009827855588924, 1997749.7491613624)
2 1.9000000000000006 7 (2, 4, 6248401.675384593, 0.1700957247677691, 0.808321736603232, 1.0094563497177258, 733576.5717032134)
0 2.0000000000000004 16 (2, 4, 17920890.097569425, 0.2570899059609757, 0.8174296995727469, 1.0082708646635523, 1941966.6035315564)
1 2.0000000000000004 5 (2, 4, 16307514.544831507, 0.23203947312535816, 0.8273302259051267, 1.0088949098777877, 1813166.5128119905)
2 2.0000000000000004 2 (2, 4, 6263621.841940145, 0.11475635171777179, 0.808321736603232, 1.008551987606203, 665636.8741940791)
0 2.1000000000000005 13 (2, 4, 17960162.746805888, 0.17661210927216073, 0.8174296995727469, 1.0075083467627906, 1769842.2588014384)
1 2.1000000000000005 5 (2, 3, 16345140.385475969, 11.292691005166366, 0.730482580696584, 1.008176365296541, 1652643.6811122294)
2 2.1000000000000005 1 (2, 4, 6276822.522865413, 0.07880769857340576, 0.808321736603232, 1.0077706675526885, 606580.4608767419)
0 2.2000000000000006 11 (2, 3, 17994437.63524152, 10.041272288274898, 0.7194136168242463, 1.006919989061064, 1619318.8885478796)
1 2.2000000000000006 2 (2, 3, 16377981.377132151, 8.630915551211535, 0.730482580696584, 1.0074588373578275, 1512235.3255666452)
2 2.2000000000000006 0 (2, 3, 6288342.143926829, 4.26717350350943, 0.7118561646579522, 1.0071706506826414, 554944.8516864269)
0 2.3000000000000007 10 (2, 3, 18024520.020207558, 7.758250761973735, 0.7194136168242463, 1.0063364429756794, 1486974.3386557356)
1 2.3000000000000007 1 (2, 3, 16406807.808549903, 6.67033942360869, 0.730482580696584, 1.006831416828166, 1388760.8622023384)
2 2.3000000000000007 0 (2, 3, 6298451.789557654, 3.2963747648126747, 0.7118561646579522, 1.0065634412899156, 509552.5599771439)
0 2.400000000000001 10 (2, 3, 18051060.75511721, 6.0563896250851235, 0.7194136168242463, 1.005823504284504, 1370029.7119637202)
1 2.400000000000001 1 (2, 3, 16432242.42201949, 5.208350297807859, 0.730482580696584, 1.0062796810679773, 1279636.1575706126)
2 2.400000000000001 0 (2, 3, 6307370.514617308, 2.5728659518542454, 0.7118561646579522, 1.006030076182735, 469447.8437051946)
0 2.500000000000001 7 (2, 3, 18074590.18885645, 4.773156164533405, 0.7194136168242463, 1.0053702576762722, 1266211.3613354408)
1 2.500000000000001 1 (2, 3, 16454792.767768912, 4.1056565183932126, 0.730482580696584, 1.0057919708777772, 1182745.9285644833)
2 2.500000000000001 0 (2, 3, 6315276.78019897, 2.0274377781113087, 0.7118561646579522, 1.0055590732778876, 433849.0266568345)
0 2.600000000000001 6 (2, 3, 18095543.76375908, 3.795254461162921, 0.7194136168242463, 1.0049678122370065, 1173645.7947115519)
1 2.600000000000001 1 (2, 3, 16474875.64656059, 3.2651150655572225, 0.730482580696584, 1.005358781013287, 1096346.4661202533)
2 2.600000000000001 0 (2, 3, 6322317.085464899, 1.6118629715557997, 0.7118561646579522, 1.005141093619891, 402112.2189998038)
0 2.700000000000001 4 (2, 3, 18114281.53715672, 3.042645391233424, 0.7194136168242463, 1.0046088684029255, 1090778.9332434058)
1 2.700000000000001 0 (2, 3, 16492835.757546933, 2.6180683836853356, 0.730482580696584, 1.0049723013619878, 1018990.8318481956)
2 2.700000000000001 0 (2, 3, 6328612.548109342, 1.2920809329204606, 0.7118561646579522, 1.0047684797944505, 373703.4645613419)
0 2.800000000000001 4 (2, 3, 18131103.208025, 2.4580749124001153, 0.7194136168242463, 1.0042873889306652, 1016313.5566623267)
1 2.800000000000001 0 (2, 3, 16508960.05735901, 2.1153849673990144, 0.730482580696584, 1.0046260678351222, 949470.8582472362)
2 2.800000000000001 0 (2, 3, 6334263.968051525, 1.0437336484617878, 0.7118561646579522, 1.0044349060144628, 348177.17528059415)
0 2.9000000000000012 3 (2, 3, 18146259.785347123, 2.000108282763343, 0.7194136168242463, 1.0039983463681816, 949160.450036012)
1 2.9000000000000012 0 (2, 3, 16523488.913478717, 1.7214960653848372, 0.730482580696584, 1.0043146940878862, 886771.821469662)
2 2.9000000000000012 0 (2, 3, 6339355.758354612, 0.8491973685314158, 0.7118561646579522, 1.0041351101150462, 325159.2983181391)
0 3.0000000000000013 2 (2, 3, 18159962.7236849, 1.6384310274564113, 0.7194136168242463, 1.0037375274104927, 888399.960305099)
1 3.0000000000000013 0 (2, 3, 16536624.838931581, 1.41037080711843, 0.730482580696584, 1.0040336634340188, 830036.7509231202)
2 3.0000000000000013 0 (2, 3, 6343959.022116123, 0.6955809267366274, 0.7118561646579522, 1.0038646863061453, 304334.07629622886)
0 3.1000000000000014 0 (2, 3, 18172391.131480563, 1.3506371483450708, 0.7194136168242463, 1.0035013799906614, 833251.5251811101)
1 3.1000000000000014 0 (2, 3, 16548539.385906167, 1.162763912737771, 0.730482580696584, 1.0037791660443536, 778538.1253321602)
2 3.1000000000000014 0 (2, 3, 6348133.979709152, 0.5733580269750035, 0.7118561646579522, 1.0036199234686751, 285433.5576554336)
damp: 3.1000000000000014
In [10]:
#Display deconvoluted images:
fig, ax = subplots(1, 3, figsize=(9,3))
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.
In [11]:
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
In [12]:
def get_angle(metadata):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return metadata["angle"]
In [13]:
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
from pyFAI.gui import jupyter
WARNING:silx.opencl.common:Unable to import pyOpenCl. Please install it from: http://pypi.python.org/pypi/pyopencl
In [14]:
gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)
In [15]:
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
In [16]:
# Display the control points freshly extracted on all 3 images:
fig,ax = subplots(1,3, figsize=(9,3))
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
WARNING:matplotlib.legend:No handles with labels found to put in legend.
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
WARNING:matplotlib.legend:No handles with labels found to put in legend.
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: 2 points
In [17]:
#Refine the model:
gonioref2d.refine2()
Cost function before refinement: 3.3673221420529115e-08
[ 2.84547605e-01 8.86540301e-02 8.93070277e-02 5.48717450e-03
5.54901577e-03 1.74619285e-02 -2.09889821e-05]
fun: 2.0154279102942266e-08
jac: array([-2.73861446e-07, 2.08617005e-07, -5.15585939e-07, -8.73413399e-08,
2.27782748e-07, -1.43998916e-06, -1.32582441e-06])
message: 'Optimization terminated successfully.'
nfev: 80
nit: 8
njev: 8
status: 0
success: True
x: array([ 2.84372950e-01, 8.86570506e-02, 8.93058252e-02, 5.48758195e-03,
5.54888418e-03, 1.74617893e-02, -2.08430380e-05])
Cost function after refinement: 2.0154279102942266e-08
GonioParam(dist=0.28437295014349545, poni1=0.08865705061324418, poni2=0.08930582517579298, rot1=0.005487581950208836, rot2=0.005548884181614612, scale1=0.017461789295730078, scale2=-2.0843037974184558e-05)
maxdelta on: dist (0) 0.2845476046521533 --> 0.28437295014349545
Out[17]:
array([ 2.84372950e-01, 8.86570506e-02, 8.93058252e-02, 5.48758195e-03,
5.54888418e-03, 1.74617893e-02, -2.08430380e-05])
In [18]:
# 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()
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.20.
File "/users/kieffer/.venv/py38/lib/python3.8/site-packages/pyFAI/multi_geometry.py", line 142, in integrate1d
res = ai.integrate1d(data, npt=npt,
MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.20.
File "<ipython-input-18-6986d308b140>", line 24, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
In [19]:
#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()
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.20.
File "<ipython-input-19-dc6d9f0ee241>", line 51, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
Out[19]:
<matplotlib.legend.Legend at 0x7f23248f3100>
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.
In [20]:
# 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
In [21]:
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)
In [22]:
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")
<ipython-input-20-445809682616>:6: RuntimeWarning: divide by zero encountered in true_divide
norm = 1/R.T.dot(numpy.ones_like(F))
<ipython-input-20-445809682616>:7: RuntimeWarning: invalid value encountered in true_divide
cor = R.T.dot(M/R.dot(F))
<ipython-input-20-445809682616>:9: RuntimeWarning: invalid value encountered in multiply
res = norm * F * cor
0 261749.75
100 213.72461
200 33.487305
300 17.201172
400 7.6657715
500 4.576172
600 4.2441406
700 3.5986328
800 2.6015625
900 1.8101807
1000 1.6519775
1100 1.4206543
1200 1.1889648
1300 0.98583984
1400 0.8175049
1500 0.68078613
1600 0.5706177
1700 0.4817505
1800 0.4093628
1900 0.3505249
2000 0.3024292
Out[22]:
Text(0.5, 1.0, 'Original')
In [23]:
deconvoluted = []
for i, img in enumerate(images):
cor = deconv_MLEM(csr, img, maxiter=2000)
deconvoluted.append(cor)
0 261749.75
<ipython-input-20-445809682616>:6: RuntimeWarning: divide by zero encountered in true_divide
norm = 1/R.T.dot(numpy.ones_like(F))
<ipython-input-20-445809682616>:7: RuntimeWarning: invalid value encountered in true_divide
cor = R.T.dot(M/R.dot(F))
<ipython-input-20-445809682616>:9: RuntimeWarning: invalid value encountered in multiply
res = norm * F * cor
100 213.72461
200 33.487305
300 17.201172
400 7.6657715
500 4.576172
600 4.2441406
700 3.5986328
800 2.6015625
900 1.8101807
1000 1.6519775
1100 1.4206543
1200 1.1889648
1300 0.98583984
1400 0.8175049
1500 0.68078613
1600 0.5706177
1700 0.4817505
1800 0.4093628
1900 0.3505249
2000 0.3024292
0 394156.12
100 159.75
200 53.876465
300 19.886719
400 9.824219
500 4.1445312
600 2.481018
700 1.7006226
800 1.2156067
900 0.89572144
1000 0.6750183
1100 0.5174103
1200 0.4018402
1300 0.31533813
1400 0.2565918
1500 0.25
1600 0.25
1700 0.25
1800 0.25
1900 0.25
2000 0.25
0 144116.5
100 69.64014
200 61.948242
300 19.260742
400 16.578613
500 16.04004
600 11.899902
700 8.777344
800 6.4873047
900 4.6914062
1000 3.3779297
1100 2.5742188
1200 2.043579
1300 1.654541
1400 1.3623047
1500 1.1380615
1600 0.9626465
1700 0.8230591
1800 0.71051025
1900 0.6185913
2000 0.54333496
In [24]:
# # 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
In [25]:
fig, ax = subplots(1, 3, figsize=(9,3))
for i,img in enumerate(deconvoluted):
jupyter.display(img, ax=ax[i], label=str(i))
In [26]:
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
In [27]:
def get_angle(metadata):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return metadata["angle"]
In [28]:
gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)
In [29]:
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
In [30]:
# Display the control points freshly extracted on all 3 images:
fig,ax = subplots(1,3, figsize=(9,3))
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
WARNING:matplotlib.legend:No handles with labels found to put in legend.
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
WARNING:matplotlib.legend:No handles with labels found to put in legend.
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
In [31]:
gonioref2d.refine2()
Cost function before refinement: 3.037350271380645e-08
[ 2.84547605e-01 8.86540301e-02 8.93070277e-02 5.48717450e-03
5.54901577e-03 1.74619285e-02 -2.09889821e-05]
fun: 2.120955314965296e-08
jac: array([ 1.29654915e-07, 1.01045028e-07, 2.88223149e-06, -7.07675797e-07,
4.29006277e-08, -2.74859755e-05, -7.14978527e-08])
message: 'Optimization terminated successfully.'
nfev: 62
nit: 6
njev: 6
status: 0
success: True
x: array([ 2.84404213e-01, 8.86563321e-02, 8.93061609e-02, 5.48644529e-03,
5.54944418e-03, 1.74618742e-02, -2.05693817e-05])
Cost function after refinement: 2.120955314965296e-08
GonioParam(dist=0.2844042128582021, poni1=0.08865633212302218, poni2=0.08930616086606516, rot1=0.005486445293758012, rot2=0.005549444178278475, scale1=0.017461874170621113, scale2=-2.056938172161214e-05)
maxdelta on: dist (0) 0.2845476046521533 --> 0.2844042128582021
Out[31]:
array([ 2.84404213e-01, 8.86563321e-02, 8.93061609e-02, 5.48644529e-03,
5.54944418e-03, 1.74618742e-02, -2.05693817e-05])
In [32]:
#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, 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)
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.20.
File "<ipython-input-32-93245475420e>", line 23, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
In [33]:
# 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="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()
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.20.
File "<ipython-input-33-6b18d6c981da>", line 7, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
Out[33]:
<matplotlib.legend.Legend at 0x7f232429cfd0>
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, ...)
In [34]:
print(f"Total execution time: {time.time()-start_time} s")
Total execution time: 186.50916409492493 s
In [ ]: