[1]:
%matplotlib widget
[2]:
import time
start_time = time.perf_counter()
import sys
print(sys.executable)
print(sys.version)
import os
os.environ["PYOPENCL_COMPILER_OUTPUT"] = "0"
/users/kieffer/.venv/py39/bin/python
3.9.5 (default, Jun  4 2021, 12:28:51)
[GCC 7.5.0]
[17]:
pix = 100e-6
shape = (1024, 1024)
npt = 1000
nimg = 1000
wl = 1e-10
I0 = 1e2
kwargs = {"npt":npt,
         "correctSolidAngle":False,
         "polarization_factor":None,
         "safe":False,
         "error_model":"Poisson",
         "method":("no", "csr", "opencl"),
         }
         # "normalization_factor": 1.0}

[25]:
import numpy
from scipy.stats import chi2 as chi2_dist
from matplotlib.pyplot import subplots
from matplotlib.colors import LogNorm
import logging
logging.basicConfig(level=logging.ERROR)
import pyFAI
print(f"pyFAI version: {pyFAI.version}")
from pyFAI.detectors import Detector
from pyFAI.azimuthalIntegrator import AzimuthalIntegrator
from pyFAI.method_registry import IntegrationMethod
from pyFAI.gui import jupyter
detector = Detector(pix, pix)
detector.shape = detector.max_shape = shape
print(detector)
flat = numpy.random.random(shape)*0.1+1
pyFAI version: 0.22.0-dev6
Detector Detector        Spline= None    PixelSize= 1.000e-04, 1.000e-04 m
[26]:
ai_init = {"dist":1.0,
           "poni1":0.0,
           "poni2":0.0,
           "rot1":-0.05,
           "rot2":+0.05,
           "rot3":0.0,
           "detector":detector,
           "wavelength":wl}
ai = AzimuthalIntegrator(**ai_init)
print(ai)
Detector Detector        Spline= None    PixelSize= 1.000e-04, 1.000e-04 m
Wavelength= 1.000000e-10m
SampleDetDist= 1.000000e+00m    PONI= 0.000000e+00, 0.000000e+00m       rot1=-0.050000  rot2= 0.050000  rot3= 0.000000 rad
DirectBeamDist= 1002.504mm      Center: x=500.417, y=501.043 pix        Tilt=4.051 deg  tiltPlanRotation= 45.036 deg
[27]:
# Generation of a "SAXS-like" curve with the shape of a lorentzian curve
unit="q_nm^-1"
q = numpy.linspace(0, ai.array_from_unit(unit=unit).max(), npt)
I = I0/(1+q**2)
fig, ax = subplots()
ax.semilogy(q, I, label="Simulated signal")
ax.set_xlabel("q ($nm^{-1}$)")
ax.set_ylabel("I (count)")
ax.set_title("SAXS-like curve with good statistical characteristics")
ax.legend()
pass
[28]:
#Reconstruction of diffusion image:
img_theo = ai.calcfrom1d(q, I, dim1_unit="q_nm^-1",
                         correctSolidAngle=True,
                         polarization_factor=None,
                         flat=flat)
kwargs["flat"] = flat
img_poisson = numpy.random.poisson(img_theo)
fig, ax = subplots(1, 2)
ax[0].imshow(img_theo, norm=LogNorm())
_=ax[0].set_title("Diffusion image (theo)")
ax[1].imshow(img_poisson, norm=LogNorm())
_=ax[1].set_title("Diffusion image (noisy)")
[33]:
factor = 2
k = kwargs.copy()
k["error_model"] = "azimuthal"
res_azim = ai.integrate1d(img_poisson, **k)
res_renorm = ai.integrate1d(img_poisson, normalization_factor=factor, **kwargs)
res_azim_renorm = ai.integrate1d(img_poisson, normalization_factor=factor, **k)
ref = ai.integrate1d(img_poisson, **kwargs)
ax = jupyter.plot1d(ref)
ax.plot(ref.radial, ref.sem, label="sem poisson unscale")
ax.plot(res_renorm.radial, res_renorm.sem*factor, "1",alpha=0.5,  label="sem poisson, scaled")
ax.plot(res_azim.radial, res_azim.sem, alpha=0.5, label="sem azimuthal, unscaled")
ax.plot(res_azim_renorm.radial, res_azim_renorm.sem*factor,"2", alpha=0.5, label="sem azimuthal, scaled")

ax.plot(ref.radial, ref.std, label="std poisson unscale")
ax.plot(res_renorm.radial, res_renorm.std*factor,"1", alpha=0.5, label="std poisson, scaled")
ax.plot(res_azim.radial, res_azim.std, alpha=0.5,label="std azimuthal, unscaled")
ax.plot(res_azim_renorm.radial, res_azim_renorm.std*factor, "2", alpha=0.5, label="std azimuthal, scaled")


ax.legend()
_=ax.set_title("Azimuthal integration")
[34]:
# factor = 2
k = kwargs.copy()
k["error_model"] = "azimuthal"
res_azim = ai.sigma_clip_ng(img_poisson, **k)
res_renorm = ai.sigma_clip_ng(img_poisson, normalization_factor=factor, **kwargs)
res_azim_renorm = ai.sigma_clip_ng(img_poisson, normalization_factor=factor, **k)
ref = ai.sigma_clip_ng(img_poisson, **kwargs)
ax = jupyter.plot1d(ref)
ax.plot(ref.radial, ref.sem, label="sem poisson unscale")
ax.plot(res_renorm.radial, res_renorm.sem*factor, alpha=0.5, label="sem poisson, scaled")
ax.plot(res_azim.radial, res_azim.sem, alpha=0.5, label="sem azimuthal, unscaled")
ax.plot(res_azim_renorm.radial, res_azim_renorm.sem*factor, alpha=0.5, label="sem azimuthal, scaled")

ax.plot(ref.radial, ref.std, label="std poisson unscale")
ax.plot(res_renorm.radial, res_renorm.std*factor, alpha=0.5, label="std poisson, scaled")
ax.plot(res_azim.radial, res_azim.std, alpha=0.5, label="std azimuthal, unscaled")
ax.plot(res_azim_renorm.radial, res_azim_renorm.std*factor, alpha=0.5, label="std azimuthal, scaled")

ax.legend()
_=ax.set_title("Sigma-clipping")
[31]:
flat
[31]:
array([[1.01637139, 1.02173552, 1.07595418, ..., 1.08856308, 1.01441441,
        1.05474796],
       [1.0415038 , 1.02462735, 1.06118086, ..., 1.05528042, 1.05411954,
        1.07374848],
       [1.07201069, 1.04516445, 1.07242547, ..., 1.00545332, 1.0417105 ,
        1.044524  ],
       ...,
       [1.08245729, 1.04522515, 1.04148592, ..., 1.07372418, 1.00892094,
        1.04526696],
       [1.05300915, 1.06319217, 1.01801696, ..., 1.00772934, 1.09840442,
        1.00121761],
       [1.03384891, 1.07645927, 1.02803138, ..., 1.073219  , 1.01400316,
        1.09369794]])
[ ]: