Pixel splitting
This notebook demonstrates the layout of pixel in polar coordinates on a small detector (5x5 pixels) to demonstrate pixel splitting and pixel reorganisation.
In a second part, it displays the effect of the splitting scheme on 2D integration.
[1]:
%matplotlib nbagg
import time
import numpy
from matplotlib.pyplot import subplots
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
start_time = time.perf_counter()
[2]:
import fabio
import pyFAI.test.utilstest
from pyFAI.gui import jupyter
import pyFAI, pyFAI.units
from pyFAI.detectors import Detector
from pyFAI.azimuthalIntegrator import AzimuthalIntegrator
from pyFAI.ext import splitPixel
[3]:
# Define a dummy detector with 1mm pixel size
det = Detector(1e-3, 1e-3, max_shape=(5,5))
print(det)
Detector Detector Spline= None PixelSize= 1.000e-03, 1.000e-03 m
[4]:
def area4(a0, a1, b0, b1,c0,c1,d0,d1):
"""
Calculate the area of the ABCD polygon with 4 with corners:
A(a0,a1)
B(b0,b1)
C(c0,c1)
D(d0,d1)
:return: area, i.e. 1/2 * (AC ^ BD)
"""
return 0.5 * (((c0 - a0) * (d1 - b1)) - ((c1 - a1) * (d0 - b0)))
[5]:
# Example of code for spotting pixel on the azimuthal discontinuity: its area has not the same sign!
chiDiscAtPi = 1
pi = numpy.pi
two_pi = 2*numpy.pi
ai = AzimuthalIntegrator(1, 2.2e-3, 2.8e-3, rot3=-0.3, detector=det)
if chiDiscAtPi:
ai.setChiDiscAtPi()
else:
ai.setChiDiscAtZero()
pos = ai.array_from_unit(typ="corner", unit="r_mm", scale=True)
a = []
s = 0
ss = 0
cnt = 0
for i0 in range(pos.shape[0]):
for i1 in range(pos.shape[1]):
p = pos[i0, i1].copy()
area = area4(*p.ravel())
area2 = None
if area>=0:
az = p[:, 1].copy()
if chiDiscAtPi:
m = numpy.where(az<0)
else:
m = numpy.where(az<pi)
az[m] = two_pi + az[m]
c1 = az.mean()
if not chiDiscAtPi and c1>two_pi:
az -= two_pi
elif chiDiscAtPi and c1>pi:
az -= two_pi
print(f"Discontinuity for pixel centered at azimuth {c1}:")
for x,y in zip(p,az):
print(x, y)
p[:, 1 ] = az
area2 = area4(*p.ravel())
print(i0, i1, area, area2)
0 0 -0.34348246455192566 None
0 1 -0.45259395241737366 None
0 2 -0.578589916229248 None
0 3 -0.5334692001342773 None
0 4 -0.4045378267765045 None
1 0 -0.41383373737335205 None
1 1 -0.6470319032669067 None
1 2 -1.1334359645843506 None
1 3 -0.8771651983261108 None
1 4 -0.5334692001342773 None
Discontinuity for pixel centered at azimuth 3.312952995300293:
[ 2.807134 -2.7702851] -2.7702851
[ 2.912044 -3.1198924] -3.1198924
[1.9697715 3.0233684] -3.2598171
[ 1.811077 -2.7309353] -2.7309353
2 0 3.0264618396759033 -0.4323282837867737
Discontinuity for pixel centered at azimuth 3.2295961380004883:
[ 1.811077 -2.7309353] -2.7309353
[1.9697715 3.0233684] -3.2598171
[1.1313709 2.6561944] -3.626991
[ 0.82462114 -2.596614 ] -2.596614
2 1 4.994504928588867 -0.7384508848190308
Discontinuity for pixel centered at azimuth 3.4415926933288574:
[ 0.82462114 -2.596614 ] -2.596614
[1.1313709 2.6561944] -3.626991
[0.82462114 1.6258177 ] -4.6573677
[ 0.28284273 -0.48539817] -0.4853983
2 2 1.7914260625839233 -0.8743038177490234
2 3 -1.1334359645843506 None
2 4 -0.578589916229248 None
Discontinuity for pixel centered at azimuth 2.9282779693603516:
[ 2.912044 -3.1198924] 3.1632931
[3.3286633 2.8702552] 2.8702552
[2.5455844 2.6561944] 2.6561944
[1.9697715 3.0233684] 3.0233684
3 0 3.8964836597442627 -0.3726010322570801
3 1 -0.5192623138427734 None
3 2 -0.7384505867958069 None
3 3 -0.6470320820808411 None
3 4 -0.45259395241737366 None
4 0 -0.30272746086120605 None
4 1 -0.37260088324546814 None
4 2 -0.4323280453681946 None
4 3 -0.4138337969779968 None
4 4 -0.34348249435424805 None
[6]:
def display_geometry(pos):
fig, ax = subplots()
patches = []
for i0 in range(pos.shape[0]):
for i1 in range(pos.shape[1]):
p = pos[i0, i1].astype("float64")
splitPixel.recenter(p, chiDiscAtPi)
p = numpy.concatenate((p, [p[0]]))
ax.plot(p[:,0], p[:,1], "--")
patches.append(Polygon(p))
p = PatchCollection(patches, alpha=0.4)
colors = numpy.linspace(0, 100, len(patches))#100 * numpy.random.rand(len(patches))
p.set_array(colors)
ax.add_collection(p)
if chiDiscAtPi:
ax.plot([0,4], [-pi, -pi])
else:
ax.plot([0,4], [two_pi, two_pi])
ax.plot([0,4], [pi, pi])
ax.plot([0,4], [0, 0])
ax.set_xlabel(unit.label)
ax.set_ylabel("Azimuthal angle (rad)")
return ax
[7]:
chiDiscAtPi = 0
unit = pyFAI.units.to_unit("r_mm")
ai = AzimuthalIntegrator(1, 2.2e-3, 2.8e-3, rot3=0.5, detector=det)
if chiDiscAtPi:
ai.setChiDiscAtPi()
low = -pi
high = pi
else:
ai.setChiDiscAtZero()
low = 0
high = two_pi
pos = ai.array_from_unit(typ="corner", unit=unit, scale=True)
ax = display_geometry(pos)
over = 0
under = 0
for i0 in range(pos.shape[0]):
for i1 in range(pos.shape[1]):
p = pos[i0, i1].copy()
area = area4(*p.ravel())
area2 = None
if area>=0:
az = p[:, 1]
if chiDiscAtPi:
m = numpy.where(az<0)
else:
m = numpy.where(az<pi)
az[m] = two_pi + az[m]
c1 = az.mean()
if not chiDiscAtPi and c1>two_pi:
az -= two_pi
elif chiDiscAtPi and c1>pi:
az -= two_pi
over += (az>high).sum()
under += (az<low).sum()
# p = numpy.concatenate((p, [p[0]]))
# ax.plot(p[:,0], p[:,1], "-")
# print(i0, i1, area)
print(f"under {low:.3f}: {under} \t Above {high:.3f}: {over}")
under 0.000: 1 Above 6.283: 3
[8]:
chiDiscAtPi = 1
pi = numpy.pi
two_pi = 2*numpy.pi
ai = AzimuthalIntegrator(1, 2.2e-3, 2.8e-3, rot3=-0.4, detector=det)
if chiDiscAtPi:
ai.setChiDiscAtPi()
low = -pi
high = pi
else:
ai.setChiDiscAtZero()
low = 0
high = two_pi
pos = ai.array_from_unit(typ="corner", unit="r_mm", scale=True)
_ = display_geometry(pos)
over = 0
under = 0
for i0 in range(pos.shape[0]):
for i1 in range(pos.shape[1]):
p = pos[i0, i1].copy()
area = area4(*p.ravel())
area2 = None
if area>=0:
az = p[:, 1]
if chiDiscAtPi:
m = numpy.where(az<0)
else:
m = numpy.where(az<pi)
az[m] = two_pi + az[m]
c1 = az.mean()
print(c1)
if c1>high:
az -= two_pi
over += (az>high).sum()
under += (az<low).sum()
# print(i0, i1, area)
print(f"under {low:.3f}: {under} \t Above {high:.3f}: {over}")
3.412953
3.329596
3.5415926
3.0282776
under -3.142: 5 Above 3.142: 1
Effect of pixel splitting on 2D integration
[9]:
img = pyFAI.test.utilstest.UtilsTest.getimage("Pilatus1M.edf")
fimg = fabio.open(img)
_ = jupyter.display(fimg.data)
[10]:
# Focus on the beam stop holder ...
poni = pyFAI.test.utilstest.UtilsTest.getimage("Pilatus1M.poni")
ai = pyFAI.load(poni)
print(ai)
ai.setChiDiscAtZero()
Detector Detector Spline= None PixelSize= 1.720e-04, 1.720e-04 m
SampleDetDist= 1.583231e+00m PONI= 3.341702e-02, 4.122778e-02m rot1=0.006487 rot2= 0.007558 rot3= 0.000000 rad
DirectBeamDist= 1583.310mm Center: x=179.981, y=263.859 pix Tilt=0.571 deg tiltPlanRotation= 130.640 deg
[11]:
kwargs = {"data":fimg.data,
"npt_rad":500,
"npt_azim":500,
"unit":"r_mm",
"dummy":-2,
"delta_dummy":2,
"azimuth_range":(150,200),
"radial_range":(0,50),
}
resn = ai.integrate2d_ng(method=("no", "histogram", "cython"), **kwargs)
resb = ai.integrate2d_ng(method=("bbox", "histogram", "cython"), **kwargs)
resp = ai.integrate2d_ng(method=("pseudo", "histogram", "cython"), **kwargs)
resf = ai.integrate2d_ng(method=("full", "histogram", "cython"), **kwargs)
fig,ax = subplots(2,2, figsize=(10,10))
jupyter.plot2d(resn, ax=ax[0,0])
ax[0,0].set_title(resn.compute_engine.split("(")[1].strip(")"))
jupyter.plot2d(resb, ax=ax[0,1])
ax[0,1].set_title(resp.compute_engine.split("(")[1].strip(")"))
jupyter.plot2d(resp, ax=ax[1,0])
ax[1,0].set_title(resp.compute_engine.split("(")[1].strip(")"))
jupyter.plot2d(resf, ax=ax[1,1])
ax[1,1].set_title(resf.compute_engine.split("(")[1].strip(")"))
pass
[12]:
# Compared performances for 2D integration ...
print("Without pixel splitting")
%timeit ai.integrate2d_ng(method=("no", "histogram", "cython"), **kwargs)
print("Bounding box pixel splitting")
%timeit ai.integrate2d_ng(method=("bbox", "histogram", "cython"), **kwargs)
print("Scalledd Bounding box pixel splitting")
%timeit ai.integrate2d_ng(method=("pseudo", "histogram", "cython"), **kwargs)
print("Full pixel splitting")
%timeit ai.integrate2d_ng(method=("full", "histogram", "cython"), **kwargs)
Without pixel splitting
11.9 ms ± 321 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Bounding box pixel splitting
25 ms ± 623 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Scalledd Bounding box pixel splitting
39.2 ms ± 1.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Full pixel splitting
111 ms ± 285 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Conclusion
This tutorial presents how pixels are located in polar space and explains why pixels on the azimuthal discontinuity requires special care. It also presents a comparison between the 4 pixel splitting schemes available in pyFAI: without splitting (no), along the bounding box (bbox), a scaled down bounding box (pseudo) and finally the splitting along the edges of each pixel (full). The corresponding runtimes are also provided.
[13]:
print(f"runtime: {time.perf_counter()-start_time:.3f}s")
runtime: 25.268s