Minor clean up and testing using a low-end GPU

This commit is contained in:
Jerome Kieffer 2012-06-14 23:32:36 +02:00
parent 22b7745cb6
commit 70a13047bf
8 changed files with 7559 additions and 5578 deletions

4
.gitignore vendored
View File

@ -25,3 +25,7 @@ pip-log.txt
#Mr Developer #Mr Developer
.mr.developer.cfg .mr.developer.cfg
#Shared libraries
*.so

View File

@ -1 +0,0 @@
#error Do not use this file, it is the result of a failed Cython compilation.

File diff suppressed because it is too large Load Diff

View File

@ -36,6 +36,8 @@ cimport numpy
import numpy import numpy
import threading import threading
import cython import cython
from stdlib cimport free, malloc
from libc.string cimport memcpy
@ -45,9 +47,11 @@ cdef class Integrator1d:
""" """
cdef ocl_xrpd1d.ocl_xrpd1D_fullsplit* cpp_integrator cdef ocl_xrpd1d.ocl_xrpd1D_fullsplit* cpp_integrator
cdef char* _devicetype cdef char* _devicetype
cdef int _nBins, _nData, _platformid, _devid, _tth_id, _dtth_id, _mask_id, _solid_angle_id cdef char* filename
cdef int _nBins, _nData, _platformid, _devid,
cdef cpp_bool _useFp64 cdef cpp_bool _useFp64
cdef float tth_min, tth_max,_tth_min, _tth_max cdef float tth_min, tth_max,_tth_min, _tth_max
cdef float* ctth_out
def __cinit__(self, filename=None): def __cinit__(self, filename=None):
""" """
@ -55,28 +59,26 @@ cdef class Integrator1d:
""" """
self._nBins = -1 self._nBins = -1
self._nData = -1 self._nData = -1
self._platformid = -1
self._devid = -1
self._useFp64 = False self._useFp64 = False
self._tth_id = 0
self._dtth_id = 0
self._mask_id = 0
self._solid_angle_id = 0
self._devicetype = "gpu" self._devicetype = "gpu"
if filename is None: if filename is None:
self.cpp_integrator = new ocl_xrpd1d.ocl_xrpd1D_fullsplit() self.cpp_integrator = new ocl_xrpd1d.ocl_xrpd1D_fullsplit()
else: else:
name = str(filename) self.filename = filename
self.cpp_integrator = new ocl_xrpd1d.ocl_xrpd1D_fullsplit(name) self.cpp_integrator = new ocl_xrpd1d.ocl_xrpd1D_fullsplit(self.filename)
def __init__(self,filename=None):
"""
Python constructor
"""
self.filename = filename
self.tth_out = None
def __dealloc__(self): def __dealloc__(self):
if self.ctth_out:
free(self.ctth_out)
del self.cpp_integrator del self.cpp_integrator
def __repr__(self):
return os.linesep.join(["Cython wrapper for ocl_xrpd1d.ocl_xrpd1D_fullsplit C++ class. Logging in %s"%self.filename,
"device: %s, platform %s device %s 64bits:%s image size: %s histogram size: %s"%(self._devicetype,self._platformid,self._devid, self._useFp64, self._nData,self._nBins),
",\t ".join(["%s: %s"%(k,v) for k,v in self.get_status().items()])])
@cython.cdivision(True) @cython.cdivision(True)
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.wraparound(False) @cython.wraparound(False)
@ -86,12 +88,15 @@ cdef class Integrator1d:
""" """
self.tth_min = lower self.tth_min = lower
self.tth_max = upper self.tth_max = upper
cdef numpy.ndarray[numpy.float32_t, ndim = 1] tth_out = numpy.empty(self._nBins,dtype=numpy.float32) # cdef numpy.ndarray[numpy.float32_t, ndim = 1] tth_out = numpy.empty(self._nBins,dtype=numpy.float32)
if self.ctth_out:
free(self.ctth_out)
self.ctth_out= <float*> malloc(self._nBins*sizeof(float))
cdef float delta = (upper - lower ) / (< float > (self._nBins)) cdef float delta = (upper - lower ) / (< float > (self._nBins))
with nogil: with nogil:
for i in range(self._nBins): for i in range(self._nBins):
tth_out[i] = self.tth_min + (0.5 +< float > i) * delta self.ctth_out[i] = <float>( self.tth_min + (0.5 +< float > i) * delta)
self.tth_out = tth_out # self.tth_out = tth_out
def getConfiguration(self, int Nimage, int Nbins, useFp64=None): def getConfiguration(self, int Nimage, int Nbins, useFp64=None):
@ -244,17 +249,19 @@ cdef class Integrator1d:
set / unset and loadTth methods have a direct impact on the execute() method. set / unset and loadTth methods have a direct impact on the execute() method.
All the rest of the methods will require at least a new configuration via configure()""" All the rest of the methods will require at least a new configuration via configure()"""
cdef int rc,i cdef int rc,i
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cimage, histogram, bins cdef numpy.ndarray[numpy.float32_t, ndim = 1] cimage, histogram, bins,tth_out
cimage = numpy.ascontiguousarray(image.ravel(),dtype=numpy.float32) cimage = numpy.ascontiguousarray(image.ravel(),dtype=numpy.float32)
outPos = numpy.zeros(self._nBins,dtype=numpy.float32) histogram = numpy.empty(self._nBins,dtype=numpy.float32)
histogram = numpy.zeros(self._nBins,dtype=numpy.float32) bins = numpy.empty(self._nBins,dtype=numpy.float32)
bins = numpy.zeros(self._nBins,dtype=numpy.float32) tth_out = numpy.empty(self._nBins,dtype=numpy.float32)
assert cimage.size == self._nData assert cimage.size == self._nData
with nogil: with nogil:
rc = self.cpp_integrator.execute(<float*> cimage.data, <float*> histogram.data, <float*> bins.data) rc = self.cpp_integrator.execute(<float*> cimage.data, <float*> histogram.data, <float*> bins.data)
if rc!=0: if rc!=0:
raise RuntimeError("OpenCL integrator failed with RC=%s"%rc) raise RuntimeError("OpenCL integrator failed with RC=%s"%rc)
return self.tth_out,histogram,bins
memcpy(tth_out.data,self.ctth_out,self._nBins*sizeof(float))
return tth_out,histogram,bins
def clean(self, int preserve_context=0): def clean(self, int preserve_context=0):
"""Free OpenCL related resources. """Free OpenCL related resources.

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf8 -*- # -*- coding: utf8 -*-
# #
# Project: Azimuthal integration # Project: Azimuthal integration
# https://forge.epn-campus.eu/projects/azimuthal # https://forge.epn-campus.eu/projects/azimuthal
# #
# File: "$Id$" # File: "$Id$"
@ -43,11 +43,11 @@ logger = logging.getLogger("pyFAI.azimuthalIntegrator")
class AzimuthalIntegrator(Geometry): class AzimuthalIntegrator(Geometry):
""" """
This class is an azimuthal integrator based on P. Boesecke's geometry and This class is an azimuthal integrator based on P. Boesecke's geometry and
histogram algorithm by Manolo S. del Rio and V.A Sole histogram algorithm by Manolo S. del Rio and V.A Sole
All geometry calculation are done in the Geometry class All geometry calculation are done in the Geometry class
""" """
def __init__(self, dist=1, poni1=0, poni2=0, rot1=0, rot2=0, rot3=0, pixel1=1, pixel2=1, splineFile=None): def __init__(self, dist=1, poni1=0, poni2=0, rot1=0, rot2=0, rot3=0, pixel1=1, pixel2=1, splineFile=None):
""" """
@ -59,7 +59,7 @@ class AzimuthalIntegrator(Geometry):
@param rot3: third rotation from sample ref to detector's ref, in radians @param rot3: third rotation from sample ref to detector's ref, in radians
@param pixel1: pixel size of the fist dimension of the detector, in meter @param pixel1: pixel size of the fist dimension of the detector, in meter
@param pixel2: pixel size of the second dimension of the detector, in meter @param pixel2: pixel size of the second dimension of the detector, in meter
@param splineFile: file containing the geometric distortion of the detector. Overrides the pixel size. @param splineFile: file containing the geometric distortion of the detector. Overrides the pixel size.
""" """
Geometry.__init__(self, dist, poni1, poni2, rot1, rot2, rot3, pixel1, pixel2, splineFile) Geometry.__init__(self, dist, poni1, poni2, rot1, rot2, rot3, pixel1, pixel2, splineFile)
self._nbPixCache = {} #key=shape, value: array self._nbPixCache = {} #key=shape, value: array
@ -80,10 +80,10 @@ class AzimuthalIntegrator(Geometry):
def makeMask(self, data, mask=None, dummy=None, delta_dummy=None, invertMask=None): def makeMask(self, data, mask=None, dummy=None, delta_dummy=None, invertMask=None):
""" """
Combines a mask Combines a mask
For the mask: 1 for good pixels, 0 for bas pixels For the mask: 1 for good pixels, 0 for bas pixels
@param data: input array of @param data: input array of
@param mask: input mask @param mask: input mask
@param dummy: value of dead pixels @param dummy: value of dead pixels
@param delta_dumy: precision of dummy pixels @param delta_dumy: precision of dummy pixels
@param invertMask: to force inversion of the input mask @param invertMask: to force inversion of the input mask
@ -133,10 +133,10 @@ class AzimuthalIntegrator(Geometry):
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: (2theta, I) in degrees @return: (2theta, I) in degrees
@rtype: 2-tuple of 1D arrays @rtype: 2-tuple of 1D arrays
@ -182,10 +182,10 @@ class AzimuthalIntegrator(Geometry):
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: (2theta, I) in degrees @return: (2theta, I) in degrees
@rtype: 2-tuple of 1D arrays @rtype: 2-tuple of 1D arrays
@ -231,13 +231,13 @@ class AzimuthalIntegrator(Geometry):
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()). @param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type chiRange: (float, float), optional @type chiRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: (2theta, I) in degrees @return: (2theta, I) in degrees
@rtype: 2-tuple of 1D arrays @rtype: 2-tuple of 1D arrays
@ -287,10 +287,10 @@ class AzimuthalIntegrator(Geometry):
def xrpd_splitPixel(self, data, nbPt, filename=None, correctSolidAngle=True, def xrpd_splitPixel(self, data, nbPt, filename=None, correctSolidAngle=True,
tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None): tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None):
""" """
Calculate the powder diffraction pattern from a set of data, an image. Calculate the powder diffraction pattern from a set of data, an image.
Cython implementation Cython implementation
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt: number of points in the output pattern @param nbPt: number of points in the output pattern
@ -298,7 +298,7 @@ class AzimuthalIntegrator(Geometry):
@param filename: file to save data in @param filename: file to save data in
@type filename: string @type filename: string
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: (2theta, I) in degrees @return: (2theta, I) in degrees
@rtype: 2-tuple of 1D arrays @rtype: 2-tuple of 1D arrays
@ -352,18 +352,18 @@ class AzimuthalIntegrator(Geometry):
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()). @param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type chiRange: (float, float), optional, disabled for now @type chiRange: (float, float), optional, disabled for now
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
OpenCL specific parameters: OpenCL specific parameters:
@param devicetype: "cpu" or "gpu" or "all" or "def" @param devicetype: "cpu" or "gpu" or "all" or "def"
@param useFp64: shall histogram be done in double precision (adviced) @param useFp64: shall histogram be done in double precision (adviced)
@param platformid: platform number @param platformid: platform number
@param deviceid: device number @param deviceid: device number
@param safe: set to false if you think your GPU is already set-up correctly (2theta, mask, solid angle...) @param safe: set to false if you think your GPU is already set-up correctly (2theta, mask, solid angle...)
@return 2theta, I, weighted histogram, unweighted histogram @return 2theta, I, weighted histogram, unweighted histogram
@ -389,7 +389,7 @@ class AzimuthalIntegrator(Geometry):
with self._ocl_sem: with self._ocl_sem:
if self._ocl is None: if self._ocl is None:
size = data.size size = data.size
fd, tmpfile = tempfile.mkstemp(".log", "pyfai-opencl") fd, tmpfile = tempfile.mkstemp(".log", "pyfai-opencl-")
os.close(fd) os.close(fd)
integr = ocl_azim.Integrator1d(tmpfile) integr = ocl_azim.Integrator1d(tmpfile)
if platformid and deviceid: if platformid and deviceid:
@ -449,19 +449,19 @@ class AzimuthalIntegrator(Geometry):
tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None): tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None):
""" """
Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image
Pure numpy implementation (VERY SLOW !!!) Pure numpy implementation (VERY SLOW !!!)
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta) @param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta)
@type nbPt: integer @type nbPt: integer
@param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi) @param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi)
@type nbPtChi: integer @type nbPtChi: integer
@param filename: file to save data in @param filename: file to save data in
@type filename: string @type filename: string
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: azimuthaly regrouped data, 2theta pos and chipos @return: azimuthaly regrouped data, 2theta pos and chipos
@rtype: 3-tuple of ndarrays @rtype: 3-tuple of ndarrays
@ -497,19 +497,19 @@ class AzimuthalIntegrator(Geometry):
tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None): tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None):
""" """
Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image
Cython implementation: fast but incaccurate Cython implementation: fast but incaccurate
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta) @param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta)
@type nbPt: integer @type nbPt: integer
@param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi) @param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi)
@type nbPtChi: integer @type nbPtChi: integer
@param filename: file to save data in @param filename: file to save data in
@type filename: string @type filename: string
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: azimuthaly regrouped data, 2theta pos and chipos @return: azimuthaly regrouped data, 2theta pos and chipos
@rtype: 3-tuple of ndarrays @rtype: 3-tuple of ndarrays
@ -550,27 +550,27 @@ class AzimuthalIntegrator(Geometry):
tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None): tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None):
""" """
Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image
Split pixels according to their coordinate and a bounding box Split pixels according to their coordinate and a bounding box
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta) @param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta)
@type nbPt: integer @type nbPt: integer
@param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi) @param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi)
@type nbPtChi: integer @type nbPtChi: integer
@param filename: file to save data in @param filename: file to save data in
@type filename: string @type filename: string
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param chiRange: The lower and upper range of the azimuthal angle. If not provided, range is simply (data.min(), data.max()). @param chiRange: The lower and upper range of the azimuthal angle. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type chiRange: (float, float), optional @type chiRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: azimuthaly regrouped data, 2theta pos. and chi pos. @return: azimuthaly regrouped data, 2theta pos. and chi pos.
@rtype: 3-tuple of ndarrays @rtype: 3-tuple of ndarrays
@ -616,27 +616,27 @@ class AzimuthalIntegrator(Geometry):
tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None): tthRange=None, chiRange=None, mask=None, dummy=None, delta_dummy=None):
""" """
Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image Calculate the 2D powder diffraction pattern (2Theta,Chi) from a set of data, an image
Split pixels according to their corner positions Split pixels according to their corner positions
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta) @param nbPt2Th: number of points in the output pattern in the Radial (horizontal) axis (2 theta)
@type nbPt: integer @type nbPt: integer
@param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi) @param nbPtChi: number of points in the output pattern along the Azimuthal (vertical) axis (chi)
@type nbPtChi: integer @type nbPtChi: integer
@param filename: file to save data in @param filename: file to save data in
@type filename: string @type filename: string
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()). @param tthRange: The lower and upper range of the 2theta. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type tthRange: (float, float), optional @type tthRange: (float, float), optional
@param chiRange: The lower and upper range of the azimuthal angle. If not provided, range is simply (data.min(), data.max()). @param chiRange: The lower and upper range of the azimuthal angle. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type chiRange: (float, float), optional @type chiRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@return: azimuthaly regrouped data, 2theta pos. and chi pos. @return: azimuthaly regrouped data, 2theta pos. and chi pos.
@rtype: 3-tuple of ndarrays @rtype: 3-tuple of ndarrays
@ -707,10 +707,10 @@ class AzimuthalIntegrator(Geometry):
qRange=None, chiRange=None, mask=None, dummy=None, qRange=None, chiRange=None, mask=None, dummy=None,
delta_dummy=None, method="bbox"): delta_dummy=None, method="bbox"):
""" """
Calculate the azimuthal integrated Saxs curve Calculate the azimuthal integrated Saxs curve
Multi algorithm implementation (tries to be bullet proof) Multi algorithm implementation (tries to be bullet proof)
@param data: 2D array from the CCD camera @param data: 2D array from the CCD camera
@type data: ndarray @type data: ndarray
@param nbPt: number of points in the output pattern @param nbPt: number of points in the output pattern
@ -718,19 +718,19 @@ class AzimuthalIntegrator(Geometry):
@param filename: file to save data to @param filename: file to save data to
@type filename: string @type filename: string
@param correctSolidAngle: if True, the data are devided by the solid angle of each pixel @param correctSolidAngle: if True, the data are devided by the solid angle of each pixel
@type correctSolidAngle: boolean @type correctSolidAngle: boolean
@param variance: array containing the variance of the data @param variance: array containing the variance of the data
@type variance: ndarray @type variance: ndarray
@param qRange: The lower and upper range of the sctter vector q. If not provided, range is simply (data.min(), data.max()). @param qRange: The lower and upper range of the sctter vector q. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type qRange: (float, float), optional @type qRange: (float, float), optional
@param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()). @param chiRange: The lower and upper range of the chi angle. If not provided, range is simply (data.min(), data.max()).
Values outside the range are ignored. Values outside the range are ignored.
@type chiRange: (float, float), optional @type chiRange: (float, float), optional
@param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels @param mask: array (same siza as image) with 0 for masked pixels, and 1 for valid pixels
@param dummy: value for dead/masked pixels @param dummy: value for dead/masked pixels
@param delta_dummy: precision for dummy value @param delta_dummy: precision for dummy value
@param method: can be "numpy", "cython", "BBox" or "splitpixel" @param method: can be "numpy", "cython", "BBox" or "splitpixel"
@return: azimuthaly regrouped data, 2theta pos. and chi pos. @return: azimuthaly regrouped data, 2theta pos. and chi pos.
@rtype: 3-tuple of ndarrays @rtype: 3-tuple of ndarrays
""" """

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf8 -*- # -*- coding: utf8 -*-
# #
# Project: Azimuthal integration # Project: Azimuthal integration
# https://forge.epn-campus.eu/projects/azimuthal # https://forge.epn-campus.eu/projects/azimuthal
# #
# File: "$Id$" # File: "$Id$"
@ -36,7 +36,7 @@ cdef extern from "math.h":
@cython.cdivision(True) @cython.cdivision(True)
cdef float getBinNr(float x0, float pos0_min, float delta) nogil: cdef float getBinNr(float x0, float pos0_min, float delta) nogil:
""" """
calculate the bin number for any point calculate the bin number for any point
param x0: current position param x0: current position
param pos0_min: position minimum param pos0_min: position minimum
param delta: bin width param delta: bin width
@ -61,48 +61,47 @@ def histoBBox1d(numpy.ndarray weights not None,
): ):
""" """
Calculates histogram of pos0 (tth) weighted by weights Calculates histogram of pos0 (tth) weighted by weights
Splitting is done on the pixel's bounding box like fit2D Splitting is done on the pixel's bounding box like fit2D
@param weights: array with intensities @param weights: array with intensities
@param pos0: 1D array with pos0: tth or q_vect @param pos0: 1D array with pos0: tth or q_vect
@param delta_pos0: 1D array with delta pos0: max center-corner distance @param delta_pos0: 1D array with delta pos0: max center-corner distance
@param pos1: 1D array with pos1: chi @param pos1: 1D array with pos1: chi
@param delta_pos1: 1D array with max pos1: max center-corner distance, unused ! @param delta_pos1: 1D array with max pos1: max center-corner distance, unused !
@param bins: number of output bins @param bins: number of output bins
@param pos0Range: minimum and maximum of the 2th range @param pos0Range: minimum and maximum of the 2th range
@param pos1Range: minimum and maximum of the chi range @param pos1Range: minimum and maximum of the chi range
@param dummy: value for bins without pixels @param dummy: value for bins without pixels
@return 2theta, I, weighted histogram, unweighted histogram @return 2theta, I, weighted histogram, unweighted histogram
""" """
cdef long size = weights.size cdef long size = weights.size
assert pos0.size == size assert pos0.size == size
assert delta_pos0.size == size assert delta_pos0.size == size
assert bins > 1 assert bins > 1
print pos0Range
cdef long bin0_max, bin0_min, bin = 0 cdef long bin0_max, bin0_min, bin = 0
cdef float data, deltaR, deltaL, deltaA,p1, epsilon = 1e-10, cdummy, ddummy cdef float data, deltaR, deltaL, deltaA,p1, epsilon = 1e-10, cdummy, ddummy
cdef float pos0_min, pos0_max, pos0_maxin, pos1_min, pos1_max, pos1_maxin, min0, max0, fbin0_min, fbin0_max cdef float pos0_min, pos0_max, pos0_maxin, pos1_min, pos1_max, pos1_maxin, min0, max0, fbin0_min, fbin0_max
cdef int checkpos1 = 0, check_mask = 0, check_dummy = 0 cdef long checkpos1 = 0, check_mask = 0, check_dummy = 0
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cdata = numpy.ascontiguousarray(weights.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] cdata = numpy.ascontiguousarray(weights.ravel(),dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0, dpos0, cpos1, dpos1,cpos0_lower, cpos0_upper cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0, dpos0, cpos1, dpos1,cpos0_lower, cpos0_upper
cpos0 = numpy.ascontiguousarray(pos0.ravel(), dtype="float32") cpos0 = numpy.ascontiguousarray(pos0.ravel(), dtype=numpy.float32)
dpos0 = numpy.ascontiguousarray(delta_pos0.ravel(), dtype="float32") dpos0 = numpy.ascontiguousarray(delta_pos0.ravel(), dtype=numpy.float32)
cdef numpy.ndarray[numpy.float64_t, ndim = 1] outData = numpy.zeros(bins, dtype="float64") cdef numpy.ndarray[numpy.float64_t, ndim = 1] outData = numpy.zeros(bins, dtype=numpy.float64)
cdef numpy.ndarray[numpy.float64_t, ndim = 1] outCount = numpy.zeros(bins, dtype="float64") cdef numpy.ndarray[numpy.float64_t, ndim = 1] outCount = numpy.zeros(bins, dtype=numpy.float64)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] outMerge = numpy.zeros(bins, dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] outMerge = numpy.zeros(bins, dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] outPos = numpy.zeros(bins, dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] outPos = numpy.zeros(bins, dtype=numpy.float32)
cdef numpy.ndarray[numpy.int8_t, ndim = 1] cmask cdef numpy.ndarray[numpy.int8_t, ndim = 1] cmask
if mask is not None: if mask is not None:
assert mask.size == size assert mask.size == size
check_mask = 1 check_mask = 1
cmask = numpy.ascontiguousarray(mask.ravel(),dtype="int8") cmask = numpy.ascontiguousarray(mask.ravel(),dtype="int8")
if (dummy is not None) and delta_dummy is not None: if (dummy is not None) and delta_dummy is not None:
check_dummy = 1 check_dummy = 1
cdummy = float(dummy) cdummy = float(dummy)
ddummy = float(delta_dummy) ddummy = float(delta_dummy)
@ -111,8 +110,8 @@ def histoBBox1d(numpy.ndarray weights not None,
else: else:
cdummy=0.0 cdummy=0.0
cpos0_lower = numpy.zeros(size, dtype="float32") cpos0_lower = numpy.zeros(size, dtype=numpy.float32)
cpos0_upper = numpy.zeros(size, dtype="float32") cpos0_upper = numpy.zeros(size, dtype=numpy.float32)
pos0_min=cpos0[0] pos0_min=cpos0[0]
pos0_max=cpos0[0] pos0_max=cpos0[0]
with nogil: with nogil:
@ -125,7 +124,7 @@ def histoBBox1d(numpy.ndarray weights not None,
pos0_max=max0 pos0_max=max0
if min0<pos0_min: if min0<pos0_min:
pos0_min=min0 pos0_min=min0
if pos0Range is not None and len(pos0Range) > 1: if pos0Range is not None and len(pos0Range) > 1:
pos0_min = min(pos0Range) pos0_min = min(pos0Range)
pos0_maxin = max(pos0Range) pos0_maxin = max(pos0Range)
@ -138,8 +137,8 @@ def histoBBox1d(numpy.ndarray weights not None,
assert pos1.size == size assert pos1.size == size
assert delta_pos1.size == size assert delta_pos1.size == size
checkpos1 = 1 checkpos1 = 1
cpos1 = numpy.ascontiguousarray(pos1.ravel(),dtype="float32") cpos1 = numpy.ascontiguousarray(pos1.ravel(),dtype=numpy.float32)
dpos1 = numpy.ascontiguousarray(delta_pos1.ravel(),dtype="float32") dpos1 = numpy.ascontiguousarray(delta_pos1.ravel(),dtype=numpy.float32)
pos1_min = min(pos1Range) pos1_min = min(pos1Range)
pos1_maxin = max(pos1Range) pos1_maxin = max(pos1Range)
pos1_max = pos1_maxin * (1 + numpy.finfo(numpy.float32).eps) pos1_max = pos1_maxin * (1 + numpy.finfo(numpy.float32).eps)
@ -153,24 +152,24 @@ def histoBBox1d(numpy.ndarray weights not None,
for idx in range(size): for idx in range(size):
if (check_mask) and cmask[idx]: if (check_mask) and cmask[idx]:
continue continue
data = cdata[idx] data = cdata[idx]
if check_dummy and fabs(data-cdummy)<ddummy: if check_dummy and fabs(data-cdummy)<ddummy:
continue continue
min0 = cpos0_lower[idx] min0 = cpos0_lower[idx]
max0 = cpos0_upper[idx] max0 = cpos0_upper[idx]
if checkpos1: if checkpos1:
if ((cpos1[idx]+dpos1[idx]) < pos1_min) or ((cpos1[idx]-dpos1[idx]) > pos1_max): if ((cpos1[idx]+dpos1[idx]) < pos1_min) or ((cpos1[idx]-dpos1[idx]) > pos1_max):
continue continue
fbin0_min = getBinNr(min0, pos0_min, delta) fbin0_min = getBinNr(min0, pos0_min, delta)
fbin0_max = getBinNr(max0, pos0_min, delta) fbin0_max = getBinNr(max0, pos0_min, delta)
bin0_min = < long > floor(fbin0_min) bin0_min = < long > floor(fbin0_min)
bin0_max = < long > floor(fbin0_max) bin0_max = < long > floor(fbin0_max)
if (bin0_max<0) or (bin0_min>=bins): if (bin0_max<0) or (bin0_min>=bins):
continue continue
if bin0_max>=bins: if bin0_max>=bins:
bin0_max=bins-1 bin0_max=bins-1
@ -224,25 +223,25 @@ def histoBBox2d(numpy.ndarray weights not None,
float dummy=0.0): float dummy=0.0):
""" """
Calculate 2D histogram of pos0(tth),pos1(chi) weighted by weights Calculate 2D histogram of pos0(tth),pos1(chi) weighted by weights
Splitting is done on the pixel's bounding box like fit2D Splitting is done on the pixel's bounding box like fit2D
@param weights: array with intensities @param weights: array with intensities
@param pos0: 1D array with pos0: tth or q_vect @param pos0: 1D array with pos0: tth or q_vect
@param delta_pos0: 1D array with delta pos0: max center-corner distance @param delta_pos0: 1D array with delta pos0: max center-corner distance
@param pos1: 1D array with pos1: chi @param pos1: 1D array with pos1: chi
@param delta_pos1: 1D array with max pos1: max center-corner distance, unused ! @param delta_pos1: 1D array with max pos1: max center-corner distance, unused !
@param bins: number of output bins (tth=100, chi=36 by default) @param bins: number of output bins (tth=100, chi=36 by default)
@param pos0Range: minimum and maximum of the 2th range @param pos0Range: minimum and maximum of the 2th range
@param pos1Range: minimum and maximum of the chi range @param pos1Range: minimum and maximum of the chi range
@param dummy: value for bins without pixels @param dummy: value for bins without pixels
@return 2theta, I, weighted histogram, unweighted histogram @return 2theta, I, weighted histogram, unweighted histogram
@return I, edges0, edges1, weighted histogram(2D), unweighted histogram (2D) @return I, edges0, edges1, weighted histogram(2D), unweighted histogram (2D)
""" """
cdef long bins0, bins1, i, j, idx cdef long bins0, bins1, i, j, idx
cdef long size = weights.size cdef long size = weights.size
assert pos0.size == size assert pos0.size == size
assert pos1.size == size assert pos1.size == size
assert delta_pos0.size == size assert delta_pos0.size == size
@ -255,32 +254,54 @@ def histoBBox2d(numpy.ndarray weights not None,
bins0 = 1 bins0 = 1
if bins1 <= 0: if bins1 <= 0:
bins1 = 1 bins1 = 1
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cdata = numpy.ascontiguousarray(weights.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] cdata = numpy.ascontiguousarray(weights.ravel(),dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0 = numpy.ascontiguousarray(pos0.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0 = numpy.ascontiguousarray(pos0.ravel(),dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] dpos0 = numpy.ascontiguousarray(delta_pos0.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] dpos0 = numpy.ascontiguousarray(delta_pos0.ravel(),dtype=numpy.float32)
# cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0_inf = cpos0 - dpos0 cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1 = numpy.ascontiguousarray(pos1.ravel(),dtype=numpy.float32)
# cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0_sup = cpos0 + dpos0 cdef numpy.ndarray[numpy.float32_t, ndim = 1] dpos1 = numpy.ascontiguousarray(delta_pos1.ravel(),dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1 = numpy.ascontiguousarray(pos1.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0_upper = numpy.zeros(size, dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] dpos1 = numpy.ascontiguousarray(delta_pos1.ravel(),dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos0_lower = numpy.zeros(size, dtype=numpy.float32)
# cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1_inf = cpos1 - dpos1 # cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1_upper = numpy.zeros(size, dtype=numpy.float32)
# cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1_sup = cpos1 + dpos1 # cdef numpy.ndarray[numpy.float32_t, ndim = 1] cpos1_lower = numpy.zeros(size, dtype=numpy.float32)
cdef numpy.ndarray[numpy.float64_t, ndim = 2] outData = numpy.zeros((bins0, bins1), dtype="float64") cdef numpy.ndarray[numpy.float64_t, ndim = 2] outData = numpy.zeros((bins0, bins1), dtype=numpy.float64)
cdef numpy.ndarray[numpy.float64_t, ndim = 2] outCount = numpy.zeros((bins0, bins1), dtype="float64") cdef numpy.ndarray[numpy.float64_t, ndim = 2] outCount = numpy.zeros((bins0, bins1), dtype=numpy.float64)
cdef numpy.ndarray[numpy.float32_t, ndim = 2] outMerge = numpy.zeros((bins0, bins1), dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 2] outMerge = numpy.zeros((bins0, bins1), dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] edges0 = numpy.zeros(bins0, dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] edges0 = numpy.zeros(bins0, dtype=numpy.float32)
cdef numpy.ndarray[numpy.float32_t, ndim = 1] edges1 = numpy.zeros(bins1, dtype="float32") cdef numpy.ndarray[numpy.float32_t, ndim = 1] edges1 = numpy.zeros(bins1, dtype=numpy.float32)
cdef float min0, max0, min1, max1, deltaR, deltaL, deltaU, deltaD, deltaA, tmp, delta0, delta1 cdef float min0, max0, min1, max1, deltaR, deltaL, deltaU, deltaD, deltaA, tmp, delta0, delta1
cdef float pos0_min, pos0_max, pos1_min, pos1_max, pos0_maxin, pos1_maxin cdef float pos0_min, pos0_max, pos1_min, pos1_max, pos0_maxin, pos1_maxin
cdef float fbin0_min, fbin0_max, fbin1_min, fbin1_max, data, epsilon = 1e-10 cdef float fbin0_min, fbin0_max, fbin1_min, fbin1_max, data, epsilon = 1e-10
cdef long bin0_max, bin0_min, bin1_max, bin1_min cdef long bin0_max, bin0_min, bin1_max, bin1_min
pos0_min=cpos0[0]
pos0_max=cpos0[0]
with nogil:
for idx in range(size):
min0 = cpos0[idx] - dpos0[idx]
max0 = cpos0[idx] + dpos0[idx]
cpos0_upper[idx]=max0
cpos0_lower[idx]=min0
if max0>pos0_max:
pos0_max=max0
if min0<pos0_min:
pos0_min=min0
# for idx in range(size):
# min1 = cpos1[idx] - dpos1[idx]
# max1 = cpos1[idx] + dpos1[idx]
# cpos1_upper[idx]=max1
# cpos1_lower[idx]=min1
# if max1>pos1_max:
# pos1_max=max1
# if min1<pos1_min:
# pos1_min=min1
if (pos0Range is not None) and (len(pos0Range) == 2): if (pos0Range is not None) and (len(pos0Range) == 2):
pos0_min = min(pos0Range) pos0_min = min(pos0Range)
pos0_maxin = max(pos0Range) pos0_maxin = max(pos0Range)
else: else:
pos0_min = cpos0.min() pos0_min = pos0_min
pos0_maxin = cpos0.max() pos0_maxin = pos0_max
if pos0_min<0: if pos0_min<0:
pos0_min=0 pos0_min=0
pos0_max = pos0_maxin * (1 + numpy.finfo(numpy.float32).eps) pos0_max = pos0_maxin * (1 + numpy.finfo(numpy.float32).eps)
@ -289,7 +310,6 @@ def histoBBox2d(numpy.ndarray weights not None,
pos1_min = min(pos1Range) pos1_min = min(pos1Range)
pos1_maxin = max(pos1Range) pos1_maxin = max(pos1Range)
else: else:
# tmp = cdelta_pos1.min()
pos1_min = cpos1.min() pos1_min = cpos1.min()
pos1_maxin = cpos1.max() pos1_maxin = cpos1.max()
pos1_max = pos1_maxin * (1 + numpy.finfo(numpy.float32).eps) pos1_max = pos1_maxin * (1 + numpy.finfo(numpy.float32).eps)
@ -305,8 +325,8 @@ def histoBBox2d(numpy.ndarray weights not None,
for idx in range(size): for idx in range(size):
data = cdata[idx] data = cdata[idx]
min0 = cpos0[idx] - dpos0[idx] min0 = cpos0_lower[idx]
max0 = cpos0[idx] + dpos0[idx] max0 = cpos0_upper[idx]
min1 = cpos1[idx] - dpos1[idx] min1 = cpos1[idx] - dpos1[idx]
max1 = cpos1[idx] + dpos1[idx] max1 = cpos1[idx] + dpos1[idx]
@ -404,7 +424,7 @@ def histoBBox2d(numpy.ndarray weights not None,
for i in range(bins0): for i in range(bins0):
for j in range(bins1): for j in range(bins1):
if outCount[i, j] > epsilon: if outCount[i, j] > epsilon:
outMerge[i, j] = outData[i, j] / outCount[i, j] outMerge[i, j] = <float> (outData[i, j] / outCount[i, j])
else: else:
outMerge[i, j] = dummy outMerge[i, j] = dummy
return outMerge.T, edges0, edges1, outData.T, outCount.T return outMerge.T, edges0, edges1, outData.T, outCount.T