Add cython transformation code for JF strixel sensors

setup.py

@@ -24,6 +24,12 @@ ext_modules = [
                             '-ftree-vectorize', '-frename-registers'],
         extra_link_args=['-fopenmp'],
     ),
+    Extension(
+        "cal_tools.jfalgs",
+        ['src/cal_tools/jfalgs.pyx'],
+        extra_compile_args=['-O3', '-march=native', '-ftree-vectorize',
+                            '-frename-registers']
+    ),
     Extension(
         "cal_tools.gotthard2.gotthard2algs",
         ["src/cal_tools/gotthard2/gotthard2algs.pyx"],

src/cal_tools/jfalgs.pyx

+
+from cython cimport boundscheck, wraparound, cdivision
+from cython.view cimport contiguous
+
+
+ctypedef fused jf_data_t:
+    unsigned short  # raw pixel data
+    float  # corrected pixel data
+    unsigned int  # mask data
+    unsigned char  # gain data
+
+
+DEF STRIXEL_Y = 86
+DEF STRIXEL_X = 1024 * 3 + 18
+
+
+def strixel_shape():
+    return STRIXEL_Y, STRIXEL_X
+
+
+@boundscheck(False)
+@wraparound(False)
+@cdivision(True)
+def strixel_transform(jf_data_t[:, :, ::contiguous] data,
+                      jf_data_t[:, :, ::contiguous] out = None):
+    """Reorder raw data to physical strixel sensor layout. """
+
+    if data.shape[1] < 256 or data.shape[2] < 256:
+        raise ValueError('Pixel shape of data may not be below (256, 256)')
+
+    if out is None:
+        import numpy as np
+        out = np.zeros((data.shape[0], STRIXEL_Y, STRIXEL_X), dtype=np.float32)
+    elif data.shape[0] > out.shape[0]:
+        raise ValueError('Cell shape of data exceeds out')
+    elif out.shape[1] < STRIXEL_Y or out.shape[2] < STRIXEL_X:
+        raise ValueError(f'Pixel shape of out may not be below '
+                         f'({STRIXEL_Y}, {STRIXEL_X})')
+
+    cdef int cell, yin, xin, xout, yout, igap
+
+    for cell in range(data.shape[0]):
+        # Normal pixels.
+        for yin in range(256):
+            yout = yin // 3
+
+            for xin in range(1024) :
+                xout = 774 * (xin // 256) + 3 * (xin % 256) + yin % 3
+                out[cell, yout, xout] = data[cell, yin, xin]
+
+        # Gap pixels.
+        for yin in range(256):
+            yout = 2 * (yin // 6)
+
+            for igap in range(3) :
+                # Left side of the gap.
+                xin = igap * 256 + 255
+                xout = igap * 774 + 765 + yin % 6
+                out[cell, yout, xout] = data[cell, yin, xin]
+                out[cell, yout+1, xout] = data[cell, yin, xin]
+
+                # Right side of the gap.
+                xin = igap * 256 + 255 + 1
+                xout = igap * 774 + 765 + 11 - yin % 6
+                out[cell, yout, xout] = data[cell, yin, xin]
+                out[cell, yout+1, xout] = data[cell, yin, xin]
+
+    return out