Skip to content
Snippets Groups Projects
Commit ab2a05d6 authored by Rafael Gort's avatar Rafael Gort
Browse files

processing 4 quadrants in parallel

parent e9214e27
No related branches found
No related tags found
1 merge request!91WIP: Dssc methods as
Hide whitespace changes
Inline Side-by-side
src/toolbox_scs/detectors/dssc.py
+ 15
− 9
View file @ ab2a05d6
......@@ -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
......
src/toolbox_scs/detectors/image_manipulation.py deleted 100644 → 0
+ 0
− 336
View file @ e9214e27
"""
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts