From 834399eba2e25f58451f89c7c54af9faad94d047 Mon Sep 17 00:00:00 2001
From: Rafael Gort <rafael.gort@xfel.eu>
Date: Wed, 26 Aug 2020 15:44:51 +0200
Subject: [PATCH] Added DSSC related geometry methods and image manipulation
routines. quick code structuring. not tested
---
src/toolbox_scs/detectors/__init__.py | 4 +-
src/toolbox_scs/detectors/dssc.py | 195 ++++++++--
src/toolbox_scs/detectors/dssc_geometry.py | 52 +++
.../detectors/image_manipulation.py | 336 ++++++++++++++++++
4 files changed, 562 insertions(+), 25 deletions(-)
create mode 100644 src/toolbox_scs/detectors/dssc_geometry.py
create mode 100644 src/toolbox_scs/detectors/image_manipulation.py
diff --git a/src/toolbox_scs/detectors/__init__.py b/src/toolbox_scs/detectors/__init__.py
index 5059545..c2998c9 100644
--- a/src/toolbox_scs/detectors/__init__.py
+++ b/src/toolbox_scs/detectors/__init__.py
@@ -10,7 +10,7 @@ from .dssc_misc import (
from .dssc_processing import (
bin_data)
from .dssc import (
- DSSCBinner, DSSCFormatter, DSSCAnalyzer)
+ DSSCBinner, DSSCFormatter, DSSC_Generator)
from .azimuthal_integrator import AzimuthalIntegrator
__all__ = (
@@ -33,7 +33,7 @@ __all__ = (
# Classes
"DSSCBinner",
"DSSCFormatter",
- "DSSCAnalyzer",
+ "DSSC_Generator",
"AzimuthalIntegrator",
# Variables
)
diff --git a/src/toolbox_scs/detectors/dssc.py b/src/toolbox_scs/detectors/dssc.py
index 189f7d7..9da1971 100644
--- a/src/toolbox_scs/detectors/dssc.py
+++ b/src/toolbox_scs/detectors/dssc.py
@@ -25,11 +25,13 @@ from .dssc_data import (
save_xarray, load_xarray, save_attributes_h5,
search_files, get_data_formatted)
from .dssc_misc import (
- load_dssc_info, get_xgm_formatted, load_geom)
+ load_dssc_info, get_xgm_formatted)
from .dssc_processing import (
bin_data, create_empty_dataset)
+from .image_manipulation import (
+ convert_hex2cart)
-__all__ = ["DSSCBinner", "DSSCFormatter", "DSSCAnalyzer"]
+__all__ = ["DSSCBinner", "DSSCFormatter", "DSSC_Generator"]
log = logging.getLogger(__name__)
@@ -257,7 +259,7 @@ class DSSCFormatter:
if os.path.exists(filepath):
try:
self._filenames = search_files(filepath)
- except ToolBoxFileError as err:
+ except ToolBoxFileError:
log.info("path did not contain any files")
else:
log.info("path did not exist")
@@ -338,30 +340,177 @@ class DSSCFormatter:
save_attributes_h5(filename, self.attributes)
-class DSSCAnalyzer:
+class DSSC_Generator:
+
+ # geomtry = None
+ tile_geodat = np.zeros((16, 2, 2), dtype='int')
+ # definitions: [[voffset, hoffset], [is_flip, axis]]
+ # axis= 0, flip up-down
+ # axis= 1, flip left-right
+ # fractional shifts not included so far
+ tile_geodat[0, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[1, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[2, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[3, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[4, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[5, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[6, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[7, :, :] = [[0, 8], [1, 1]]
+ tile_geodat[8, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[9, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[10, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[11, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[12, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[13, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[14, :, :] = [[0, 8], [1, 0]]
+ tile_geodat[15, :, :] = [[0, 8], [1, 0]]
+
+ module_geodat = np.zeros((16, 2), dtype='int')
+ module_geodat[0, :] = [0, 0]
+ module_geodat[1, :] = [19, 0]
+ module_geodat[2, :] = [18, 0]
+ module_geodat[3, :] = [19, 0]
+
+ module_geodat[4, :] = [0, 0]
+ module_geodat[5, :] = [19, 0]
+ module_geodat[6, :] = [19, 0]
+ module_geodat[7, :] = [18, 0]
+
+ module_geodat[11, :] = [0, 0]
+ module_geodat[10, :] = [19, 0]
+ module_geodat[9, :] = [18, 0]
+ module_geodat[8, :] = [19, 0]
+
+ module_geodat[15, :] = [0, 0]
+ module_geodat[14, :] = [19, 0]
+ module_geodat[13, :] = [19, 0]
+ module_geodat[12, :] = [18, 0]
+
+ quad_geodat = np.zeros((4, 2), dtype='int')
+ quad_geodat[0, :] = [10, 17] # [quadrant 1, insert, gap]
+ quad_geodat[1, :] = [8, 26] # [quadrant 2, insert, gap]
+ quad_geodat[2, :] = [-9, 19] # [quadrant 3, insert, gap]
+ quad_geodat[3, :] = [-9, 25] # [quadrant 4, insert, gap]
+
def __init__(self):
+ pass
+
+ def read_geometry(self, filename):
+ pass
+
+ def save_geometry(self, filename):
+ pass
+
+ def full_detector(self, dssc_data):
+ print('assemble full detector: ')
+ # determine dimension of full detector
+ [qdimy, qdimx] = q1.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]]
+ ) + 2 * qdimy
+ fdimx = np.max(
+ [np.abs(self.quad_geodat[2, 0]) + self.quad_geodat[2, 1],
+ np.abs(self.quad_geodat[0, 0]) + self.quad_geodat[0, 1]]
+ ) + 2 * qdimx
+
+ # assemble full detector
+ full_dssc = np.zeros((fdimy, fdimx))
+ full_dssc[self.quad_geodat[1, 0]:qdimy+self.quad_geodat[1, 0], 0:qdimx ] \
+ = self.quadrant(2, dssc_data)
+ 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)#
+
+ 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)
+
+ 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)
+
+ return full_dssc
+
+ def quadrant(self, q, dssc_data):
+ """
+ q: quadrant index 0-3
+ dssc_data: (15, 128, 512)
+ """
+ print("collecting quadrant: Q" + str(q))
+ if q == 1:
+ module0 = self.convert_module(0, dssc_data[0, :, :])
+ mdimx = module0.shape[1]
+ qu = [module0]
+ for n in range(1, 4):
+ qu.append(np.zeros((self.module_geodat[n, 0], mdimx), dtype='float32'))
+ qu.append(self.convert_module(n, dssc_data[n, :, :]))
+ elif q == 2:
+ module0 = self.convert_module(4, dssc_data[4, :, :])
+ mdimx = module0.shape[1]
+ qu = [module0]
+ for n in range(5, 8):
+ qu.append(np.zeros((self.module_geodat[n, 0], mdimx), dtype='float32'))
+ qu.append(self.convert_module(n, dssc_data[n, :, :]))
+ elif q == 3:
+ module0 = self.convert_module(11, dssc_data[11, :, :])
+ mdimx = module0.shape[1]
+ qu = [module0]
+ for n in range(10, 7, -1):
+ qu.append(np.zeros((self.module_geodat[n, 0], mdimx), dtype='float32'))
+ qu.append(self.convert_module(n, dssc_data[n, :, :]))
+ elif q == 4:
+ module0 = self.convert_module(15, dssc_data[15, :, :])
+ mdimx = module0.shape[1]
+ qu = [module0]
+ for n in range(14, 11, -1):
+ qu.append(np.zeros((self.module_geodat[n, 0], mdimx), dtype='float32'))
+ qu.append(self.convert_module(n, dssc_data[n, :, :]))
+
+ return np.vstack(qu)
+
+ def convert_module(self, m, module_hex):
"""
- Placeholder for future class handling image operation/manipulation
- (the actual data evaluation).
+ convert hex to cartesian
+ add gap according to geometry file
+ flip according to geometry file
+ returns the module in cartesian coordinates with the correct gap between tiles
+ """
+ print("processing module: " + str(m))
+ tile_hex = self.split_module(module_hex)
+ offset, flip_module = self._get_tile_geodat(m)
- Parameters
- ----------
- ToDo ...
+ module_cart = np.zeros((120 + offset[0], 2*276 + offset[1]), dtype='float32')
+ module_cart[0:120, 0:276] = convert_hex2cart(tile_hex[0, :, :])
+ module_cart[0:120, 276 + offset[1]:552 + offset[1]] = convert_hex2cart(tile_hex[1, :, :])
- Returns
- -------
- ToDo ...
+ if flip_module[0]:
+ module_cart = np.flip(module_cart, axis=flip_module[1])
+ #print('flipped module: ' + str(flip_module[1]))
- Raises
- ------
- ToDo....
+ return module_cart
- Example
- -------
- """
- # self.distance = distance
- # self.pxpitchh = 236
- # self.pxpitchv = 204
- # self.aspect = self.pxpitchv/self.pxpitchh
- # .....
+ def split_module(self, module_hex):
+ # return two tiles
+ tile_hex = np.zeros((2, 128, 256), dtype='float32')
+ tile_hex[0, :, :] = module_hex[:, 0:256]
+ tile_hex[1, :, :] = module_hex[:, 256:512]
+ return tile_hex
+
+ def _get_tile_geodat(self, m):
+ return self.tile_geodat[m, 0, :], self.tile_geodat[m, 1, :]
+
+ def _default_geometry(self):
pass
+
+# for all modules do
+# 1) fix pedestals by ASIC mask and add/subtract pedestal offset
+# 2) split module into its tiles
+# 3) convert tiles from hex to cart
+# 4) reposition tiles according to offsets (from geometry file)
+# 5) figure out orientation and flip up-down or left-right (flag in geometry file)
diff --git a/src/toolbox_scs/detectors/dssc_geometry.py b/src/toolbox_scs/detectors/dssc_geometry.py
new file mode 100644
index 0000000..af14223
--- /dev/null
+++ b/src/toolbox_scs/detectors/dssc_geometry.py
@@ -0,0 +1,52 @@
+"""
+ DSSC-geometry related sub-routines.
+
+ comment: contributions should comply with pep8 code structure guidelines.
+"""
+import logging
+
+import numpy as np
+import skimage.util as sku
+import skimage.transform as skt
+
+log = logging.getLogger(__name__)
+
+
+def slice_module(module):
+ tile = np.zeros((2, 128, 256), dtype='float32')
+ tile[0, :, :] = module[:, 0:256]
+ tile[1, :, :] = module[:, 256:512]
+ print(tile.shape)
+ return tile
+
+
+def data_from_symmetry(data, center=[584.0, 544.0], threshold=0.0):
+ """
+ Parameters
+ ----------
+ data: np.ndarray
+ 2-dimensional numpy-array containing the raw image data.
+ center: list
+ a list with two values indicating the physical center, i.e. the
+ position of the beam.
+ threshold: float
+ The lower magnitude limit, below which data from the image should
+ be replaced with data from the rotated image.
+
+ Returns
+ -------
+ data_rotated: np.ndarray
+ A 2-dimensional numpy-array containing data from the rotated image,
+ where it is missing in the original one.
+ """
+ image = np.nan_to_num(sku.img_as_float(data))
+ shape = np.shape(image)
+
+ shifts = [(shape[i] - 2*(shape[i]/2.0 - center[i])) for i in range(2)]
+
+ tform = skt.EuclideanTransform(rotation=np.pi, translation=shifts[::-1])
+ data_rotated = skt.warp(image, tform)
+
+ mask = (image <= threshold)
+ data_rotated[~mask] = 0.0
+ return data_rotated
diff --git a/src/toolbox_scs/detectors/image_manipulation.py b/src/toolbox_scs/detectors/image_manipulation.py
new file mode 100644
index 0000000..e4ecc9a
--- /dev/null
+++ b/src/toolbox_scs/detectors/image_manipulation.py
@@ -0,0 +1,336 @@
+"""
+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