diff --git a/setup.py b/setup.py index aae7823b411bdbf59ab410e4e3b2d4d4fe510821..f5c28b08850f7c7bce53ae148976a5b76322043a 100644 --- a/setup.py +++ b/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"], diff --git a/src/cal_tools/jfalgs.pyx b/src/cal_tools/jfalgs.pyx new file mode 100644 index 0000000000000000000000000000000000000000..cc39b3e5a84f8cdead19467f1b04d55d1f739319 --- /dev/null +++ b/src/cal_tools/jfalgs.pyx @@ -0,0 +1,68 @@ + +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