mirror of https://github.com/silx-kit/pyFAI.git
work on refactoring distortion correction
This commit is contained in:
Jerome Kieffer
parent
c6b4a4a21d
commit
5fe5ff7321
|
|
@ -32,6 +32,9 @@ Files: version.py
|
|||
pyFAI/diffmap.py
|
||||
pyFAI/ocl_sort.py
|
||||
pyFAI/decorators.py
|
||||
pyFAI/distortion.py
|
||||
pyFAI/ext/_distortion.py
|
||||
pyFAI/test/test_distortion.py
|
||||
run_tests.py
|
||||
Copyright: 2015-2016 European Synchrotron Radiation Facility
|
||||
License: MIT/X11 (BSD like)
|
||||
|
|
|
|||
|
|
@ -4,30 +4,33 @@
|
|||
# Project: Azimuthal integration
|
||||
# https://github.com/pyFAI/pyFAI
|
||||
#
|
||||
# Copyright (C) European Synchrotron Radiation Facility, Grenoble, France
|
||||
#
|
||||
# Principal author: Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
|
||||
#
|
||||
# This program is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# This program is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# GNU General Public License for more details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
# Copyright 2013-2016 (C) European Synchrotron Radiation Facility, Grenoble, France
|
||||
#
|
||||
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
# of this software and associated documentation files (the "Software"), to deal
|
||||
# in the Software without restriction, including without limitation the rights
|
||||
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
# copies of the Software, and to permit persons to whom the Software is
|
||||
# furnished to do so, subject to the following conditions:
|
||||
# .
|
||||
# The above copyright notice and this permission notice shall be included in
|
||||
# all copies or substantial portions of the Software.
|
||||
# .
|
||||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
# THE SOFTWARE.
|
||||
|
||||
from __future__ import absolute_import, print_function, division, with_statement
|
||||
|
||||
__author__ = "Jérôme Kieffer"
|
||||
__contact__ = "Jerome.Kieffer@ESRF.eu"
|
||||
__license__ = "GPLv3+"
|
||||
__license__ = "MIT"
|
||||
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
|
||||
__date__ = "02/05/2016"
|
||||
__date__ = "03/05/2016"
|
||||
__status__ = "development"
|
||||
|
||||
import logging
|
||||
|
|
@ -61,38 +64,45 @@ class Distortion(object):
|
|||
|
||||
New version compatible both with CSR and LUT...
|
||||
"""
|
||||
def __init__(self, detector="detector", shape=None, method="LUT", device=None, workgroup=8):
|
||||
def __init__(self, detector="detector", shape=None, resize=False, empty=0,
|
||||
mask=None, method="LUT", device=None, workgroup=8):
|
||||
"""
|
||||
@param detector: detector instance or detector name
|
||||
@param shape: shape of the output image
|
||||
@param resize: allow the output shape to be different from the input shape
|
||||
@param empty: value to be given for empty bins
|
||||
@param method: "lut" or "csr", the former is faster
|
||||
@param device: Name of the device: None for OpenMP, "cpu" or "gpu" or the id of the OpenCL device a 2-tuple of integer
|
||||
@param workgroup: workgroup size for CSR on OpenCL
|
||||
"""
|
||||
self._shape_out = None
|
||||
if isinstance(detector, six.string_types):
|
||||
self.detector = detectors.detector_factory(detector)
|
||||
else: # we assume it is a Detector instance
|
||||
self.detector = detector
|
||||
if shape is not None:
|
||||
self.shape = tuple([int(i) for i in shape])
|
||||
self.shape_in = self.detector.shape
|
||||
if mask is not None:
|
||||
self.mask = numpy.ascontiguousarray(mask, numpy.int8)
|
||||
else:
|
||||
inshape = self.detector.shape
|
||||
corner_pos = self.self.detector.get_pixel_corners()
|
||||
corner_pos.shape = (-1, 3)
|
||||
pos0_min, pos1_min, pos2_min = corner_pos.min(axis=0)
|
||||
pos0_max, pos1_max, pos2_max = corner_pos.max(axis=0)
|
||||
# z is coord 0
|
||||
self.mask = numpy.ascontiguousarray(self.detector.mask, numpy.int8)
|
||||
self.resize = resize
|
||||
if shape is not None:
|
||||
self._shape_out = tuple([int(i) for i in shape])
|
||||
elif not self.resize:
|
||||
if self.detector.shape is not None:
|
||||
self._shape_out = self.detector.shape
|
||||
else:
|
||||
raise RuntimeError("You need to provide either the detector or its shape")
|
||||
|
||||
self.out_shape = (int(ceil((pos1_max - pos1_min) / self.detector.pixel1)),
|
||||
int(ceil((pos2_max - pos2_min) / self.detector.pixel2)))
|
||||
self._sem = threading.Semaphore()
|
||||
|
||||
self.bin_size = None
|
||||
self.max_size = None
|
||||
self.pos = None
|
||||
self.lut = None
|
||||
self.delta0 = self.delta1 = None # max size of an pixel on a regular grid ...
|
||||
self.delta1 = self.delta2 = None # max size of an pixel on a regular grid ...
|
||||
self.offset1 = self.offset2 = 0 # position of the first bin
|
||||
self.integrator = None
|
||||
self.empty = empty # "dummy" value for empty bins
|
||||
if not method:
|
||||
self.method = "lut"
|
||||
else:
|
||||
|
|
@ -107,19 +117,21 @@ class Distortion(object):
|
|||
return os.linesep.join(["Distortion correction %s on device %s for detector shape %s:" % (self.method, self.device, self.shape),
|
||||
self.detector.__repr__()])
|
||||
|
||||
def reset(self, method=None, device=None, workgroup=None):
|
||||
def reset(self, method=None, device=None, workgroup=None, prepare=True):
|
||||
"""
|
||||
reset the distortion correction and re-calculate the look-up table
|
||||
|
||||
@param method: can be "lut" or "csr", "lut" looks faster
|
||||
@param device: can be None, "cpu" or "gpu" or the id as a 2-tuple of integer
|
||||
@param worgroup: enforce the workgroup size for CSR.
|
||||
@param prepare: set to false to only reset and not re-initialize
|
||||
"""
|
||||
with self._sem:
|
||||
self.max_size = None
|
||||
self.pos = None
|
||||
self.lut = None
|
||||
self.delta0 = self.delta1 = None
|
||||
self.delta1 = self.delta2 = None
|
||||
self.offset1 = self.offset2 = 0
|
||||
self.integrator = None
|
||||
if method is not None:
|
||||
self.method = method.lower()
|
||||
|
|
@ -127,38 +139,73 @@ class Distortion(object):
|
|||
self.device = device
|
||||
if workgroup is not None:
|
||||
self.workgroup = int(workgroup)
|
||||
if prepare:
|
||||
self.calc_init()
|
||||
|
||||
self.calc_init()
|
||||
@property
|
||||
def shape_out(self):
|
||||
"""
|
||||
Calculate/cache the output shape
|
||||
|
||||
@return output shape
|
||||
"""
|
||||
if self._shape_out is None:
|
||||
self.calc_pos()
|
||||
return self._shape_out
|
||||
|
||||
@timeit
|
||||
def calc_pos(self):
|
||||
"""Calculate the pixel position on the regular grid
|
||||
|
||||
@return: pixel corner positions (in pixel units) on the regular grid
|
||||
@rtyep: ndarray of shape (nrow, ncol, 4, 2)
|
||||
"""
|
||||
if self.delta1 is None:
|
||||
with self._sem:
|
||||
if self.delta1 is None:
|
||||
pos_corners = numpy.empty((self.shape[0] + 1, self.shape[1] + 1, 2), dtype=numpy.float64)
|
||||
d1 = numpy.outer(numpy.arange(self.shape[0] + 1, dtype=numpy.float64), numpy.ones(self.shape[1] + 1, dtype=numpy.float64)) - 0.5
|
||||
d2 = numpy.outer(numpy.ones(self.shape[0] + 1, dtype=numpy.float64), numpy.arange(self.shape[1] + 1, dtype=numpy.float64)) - 0.5
|
||||
pos_corners[:, :, 0], pos_corners[:, :, 1] = self.detector.calc_cartesian_positions(d1, d2)[:2]
|
||||
pos_corners[:, :, 0] /= self.detector.pixel1
|
||||
pos_corners[:, :, 1] /= self.detector.pixel2
|
||||
pos = numpy.empty((self.shape[0], self.shape[1], 4, 2), dtype=numpy.float32)
|
||||
pos[:, :, 0, :] = pos_corners[:-1, :-1]
|
||||
pos[:, :, 1, :] = pos_corners[:-1, 1: ]
|
||||
pos[:, :, 2, :] = pos_corners[1: , 1: ]
|
||||
pos[:, :, 3, :] = pos_corners[1: , :-1]
|
||||
self.pos = pos
|
||||
self.delta0 = int((numpy.ceil(pos_corners[1:, :, 0]) - numpy.floor(pos_corners[:-1, :, 0])).max())
|
||||
self.delta1 = int((numpy.ceil(pos_corners[:, 1:, 1]) - numpy.floor(pos_corners[:, :-1, 1])).max())
|
||||
pixel_size = numpy.array([self.detector.pixel1, self.detector.pixel2])
|
||||
# make it a 4D array
|
||||
pixel_size.shape = 1, 1, 1, 2
|
||||
pixel_size.strides = 0, 0, 0, pixel_size.strides[-1]
|
||||
self.pos = self.detector.get_pixel_corners()[..., 1:] / pixel_size
|
||||
if self._shape_out is None:
|
||||
# if defined, it is probably because resize=False
|
||||
corner_pos = self.pos.view()
|
||||
corner_pos.shape = -1, 2
|
||||
pos1_min, pos2_min = corner_pos.min(axis=0)
|
||||
pos1_max, pos2_max = corner_pos.max(axis=0)
|
||||
self._shape_out = (int(ceil(pos1_max - pos1_min)),
|
||||
int(ceil(pos2_max - pos2_min)))
|
||||
self.offset1, self.offset2 = pos1_min, pos2_min
|
||||
pixel_delta = self.pos.view()
|
||||
pixel_delta.shape = -1, 4, 2
|
||||
self.delta1, self.delta2 = (pixel_delta.max(axis=1) - pixel_delta.min(axis=1)).max(axis=0)
|
||||
# pos_corners = numpy.empty((self.shape[0] + 1, self.shape[1] + 1, 2), dtype=numpy.float64)
|
||||
# d1 = numpy.outer(numpy.arange(self.shape[0] + 1, dtype=numpy.float64), numpy.ones(self.shape[1] + 1, dtype=numpy.float64)) - 0.5
|
||||
# d2 = numpy.outer(numpy.ones(self.shape[0] + 1, dtype=numpy.float64), numpy.arange(self.shape[1] + 1, dtype=numpy.float64)) - 0.5
|
||||
# pos_corners[:, :, 0], pos_corners[:, :, 1] = self.detector.calc_cartesian_positions(d1, d2)[:2]
|
||||
# pos_corners[:, :, 0] /= self.detector.pixel1
|
||||
# pos_corners[:, :, 1] /= self.detector.pixel2
|
||||
# pos = numpy.empty((self.shape[0], self.shape[1], 4, 2), dtype=numpy.float32)
|
||||
# pos[:, :, 0, :] = pos_corners[:-1, :-1]
|
||||
# pos[:, :, 1, :] = pos_corners[:-1, 1: ]
|
||||
# pos[:, :, 2, :] = pos_corners[1: , 1: ]
|
||||
# pos[:, :, 3, :] = pos_corners[1: , :-1]
|
||||
# self.pos = pos
|
||||
# self.delta0 = int((numpy.ceil(pos_corners[1:, :, 0]) - numpy.floor(pos_corners[:-1, :, 0])).max())
|
||||
# self.delta1 = int((numpy.ceil(pos_corners[:, 1:, 1]) - numpy.floor(pos_corners[:, :-1, 1])).max())
|
||||
return self.pos
|
||||
|
||||
@timeit
|
||||
def calc_size(self):
|
||||
"""
|
||||
def calc_size(self, use_cython=True):
|
||||
"""Calculate the number of pixels falling into every single bin and
|
||||
|
||||
@return: max of pixel falling into a single bin
|
||||
|
||||
Considering the "half-CCD" spline from ID11 which describes a (1025,2048) detector,
|
||||
the physical location of pixels should go from:
|
||||
[-17.48634 : 1027.0543, -22.768829 : 2028.3689]
|
||||
We chose to discard pixels falling outside the [0:1025,0:2048] range with a lose of intensity
|
||||
|
||||
"""
|
||||
if self.pos is None:
|
||||
pos = self.calc_pos()
|
||||
|
|
@ -167,32 +214,33 @@ class Distortion(object):
|
|||
if self.max_size is None:
|
||||
with self._sem:
|
||||
if self.max_size is None:
|
||||
if _distortion:
|
||||
self.bin_size = _distortion.calc_size(self.pos, self.shape, self.detector.mask)
|
||||
if _distortion and use_cython:
|
||||
self.bin_size = _distortion.calc_size(self.pos, self._shape_out, self.mask, (self.offset1, self.offset2))
|
||||
else:
|
||||
mask = self.detector.mask
|
||||
pos0min = numpy.floor(pos[:, :, :, 0].min(axis=-1)).astype(numpy.int32).clip(0, self.shape[0])
|
||||
pos1min = numpy.floor(pos[:, :, :, 1].min(axis=-1)).astype(numpy.int32).clip(0, self.shape[1])
|
||||
pos0max = (numpy.ceil(pos[:, :, :, 0].max(axis=-1)).astype(numpy.int32) + 1).clip(0, self.shape[0])
|
||||
pos1max = (numpy.ceil(pos[:, :, :, 1].max(axis=-1)).astype(numpy.int32) + 1).clip(0, self.shape[1])
|
||||
self.bin_size = numpy.zeros(self.shape, dtype=numpy.int32)
|
||||
max0 = 0
|
||||
max1 = 0
|
||||
for i in range(self.shape[0]):
|
||||
for j in range(self.shape[1]):
|
||||
mask = self.mask
|
||||
pos0min = (numpy.floor(pos[:, :, :, 0].min(axis=-1) - self.offset1).astype(numpy.int32)).clip(0, self._shape_out[0])
|
||||
pos1min = (numpy.floor(pos[:, :, :, 1].min(axis=-1) - self.offset2).astype(numpy.int32)).clip(0, self._shape_out[1])
|
||||
pos0max = (numpy.ceil(pos[:, :, :, 0].max(axis=-1) - self.offset1 + 1).astype(numpy.int32)).clip(0, self._shape_out[0])
|
||||
pos1max = (numpy.ceil(pos[:, :, :, 1].max(axis=-1) - self.offset2 + 1).astype(numpy.int32)).clip(0, self._shape_out[1])
|
||||
self.bin_size = numpy.zeros(self._shape_out, dtype=numpy.int32)
|
||||
# max0 = 0
|
||||
# max1 = 0
|
||||
for i in range(self.shape_in[0]):
|
||||
for j in range(self.shape_in[1]):
|
||||
if (mask is not None) and mask[i, j]:
|
||||
continue
|
||||
if (pos0max[i, j] - pos0min[i, j]) > max0:
|
||||
old = max0
|
||||
max0 = pos0max[i, j] - pos0min[i, j]
|
||||
logger.debug(old, "new max0", max0, i, j)
|
||||
if (pos1max[i, j] - pos1min[i, j]) > max1:
|
||||
old = max1
|
||||
max1 = pos1max[i, j] - pos1min[i, j]
|
||||
logger.debug(old, "new max1", max1, i, j)
|
||||
# if (pos0max[i, j] - pos0min[i, j]) > max0:
|
||||
# old = max0
|
||||
# max0 = pos0max[i, j] - pos0min[i, j]
|
||||
# logger.debug(old, "new max0", max0, i, j)
|
||||
# if (pos1max[i, j] - pos1min[i, j]) > max1:
|
||||
# old = max1
|
||||
# max1 = pos1max[i, j] - pos1min[i, j]
|
||||
# logger.debug(old, "new max1", max1, i, j)
|
||||
|
||||
self.bin_size[pos0min[i, j]:pos0max[i, j], pos1min[i, j]:pos1max[i, j]] += 1
|
||||
self.max_size = self.bin_size.max()
|
||||
return self.bin_size
|
||||
|
||||
def calc_init(self):
|
||||
"""
|
||||
|
|
@ -228,16 +276,16 @@ class Distortion(object):
|
|||
mask = self.detector.mask
|
||||
if _distortion:
|
||||
if self.method == "lut":
|
||||
self.lut = _distortion.calc_LUT(self.pos, self.shape, self.bin_size, max_pixel_size=(self.delta0, self.delta1))
|
||||
self.lut = _distortion.calc_LUT(self.pos, self.shape, self.bin_size, max_pixel_size=(self.delta1, self.delta2))
|
||||
else:
|
||||
self.lut = _distortion.calc_CSR(self.pos, self.shape, self.bin_size, max_pixel_size=(self.delta0, self.delta1))
|
||||
self.lut = _distortion.calc_CSR(self.pos, self.shape, self.bin_size, max_pixel_size=(self.delta1, self.delta2))
|
||||
else:
|
||||
lut = numpy.recarray(shape=(self.shape[0], self.shape[1], self.max_size), dtype=[("idx", numpy.uint32), ("coef", numpy.float32)])
|
||||
lut[:, :, :].idx = 0
|
||||
lut[:, :, :].coef = 0.0
|
||||
outMax = numpy.zeros(self.shape, dtype=numpy.uint32)
|
||||
idx = 0
|
||||
buffer = numpy.empty((self.delta0, self.delta1))
|
||||
buffer = numpy.empty((self.delta1, self.delta2))
|
||||
quad = Quad(buffer)
|
||||
for i in range(self.shape[0]):
|
||||
for j in range(self.shape[1]):
|
||||
|
|
@ -477,6 +525,7 @@ class Quad(object):
|
|||
return 0.5 * (K2 - K1) * (self.pCD * (K2 + K1) + 2 * self.cCD)
|
||||
else:
|
||||
return 0
|
||||
|
||||
def calc_area_DA(self, L1, L2):
|
||||
|
||||
if numpy.isfinite(self.pDA):
|
||||
|
|
|
|||
21252
pyFAI/ext/_distortion.c
21252
pyFAI/ext/_distortion.c
File diff suppressed because it is too large
Load Diff
|
|
@ -4,26 +4,30 @@
|
|||
# Project: Fast Azimuthal integration
|
||||
# https://github.com/pyFAI/pyFAI
|
||||
#
|
||||
# Copyright (C) European Synchrotron Radiation Facility, Grenoble, France
|
||||
# Copyright (C) 2013-2016 European Synchrotron Radiation Facility, Grenoble, France
|
||||
#
|
||||
# Principal author: Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
|
||||
#
|
||||
# This program is free software: you can redistribute it and/or modify
|
||||
# it under the terms of the GNU General Public License as published by
|
||||
# the Free Software Foundation, either version 3 of the License, or
|
||||
# (at your option) any later version.
|
||||
#
|
||||
# This program is distributed in the hope that it will be useful,
|
||||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
# GNU General Public License for more details.
|
||||
#
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
#
|
||||
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
# of this software and associated documentation files (the "Software"), to deal
|
||||
# in the Software without restriction, including without limitation the rights
|
||||
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
# copies of the Software, and to permit persons to whom the Software is
|
||||
# furnished to do so, subject to the following conditions:
|
||||
# .
|
||||
# The above copyright notice and this permission notice shall be included in
|
||||
# all copies or substantial portions of the Software.
|
||||
# .
|
||||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
# THE SOFTWARE.
|
||||
|
||||
__author__ = "Jerome Kieffer"
|
||||
__license__ = "GPLv3+"
|
||||
__license__ = "MIT"
|
||||
__date__ = "03/05/2016"
|
||||
__copyright__ = "2011-2016, ESRF"
|
||||
__contact__ = "jerome.kieffer@esrf.fr"
|
||||
|
|
@ -31,7 +35,7 @@ __contact__ = "jerome.kieffer@esrf.fr"
|
|||
import cython
|
||||
cimport numpy
|
||||
import numpy
|
||||
from cython cimport view
|
||||
from cython cimport view, floating
|
||||
from cython.parallel import prange
|
||||
from cpython.ref cimport PyObject, Py_XDECREF
|
||||
from libc.string cimport memset, memcpy
|
||||
|
|
@ -212,584 +216,6 @@ cdef inline void integrate(float[:, :] box, float start, float stop, float slope
|
|||
AA -= dA
|
||||
h += 1
|
||||
|
||||
cdef class Quad:
|
||||
"""
|
||||
Basic quadrilatere object
|
||||
|
||||
|
|
||||
|
|
||||
| xxxxxA
|
||||
| xxxxxxxI'xxxxxxxx x
|
||||
xxxxxxxxIxxxxxx | x
|
||||
Bxxxxxxxxxxxx | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x | | x
|
||||
x O| P A' x
|
||||
-----------------J------------------+--------------------------------L-----------------------
|
||||
x | x
|
||||
x | x
|
||||
x | x
|
||||
x | xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxD
|
||||
CxxxxxxxxxxxxxxxxxKxxxxx
|
||||
|
|
||||
|
|
||||
|
|
||||
|
|
||||
"""
|
||||
cdef float[:, :] box
|
||||
cdef float A0, A1, B0, B1, C0, C1, D0, D1, pAB, pBC, pCD, pDA, cAB, cBC, cCD, cDA, area
|
||||
cdef int offset0, offset1, box_size0, box_size1
|
||||
cdef bint has_area, has_slope
|
||||
|
||||
def __cinit__(self, float[:, :] buffer):
|
||||
self.box = buffer
|
||||
self.A0 = self.A1 = 0
|
||||
self.B0 = self.B1 = 0
|
||||
self.C0 = self.C1 = 0
|
||||
self.D0 = self.D1 = 0
|
||||
self.offset0 = self.offset1 = 0
|
||||
self.box_size0 = self.box_size1 = 0
|
||||
self.pAB = self.pBC = self.pCD = self.pDA = 0
|
||||
self.cAB = self.cBC = self.cCD = self.cDA = 0
|
||||
self.area = 0
|
||||
self.has_area = 0
|
||||
self.has_slope = 0
|
||||
|
||||
def reinit(self, A0, A1, B0, B1, C0, C1, D0, D1):
|
||||
self.box[:, :] = 0.0
|
||||
self.A0 = A0
|
||||
self.A1 = A1
|
||||
self.B0 = B0
|
||||
self.B1 = B1
|
||||
self.C0 = C0
|
||||
self.C1 = C1
|
||||
self.D0 = D0
|
||||
self.D1 = D1
|
||||
self.offset0 = (<int> floor(min(self.A0, self.B0, self.C0, self.D0)))
|
||||
self.offset1 = (<int> floor(min(self.A1, self.B1, self.C1, self.D1)))
|
||||
self.box_size0 = (<int> ceil(max(self.A0, self.B0, self.C0, self.D0))) - self.offset0
|
||||
self.box_size1 = (<int> ceil(max(self.A1, self.B1, self.C1, self.D1))) - self.offset1
|
||||
self.A0 -= self.offset0
|
||||
self.A1 -= self.offset1
|
||||
self.B0 -= self.offset0
|
||||
self.B1 -= self.offset1
|
||||
self.C0 -= self.offset0
|
||||
self.C1 -= self.offset1
|
||||
self.D0 -= self.offset0
|
||||
self.D1 -= self.offset1
|
||||
self.pAB = self.pBC = self.pCD = self.pDA = 0
|
||||
self.cAB = self.cBC = self.cCD = self.cDA = 0
|
||||
self.area = 0
|
||||
self.has_area = 0
|
||||
self.has_slope = 0
|
||||
|
||||
def __repr__(self):
|
||||
res = ["offset %i,%i size %i, %i" % (self.offset0, self.offset1, self.box_size0, self.box_size1), "box:"]
|
||||
for i in range(self.box_size0):
|
||||
line = ""
|
||||
for j in range(self.box_size1):
|
||||
line += "\t%.3f"%self.box[i,j]
|
||||
res.append(line)
|
||||
return os.linesep.join(res)
|
||||
|
||||
cpdef float get_box(self, int i, int j):
|
||||
return self.box[i,j]
|
||||
cpdef int get_offset0(self):
|
||||
return self.offset0
|
||||
cpdef int get_offset1(self):
|
||||
return self.offset1
|
||||
cpdef int get_box_size0(self):
|
||||
return self.box_size0
|
||||
cpdef int get_box_size1(self):
|
||||
return self.box_size1
|
||||
|
||||
cpdef init_slope(self):
|
||||
if not self.has_slope:
|
||||
if self.B0 != self.A0:
|
||||
self.pAB = (self.B1 - self.A1) / (self.B0 - self.A0)
|
||||
self.cAB = self.A1 - self.pAB * self.A0
|
||||
if self.C0 != self.B0:
|
||||
self.pBC = (self.C1 - self.B1) / (self.C0 - self.B0)
|
||||
self.cBC = self.B1 - self.pBC * self.B0
|
||||
if self.D0 != self.C0:
|
||||
self.pCD = (self.D1 - self.C1) / (self.D0 - self.C0)
|
||||
self.cCD = self.C1 - self.pCD * self.C0
|
||||
if self.A0 != self.D0:
|
||||
self.pDA = (self.A1 - self.D1) / (self.A0 - self.D0)
|
||||
self.cDA = self.D1 - self.pDA * self.D0
|
||||
self.has_slope = 1
|
||||
|
||||
cpdef float calc_area_AB(self, float I1, float I2):
|
||||
if self.B0 != self.A0:
|
||||
return 0.5 * (I2 - I1) * (self.pAB * (I2 + I1) + 2 * self.cAB)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
cpdef float calc_area_BC(self, float J1, float J2):
|
||||
if self.B0 != self.C0:
|
||||
return 0.5 * (J2 - J1) * (self.pBC * (J1 + J2) + 2 * self.cBC)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
cpdef float calc_area_CD(self, float K1, float K2):
|
||||
if self.C0 != self.D0:
|
||||
return 0.5 * (K2 - K1) * (self.pCD * (K2 + K1) + 2 * self.cCD)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
cpdef float calc_area_DA(self, float L1, float L2):
|
||||
if self.D0 != self.A0:
|
||||
return 0.5 * (L2 - L1) * (self.pDA * (L1 + L2) + 2 * self.cDA)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
cpdef float calc_area(self):
|
||||
if not self.has_area:
|
||||
self.area = 0.5 * ((self.C0 - self.A0) * (self.D1 - self.B1) - (self.C1 - self.A1) * (self.D0 - self.B0))
|
||||
self.has_area = 1
|
||||
return self.area
|
||||
|
||||
def populate_box(self):
|
||||
cdef int i0, i1
|
||||
cdef float area, value
|
||||
if not self.has_slope:
|
||||
self.init_slope()
|
||||
integrate(self.box, self.B0, self.A0, self.pAB, self.cAB)
|
||||
integrate(self.box, self.A0, self.D0, self.pDA, self.cDA)
|
||||
integrate(self.box, self.D0, self.C0, self.pCD, self.cCD)
|
||||
integrate(self.box, self.C0, self.B0, self.pBC, self.cBC)
|
||||
area = self.calc_area()
|
||||
for i0 in range(self.box_size0):
|
||||
for i1 in range(self.box_size1):
|
||||
value = self.box[i0, i1] / area
|
||||
self.box[i0, i1] = value
|
||||
if value < 0.0:
|
||||
print(self.box)
|
||||
self.box[:, :] = 0
|
||||
print("AB")
|
||||
self.integrate(self.B0, self.A0, self.pAB, self.cAB)
|
||||
print(self.box)
|
||||
self.box[:, :] = 0
|
||||
print("DA")
|
||||
self.integrate(self.A0, self.D0, self.pDA, self.cDA)
|
||||
print(self.box)
|
||||
self.box[:, :] = 0
|
||||
print("CD")
|
||||
self.integrate(self.D0, self.C0, self.pCD, self.cCD)
|
||||
print(self.box)
|
||||
self.box[:, :] = 0
|
||||
print("BC")
|
||||
self.integrate(self.C0, self.B0, self.pBC, self.cBC)
|
||||
print(self.box)
|
||||
print(self)
|
||||
raise RuntimeError()
|
||||
|
||||
def integrate(self, float start, float stop, float slope, float intercept):
|
||||
cdef int i, h = 0
|
||||
cdef float P, dP, A, AA, dA, sign
|
||||
# print(start, stop, calc_area(start, stop)
|
||||
if start < stop: # positive contribution
|
||||
P = ceil(start)
|
||||
dP = P - start
|
||||
# print("Integrate", start, P, stop, calc_area(start, stop)
|
||||
if P > stop: # start and stop are in the same unit
|
||||
A = calc_area(start, stop, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
dA = (stop - start) # always positive
|
||||
# print(AA, sign, dA
|
||||
h = 0
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(start)), h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
else:
|
||||
if dP > 0:
|
||||
A = calc_area(start, P, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
h = 0
|
||||
dA = dP
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(P)) - 1, h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
# subsection P1->Pn
|
||||
for i in range((<int> floor(P)), (<int> floor(stop))):
|
||||
A = calc_area(i, i + 1, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
|
||||
h = 0
|
||||
dA = 1.0
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[i, h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
# Section Pn->B
|
||||
P = floor(stop)
|
||||
dP = stop - P
|
||||
if dP > 0:
|
||||
A = calc_area(P, stop, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
h = 0
|
||||
dA = fabs(dP)
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(P)), h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
elif start > stop: # negative contribution. Nota is start=stop: no contribution
|
||||
P = floor(start)
|
||||
if stop > P: # start and stop are in the same unit
|
||||
A = calc_area(start, stop, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
dA = (start - stop) # always positive
|
||||
h = 0
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(start)), h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
else:
|
||||
dP = P - start
|
||||
if dP < 0:
|
||||
A = calc_area(start, P, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
h = 0
|
||||
dA = fabs(dP)
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(P)) , h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
# subsection P1->Pn
|
||||
for i in range((<int> start), (<int> ceil(stop)), -1):
|
||||
A = calc_area(i, i - 1, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
h = 0
|
||||
dA = 1
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[i - 1, h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
# Section Pn->B
|
||||
P = ceil(stop)
|
||||
dP = stop - P
|
||||
if dP < 0:
|
||||
A = calc_area(P, stop, slope, intercept)
|
||||
if A != 0:
|
||||
AA = fabs(A)
|
||||
sign = A / AA
|
||||
h = 0
|
||||
dA = fabs(dP)
|
||||
while AA > 0:
|
||||
if dA > AA:
|
||||
dA = AA
|
||||
AA = -1
|
||||
self.box[(<int> floor(stop)), h] += sign * dA
|
||||
AA -= dA
|
||||
h += 1
|
||||
|
||||
|
||||
class Distortion(object):
|
||||
"""
|
||||
|
||||
This class applies a distortion correction on an image.
|
||||
|
||||
It is also able to apply an inversion of the correction.
|
||||
|
||||
"""
|
||||
def __init__(self, detector="detector", shape=None):
|
||||
"""
|
||||
@param detector: detector instance or detector name
|
||||
"""
|
||||
if isinstance(detector, six.string_types):
|
||||
self.detector = detector_factory(detector)
|
||||
else: # we assume it is a Detector instance
|
||||
self.detector = detector
|
||||
if shape:
|
||||
self.shape = shape
|
||||
elif "max_shape" in dir(self.detector):
|
||||
self.shape = self.detector.max_shape
|
||||
self.shape = tuple([int(i) for i in self.shape])
|
||||
self._sem = threading.Semaphore()
|
||||
self.lut_size = None
|
||||
self.pos = None
|
||||
self.LUT = None
|
||||
self.delta0 = self.delta1 = None # max size of an pixel on a regular grid ...
|
||||
|
||||
def __repr__(self):
|
||||
return os.linesep.join(["Distortion correction for detector:",
|
||||
self.detector.__repr__()])
|
||||
|
||||
def calc_pos(self):
|
||||
if self.pos is None:
|
||||
with self._sem:
|
||||
if self.pos is None:
|
||||
pos_corners = numpy.empty((self.shape[0] + 1, self.shape[1] + 1, 2), dtype=numpy.float64)
|
||||
d1 = expand2d(numpy.arange(self.shape[0] + 1.0), self.shape[1] + 1, False) - 0.5
|
||||
d2 = expand2d(numpy.arange(self.shape[1] + 1.0), self.shape[0] + 1, True) - 0.5
|
||||
p = self.detector.calc_cartesian_positions(d1, d2)
|
||||
if p[-1] is not None:
|
||||
logger.warning("makes little sense to correct for distortion non-flat detectors: %s"%self.detector)
|
||||
pos_corners[:, :, 0], pos_corners[:, :, 1] = p[:2]
|
||||
pos_corners[:, :, 0] /= self.detector.pixel1
|
||||
pos_corners[:, :, 1] /= self.detector.pixel2
|
||||
pos = numpy.empty((self.shape[0], self.shape[1], 4, 2), dtype=numpy.float32)
|
||||
pos[:, :, 0, :] = pos_corners[:-1, :-1]
|
||||
pos[:, :, 1, :] = pos_corners[:-1, 1:]
|
||||
pos[:, :, 2, :] = pos_corners[1:, 1:]
|
||||
pos[:, :, 3, :] = pos_corners[1:, :-1]
|
||||
self.pos = pos
|
||||
self.delta0 = int((numpy.ceil(pos_corners[1:, :, 0]) - numpy.floor(pos_corners[:-1, :, 0])).max())
|
||||
self.delta1 = int((numpy.ceil(pos_corners[:, 1:, 1]) - numpy.floor(pos_corners[:, :-1, 1])).max())
|
||||
return self.pos
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
def calc_LUT_size(self):
|
||||
"""
|
||||
Considering the "half-CCD" spline from ID11 which describes a (1025,2048) detector,
|
||||
the physical location of pixels should go from:
|
||||
[-17.48634 : 1027.0543, -22.768829 : 2028.3689]
|
||||
We chose to discard pixels falling outside the [0:1025,0:2048] range with a lose of intensity
|
||||
|
||||
We keep self.pos: pos_corners will not be compatible with systems showing non adjacent pixels (like some xpads)
|
||||
|
||||
"""
|
||||
cdef int i, j, k, l, shape0, shape1
|
||||
cdef numpy.ndarray[numpy.float32_t, ndim = 4] pos
|
||||
cdef int[:, :] pos0min, pos1min, pos0max, pos1max
|
||||
cdef numpy.ndarray[numpy.int32_t, ndim = 2] lut_size
|
||||
if self.pos is None:
|
||||
pos = self.calc_pos()
|
||||
else:
|
||||
pos = self.pos
|
||||
if self.lut_size is None:
|
||||
with self._sem:
|
||||
if self.lut_size is None:
|
||||
shape0, shape1 = self.shape
|
||||
pos0min = numpy.floor(pos[:, :, :, 0].min(axis=-1)).astype(numpy.int32).clip(0, self.shape[0])
|
||||
pos1min = numpy.floor(pos[:, :, :, 1].min(axis=-1)).astype(numpy.int32).clip(0, self.shape[1])
|
||||
pos0max = (numpy.ceil(pos[:, :, :, 0].max(axis=-1)).astype(numpy.int32) + 1).clip(0, self.shape[0])
|
||||
pos1max = (numpy.ceil(pos[:, :, :, 1].max(axis=-1)).astype(numpy.int32) + 1).clip(0, self.shape[1])
|
||||
lut_size = numpy.zeros(self.shape, dtype=numpy.int32)
|
||||
with nogil:
|
||||
for i in range(shape0):
|
||||
for j in range(shape1):
|
||||
for k in range(pos0min[i, j], pos0max[i, j]):
|
||||
for l in range(pos1min[i, j], pos1max[i, j]):
|
||||
lut_size[k, l] += 1
|
||||
self.lut_size = lut_size.max()
|
||||
return lut_size
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
@cython.cdivision(True)
|
||||
def calc_LUT(self):
|
||||
cdef:
|
||||
int i, j, ms, ml, ns, nl, shape0, shape1, delta0, delta1, buffer_size, i0, i1, size
|
||||
int offset0, offset1, box_size0, box_size1
|
||||
numpy.int32_t k, idx = 0
|
||||
float A0, A1, B0, B1, C0, C1, D0, D1, pAB, pBC, pCD, pDA, cAB, cBC, cCD, cDA, area, value
|
||||
float[:, :, :, :] pos
|
||||
numpy.ndarray[lut_point, ndim = 3] lut
|
||||
numpy.ndarray[numpy.int32_t, ndim = 2] outMax = numpy.zeros(self.shape, dtype=numpy.int32)
|
||||
float[:, :] buffer
|
||||
shape0, shape1 = self.shape
|
||||
|
||||
if self.lut_size is None:
|
||||
self.calc_LUT_size()
|
||||
if self.LUT is None:
|
||||
with self._sem:
|
||||
if self.LUT is None:
|
||||
pos = self.pos
|
||||
if (self.lut_size == 0): #fix 271
|
||||
raise RuntimeError("The look-up table has dimension (0) which is a non-sense."
|
||||
+ "Did you mask out all pixel or is your image out of the geometry range ?")
|
||||
lut = numpy.recarray(shape=(self.shape[0], self.shape[1], self.lut_size), dtype=[("idx", numpy.int32), ("coef", numpy.float32)])
|
||||
size = self.shape[0] * self.shape[1] * self.lut_size * sizeof(lut_point)
|
||||
memset(&lut[0, 0, 0], 0, size)
|
||||
logger.info("LUT shape: (%i,%i,%i) %.3f MByte" % (lut.shape[0], lut.shape[1], lut.shape[2], size / 1.0e6))
|
||||
buffer = numpy.empty((self.delta0, self.delta1), dtype=numpy.float32)
|
||||
buffer_size = self.delta0 * self.delta1 * sizeof(float)
|
||||
logger.info("Max pixel size: %ix%i; Max source pixel in target: %i" % (buffer.shape[1], buffer.shape[0], self.lut_size))
|
||||
with nogil:
|
||||
# i,j, idx are indexes of the raw image uncorrected
|
||||
for i in range(shape0):
|
||||
for j in range(shape1):
|
||||
# reinit of buffer
|
||||
buffer[:, :] = 0
|
||||
A0 = pos[i, j, 0, 0]
|
||||
A1 = pos[i, j, 0, 1]
|
||||
B0 = pos[i, j, 1, 0]
|
||||
B1 = pos[i, j, 1, 1]
|
||||
C0 = pos[i, j, 2, 0]
|
||||
C1 = pos[i, j, 2, 1]
|
||||
D0 = pos[i, j, 3, 0]
|
||||
D1 = pos[i, j, 3, 1]
|
||||
offset0 = (<int> floor(min(A0, B0, C0, D0)))
|
||||
offset1 = (<int> floor(min(A1, B1, C1, D1)))
|
||||
box_size0 = (<int> ceil(max(A0, B0, C0, D0))) - offset0
|
||||
box_size1 = (<int> ceil(max(A1, B1, C1, D1))) - offset1
|
||||
A0 -= <float> offset0
|
||||
A1 -= <float> offset1
|
||||
B0 -= <float> offset0
|
||||
B1 -= <float> offset1
|
||||
C0 -= <float> offset0
|
||||
C1 -= <float> offset1
|
||||
D0 -= <float> offset0
|
||||
D1 -= <float> offset1
|
||||
if B0 != A0:
|
||||
pAB = (B1 - A1) / (B0 - A0)
|
||||
cAB = A1 - pAB * A0
|
||||
else:
|
||||
pAB = cAB = 0.0
|
||||
if C0 != B0:
|
||||
pBC = (C1 - B1) / (C0 - B0)
|
||||
cBC = B1 - pBC * B0
|
||||
else:
|
||||
pBC = cBC = 0.0
|
||||
if D0 != C0:
|
||||
pCD = (D1 - C1) / (D0 - C0)
|
||||
cCD = C1 - pCD * C0
|
||||
else:
|
||||
pCD = cCD = 0.0
|
||||
if A0 != D0:
|
||||
pDA = (A1 - D1) / (A0 - D0)
|
||||
cDA = D1 - pDA * D0
|
||||
else:
|
||||
pDA = cDA = 0.0
|
||||
integrate(buffer, B0, A0, pAB, cAB)
|
||||
integrate(buffer, A0, D0, pDA, cDA)
|
||||
integrate(buffer, D0, C0, pCD, cCD)
|
||||
integrate(buffer, C0, B0, pBC, cBC)
|
||||
area = 0.5 * ((C0 - A0) * (D1 - B1) - (C1 - A1) * (D0 - B0))
|
||||
for ms in range(box_size0):
|
||||
ml = ms + offset0
|
||||
if ml < 0 or ml >= shape0:
|
||||
continue
|
||||
for ns in range(box_size1):
|
||||
# ms,ns are indexes of the corrected image in short form, ml & nl are the same
|
||||
nl = ns + offset1
|
||||
if nl < 0 or nl >= shape1:
|
||||
continue
|
||||
value = buffer[ms, ns] / area
|
||||
if value <= 0:
|
||||
continue
|
||||
k = outMax[ml, nl]
|
||||
lut[ml, nl, k].idx = idx
|
||||
lut[ml, nl, k].coef = value
|
||||
outMax[ml, nl] = k + 1
|
||||
idx += 1
|
||||
self.LUT = lut.reshape(self.shape[0] * self.shape[1], self.lut_size)
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
def correct(self, image):
|
||||
"""
|
||||
Correct an image based on the look-up table calculated ...
|
||||
|
||||
@param image: 2D-array with the image
|
||||
@return: corrected 2D image
|
||||
"""
|
||||
cdef:
|
||||
int i, j, lshape0, lshape1, idx, size
|
||||
float coef
|
||||
lut_point[:, :] LUT
|
||||
float[:] lout, lin
|
||||
if self.LUT is None:
|
||||
self.calc_LUT()
|
||||
LUT = self.LUT
|
||||
lshape0 = LUT.shape[0]
|
||||
lshape1 = LUT.shape[1]
|
||||
img_shape = image.shape
|
||||
if (img_shape[0] < self.shape[0]) or (img_shape[1] < self.shape[1]):
|
||||
new_image = numpy.zeros(self.shape, dtype=numpy.float32)
|
||||
new_image[:img_shape[0], :img_shape[1]] = image
|
||||
image = new_image
|
||||
logger.warning("Patching image as image is %ix%i and spline is %ix%i" % (img_shape[1], img_shape[0], self.shape[1], self.shape[0]))
|
||||
|
||||
out = numpy.zeros(self.shape, dtype=numpy.float32)
|
||||
lout = out.ravel()
|
||||
lin = numpy.ascontiguousarray(image.ravel(), dtype=numpy.float32)
|
||||
size = lin.size
|
||||
for i in prange(lshape0, nogil=True, schedule="static"):
|
||||
for j in range(lshape1):
|
||||
idx = LUT[i, j].idx
|
||||
coef = LUT[i, j].coef
|
||||
if coef <= 0:
|
||||
continue
|
||||
if idx >= size:
|
||||
with gil:
|
||||
logger.warning("Accessing %i >= %i !!!" % (idx, size))
|
||||
continue
|
||||
lout[i] += lin[idx] * coef
|
||||
return out[:img_shape[0], :img_shape[1]]
|
||||
|
||||
@timeit
|
||||
def uncorrect(self, image):
|
||||
"""
|
||||
Take an image which has been corrected and transform it into it's raw (with loss of information)
|
||||
|
||||
@param image: 2D-array with the image
|
||||
@return: uncorrected 2D image and a mask (pixels in raw image
|
||||
"""
|
||||
if self.LUT is None:
|
||||
self.calc_LUT()
|
||||
out = numpy.zeros(self.shape, dtype=numpy.float32)
|
||||
mask = numpy.zeros(self.shape, dtype=numpy.int8)
|
||||
lmask = mask.ravel()
|
||||
lout = out.ravel()
|
||||
lin = image.ravel()
|
||||
tot = self.LUT.coef.sum(axis=-1)
|
||||
for idx in range(self.LUT.shape[0]):
|
||||
t = tot[idx]
|
||||
if t <= 0:
|
||||
lmask[idx] = 1
|
||||
continue
|
||||
val = lin[idx] / t
|
||||
lout[self.LUT[idx].idx] += val * self.LUT[idx].coef
|
||||
return out, mask
|
||||
|
||||
################################################################################
|
||||
# Functions used in python classes from PyFAI.distortion
|
||||
################################################################################
|
||||
|
|
@ -797,42 +223,53 @@ class Distortion(object):
|
|||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
def calc_size(float[:, :, :, :] pos not None, shape, numpy.int8_t[:, :] mask=None):
|
||||
def calc_size(floating[:, :, :, ::1] pos not None,
|
||||
shape,
|
||||
numpy.int8_t[:, ::1] mask=None,
|
||||
offset=None):
|
||||
"""
|
||||
Calculate the number of items per output pixel
|
||||
|
||||
@param pos: 4D array with position in space
|
||||
@param shape: shape of the output array
|
||||
@param mask: input data mask
|
||||
@param offset: 2-tuple of float with the minimal index of
|
||||
@return: number of input element per output elements
|
||||
"""
|
||||
cdef:
|
||||
int i, j, k, l, shape0, shape1, min0, min1, max0, max1
|
||||
int i, j, k, l, shape_out0, shape_out1, shape_in0, shape_in1, min0, min1, max0, max1
|
||||
numpy.ndarray[numpy.int32_t, ndim = 2] lut_size = numpy.zeros(shape, dtype=numpy.int32)
|
||||
float A0, A1, B0, B1, C0, C1, D0, D1
|
||||
bint do_mask = mask is None
|
||||
shape0, shape1 = shape
|
||||
|
||||
if do_mask:
|
||||
assert mask.shape[0] == shape0
|
||||
assert mask.shape[1] == shape1
|
||||
|
||||
float A0, A1, B0, B1, C0, C1, D0, D1, offset0, offset1
|
||||
bint do_mask = mask is not None
|
||||
|
||||
shape_in0, shape_in1 = pos.shape[0], pos.shape[1]
|
||||
shape_out0, shape_out1 = shape
|
||||
|
||||
if do_mask and ((mask.shape[0] != shape_in0) or (mask.shape[1] != shape_in1)):
|
||||
err = 'Mismatch between shape of detector (%s, %s) and shape of mask (%s, %s)' % (shape_in0, shape_in1, mask.shape[0], mask.shape[1])
|
||||
logger.error(err)
|
||||
raise RuntimeError(err)
|
||||
|
||||
if offset is not None:
|
||||
offset0, offset1 = offset
|
||||
|
||||
with nogil:
|
||||
for i in range(shape0):
|
||||
for j in range(shape1):
|
||||
for i in range(shape_in0):
|
||||
for j in range(shape_in1):
|
||||
if do_mask and mask[i, j]:
|
||||
continue
|
||||
A0 = pos[i, j, 0, 0]
|
||||
A1 = pos[i, j, 0, 1]
|
||||
B0 = pos[i, j, 1, 0]
|
||||
B1 = pos[i, j, 1, 1]
|
||||
C0 = pos[i, j, 2, 0]
|
||||
C1 = pos[i, j, 2, 1]
|
||||
D0 = pos[i, j, 3, 0]
|
||||
D1 = pos[i, j, 3, 1]
|
||||
min0 = clip(<int> floor(min(A0, B0, C0, D0)), 0, shape0)
|
||||
min1 = clip(<int> floor(min(A1, B1, C1, D1)), 0, shape1)
|
||||
max0 = clip(<int> ceil(max(A0, B0, C0, D0)) + 1, 0, shape0)
|
||||
max1 = clip(<int> ceil(max(A1, B1, C1, D1)) + 1, 0, shape1)
|
||||
A0 = pos[i, j, 0, 0] - offset0
|
||||
A1 = pos[i, j, 0, 1] - offset1
|
||||
B0 = pos[i, j, 1, 0] - offset0
|
||||
B1 = pos[i, j, 1, 1] - offset1
|
||||
C0 = pos[i, j, 2, 0] - offset0
|
||||
C1 = pos[i, j, 2, 1] - offset1
|
||||
D0 = pos[i, j, 3, 0] - offset0
|
||||
D1 = pos[i, j, 3, 1] - offset1
|
||||
min0 = clip(<int> floor(min(A0, B0, C0, D0)), 0, shape_out0)
|
||||
min1 = clip(<int> floor(min(A1, B1, C1, D1)), 0, shape_out1)
|
||||
max0 = clip(<int> ceil(max(A0, B0, C0, D0)) + 1, 0, shape_out0)
|
||||
max1 = clip(<int> ceil(max(A1, B1, C1, D1)) + 1, 0, shape_out1)
|
||||
for k in range(min0, max0):
|
||||
for l in range(min1, max1):
|
||||
lut_size[k, l] += 1
|
||||
|
|
|
|||
|
|
@ -4,27 +4,27 @@
|
|||
# Project: Azimuthal integration
|
||||
# https://github.com/pyFAI/pyFAI
|
||||
#
|
||||
# Copyright (C) 2015 European Synchrotron Radiation Facility, Grenoble, France
|
||||
# Copyright (C) 2013-2015 European Synchrotron Radiation Facility, Grenoble, France
|
||||
#
|
||||
# Principal author: Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
|
||||
#
|
||||
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
# of this software and associated documentation files (the "Software"), to deal
|
||||
# in the Software without restriction, including without limitation the rights
|
||||
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
# copies of the Software, and to permit persons to whom the Software is
|
||||
# furnished to do so, subject to the following conditions:
|
||||
#
|
||||
# The above copyright notice and this permission notice shall be included in
|
||||
# all copies or substantial portions of the Software.
|
||||
#
|
||||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
# THE SOFTWARE.
|
||||
# Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
# of this software and associated documentation files (the "Software"), to deal
|
||||
# in the Software without restriction, including without limitation the rights
|
||||
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
# copies of the Software, and to permit persons to whom the Software is
|
||||
# furnished to do so, subject to the following conditions:
|
||||
# .
|
||||
# The above copyright notice and this permission notice shall be included in
|
||||
# all copies or substantial portions of the Software.
|
||||
# .
|
||||
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
||||
# THE SOFTWARE.
|
||||
|
||||
from __future__ import absolute_import, division, print_function
|
||||
|
||||
|
|
@ -33,7 +33,7 @@ __author__ = "Jérôme Kieffer"
|
|||
__contact__ = "Jerome.Kieffer@ESRF.eu"
|
||||
__license__ = "MIT"
|
||||
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
|
||||
__date__ = "29/01/2016"
|
||||
__date__ = "03/05/2016"
|
||||
|
||||
|
||||
import unittest
|
||||
|
|
@ -42,7 +42,7 @@ import fabio
|
|||
from .utilstest import UtilsTest, getLogger
|
||||
logger = getLogger(__file__)
|
||||
from .. import detectors
|
||||
from ..ext import _distortion
|
||||
from .. import distortion
|
||||
|
||||
|
||||
class TestHalfCCD(unittest.TestCase):
|
||||
|
|
@ -57,7 +57,8 @@ class TestHalfCCD(unittest.TestCase):
|
|||
self.halfFrelon = UtilsTest.getimage(self.__class__.halfFrelon)
|
||||
self.splineFile = UtilsTest.getimage(self.__class__.splineFile)
|
||||
self.det = detectors.FReLoN(self.splineFile)
|
||||
self.dis = _distortion.Distortion(self.det)
|
||||
self.dis = distortion.Distortion(self.det, self.det.shape, resize=False,
|
||||
mask=numpy.zeros(self.det.shape, "int8"))
|
||||
self.fit2d = fabio.open(self.fit2dFile).data
|
||||
self.raw = fabio.open(self.halfFrelon).data
|
||||
|
||||
|
|
@ -65,6 +66,13 @@ class TestHalfCCD(unittest.TestCase):
|
|||
unittest.TestCase.tearDown(self)
|
||||
self.fit2dFile = self.halfFrelon = self.splineFile = self.det = self.dis = self.fit2d = self.raw = None
|
||||
|
||||
def test_size(self):
|
||||
self.dis.reset(prepare=False)
|
||||
ny = self.dis.calc_size(False)
|
||||
self.dis.reset(prepare)
|
||||
cy = self.dis.calc_size(True)
|
||||
self.assertEqual(abs(ny - cy).max(), 0, "equivalence of the cython and numpy model")
|
||||
|
||||
def test_vs_fit2d(self):
|
||||
"""
|
||||
Compare spline correction vs fit2d's code
|
||||
|
|
@ -93,6 +101,7 @@ class TestHalfCCD(unittest.TestCase):
|
|||
|
||||
def suite():
|
||||
testsuite = unittest.TestSuite()
|
||||
testsuite.addTest(TestHalfCCD("test_size"))
|
||||
testsuite.addTest(TestHalfCCD("test_vs_fit2d"))
|
||||
# testsuite.addTest(test_azim_halfFrelon("test_numpy_vs_fit2d"))
|
||||
# testsuite.addTest(test_azim_halfFrelon("test_cythonSP_vs_fit2d"))
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
4
setup.py
4
setup.py
|
|
@ -28,7 +28,7 @@ __author__ = "Jerome Kieffer"
|
|||
__contact__ = "Jerome.Kieffer@ESRF.eu"
|
||||
__license__ = "GPLv3+"
|
||||
__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
|
||||
__date__ = "11/04/2016"
|
||||
__date__ = "03/05/2016"
|
||||
__status__ = "stable"
|
||||
|
||||
install_warning = True
|
||||
|
|
@ -204,7 +204,7 @@ ext_modules = [
|
|||
Extension('relabel'),
|
||||
Extension("bilinear", can_use_openmp=True),
|
||||
Extension('_distortion', can_use_openmp=True),
|
||||
Extension('_distortionCSR', can_use_openmp=True),
|
||||
# Extension('_distortionCSR', can_use_openmp=True),
|
||||
Extension('_bispev', can_use_openmp=True),
|
||||
Extension('_convolution', can_use_openmp=True),
|
||||
Extension('_blob'),
|
||||
|
|
|
|||
Loading…
Reference in New Issue