From ab2a05d62f1c5433551cb4126b67b11e1e493d43 Mon Sep 17 00:00:00 2001
From: Rafael Gort <rafael.gort@xfel.eu>
Date: Thu, 17 Sep 2020 13:33:28 +0200
Subject: [PATCH] processing 4 quadrants in parallel
---
src/toolbox_scs/detectors/dssc.py | 24 +-
.../detectors/image_manipulation.py | 336 ------------------
2 files changed, 15 insertions(+), 345 deletions(-)
delete mode 100644 src/toolbox_scs/detectors/image_manipulation.py
diff --git a/src/toolbox_scs/detectors/dssc.py b/src/toolbox_scs/detectors/dssc.py
index 213940c..9ee7b47 100644
--- a/src/toolbox_scs/detectors/dssc.py
+++ b/src/toolbox_scs/detectors/dssc.py
@@ -389,9 +389,11 @@ class DSSC_Generator:
def full_detector(self, dssc_data):
print('assemble full detector: ')
- # determine dimension of full detector
- q2 = self.quadrant(2, dssc_data)
- [qdimy, qdimx] = q2.shape
+ # assemble full detector
+ q_data = joblib.Parallel(n_jobs=4)(joblib.delayed(self.quadrant)(i, dssc_data) for i in [1,2,3,4])
+
+ # determine dimensions
+ [qdimy, qdimx] = q_data[0].shape
fdimy = np.max(
[np.abs(self.quad_geodat[1, 0]) + self.quad_geodat[1, 1],
np.abs(self.quad_geodat[3, 0]) + self.quad_geodat[3, 1]]
@@ -401,26 +403,30 @@ class DSSC_Generator:
np.abs(self.quad_geodat[0, 0]) + self.quad_geodat[0, 1]]
) + 2 * qdimx
- # assemble full detector
full_dssc = np.zeros((fdimy, fdimx))
+
+ # quadrant 2
full_dssc[self.quad_geodat[1, 0]:qdimy+self.quad_geodat[1, 0], 0:qdimx ] \
- = self.quadrant(2, dssc_data)
+ = q_data[1]
+
+ # quadrant 1
full_dssc[qdimy+self.quad_geodat[1, 0]+self.quad_geodat[1, 1]:2*qdimy+self.quad_geodat[1, 0]+self.quad_geodat[1, 1],
self.quad_geodat[0, 0]:self.quad_geodat[0, 0]+qdimx] \
- = self.quadrant(1, dssc_data)#
+ = q_data[0]
+ # quadrant 4
dummy1 = [qdimy + self.quad_geodat[1, 0] + self.quad_geodat[1, 1] + self.quad_geodat[3, 0],
2 * qdimy + self.quad_geodat[1, 0] + self.quad_geodat[1, 1] + self.quad_geodat[3, 0]]
dummy2 = [self.quad_geodat[0, 0] + self.quad_geodat[0, 1] + qdimx,
self.quad_geodat[0, 0] + self.quad_geodat[0, 1] + 2*qdimx]
- full_dssc[dummy1[0]:dummy1[1], dummy2[0]:dummy2[1]] = self.quadrant(4, dssc_data)
+ full_dssc[dummy1[0]:dummy1[1], dummy2[0]:dummy2[1]] = q_data[3]
+ # quadrant 3
dummy3 = [self.quad_geodat[1, 0] + self.quad_geodat[1, 1] + self.quad_geodat[3, 0] - self.quad_geodat[3, 1],
qdimy+self.quad_geodat[1, 0] + self.quad_geodat[1, 1] + self.quad_geodat[3, 0] - self.quad_geodat[3, 1]]
dummy4 = [self.quad_geodat[0, 0] + self.quad_geodat[0, 1] + qdimx + self.quad_geodat[2, 0],
self.quad_geodat[0, 0] + self.quad_geodat[0, 1] + 2*qdimx + self.quad_geodat[2, 0]]
-
- full_dssc[dummy3[0]:dummy3[1],dummy4[0]:dummy4[1]] = self.quadrant(3, dssc_data)
+ full_dssc[dummy3[0]:dummy3[1],dummy4[0]:dummy4[1]] = q_data[2]
return full_dssc
diff --git a/src/toolbox_scs/detectors/image_manipulation.py b/src/toolbox_scs/detectors/image_manipulation.py
deleted file mode 100644
index e4ecc9a..0000000
--- a/src/toolbox_scs/detectors/image_manipulation.py
+++ /dev/null
@@ -1,336 +0,0 @@
-"""
-Collection of Filters for image processing:
-Thiran filters for delay shifters. The filter is used for image shearing.
-"""
-import logging
-
-import numpy as np
-import math as math
-from skimage.util import pad as skpad
-from skimage import data as skdata
-
-log = logging.getLogger(__name__)
-
-
-def n_choose_k(n, k):
- """
- :n choose k problem: returns the number of ways to choose a subset k from a
- set of n elements this is eqivalent to:
- factorial(n)/factorial(k)/factorial(n-k)
-
- written by A. Scherz
- """
- return math.factorial(n) / math.factorial(k) / math.factorial(n - k)
-
-
-def thiran_type_I(order=3):
- return Thiran(1, order)
-
-
-def thiran_type_II(order=2):
- return Thiran(2, order)
-
-
-class Thiran:
- """
- there are two types of 1D fractional filters based on Thiran filters. both
- types can have an order N=1 to inf. some types, provide simple algorithms
- for applying the delay to a signal. for type II and orders larger than N.
- he transfer and impulse response fct have to be calculated and are more
- computational expensive
- The generator defaults the Thiran filter of type I and order N=3
- For more details see of the filters and algorithms,
- see pg. 685, IEEE Trans. on Im. Proc. 17, 679 (2008)
-
- written by A. Scherz, XFEL, 2012-11-26
- updated by A. Scherz, XFEL, 2019-02-11 for python
- added to SpeckleToolbox by A. Scherz, XFEL, 2020-08-06
- """
-
- type_I = staticmethod(thiran_type_I)
-
- type_II = staticmethod(thiran_type_II)
-
- def __init__(self, filter_type=1, order=3, delay=0):
- self._type = filter_type
- self._n = order
- self._delay = delay
- self._b = np.zeros(order)
- self._a_minus_delay = 0
- self._a_plus_delay = 0
-
- def __str__(self):
- print('Thiran 1D Fractional filter')
- print(' Type : ' + self._type.__str__())
- print(' Order : ' + self._n.__str__())
- print(' Delay : ' + self._delay.__str__())
-
- return ''
-
- def apply_delay(self, signal):
- # periodic boundaries are important. add padding of at least the order
- # N left and right
- # this should be applied before calling this method for type I
- shifted_signal = np.array(signal)
-
- if self._type == 1: # type I Thiran
- for j in range(len(shifted_signal) - self._n - 1, self._n - 1, -1):
- for k in range(1, self._n + 1):
- shifted_signal[j] += self._b[k - 1] \
- * (shifted_signal[j - k] - shifted_signal[j + k])
- else: # type II of order N = 2
- for j in range(self._n - 1, len(shifted_signal) - self._n - 1, 1):
- shifted_signal[j] += self._a_plus_delay \
- * (shifted_signal[j - 1] - shifted_signal[j + 1])
- for j in range(len(shifted_signal) - self._n - 1, self._n - 1, -1):
- shifted_signal[j] += self._a_minus_delay \
- * (shifted_signal[j + 1] - shifted_signal[j - 1])
-
- return shifted_signal # shifted signal
-
- def delay(self, tau):
- if (tau < -0.5) or (tau > 0.5):
- raise ValueError('WARNING out of range: delay: 0.0 - 0.5')
- self._delay = tau
- # update coefficients for delay
- self._coefficients()
- return self
-
- def _coefficients(self):
- for k in range(1, self._n + 1):
- p_k = 1
- for l in range(self._n + 1):
- p_k *= (self._delay - l) / (self._delay - l - k)
- # n_choose_k function is factorial(n)/factorial(k)/factorial(n-k)
- self._b[k - 1] = (-1) ** k * p_k * n_choose_k(self._n, k)
- self._a_plus_delay = \
- (-4 + self._delay**2 + np.sqrt(12 - 3*self._delay**2))
- self._a_minus_delay = self._a_plus_delay
- self._a_plus_delay /= (3*self._delay + 2 + self._delay**2)
- self._a_minus_delay /= (-3*self._delay + 2 + self._delay**2)
-
-
-def is_even(num):
- """
- is num an even number?
- :param num: integer
- :return: boolean True = num is even number
- """
- return not (num & 1)
-
-
-def next_power_of_2(x=1):
- """
- Returns the next higher power of 2 compared to x
- :param x: current dimension
- :return: next higher power of 2
- """
- return 1 << (x - 1).bit_length()
-
-
-def get_im_center(im: np.array) -> (int, int):
- """
- the center of image is the DC components with respect to discrete Fourier
- transformation
- :returns center pos (int):
- """
- return int(im.shape[0] / 2), int(im.shape[1] / 2)
-
-
-def pad4fft(im: np.array = None,
- dim: [None, bool, int] = None,
- fill: [None, str, float] = None) -> np.array:
- """
- pad image using dim parameter and pad image with given fill
- :param im:
- :param dim: None == up to next power of 2, True/False == oversampling, or
- int with given dimension
- :param fill: None == zeros, str == special format see skimage.pad, float ==
- fill value
- :return: object itself
- """
- if im is not None:
- im = np.array(im).squeeze()
- else:
- im = skdata.lena()
- #print('image size is ' + str(im.shape))
- #print('dimension is supposed to be ' + str(dim))
-
- if dim is None:
- dim = max((next_power_of_2(im.shape[0]), next_power_of_2(im.shape[1])))
- dim = (dim, dim)
- elif type(dim) is bool:
- oversampling = dim
- print('is boolean option')
- dim = max((next_power_of_2(im.shape[0]), next_power_of_2(im.shape[1])))
- dim *= 2 if oversampling else 1 # another factor 2 for oversampling
- dim = (dim, dim)
- elif (dim[0] < im.shape[0]) and (dim[1] < im.shape[1]):
- return im
-
- #print('dimension is ' + str(dim))
-
- w0 = ((dim[0] - im.shape[0]) // 2)
- w1 = ((dim[1] - im.shape[1]) // 2)
-
- width0 = (w0, w0) if is_even(im.shape[0]) else (w0, w0 + 1)
- width1 = (w1, w1) if is_even(im.shape[1]) else (w1, w1 + 1)
-
- if fill is None:
- im = skpad(im, (width0, width1), mode='constant', constant_values=0)
- elif type(fill) is str:
- im = skpad(im, (width0, width1), mode=fill)
- else:
- im = skpad(im, (width0, width1), mode='constant', constant_values=fill)
-
- return im
-
-
-def shear_im(im: np.array, delay: float,
- im_axis: int = 1,
- thiran: Thiran = thiran_type_I()) -> np.array:
- """
- : direction (str): 1=='hor' or 0=='ver', required
- : thiran: Thiran type 1 of order 3 is default (default)
- """
- im = np.float32(im)
- cx, cy = get_im_center(im)
- # print(im.shape)
- # print(cx, cy)
- if im_axis == 0:
- for k in range(2 * cy):
- shift = delay * (k - cy)
- # full pixel shifts
- im[:, k] = np.roll(im[:, k], round(shift))
- # now shift a pixel fraction
- im[:, k] = thiran.delay(shift-round(shift)).apply_delay(im[:, k])
- # thiran_filter(im[:, k], thiran(order, shift - round(shift)))
- # print('apply vertical shear')
- elif im_axis == 1:
- for k in range(2 * cx):
- shift = delay * (k - cx)
- # full pixel shifts
- im[k, :] = np.roll(im[k, :], round(shift))
- # now shift a pixel fraction: shift between +/- [0-0.5]
- im[k, :] = thiran.delay(shift - round(shift)).apply_delay(im[k, :])
- # thiran_filter(im[k, :], thiran(order, shift - round(shift)))
- #print('apply horizontal shear')
- else:
- # TODO: throw exception
- print('direction undefined')
-
- return im
-
-
-def convert_cart2hex(im, thiran=thiran_type_II()):
- """
- [ imH ] = ConversionCart2Hex( imC, orderN, showImage )
- imC: image in cartesian coord.
- thiran: Thiran filter type 1 of order 3 is default
- show: optional, if true an image representation in hex is shown.
- note: show hex image may take considerable amount of time for larger
- images!)
- Nov.28, 2012, A. Scherz
- """
- # first step: zero pad the image for extending boundaries
-
- im_hex = pad4fft(im, dim=[2*im.shape[0], 2*im.shape[1]], fill='reflect')
- cx, cy = get_im_center(im_hex)
-
- sqrt3 = math.sqrt(3)
- shear1 = sqrt3 - math.sqrt(6 / sqrt3)
- shear2 = math.sqrt(sqrt3 / 6)
- shear3 = 2 - math.sqrt(6 / sqrt3)
- # this is conversion from cartesian to hexagonal lattice
- im_hex = shear_im(im_hex, -shear1, im_axis=0, thiran=thiran)
- im_hex = shear_im(im_hex, -shear2, im_axis=1, thiran=thiran)
- im_hex = shear_im(im_hex, -shear3, im_axis=0, thiran=thiran)
-
- # reorder of pixels required
- height = int(np.round(cart2hex([im.shape[0], 0]))[0])
- width = int(np.round(cart2hex([0, im.shape[1]]))[1])
- y0 = height % 2
-
- # print(height, width, y0)
- for y in range(-height // 2, height // 2 + y0 + 1):
- im_hex[y - y0 + cx, :] = np.roll(im_hex[y - y0 + cx, :],
- math.floor((y - y0) / 2))
-
- # now crop image
- im_hex = im_hex[-height // 2 + cx: height // 2 + cx + y0,
- -width // 2 + cy: width // 2 + cy + 1]
-
- # TODO check height and width
- # implement show
-
- return im_hex
-
-
-def convert_hex2cart(im_hex, thiran=thiran_type_II()):
- """
- [ imH ] = ConversionCart2Hex( imC, orderN, showImage )
- imC: image in cartesian coord.
- thiran: Thiran filter type 1 of order 3 is default
- show: optional, if true an image representation in hex is shown.
- note: show hex image may take considerable amount of time for larger
- images!)
- Nov.28, 2012, A. Scherz
- """
- # first step: zero pad the image for extending boundaries
-
- im = pad4fft(im_hex,
- dim=[2*im_hex.shape[0], 2*im_hex.shape[1]],
- fill='reflect')
- cx, cy = get_im_center(im)
-
- # need to skew image to hexagonal shape before apply shear operations
- height = im_hex.shape[0]
- y0 = height % 2
- for y in range(- height // 2, height // 2 + y0 + 1):
- im[y - y0 + cx, :] = np.roll(im[y - y0 + cx, :],
- - math.floor((y + y0) / 2))
-
- sqrt3 = math.sqrt(3)
- shear1 = sqrt3 - math.sqrt(6 / sqrt3)
- shear2 = math.sqrt(sqrt3 / 6)
- shear3 = 2 - math.sqrt(6 / sqrt3)
- # this is conversion from cartesian to hexagonal lattice
- im = shear_im(im, shear3, im_axis=0, thiran=thiran)
- im = shear_im(im, shear2, im_axis=1, thiran=thiran)
- im = shear_im(im, shear1, im_axis=0, thiran=thiran)
-
- # the image should now form a rectangular box of following size:
- height = int(np.round(hex2cart([im_hex.shape[0], 0]))[0])
- width = int(np.round(hex2cart([0, im_hex.shape[1]]))[1])
- y0 = height % 2
-
-# print(height, width, y0)
-# for y in range(-height // 2, height // 2 + y0 + 1):
-# im_hex[y - y0 + cx, :] = np.roll(im_hex[y - y0 + cx, :],
-# math.floor((y - y0) / 2))
-
- # now crop image
- im = im[-height // 2 + cx: height // 2 + cx + y0,
- -width // 2 + cy: width // 2 + cy + 1]
-
- # TODO check height and width
- # implement show
-
- return im
-
-
-def cart2hex(r_cart):
- # TODO: this function can eventually move somewhere else
- sqrt3 = math.sqrt(3)
- # conversion matrix
- rhex = math.sqrt(2 / sqrt3) * np.array([[1, 1 / 2], [0, sqrt3 / 2]])
-
- return np.dot(rhex, r_cart)
-
-
-def hex2cart(r_hex):
- sqrt3 = math.sqrt(3)
- rhex = math.sqrt(2 / sqrt3) * np.array([[1, 1 / 2], [0, sqrt3 / 2]])
- rcart = np.linalg.inv(rhex)
-
- return np.dot(rcart, r_hex)
--
GitLab