mask problem with numpy method

[picca]

Hello, while playing with the integrate1d, I was using method'"numpy" with a big mask.

I masked almost all the image (keeping only the 120 first lines) (Xpad_flat) but I get a result with plenty of noise when nothing should be computed. With the lut method it works well

[kif]

You spotted a hidden/forgotten feature !

Some people consider the mask should be -1 for masked pixels, other the people it should be zero like you multiply the given value by zero.

So the original version of pyFAI (in pure Python, a very long time ago) tried to guess if the mask was direct or if it needed to be inverted. The criteria used was the coverage of the mask. No scientist would ever use only a tiny fraction of an expensive detector. So if more then half the detector is masked out, it is obviously the mask which is wrong and needs to be inverted.

That's for the historical perspective. I agree it is a bug/nasty feature. Now pyFAI CAN enforce scientists to change their point of view and this feature should be removed.

For you information, you can use look-up-table based integration with the "nosplit" option to get the same result as numpy.

Cheers,

Jerome

[picca]

Thanks for the info ;)

[Miradin]

Hello Jerome,

No scientist would ever use only a tiny fraction of an expensive detector. So if more then half the detector is masked out, it is obviously the mask which is wrong and needs to be inverted.

I would not agree. Sometimes there is a need to mask more than half of the detector in order to catch a feature. Samples with strong anisotropy as an example. Maybe just warning text would be enough without actual mask invertion?

Nevertheless, I feel like mask is not working correctly. To make it more pronounce I have commented out line in azimuthalIntegrator.py about mask invertion (#numpy.logical_not(mask, mask)) and wrote following script:

import pyFAI
import numpy as np
import matplotlib.pyplot as plt

centx = 243.5
centy = 309.5
sdd =  1 #sample-to-detector distance in meters
num_points = 400 #number of points
units = "q_A^-1" 
wavelength = 1.54e-10 #wavelength in meters
pxsize = 0.172
dete = pyFAI.detectors.Pilatus300k()
poni1 = centx*pxsize*1e-3
poni2 = centy*pxsize*1e-3
data = np.zeros((487, 619), dtype=np.float)
data[200:300,0:619] = 50
#lets mask out everything which is not 50 using logic: valid = 1, invalid = 0
mask1 = np.zeros((487, 619), dtype=np.float)
mask1[data>0] = 1
mask1 = np.logical_not(mask1).astype(float)

ai_long = pyFAI.AzimuthalIntegrator(dist=sdd, poni1=poni1, poni2=poni2, detector=dete, wavelength=wavelength)
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask1, method="nosplit", unit=units)
plt.plot(q, I, 'o')

#lets try mask with logic: valid = 0, invalid <0
mask2 = np.zeros((487, 619), dtype=np.float)
mask2[data<50] = -2
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask2, method="nosplit", unit=units)
plt.plot(q, I, 'o')

Theoretically, I should get a constant line with "intensity" 50. But this is what I got instead:

I have tried "BBox" and "splitpixel" methods as well, but all of them behaves like this.

I was trying to apply mask with NaNs, but failed. Any chance to solve this?

With best regards, Alexander

[vallsv]

A proposal could be to use an object instead of a numpy array for masks (compatible with numpy arrays). It could be provided to our functions to bypass the aumatic normalization apply on the masks, without breaking the API. Then you could explicitly describe which bits==1 or bits==0 are masking your data.

mask2 = np.zeros((487, 619), dtype=np.float)
mask2[data<50] = -2
mymask2 = pyFAI.mask(mask2, masking_bits=1)
ai.integrate1d(data, num_points, error_model="poisson", mask=mymask2)

Magic things are never a good idea, but we can't break our API.

[Miradin]

Hello Valentin, sorry to say, I couldn't find pyFAI.mask function. Search through the whole source by phrase "masking_bits" also gave 0 result. Using pyFAI version 0.15.0.

[vallsv]

Oh sorry. I was talking about an hypothetical way to patch our code to have something which could work. Not about something which is available.

[Miradin]

Well, as far as it is not available yet, I am thinking of way around.

1. make integration using a mask
2. make integration of the mask with same binning to get the weight of each point
3. multiply...

[vallsv]

We have some other works to do around, but if you could wait 2 weeks, we could find a solution. I have no idea for a workaround.

[Miradin]

Ok, here is a code that works on synthetic data. I will test it later on some real.

import pyFAI
import numpy as np
import matplotlib.pyplot as plt

centx = 243.5
centy = 309.5
sdd =  1 #sample-to-detector distance in meters
num_points = 400 #number of points
units = "q_A^-1" 
wavelength = 1.54e-10 #wavelength in meters
pxsize = 0.172
dete = pyFAI.detectors.Pilatus300k()
poni1 = centx*pxsize*1e-3
poni2 = centy*pxsize*1e-3
data = np.zeros((487, 619), dtype=np.float)
data[200:300,0:619] = 50
#plt.imshow(data)
mask1 = np.zeros((487, 619), dtype=np.float)
mask1[data>0] = 1
mask1 = np.logical_not(mask1).astype(float)

ai_long = pyFAI.AzimuthalIntegrator(dist=sdd, poni1=poni1, poni2=poni2, detector=dete, wavelength=wavelength)

q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask1, method="nosplit", unit=units)
plt.plot(q, I, '*', label="usual masking")

data2 = np.logical_not(mask1).astype(float)
q2, I2, sigma = ai_long.integrate1d(data2, num_points, error_model="poisson", mask=mask1, method="nosplit", unit=units)

I3 = I/I2
plt.plot(q2, I3, '.', label="corrected")
plt.legend()

As I have mentioned before, to make this script to work, I have commented a line in azimuthalIntegrator.py to skip automated mask invertion.

[kif]

On Thu, 14 Jun 2018 04:22:35 -0700 Miradin notifications@github.com wrote:

Hello Jerome,

No scientist would ever use only a tiny fraction of an expensive detector. So if more then half the detector is masked out, it is obviously the mask which is wrong and needs to be inverted.

I would not agree. Sometimes there is a need to mask more than half of the detector in order to catch a feature. Samples with strong anisotropy as an example. Maybe just warning text would be enough without actual mask invertion?

Nevertheless, I feel like mask is not working correctly. To make it more pronounce I have commented out line in azimuthalIntegrator.py about mask invertion (#numpy.logical_not(mask, mask)) and wrote following script:

import pyFAI
import numpy as np
import matplotlib.pyplot as plt

centx = 243.5
centy = 309.5
sdd =  1 #sample-to-detector distance in meters
num_points = 400 #number of points
units = "q_A^-1" 
wavelength = 1.54e-10 #wavelength in meters
pxsize = 0.172
dete = pyFAI.detectors.Pilatus300k()
poni1 = centx*pxsize*1e-3
poni2 = centy*pxsize*1e-3
data = np.zeros((487, 619), dtype=np.float)
data[200:300,0:619] = 50
#lets mask out everything which is not 50 using logic: valid = 1, invalid = 0
mask1 = np.zeros((487, 619), dtype=np.float)
mask1[data>0] = 1  
mask1 = np.logical_not(mask1).astype(float)

ai_long = pyFAI.AzimuthalIntegrator(dist=sdd, poni1=poni1, poni2=poni2, detector=dete, wavelength=wavelength)
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask1, method="nosplit", unit=units)
plt.plot(q, I, 'o')

#lets try mask with logic: valid = 0, invalid <0
mask2 = np.zeros((487, 619), dtype=np.float)
mask2[data<50] = -2
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask2, method="nosplit", unit=units)
plt.plot(q, I, 'o')

Theoretically, I should get a constant line with "intensity" 50. But this is what I got instead:

I have tried "BBox" and "splitpixel" methods as well, but all of them behaves like this.

I was trying to apply mask with NaNs, but failed. Any chance to solve this?

Hi Alexander,

I think you are mixing many things ...

Apparently your detector and your dataset have not the same shape (they are transposed). Did you see the image about the axis in pyFAI ? https://pyfai.readthedocs.io/en/latest/geometry.html#detector-position You should double check you beam-center position !

You wish to obtain a "flat" curve at 50 ?

The effect you are seeing has nothing to do with the mask but results of the correction for solid-angle. Disable the solid-angle correction and your integrated theoretical curve will be flat. option: correctSolidAngle=False

For inhibiting the mask "auto-inversion", there is a work-around: None of the "advanced" integration methods suffer from this bug as it is related to the histogram of numpy. You can use "csr_nosplit" as method: the results are exactly the same and the memory footprint is not that bigger on a Pilatus 300k detector.

See the attached document where I corrected the beam center, the image sizes, import which will become invalid in 0.16, ...

Cheers, -- Jérôme Kieffer tel +33 476 882 445

[kif]

%matplotlib nbagg

import pyFAI
import numpy as np
import matplotlib.pyplot as plt
from pyFAI import detectors
from pyFAI.azimuthalIntegrator import AzimuthalIntegrator

/mntdirect/_scisoft/users/jupyter/jupy35/lib/python3.5/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

centx = 243.5
centy = 309.5
sdd =  1 #sample-to-detector distance in meters
num_points = 400 #number of points
units = "q_A^-1" 
wavelength = 1.54e-10 #wavelength in meters
pxsize = 0.172
dete = pyFAI.detectors.Pilatus300k()
poni2 = centx*pxsize*1e-3
poni1 = centy*pxsize*1e-3
data = np.zeros(dete.shape, dtype=np.float)
data[:, 200:300] = 50

print(data.shape, dete.shape)
plt.figure()
plt.imshow(data)
plt.show()

(619, 487) (619, 487)

<IPython.core.display.Javascript object>

plt.figure()
#mask with logic: valid = 0, invalid=1
mask1 = np.ones(dete.shape, dtype="int8")
mask1[data>0] = 0
ai_long = AzimuthalIntegrator(dist=sdd, poni1=poni1, poni2=poni2, detector=dete, wavelength=wavelength)
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", mask=mask1, 
                                  method="csr_nosplit", unit=units, correctSolidAngle=False)
plt.plot(q, I, '-r', label="mask1")

ai_long.reset()
dete.mask = np.logical_or(mask1, dete.mask)
q, I, sigma = ai_long.integrate1d(data, num_points, error_model="poisson", #mask=mask1, 
                                  method="csr_nosplit", unit=units, correctSolidAngle=False)
plt.plot(q, I, 'ob', label="mask in detector")
plt.legend()
plt.show()

<IPython.core.display.Javascript object>

WARNING:pyFAI.ext.splitBBoxCSR:Pixel splitting desactivated !
/mntdirect/_scisoft/users/jupyter/jupy35/lib/python3.5/site-packages/pyFAI/azimuthalIntegrator.py:2544: RuntimeWarning: invalid value encountered in true_divide
  sigma = numpy.sqrt(a) / (b * normalization_factor)
WARNING:pyFAI.ext.splitBBoxCSR:Pixel splitting desactivated !

ai_long.integrate1d?

[Miradin]

Hello Jerome, you are absolutely right about transposed shape I used in example. For that generated data, beam center coordinates could be anywhere inside the detector. Correction of data/mask shape has no effect on result, but option about solid angle correction did the trick. That's an artificial data and I completely forgot about solid angle correction implemented in pyFAI, and was surprised about values exceeds expected value. I've never tried to use that option, as it is "on" by default (like it should be). Thank you for pointing that for me. Have a nice weekend, Alexander

[vallsv]

Hi, the PR above introduce a class argument to switch off the mask inversion. If there not much problems it could be merged soon.

We think about removing this "mask inversion" feature, but in the mid time, after the merge, you could disable it on your side like that:

# Globally
pyFAI.AzimuthalIntegrator.USE_LEGACY_MASK_NORMALIZATION = False

# Or per integrators
ai = pyFAI.AzimuthalIntegrator(...blahblah...)
ai.USE_LEGACY_MASK_NORMALIZATION = False