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.
[1]:
%matplotlib nbagg
[2]:
import os
import time
start_time = time.perf_counter()
[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
[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
[5]:
from scipy.sparse import save_npz, load_npz, csr_matrix, csc_matrix, linalg
#Saved in notebook called "raytracing"
csr = load_npz("csr.npz")
[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))
# Detector: PILATUS 2M, 24-0111
# 2016-09-25T21:17:04.826413
# Pixel_size 172e-6 m x 172e-6 m
# Silicon sensor, thickness 0.000450 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.300)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Start_angle -10.00 deg.
# Angle_increment 20.00 deg.
# Omega 0.00 deg.
# Omega_increment 0.00 deg.
# Phi -10.00 deg.
# Phi_increment 20.00 deg.
# Kappa 0.0 deg.
# Oscillation_axis PHI
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0 counts
# Temperature 0.00 K
# Blower 0.0 C
# Lakeshore 0.00 K
# Detector: PILATUS 2M 24-0111
# Silicon sensor, thickness 0.00045 m
# Pixel_size 0.000172 m x 0.000172 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.200)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0.0
# Start_angle -10.0 deg.
# Angle_increment 20.0 deg.
# Kappa 0.0 deg.
# Phi -10.0 deg.
# Phi_increment 20.0 deg.
# Omega 0.0 deg.
# Omega_increment 0.0 deg.
# Oscillation_axis PHI
# Detector: PILATUS 2M, 24-0111
# 2016-09-25T21:18:25.326313
# Pixel_size 172e-6 m x 172e-6 m
# Silicon sensor, thickness 0.000450 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.300)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Start_angle -10.00 deg.
# Angle_increment 20.00 deg.
# Omega 0.00 deg.
# Omega_increment 0.00 deg.
# Phi -10.00 deg.
# Phi_increment 20.00 deg.
# Kappa 0.0 deg.
# Oscillation_axis PHI
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0 counts
# Temperature 0.00 K
# Blower 0.0 C
# Lakeshore 0.00 K
# Detector: PILATUS 2M 24-0111
# Silicon sensor, thickness 0.00045 m
# Pixel_size 0.000172 m x 0.000172 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.200)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0.0
# Start_angle -10.0 deg.
# Angle_increment 20.0 deg.
# Kappa 0.0 deg.
# Phi -10.0 deg.
# Phi_increment 20.0 deg.
# Omega 0.0 deg.
# Omega_increment 0.0 deg.
# Oscillation_axis PHI
# Detector: PILATUS 2M, 24-0111
# 2016-09-25T21:19:54.326067
# Pixel_size 172e-6 m x 172e-6 m
# Silicon sensor, thickness 0.000450 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.300)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Start_angle -10.00 deg.
# Angle_increment 20.00 deg.
# Omega 0.00 deg.
# Omega_increment 0.00 deg.
# Phi -10.00 deg.
# Phi_increment 20.00 deg.
# Kappa 0.0 deg.
# Oscillation_axis PHI
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0 counts
# Temperature 0.00 K
# Blower 0.0 C
# Lakeshore 0.00 K
# Detector: PILATUS 2M 24-0111
# Silicon sensor, thickness 0.00045 m
# Pixel_size 0.000172 m x 0.000172 m
# Exposure_time 19.9977 s
# Exposure_period 19.9977 s
# Tau = 0.0 s
# Count_cutoff 1053499 counts
# Threshold_setting: 12000.0 eV
# Gain_setting: low gain (vrf = -0.200)
# N_excluded_pixels = 32
# Excluded_pixels: badpix_mask.tif
# Flat_field: FF_p2m0111_E26000_T13000_vrf_m0p30.tif
# Trim_file: p2m0111_E26000_T13000_vrf_m0p30.bin
# Image_path: /ramdisk/
# Wavelength 0.7 A
# Detector_distance 0.157 m
# Detector_Voffset 0.0 m
# Beam_xy (542.41, 732.4) pixels
# Flux 0.0
# Start_angle -10.0 deg.
# Angle_increment 20.0 deg.
# Kappa 0.0 deg.
# Phi -10.0 deg.
# Phi_increment 20.0 deg.
# Omega 0.0 deg.
# Omega_increment 0.0 deg.
# Oscillation_axis PHI
[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.40802235349464, 11.10922455979307, 1.824532067391076, 3.948079607055263, 44886732.12361038)
(1, 22, 71.1619233914243, 12.075457104586718, 1.8280417920600356, 3.91828178168281, 38825920.361508764)
(1, 23, 30.219711117658644, 5.206798871816012, 1.862128839056976, 3.9048464388452633, 16554075.29726397)
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.
[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)
--- Logging error ---
Traceback (most recent call last):
File "/usr/lib/python3.9/logging/__init__.py", line 1082, in emit
stream.write(msg + self.terminator)
OSError: [Errno 5] Input/output error
Call stack:
File "/usr/lib/python3.9/runpy.py", line 197, in _run_module_as_main
return _run_code(code, main_globals, None,
File "/usr/lib/python3.9/runpy.py", line 87, in _run_code
exec(code, run_globals)
File "/usr/lib/python3/dist-packages/ipykernel_launcher.py", line 16, in <module>
app.launch_new_instance()
File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 845, in launch_instance
app.start()
File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 612, in start
self.io_loop.start()
File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 199, in start
self.asyncio_loop.run_forever()
File "/usr/lib/python3.9/asyncio/base_events.py", line 596, in run_forever
self._run_once()
File "/usr/lib/python3.9/asyncio/base_events.py", line 1890, in _run_once
handle._run()
File "/usr/lib/python3.9/asyncio/events.py", line 80, in _run
self._context.run(self._callback, *self._args)
File "/usr/lib/python3/dist-packages/tornado/ioloop.py", line 688, in <lambda>
lambda f: self._run_callback(functools.partial(callback, future))
File "/usr/lib/python3/dist-packages/tornado/ioloop.py", line 741, in _run_callback
ret = callback()
File "/usr/lib/python3/dist-packages/tornado/gen.py", line 814, in inner
self.ctx_run(self.run)
File "/usr/lib/python3/dist-packages/tornado/gen.py", line 775, in run
yielded = self.gen.send(value)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 362, in process_one
yield gen.maybe_future(dispatch(*args))
File "/usr/lib/python3/dist-packages/tornado/gen.py", line 234, in wrapper
yielded = ctx_run(next, result)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 265, in dispatch_shell
yield gen.maybe_future(handler(stream, idents, msg))
File "/usr/lib/python3/dist-packages/tornado/gen.py", line 234, in wrapper
yielded = ctx_run(next, result)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 540, in execute_request
self.do_execute(
File "/usr/lib/python3/dist-packages/tornado/gen.py", line 234, in wrapper
yielded = ctx_run(next, result)
File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 302, in do_execute
res = shell.run_cell(code, store_history=store_history, silent=silent)
File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 539, in run_cell
return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2877, in run_cell
result = self._run_cell(
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2923, in _run_cell
return runner(coro)
File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 68, in _pseudo_sync_runner
coro.send(None)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3146, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3338, in run_ast_nodes
if (await self.run_code(code, result, async_=asy)):
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3418, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "<ipython-input-8-c9827655079d>", line 18, in <module>
display(interactive_plot)
File "/usr/lib/python3/dist-packages/IPython/core/display.py", line 313, in display
format_dict, md_dict = format(obj, include=include, exclude=exclude)
File "/usr/lib/python3/dist-packages/IPython/core/formatters.py", line 146, in format
if self.ipython_display_formatter(obj):
File "<decorator-gen-4>", line 2, in __call__
File "/usr/lib/python3/dist-packages/IPython/core/formatters.py", line 224, in catch_format_error
r = method(self, *args, **kwargs)
File "/usr/lib/python3/dist-packages/IPython/core/formatters.py", line 918, in __call__
method()
File "/usr/lib/python3/dist-packages/ipywidgets/widgets/widget.py", line 504, in _ipython_display_
loud_error('Widget Javascript not detected. It may not be '
File "/usr/lib/python3/dist-packages/ipywidgets/widgets/widget.py", line 494, in loud_error
self.log.warn(message)
Message: 'Widget Javascript not detected. It may not be installed or enabled properly.'
Arguments: ()
Widget Javascript not detected. It may not be installed or enabled properly.
[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 12873 (2, 19, 4317963.467172676, 5.796222507131833, 1.6995853423249883, 3.3573597881674195, 41659950.47115388)
1 0.1 12250 (2, 19, 3751081.3232033723, 6.130294950181751, 1.7033086555221173, 3.4429611365152657, 36395298.00027828)
2 0.1 17393 (2, 20, 1574328.6101574465, 2.2157599000276864, 1.7393055159095367, 3.4852624162490127, 15078118.68779698)
0 0.2 11688 (2, 14, 7881767.948629593, 5.131233708982173, 1.4661572073658458, 2.702575946099344, 34991606.27574904)
1 0.2 10719 (2, 14, 6919721.86173985, 4.713794281607347, 1.4712255190981613, 2.725752576968081, 31156874.172323484)
2 0.2 14893 (2, 14, 2832176.59436549, 2.766164129458951, 1.465294038864149, 2.718024167657065, 12392353.535399606)
0 0.30000000000000004 9574 (2, 11, 10545525.37060314, 3.987369950737397, 1.3061977941858618, 2.263003702847592, 28202640.346715998)
1 0.30000000000000004 9043 (2, 10, 9349279.678903326, 11.468099752303793, 1.254986069286585, 2.228134635879879, 25507028.269693363)
2 0.30000000000000004 12141 (2, 11, 3753059.6243482213, 1.942482551631614, 1.304630392503682, 2.2296885229250427, 9872793.454374196)
0 0.4 8081 (2, 9, 12462825.01113654, 3.7306288707511173, 1.1903981931430425, 1.8837027621149782, 22441478.459732875)
1 0.4 7392 (2, 9, 11129232.747765236, 3.303110359534544, 1.1940946670588206, 1.881253238023583, 20507006.502665073)
2 0.4 9534 (2, 9, 4409327.444722604, 1.8114627013704527, 1.184041457906668, 1.8821012710573697, 7805486.6059175)
0 0.5 6498 (2, 7, 13830862.81673486, 12.863992665307908, 1.0566135480588001, 1.6323312423999208, 17878127.66292845)
1 0.5 5698 (2, 7, 12413159.338361552, 11.031231743310872, 1.061501707189753, 1.631137536109166, 16444501.716158308)
2 0.5 7015 (2, 8, 4874986.939884031, 0.9782623338632082, 1.1182593692324785, 1.627364688246967, 6193292.096289256)
0 0.6 4912 (2, 7, 14813190.14485558, 2.3120004831212486, 1.0566135480588001, 1.4278850455974332, 14360843.105051994)
1 0.6 4047 (2, 7, 13341240.24911563, 1.9907574085941135, 1.061501707189753, 1.4272524058185623, 13265388.65584298)
2 0.6 4768 (2, 7, 5208157.30400008, 1.0304064396605574, 1.0493375936706058, 1.4259276946031396, 4961402.471124767)
0 0.7 3451 (2, 6, 15528525.78473201, 4.573191036769078, 0.9838980775844784, 1.2647141424317043, 11668856.571398139)
1 0.7 2691 (2, 6, 14019916.159133047, 3.997438608193384, 0.9892525103201936, 1.2676716747197634, 10809339.462154008)
2 0.7 3040 (2, 6, 5450171.102339072, 2.021326306258463, 0.9757893362695215, 1.2647069447895245, 4023716.3339958256)
0 0.7999999999999999 2260 (2, 6, 16058574.585222537, 1.3011586432387499, 0.9838980775844784, 1.133941220345405, 9601236.8081428)
1 0.7999999999999999 1704 (2, 6, 14524191.759475194, 1.1401880783757614, 0.9892525103201936, 1.1360396528326577, 8911491.809645351)
2 0.7999999999999999 1833 (2, 6, 5629177.748824183, 0.571527814501679, 0.9757893362695215, 1.1339626600715018, 3306188.409147119)
0 0.8999999999999999 1351 (2, 5, 16458568.61127798, 5.797585932020075, 0.9053748607917161, 1.0388487654866818, 7999006.764369947)
1 0.8999999999999999 1045 (2, 5, 14905457.245277055, 5.192896410939718, 0.9127146572917646, 1.0433171669803352, 7434801.131982487)
2 0.8999999999999999 1107 (2, 5, 5764084.152363339, 2.593635270631173, 0.8966978064226785, 1.038967904747259, 2751636.4955824465)
0 0.9999999999999999 761 (2, 5, 16765847.53074498, 2.368954793955956, 0.9053748607917161, 1.0317602144326814, 6743280.737938973)
1 0.9999999999999999 634 (2, 5, 15198742.813141888, 2.126342788466719, 0.9127146572917646, 1.035460232111192, 6274138.437075575)
2 0.9999999999999999 668 (2, 5, 5867616.952693182, 1.0563363540688993, 0.8966978064226785, 1.0321485740970646, 2317861.2879105886)
0 1.0999999999999999 428 (2, 5, 17005892.485248163, 1.0309911501986544, 0.9053748607917161, 1.026418076158437, 5747090.77170436)
1 1.0999999999999999 349 (2, 5, 15428082.206907438, 0.9269451015375582, 0.9127146572917646, 1.029529185179549, 5351440.882159874)
2 1.0999999999999999 379 (2, 5, 5948434.108906353, 0.4585948285766914, 0.8966978064226785, 1.0269317855060405, 1974247.1718553898)
0 1.2 246 (2, 4, 17196335.949037977, 11.644564244075315, 0.8173110182940442, 1.022475811785219, 4947158.669781811)
1 1.2 200 (2, 4, 15610167.423623178, 10.444173218311086, 0.8273636539071435, 1.0239866977575625, 4609365.093863158)
2 1.2 233 (2, 5, 6012513.016835232, 0.21073941706814514, 0.8966978064226785, 1.0228641254225557, 1698642.681363863)
0 1.3 149 (2, 4, 17349574.26837475, 6.524931878893561, 0.8173110182940442, 1.0192354250488684, 4297295.544389076)
1 1.3 108 (2, 4, 15756764.287056616, 5.860557351764658, 0.8273636539071435, 1.0205432373813297, 4005782.1565249437)
2 1.3 150 (2, 4, 6064048.678913335, 2.9262043382263943, 0.8083995272547733, 1.0196378880051784, 1474943.4334682485)
0 1.4000000000000001 89 (2, 4, 17474466.244284548, 3.788753992733958, 0.8173110182940442, 1.016641970474952, 3763527.199018392)
1 1.4000000000000001 74 (2, 4, 15876297.184589313, 3.4068850587685646, 0.8273636539071435, 1.0177842513265012, 3509561.2700879867)
2 1.4000000000000001 99 (2, 4, 6106034.884437361, 1.697291202662571, 0.8083995272547733, 1.0170398078814604, 1291339.2257554878)
0 1.5000000000000002 56 (2, 4, 17577445.187740076, 2.2712968037052828, 0.8173110182940442, 1.014535672277159, 3320615.8589424356)
1 1.5000000000000002 44 (2, 4, 15974892.867314717, 2.044302573195794, 0.8273636539071435, 1.0155413920192156, 3097498.151464875)
2 1.5000000000000002 59 (2, 4, 6140643.388204399, 1.0166026089939373, 0.8083995272547733, 1.014919031080728, 1139076.905576041)
0 1.6000000000000003 41 (2, 4, 17663255.443435524, 1.4012098388473304, 0.8173110182940442, 1.0128027195968488, 2949608.9005856235)
1 1.6000000000000003 27 (2, 4, 16057074.708724981, 1.2621643559053994, 0.8273636539071435, 1.0136945630474439, 2752122.216704633)
2 1.6000000000000003 39 (2, 4, 6169474.387042507, 0.6267066869810065, 0.8083995272547733, 1.0131667883864373, 1011595.0003661498)
0 1.7000000000000004 32 (2, 4, 17735446.964808192, 0.8870736721959765, 0.8173110182940442, 1.0113604223600619, 2636114.4171607345)
1 1.7000000000000004 17 (2, 4, 16126230.348298585, 0.7995753771982965, 0.8273636539071435, 1.0121671030568424, 2460140.2694889517)
2 1.7000000000000004 28 (2, 4, 6193724.373025578, 0.39651233609496156, 0.8083995272547733, 1.011703253260331, 903918.1530269532)
0 1.8000000000000005 23 (2, 4, 17796712.32564079, 0.5748774391160818, 0.8173110182940442, 1.0101476169316441, 2369083.009259818)
1 1.8000000000000005 10 (2, 4, 16184931.098959597, 0.5184636865151532, 0.8273636539071435, 1.010895426586428, 2211330.1248646826)
2 1.8000000000000005 13 (2, 4, 6214300.327047811, 0.25683218763087795, 0.8083995272547733, 1.0104688792605792, 812231.0093099748)
0 1.9000000000000006 17 (2, 4, 17849119.78546627, 0.38055384693645217, 0.8173110182940442, 1.0091183238445998, 2139940.8546512546)
1 1.9000000000000006 6 (2, 4, 16235153.239121253, 0.34337385233492296, 0.8273636539071435, 1.0098112990273662, 1997749.58149476)
2 1.9000000000000006 7 (2, 4, 6231898.588182204, 0.16994197458397006, 0.8083995272547733, 1.0094185807516582, 733575.8981295343)
0 2.0000000000000004 16 (2, 4, 17894276.82018337, 0.25684050390750757, 0.8173110182940442, 1.0082374671711953, 1941967.1739861432)
1 2.0000000000000004 5 (2, 4, 16278433.50550687, 0.23184287568691414, 0.8273636539071435, 1.0088799092022736, 1813166.359277846)
2 2.0000000000000004 2 (2, 4, 6247060.133267298, 0.11465298504683204, 0.8083995272547733, 1.0085177520827633, 665636.2602974311)
0 2.1000000000000005 13 (2, 4, 17933446.47460093, 0.17644006041373272, 0.8173110182940442, 1.0074808608764898, 1769842.7805640502)
1 2.1000000000000005 5 (2, 3, 16315979.85479131, 11.281501689686154, 0.7305030202731361, 1.0081627909115904, 1652643.5407427126)
2 2.1000000000000005 1 (2, 4, 6260209.873521866, 0.07873692804819597, 0.8083995272547733, 1.0077394974264007, 606579.8992925673)
0 2.2000000000000006 11 (2, 3, 17967631.255237713, 10.02570219965494, 0.7192385713016005, 1.0068973030058297, 1619319.3782021266)
1 2.2000000000000006 2 (2, 3, 16348751.400854087, 8.6223687563185, 0.7305030202731361, 1.0074464266933758, 1512235.1961814666)
2 2.2000000000000006 0 (2, 3, 6271684.967091701, 4.262651422225273, 0.7118466192253221, 1.0071422918332718, 554944.3359663769)
0 2.3000000000000007 10 (2, 3, 17997634.385219257, 7.7461890955557084, 0.7192385713016005, 1.0063156202403392, 1486974.7872277477)
1 2.3000000000000007 1 (2, 3, 16377516.827369843, 6.663737547382244, 0.7305030202731361, 1.0068200278865433, 1388760.7426374236)
2 2.3000000000000007 0 (2, 3, 6281755.4785107225, 3.2928809427097185, 0.7118466192253221, 1.0065374097608915, 509552.08513335115)
0 2.400000000000001 10 (2, 3, 18024105.064457696, 6.046951996429845, 0.7192385713016005, 1.0058043267158263, 1370030.1246592777)
1 2.400000000000001 1 (2, 3, 16402897.576303367, 5.203197772851196, 0.7305030202731361, 1.0062691938343942, 1279636.0468112347)
2 2.400000000000001 0 (2, 3, 6290639.635306549, 2.5701385975888673, 0.7118466192253221, 1.0060060996871103, 469447.40516314365)
0 2.500000000000001 7 (2, 3, 18047572.28778381, 4.765702922172293, 0.7192385713016005, 1.0053525395601839, 1266211.7424624844)
1 2.500000000000001 1 (2, 3, 16425400.136154233, 4.101596516720589, 0.7305030202731361, 1.0057822832624919, 1182745.8257487365)
2 2.500000000000001 0 (2, 3, 6298515.222280575, 2.0252883312218133, 0.7118466192253221, 1.005536920003507, 433848.62049401563)
0 2.600000000000001 6 (2, 3, 18068470.3809244, 3.789317382926453, 0.7192385713016005, 1.0049513945432718, 1173646.1478311655)
1 2.600000000000001 1 (2, 3, 16445440.43439186, 3.261887422315954, 0.7305030202731361, 1.005349805653019, 1096346.3704464126)
2 2.600000000000001 0 (2, 3, 6305528.181777495, 1.6101539118839314, 0.7118466192253221, 1.005120565046595, 402111.8418245866)
0 2.700000000000001 4 (2, 3, 18087158.474343535, 3.037877886887301, 0.7192385713016005, 1.0045936142014424, 1090779.2614189724)
1 2.700000000000001 0 (2, 3, 16463362.446917284, 2.615481193390746, 0.7305030202731361, 1.0049639630772658, 1018990.7426285609)
2 2.700000000000001 0 (2, 3, 6311799.169936641, 1.2907107947294874, 0.7118466192253221, 1.00474940497514, 373703.11342852906)
0 2.800000000000001 4 (2, 3, 18103935.49287357, 2.4542177326042034, 0.7192385713016005, 1.0042731797129323, 1016313.8624672747)
1 2.800000000000001 0 (2, 3, 16479452.527503598, 2.1132951337113703, 0.7305030202731361, 1.00461830159302, 949470.7748642057)
2 2.800000000000001 0 (2, 3, 6317428.601653959, 1.0426267546675387, 0.7118466192253221, 1.0044171369863268, 348176.8476259038)
0 2.9000000000000012 3 (2, 3, 18119051.795134783, 1.9969656036567343, 0.7192385713016005, 1.0039850790738771, 949160.7357113308)
1 2.9000000000000012 0 (2, 3, 16493950.538080685, 1.7197958032303933, 0.7305030202731361, 1.0043074433575225, 886771.743394762)
2 2.9000000000000012 0 (2, 3, 6322500.566891307, 0.848296704148417, 0.7118466192253221, 1.0041185181907468, 325158.9919001349)
0 3.0000000000000013 2 (2, 3, 18132718.286304507, 1.6358535712353708, 0.7192385713016005, 1.003725112007631, 888400.2277912328)
1 3.0000000000000013 0 (2, 3, 16507058.565306526, 1.4089781580763365, 0.7305030202731361, 1.0040268788498157, 830036.6776785185)
2 3.0000000000000013 0 (2, 3, 6327085.8960832665, 0.6948431296752139, 0.7118466192253221, 1.003849159055634, 304333.7891439648)
0 3.1000000000000014 0 (2, 3, 18145113.60832791, 1.3485101354897395, 0.7192385713016005, 1.0034897374936211, 833251.7761661737)
1 3.1000000000000014 0 (2, 3, 16518947.799855994, 1.1616160022114215, 0.7305030202731361, 1.0037728043123038, 778538.0564906027)
2 3.1000000000000014 0 (2, 3, 6331244.57828545, 0.5727498258681804, 0.7118466192253221, 1.0036053622471368, 285433.2880282338)
damp: 3.1000000000000014
[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.
[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
[12]:
def get_angle(metadata):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return metadata["angle"]
[13]:
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
from pyFAI.gui import jupyter
[14]:
gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)
[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
[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
[17]:
#Refine the model:
gonioref2d.refine2()
Cost function before refinement: 3.2988422084873674e-08
[ 2.84547625e-01 8.86548910e-02 8.93078327e-02 5.48988055e-03
5.54612858e-03 1.74619292e-02 -2.09754960e-05]
fun: 2.0278488020441964e-08
jac: array([ 1.56336722e-07, 1.64632139e-06, 4.24729213e-06, -1.24910174e-06,
5.95437752e-07, -3.71280537e-05, -3.50152418e-07])
message: 'Optimization terminated successfully'
nfev: 49
nit: 5
njev: 5
status: 0
success: True
x: array([ 2.84379519e-01, 8.86590290e-02, 8.93059801e-02, 5.48964897e-03,
5.54684765e-03, 1.74617422e-02, -2.06534646e-05])
Cost function after refinement: 2.0278488020441964e-08
GonioParam(dist=0.2843795190151917, poni1=0.08865902895078956, poni2=0.08930598013739366, rot1=0.00548964896762195, rot2=0.005546847645358473, scale1=0.017461742189945135, scale2=-2.065346464017796e-05)
maxdelta on: dist (0) 0.2845476247809866 --> 0.2843795190151917
[17]:
array([ 2.84379519e-01, 8.86590290e-02, 8.93059801e-02, 5.48964897e-03,
5.54684765e-03, 1.74617422e-02, -2.06534646e-05])
[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()
MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)
[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()
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.
[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
[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)
[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")
pass
<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 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.807373
1000 1.6512451
1100 1.4208984
1200 1.1898193
1300 0.98669434
1400 0.8183594
1500 0.68151855
1600 0.5710449
1700 0.4821167
1800 0.40985107
1900 0.35058594
2000 0.3024292
[23]:
deconvoluted = []
for i, img in enumerate(images):
cor = deconv_MLEM(csr, img, maxiter=2000)
deconvoluted.append(cor)
<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 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.807373
1000 1.6512451
1100 1.4208984
1200 1.1898193
1300 0.98669434
1400 0.8183594
1500 0.68151855
1600 0.5710449
1700 0.4821167
1800 0.40985107
1900 0.35058594
2000 0.3024292
0 394157.62
100 159.6875
200 53.871094
300 19.878906
400 9.826172
500 4.1464844
600 2.4799194
700 1.6993408
800 1.2145691
900 0.89471436
1000 0.67407227
1100 0.5165405
1200 0.401062
1300 0.31459045
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.4853516
900 4.690918
1000 3.3764648
1100 2.5717773
1200 2.041504
1300 1.652832
1400 1.3607178
1500 1.1367188
1600 0.9614258
1700 0.8220215
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
[24]:
fig, ax = subplots(1, 3, figsize=(9,3))
for i,img in enumerate(deconvoluted):
jupyter.display(img, ax=ax[i], label=str(i))
[25]:
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
[26]:
def get_angle(metadata):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return metadata["angle"]
[27]:
gonioref2d = GoniometerRefinement.sload("id28.json", pos_function=get_angle)
[28]:
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
[29]:
# 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
[30]:
gonioref2d.refine2()
Cost function before refinement: 3.108796292554342e-08
[ 2.84547625e-01 8.86548910e-02 8.93078327e-02 5.48988055e-03
5.54612858e-03 1.74619292e-02 -2.09754960e-05]
fun: 2.1568343080973774e-08
jac: array([ 2.16702984e-07, 1.33092694e-06, 4.22263922e-06, -9.51162108e-07,
4.44783629e-07, -3.57616307e-05, -8.46252257e-08])
message: 'Optimization terminated successfully'
nfev: 49
nit: 5
njev: 5
status: 0
success: True
x: array([ 2.84401353e-01, 8.86585971e-02, 8.93077083e-02, 5.48868130e-03,
5.54699445e-03, 1.74619491e-02, -2.08352886e-05])
Cost function after refinement: 2.1568343080973774e-08
GonioParam(dist=0.2844013527629382, poni1=0.08865859710300189, poni2=0.08930770829263936, rot1=0.005488681303419845, rot2=0.005546994446350095, scale1=0.017461949080003036, scale2=-2.0835288605712943e-05)
maxdelta on: dist (0) 0.2845476247809866 --> 0.2844013527629382
[30]:
array([ 2.84401353e-01, 8.86585971e-02, 8.93077083e-02, 5.48868130e-03,
5.54699445e-03, 1.74619491e-02, -2.08352886e-05])
[31]:
#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)
[32]:
# 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, …)
[33]:
print(f"Total execution time: {time.perf_counter()-start_time:.3f} s")
Total execution time: 516.321 s