init LocCorr - grasshopper CMOS & pusirobot motors based corrector

LocCorr/CMakeLists.txt

@@ -0,0 +1,76 @@
+cmake_minimum_required(VERSION 3.0)
+set(PROJ loccorr)
+set(MINOR_VERSION "1")
+set(MID_VERSION "0")
+set(MAJOR_VERSION "0")
+set(VERSION "${MAJOR_VERSION}.${MID_VERSION}.${MINOR_VERSION}")
+
+project(${PROJ} VERSION ${PROJ_VERSION} LANGUAGES C)
+
+# default flags
+set(CMAKE_C_FLAGS "-O3 -std=gnu99 -march=native")
+
+set(CMAKE_COLOR_MAKEFILE ON)
+
+# here is one of two variants: all .c in directory or .c files in list
+aux_source_directory(${CMAKE_CURRENT_SOURCE_DIR} SOURCES)
+
+# cmake -DEBUG=1 -> debugging
+if(DEFINED EBUG AND EBUG EQUAL 1)
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wextra -Wall -Werror -W")
+    set(CMAKE_BUILD_TYPE DEBUG)
+    set(CMAKE_VERBOSE_MAKEFILE "ON")
+    add_definitions(-DEBUG)
+else()
+    set(CMAKE_BUILD_TYPE RELEASE)
+    set(CMAKE_VERBOSE_MAKEFILE "ON")
+endif()
+
+###### pkgconfig ######
+# pkg-config modules (for pkg-check-modules)
+set(MODULES usefull_macros cfitsio)
+
+# find packages:
+SET(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR})
+find_package(PkgConfig REQUIRED)
+find_package(FLYCAP    REQUIRED)
+pkg_check_modules(${PROJ} REQUIRED ${MODULES})
+
+include(FindOpenMP)
+if(OPENMP_FOUND)
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+    add_definitions(-DOMP_FOUND)
+endif()
+if(NOT DEFINED PROCESSOR_COUNT)
+    set(PROCESSOR_COUNT 2) # by default 2 cores
+    execute_process(COMMAND getconf _NPROCESSORS_ONLN OUTPUT_VARIABLE PROCESSOR_COUNT OUTPUT_STRIP_TRAILING_WHITESPACE)
+endif()
+message("In multithreaded operations will use ${PROCESSOR_COUNT} threads")
+
+
+# change wrong behaviour with install prefix
+if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT AND CMAKE_INSTALL_PREFIX MATCHES "/usr/local")
+else()
+    message("Change default install path to /usr/local")
+    set(CMAKE_INSTALL_PREFIX "/usr/local")
+endif()
+message("Install dir prefix: ${CMAKE_INSTALL_PREFIX}")
+
+# exe file
+add_executable(${PROJ} ${SOURCES})
+# -I
+include_directories(${${PROJ}_INCLUDE_DIRS} ${FLYCAP_INCLUDE_DIRS})
+# -L
+link_directories(${${PROJ}_LIBRARY_DIRS})
+# -D
+add_definitions(${CFLAGS} -DLOCALEDIR=\"${LOCALEDIR}\"
+    -DPACKAGE_VERSION=\"${VERSION}\" -DGETTEXT_PACKAGE=\"${PROJ}\"
+    -DMINOR_VERSION=\"${MINOR_VERSION}\" -DMID_VERSION=\"${MID_VERSION}\"
+    -DMAJOR_VERSION=\"${MAJOR_VESION}\" -DTHREAD_NUMBER=${PROCESSOR_COUNT})
+
+# -l
+target_link_libraries(${PROJ} ${${PROJ}_LIBRARIES} ${FLYCAP_LIBRARIES} -lm)
+
+# Installation of the program
+INSTALL(TARGETS ${PROJ} DESTINATION "bin")

LocCorr/FindFLYCAP.cmake

@@ -0,0 +1,40 @@
+# - Try to find libflycapture
+# Once done this will define
+#
+#  FLYCAP_FOUND - system has libflycapture
+#  FLYCAP_INCLUDE_DIR - include directory
+#  FLYCAP_LIBRARIES - Link these to use libflycapture
+
+# Copyright (c) 2021, Edward V. Emelianov <edward.emelianoff@gmail.com>
+#
+# Redistribution and use is allowed according to the terms of the GPLv2/GPLv3.
+
+include(GNUInstallDirs)
+
+find_path(FLYCAP_INCLUDE_DIR  FlyCapture2.h
+    PATH_SUFFIXES libflycapture flycapture
+    PATHS /usr/include /usr/local/include /opt/include /opt/local/include
+)
+#find_path(FLYCAP_LIBRARY_DIR  libflycapture.so
+#    PATHS /lib /lib64 /usr/lib /usr/lib64 /opt/lib /opt/lib64 /usr/local/lib /usr/local/lib64
+#)
+find_library(FLYCAP_LIBRARY NAMES flycapture
+    PATHS /lib /lib64 /usr/lib /usr/lib64 /opt/lib /opt/lib64 /usr/local/lib /usr/local/lib64
+)
+find_library(FLYCAP_LIBRARYC NAMES flycapture-c
+    PATHS /lib /lib64 /usr/lib /usr/lib64 /opt/lib /opt/lib64 /usr/local/lib /usr/local/lib64
+)
+
+find_package_handle_standard_args(FLYCAP DEFAULT_MSG FLYCAP_INCLUDE_DIR FLYCAP_LIBRARY FLYCAP_LIBRARYC)
+
+if(FLYCAP_FOUND)
+    set(FLYCAP_INCLUDE_DIRS ${FLYCAP_INCLUDE_DIR})
+    set(FLYCAP_LIBRARIES ${FLYCAP_LIBRARY} ${FLYCAP_LIBRARYC})
+##    set(FLYCAP_LIBRARY_DIRS ${FLYCAP_LIBRARY_DIR})
+#    message("FLYCAP include dir = ${FLYCAP_INCLUDE_DIRS}")
+#    message("FLYCAP lib = ${FLYCAP_LIBRARIES}")
+#    message("FLYCAP libdir = ${FLYCAP_LIBRARY_DIRS}")
+    mark_as_advanced(FLYCAP_INCLUDE_DIRS FLYCAP_LIBRARIES FLYCAP_LIBRARY_DIRS)
+else()
+    message("FLYCAP not found")
+endif(FLYCAP_FOUND)

LocCorr/Readme.md

@@ -0,0 +1,3 @@
+# Local corrector
+
+Takes fresh images from single file or directory using inotify, find objects and calculate their centroids.

LocCorr/binmorph.c

@@ -0,0 +1,571 @@
+/*
+ * binmorph.c - functions for morphological operations on binary image
+ *
+ * Copyright 2015 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h> // memcpy
+#include <math.h>
+#include <sys/time.h>
+#include <usefull_macros.h>
+#include "binmorph.h"
+#include "imagefile.h"
+
+// global arrays for erosion/dilation masks
+static uint8_t *ER = NULL, *DIL = NULL;
+
+/*
+ * =================== AUXILIARY FUNCTIONS ===================>
+ */
+
+/**
+ * This function inits masks arrays for erosion and dilation
+ * You may call it yourself or it will be called when one of
+ * `erosion` or `dilation` functions will be ran first time
+ */
+static void morph_init(){
+    if(ER) return;
+    ER = MALLOC(uint8_t, 256);
+    DIL = MALLOC(uint8_t, 256);
+    OMP_FOR()
+    for(int i = 0; i < 256; i++){
+        ER[i]  = i & ((i << 1) | 1) & ((i >> 1) | (0x80)); // don't forget that << and >> set borders to zero
+        DIL[i] = i | (i << 1) | (i >> 1);
+    }
+}
+
+/*
+ * =================== MORPHOLOGICAL OPERATIONS ===================>
+ */
+
+/**
+ * Remove all non-4-connected pixels
+ * @param image (i) - input image
+ * @param W, H      - size of binarized image (in pixels)
+ * @return allocated memory area with converted input image
+ */
+uint8_t *filter4(uint8_t *image, int W, int H){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    uint8_t *ret = MALLOC(uint8_t, W*H);
+    int W0 = (W + 7) / 8; // width in bytes
+    int w = W0-1, h = H-1;
+    {int y = 0;
+    // top of image, y = 0
+    #define IM_UP
+    #include "fc_filter.h"
+    #undef IM_UP
+    }
+    {
+    // mid of image, y = 1..h-1
+    #include "fc_filter.h"
+    }
+    {
+    // image bottom, y = h
+    int y = h;
+    #define IM_DOWN
+    #include "fc_filter.h"
+    #undef IM_DOWN
+    }
+    return ret;
+}
+
+/**
+ * Remove all non-8-connected pixels (single points)
+ * @param image (i) - input image
+ * @param W, H      - size of binarized image (in pixels)
+ * @return allocated memory area with converted input image
+ */
+uint8_t *filter8(uint8_t *image, int W, int H){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    uint8_t *ret = MALLOC(uint8_t, W*H);
+    int W0 = (W + 7) / 8; // width in bytes
+    int w = W0-1, h = H-1;
+    {int y = 0;
+    #define IM_UP
+    #include "ec_filter.h"
+    #undef IM_UP
+    }
+    {
+    #include "ec_filter.h"
+    }
+    {
+    int y = h;
+    #define IM_DOWN
+    #include "ec_filter.h"
+    #undef IM_DOWN
+    }
+    return ret;
+}
+
+/**
+ * Make morphological operation of dilation
+ * @param image (i) - input image
+ * @param W, H      - size of image (pixels)
+ * @return allocated memory area with dilation of input image
+ */
+uint8_t *dilation(uint8_t *image, int W, int H){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    int W0 = (W + 7) / 8; // width in bytes
+    int w = W0-1, h = H-1, rest = 7 - (W - w*8);
+    uint8_t lastmask = ~(1<<rest);
+    if(!DIL) morph_init();
+    uint8_t *ret = MALLOC(uint8_t, W0*H);
+    {int y = 0;
+    // top of image, y = 0
+    #define IM_UP
+    #include "dilation.h"
+    #undef IM_UP
+    }
+    {
+    // mid of image, y = 1..h-1
+    #include "dilation.h"
+    }
+    {int y = h;
+    // image bottom, y = h
+    #define IM_DOWN
+    #include "dilation.h"
+    #undef IM_DOWN
+    }
+    return ret;
+}
+
+/**
+ * Make morphological operation of erosion by cross 3x3 pixels
+ * @param image (i) - input image
+ * @param W, H      - size of image (in pixels)
+ * @return allocated memory area with erosion of input image
+ */
+uint8_t *erosion(uint8_t *image, int W, int H){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    if(!ER) morph_init();
+    int W0 = (W + 7) / 8; // width in bytes
+    int w = W0-1, h = H-1, rest = 8 - (W - w*8);
+    uint8_t lastmask = ~(1<<rest);
+    uint8_t *ret = MALLOC(uint8_t, W0*H);
+    DBG("rest=%d, mask:0x%x", rest, lastmask);
+    OMP_FOR()
+    for(int y = 1; y < h; y++){ // reset first & last rows of image
+        uint8_t *iptr = &image[W0*y];
+        uint8_t *optr = &ret[W0*y];
+        register uint8_t p;
+        // x == 0, reset left bit
+        p = ER[*iptr] & *iptr & iptr[-W0] & iptr[W0];
+        if(!(*(++iptr) & 0x80)) p &= 0xfe;
+        *optr++ = p & 0x7f;
+        for(int x = 1; x < w; x++, iptr++, optr++){
+            p = ER[*iptr] & *iptr & iptr[-W0] & iptr[W0];
+            if(!(iptr[-1] & 1)) p &= 0x7f;
+            if(!(iptr[1] & 0x80)) p &= 0xfe;
+            *optr = p;
+        }
+        // x == w, reset right bit
+        p = ER[*iptr] & *iptr & iptr[-W0] & iptr[W0];
+        if(!(iptr[-1] & 1)) p &= 0x7f;
+        *optr = p & lastmask;
+    }
+    return ret;
+}
+
+// Make erosion N times
+uint8_t *erosionN(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < 1 || H < 1) return NULL;
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1){
+        uint8_t *copy = MALLOC(uint8_t, W*H);
+        memcpy(copy, image, W*H);
+        return copy;
+    }
+    uint8_t *cur = image, *next = NULL;
+    for(int i = 0; i < N; ++i){
+        next = erosion(cur, W, H);
+        if(cur != image) FREE(cur);
+        cur = next;
+    }
+    return next;
+}
+// Make dilation N times
+uint8_t *dilationN(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < 1 || H < 1) return NULL;
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1){
+        uint8_t *copy = MALLOC(uint8_t, W*H);
+        memcpy(copy, image, W*H);
+        return copy;
+    }
+    uint8_t *cur = image, *next = NULL;
+    for(int i = 0; i < N; ++i){
+        next = dilation(cur, W, H);
+        if(cur != image) FREE(cur);
+        cur = next;
+    }
+    return next;
+}
+
+// Ntimes opening
+uint8_t *openingN(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
+    uint8_t *er = erosionN(image, W, H, N);
+    uint8_t *op = dilationN(er, W, H, N);
+    FREE(er);
+    return op;
+}
+
+// Ntimes closing
+uint8_t *closingN(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
+    uint8_t *di = dilationN(image, W, H, N);
+    uint8_t *cl = erosionN(di, W, H, N);
+    FREE(di);
+    return cl;
+}
+
+// top hat operation: image - opening(image)
+uint8_t *topHat(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
+    uint8_t *op = openingN(image, W, H, N);
+    int W0 = (W + 7) / 8; // width in bytes
+    int wh = W0 * H;
+    OMP_FOR()
+    for(int i = 0; i < wh; ++i)
+        op[i] = image[i] & (~op[i]);
+    return op;
+}
+
+// bottom hat operation: closing(image) - image
+uint8_t *botHat(uint8_t *image, int W, int H, int N){
+    FNAME();
+    if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
+    uint8_t *op = closingN(image, W, H, N);
+    int W0 = (W + 7) / 8; // width in bytes
+    int wh = W0 * H;
+    OMP_FOR()
+    for(int i = 0; i < wh; ++i)
+        op[i] &= ~image[i];
+    return op;
+}
+
+/*
+ * <=================== MORPHOLOGICAL OPERATIONS ===================
+ */
+
+/*
+ * =================== LOGICAL OPERATIONS ===================>
+ */
+#if 0
+/**
+ * Logical AND of two images
+ * @param im1, im2 (i) - two images
+ * @param W, H         - their size (of course, equal for both images)
+ * @return allocated memory area with   image = (im1 AND im2)
+ */
+uint8_t *imand(uint8_t *im1, uint8_t *im2, int W, int H){
+    uint8_t *ret = MALLOC(uint8_t, W*H);
+    int y;
+    OMP_FOR()
+    for(y = 0; y < H; y++){
+        int x, S = y*W;
+        uint8_t *rptr = &ret[S], *p1 = &im1[S], *p2 = &im2[S];
+        for(x = 0; x < W; x++)
+            *rptr++ = *p1++ & *p2++;
+    }
+    return ret;
+}
+
+/**
+ * Substitute image 2 from image 1: reset to zero all bits of image 1 which set to 1 on image 2
+ * @param im1, im2 (i) - two images
+ * @param W, H         - their size (of course, equal for both images)
+ * @return allocated memory area with    image = (im1 AND (!im2))
+ */
+uint8_t *substim(uint8_t *im1, uint8_t *im2, int W, int H){
+    uint8_t *ret = MALLOC(uint8_t, W*H);
+    int y;
+    OMP_FOR()
+    for(y = 0; y < H; y++){
+        int x, S = y*W;
+        uint8_t *rptr = &ret[S], *p1 = &im1[S], *p2 = &im2[S];
+        for(x = 0; x < W; x++)
+            *rptr++ = *p1++ & (~*p2++);
+    }
+    return ret;
+}
+#endif
+/*
+ * <=================== LOGICAL OPERATIONS ===================
+ */
+
+/*
+ * =================== CONNECTED COMPONENTS LABELING ===================>
+ */
+
+//#define TESTMSGS
+#ifdef TESTMSGS
+#define TEST(...) printf(__VA_ARGS__)
+#else
+#define TEST(...)
+#endif
+
+/**
+ * label 4-connected components on image
+ * (slow algorythm, but easy to parallel)
+ *
+ * @param I (i)    - image ("packed")
+ * @param W,H      - size of the image (W - width in pixels)
+ * @param CC (o)   - connected components boxes
+ * @return an array of labeled components
+ */
+size_t *cclabel4(uint8_t *Img, int W, int H, ConnComps **CC){
+    size_t *assoc;
+    // check table and rename all "oldval" into "newval"
+    inline void remark(size_t newval, size_t oldval){
+        TEST("\tnew = %zd, old=%zd; ", newval, oldval);
+        // find the least values
+        do{newval = assoc[newval];}while(assoc[newval] != newval);
+        do{oldval = assoc[oldval];}while(assoc[oldval] != oldval);
+        TEST("\trealnew = %zd, realold=%zd ", newval, oldval);
+        // now change larger value to smaller
+        if(newval > oldval){
+            assoc[newval] = oldval;
+            TEST("change %zd to %zd\n", newval, oldval);
+        }else{
+            assoc[oldval] = newval;
+            TEST("change %zd to %zd\n", oldval, newval);
+        }
+    }
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    uint8_t *f = filter4(Img, W, H); // remove all non 4-connected pixels
+    DBG("convert to size_t");
+    size_t *labels = bin2ST(f, W, H);
+    FREE(f);
+    DBG("Calculate");
+    size_t Nmax = W*H/4; // max number of 4-connected labels
+    assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array
+    size_t last_assoc_idx = 1; // last index filled in assoc array
+    for(int y = 0; y < H; ++y){
+        bool found = false;
+        size_t *ptr = &labels[y*W];
+        size_t curmark = 0; // mark of pixel to the left
+        for(int x = 0; x < W; ++x, ++ptr){
+            if(!*ptr){found = false; continue;} // empty pixel
+            size_t U = (y) ? ptr[-W] : 0; // upper mark
+            if(found){ // there's a pixel to the left
+                if(U && U != curmark){ // meet old mark -> remark one of them in assoc[]
+                    TEST("(%d, %d): remark %zd --> %zd\n", x, y, U, curmark);
+                    remark(U, curmark);
+                    curmark = U; // change curmark to upper mark (to reduce further checks)
+                }
+            }else{ // new mark -> change curmark
+                found = true;
+                if(U) curmark = U; // current mark will copy upper value
+                else{ // current mark is new value
+                    curmark = last_assoc_idx++;
+                    assoc[curmark] = curmark;
+                    TEST("(%d, %d): new mark=%zd\n", x, y, curmark);
+                }
+            }
+            *ptr = curmark;
+        }
+    }
+    size_t *indexes = MALLOC(size_t, last_assoc_idx); // new indexes
+    size_t cidx = 1;
+    TEST("\n\n\nRebuild indexes\n\n");
+    for(size_t i = 1; i < last_assoc_idx; ++i){
+        TEST("%zd\t%zd ",i,assoc[i]);
+        // search new index
+        register size_t realnew = i, newval = 0;
+        do{
+            realnew = assoc[realnew];
+            TEST("->%zd ", realnew);
+            if(indexes[realnew]){ // find least index
+                newval = indexes[realnew];
+                TEST("real: %zd ", newval);
+                break;
+            }
+        }while(assoc[realnew] != realnew);
+        if(newval){
+            TEST(" ==> %zd\n", newval);
+            indexes[i] = newval;
+            continue;
+        }
+        TEST("new index %zd\n", cidx);
+        // enter next label
+        indexes[i] = cidx++;
+    }
+    // cidx now is amount of detected objects + 1 - size of output array (0th idx is not used)
+    DBG("amount after rebuild: %zd", cidx-1);
+    #ifdef TESTMSGS
+    printf("\n\n\nI\tASS[I]\tIDX[I]\n");
+    for(size_t i = 1; i < last_assoc_idx; ++i)
+        printf("%zd\t%zd\t%zd\n",i,assoc[i],indexes[i]);
+    #endif
+    Box *boxes = MALLOC(Box, cidx);
+    OMP_FOR()
+    for(size_t i = 1; i < cidx; ++i){ // init borders
+        boxes[i].xmin = W;
+        boxes[i].ymin = H;
+    }
+    //nonOMP: 36.5ms
+/* linear:
+    for(int y = 0; y < H; ++y){
+        size_t *lptr = &labels[y*W];
+        for(int x = 0; x < W; ++x, ++lptr){
+            if(!*lptr) continue;
+            register size_t mark = indexes[*lptr];
+            *lptr = mark;
+            Box *b = &boxes[mark];
+            ++b->area;
+            if(b->xmax < x) b->xmax = x;
+            if(b->xmin > x) b->xmin = x;
+            if(b->ymax < y) b->ymax = y;
+            if(b->ymin > y) b->ymin = y;
+        }
+    }*/
+// parallel:
+#pragma omp parallel shared(boxes)
+    {
+        Box *l_boxes = MALLOC(Box, cidx);
+        for(size_t i = 1; i < cidx; ++i){ // init borders
+            l_boxes[i].xmin = W;
+            l_boxes[i].ymin = H;
+        }
+        #pragma omp for nowait
+    for(int y = 0; y < H; ++y){
+        size_t *lptr = &labels[y*W];
+        for(int x = 0; x < W; ++x, ++lptr){
+            if(!*lptr) continue;
+            register size_t mark = indexes[*lptr];
+            *lptr = mark;
+            Box *b = &l_boxes[mark];
+            ++b->area;
+            if(b->xmax < x) b->xmax = x;
+            if(b->xmin > x) b->xmin = x;
+            if(b->ymax < y) b->ymax = y;
+            if(b->ymin > y) b->ymin = y;
+        }
+    }
+    #pragma omp critical
+    for(size_t i = 1; i < cidx; ++i){
+        Box *ob = &boxes[i], *ib = &l_boxes[i];
+        if(ob->xmax < ib->xmax) ob->xmax = ib->xmax;
+        if(ob->xmin > ib->xmin) ob->xmin = ib->xmin;
+        if(ob->ymax < ib->ymax) ob->ymax = ib->ymax;
+        if(ob->ymin > ib->ymin) ob->ymin = ib->ymin;
+        ob->area += ib->area;
+    }
+    FREE(l_boxes);
+    }
+    FREE(assoc);
+    FREE(indexes);
+#ifdef TESTMSGS
+    for(size_t i = 1; i < cidx; ++i){
+        printf("%8zd\t%6d\t(%4d..%4d, %4d..%4d)\t%.2f\n", i, boxes[i].area,
+               boxes[i].xmin, boxes[i].xmax, boxes[i].ymin, boxes[i].ymax,
+               (1.+boxes[i].xmax-boxes[i].xmin)/(1.+boxes[i].ymax-boxes[i].ymin));
+    }printf("\n\n");
+#endif
+    if(CC){
+        *CC = MALLOC(ConnComps, 1);
+        (*CC)->Nobj = cidx; (*CC)->boxes = boxes;
+    }else{
+        FREE(boxes);
+    }
+    return labels;
+}
+
+#if 0
+// label 8-connected components, look cclabel4
+size_t *cclabel8(size_t *labels, int W, int H, size_t *Nobj){
+    if(W < MINWIDTH || H < MINHEIGHT) return NULL;
+    #define LABEL_8
+    #include "cclabling.h"
+    #undef LABEL_8
+    return labels;
+}
+#endif
+
+
+#if 0
+/**
+ * Make connection-component labeling
+ * output image vould have uint16_t data
+ * @param img (i)       - input image
+ * @param threshold (i) - threshold level in value of dynamic range (0,1)
+ * @param Nobj (o)      - amount of object found (or NULL if not needed)
+ */
+IMAGE *cclabel4(IMAGE *img, double threshold, size_t *Nobj){
+    /*if(N != 4 || N != 8){
+        ERRX(_("Can work only for 4- or 8-connected components"));
+    }*/
+    double thrval;
+    uint16_t *binary = binarize(img, threshold, &thrval);
+    if(!binary) return NULL;
+    int W_0;
+    uint8_t *Ima = u16tochar(binary, img->width, img->height, &W_0);
+    FREE(binary);
+    size_t N;
+    uint16_t *dat = _cclabel4(Ima, img->width, img->height, W_0, &N);
+    if(Nobj) *Nobj = N;
+    FREE(Ima);
+    IMAGE *ret = buildFITSfromdat(img->height, img->width, SHORT_IMG, (uint8_t*)dat);
+    FREE(dat);
+    char buf[80];
+    snprintf(buf, 80, "COMMENT found %zd 4-connected components, threshold value %g",
+        N, (double)thrval);
+    list_add_record(&ret->keylist, buf);
+    snprintf(buf, 80, "COMMENT    (%g%% fromdata range%s)", fabs(threshold)*100.,
+        (threshold < 0.) ? ", inverted" : "");
+    list_add_record(&ret->keylist, buf);
+    return ret;
+}
+
+IMAGE *cclabel8(IMAGE *img, double threshold, size_t *Nobj){
+    double thrval;
+    uint16_t *binary = binarize(img, threshold, &thrval);
+    if(!binary) return NULL;
+    size_t N;
+    _cclabel8(binary, img->width, img->height, &N);
+    if(Nobj) *Nobj = N;
+    IMAGE *ret = buildFITSfromdat(img->height, img->width, SHORT_IMG, (uint8_t*)binary);
+    FREE(binary);
+    char buf[80];
+    snprintf(buf, 80, "COMMENT found %zd 8-connected components, threshold value %g",
+        N, (double)thrval);
+    list_add_record(&ret->keylist, buf);
+    snprintf(buf, 80, "COMMENT    (%g%% fromdata range%s)", fabs(threshold)*100.,
+        (threshold < 0.) ? ", inverted" : "");
+    list_add_record(&ret->keylist, buf);
+    return ret;
+}
+#endif
+/*
+ * <=================== CONNECTED COMPONENTS LABELING ===================
+ */
+
+
+/*
+ * <=================== template ===================>
+ */

LocCorr/binmorph.h

@@ -0,0 +1,85 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#ifndef BINMORPH_H__
+#define BINMORPH_H__
+
+#include <stdbool.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include "fits.h"
+
+// minimal image size for morph operations
+#define MINWIDTH    (9)
+#define MINHEIGHT   (3)
+
+// simple box with given borders
+typedef struct{
+    uint16_t xmin;
+    uint16_t xmax;
+    uint16_t ymin;
+    uint16_t ymax;
+    uint32_t area; // total amount of object pixels inside the box
+} Box;
+
+typedef struct{
+    size_t Nobj;
+    Box *boxes;
+} ConnComps;
+
+
+// morphological operations:
+uint8_t *dilation(uint8_t *image, int W, int H);
+uint8_t *dilationN(uint8_t *image, int W, int H, int N);
+uint8_t *erosion(uint8_t *image, int W, int H);
+uint8_t *erosionN(uint8_t *image, int W, int H, int N);
+uint8_t *openingN(uint8_t *image, int W, int H, int N);
+uint8_t *closingN(uint8_t *image, int W, int H, int N);
+uint8_t *topHat(uint8_t *image, int W, int H, int N);
+uint8_t *botHat(uint8_t *image, int W, int H, int N);
+
+// clear non 4-connected pixels
+uint8_t *filter4(uint8_t *image, int W, int H);
+// clear single pixels
+uint8_t *filter8(uint8_t *image, int W, int H);
+
+size_t *cclabel4(uint8_t *Img, int W, int H, ConnComps **CC);
+size_t *cclabel8(uint8_t *Img, int W, int H, size_t *Nobj);
+
+#if 0
+// logical operations
+uint8_t *imand(uint8_t *im1, uint8_t *im2, int W, int H);
+uint8_t *substim(uint8_t *im1, uint8_t *im2, int W, int H);
+/*
+// conncomp
+// this is a box structure containing one object; data is aligned by original image bytes!
+typedef struct {
+    uint8_t *data; // pattern data in "packed" format
+    int x,   // x coordinate of LU-pixel of box in "unpacked" image (by pixels)
+    y,       // y -//-
+    x_0;     // x coordinate in "packed" image (morph operations should work with it)
+    size_t N;// number of component, starting from 1
+} CCbox;
+*/
+
+IMAGE *cclabel4(IMAGE *I, double threshold, size_t *Nobj);
+IMAGE *cclabel8(IMAGE *I, double threshold, size_t *Nobj);
+#endif
+
+#endif // BINMORPH_H__

LocCorr/cclabling.h

@@ -0,0 +1,136 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+//#define TESTMSGS
+
+#ifdef TESTMSGS
+#define TEST(...) printf(__VA_ARGS__)
+#else
+#define TEST(...)
+#endif
+
+size_t Nmax = W*H/4; // max number of 4-connected labels
+size_t *assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array
+size_t last_assoc_idx = 1; // last index filled in assoc array
+// was: 35
+// check table and rename all "oldval" into "newval"
+inline void remark(size_t newval, size_t oldval){
+    TEST("\tnew = %zd, old=%zd; ", newval, oldval);
+    // find the least values
+    do{newval = assoc[newval];}while(assoc[newval] != newval);
+    do{oldval = assoc[oldval];}while(assoc[oldval] != oldval);
+    TEST("\trealnew = %zd, realold=%zd ", newval, oldval);
+    // now change larger value to smaller
+    if(newval > oldval){
+        assoc[newval] = oldval;
+        TEST("change %zd to %zd\n", newval, oldval);
+    }else{
+        assoc[oldval] = newval;
+        TEST("change %zd to %zd\n", oldval, newval);
+    }
+}
+
+for(int y = 0; y < H; ++y){
+    bool found = false;
+    size_t *ptr = &labels[y*W];
+    size_t curmark = 0; // mark of pixel to the left
+    for(int x = 0; x < W; ++x, ++ptr){
+        if(!*ptr){found = false; continue;} // empty pixel
+        size_t U = (y) ? ptr[-W] : 0; // upper mark
+        if(found){ // there's a pixel to the left
+            if(U && U != curmark){ // meet old mark -> remark one of them in assoc[]
+                TEST("(%d, %d): remark %zd --> %zd\n", x, y, U, curmark);
+                remark(U, curmark);
+                curmark = U; // change curmark to upper mark (to reduce further checks)
+            }
+        }else{ // new mark -> change curmark
+            found = true;
+            if(U) curmark = U; // current mark will copy upper value
+            else{ // current mark is new value
+                curmark = last_assoc_idx++;
+                assoc[curmark] = curmark;
+                TEST("(%d, %d): new mark=%zd\n", x, y, curmark);
+            }
+        }
+        *ptr = curmark;
+    }
+}
+size_t *indexes = MALLOC(size_t, last_assoc_idx); // new indexes
+size_t cidx = 1;
+TEST("\n\n\nRebuild indexes\n\n");
+for(size_t i = 1; i < last_assoc_idx; ++i){
+    TEST("%zd\t%zd ",i,assoc[i]);
+    // search new index
+    register size_t realnew = i, newval = 0;
+    do{
+        realnew = assoc[realnew];
+        TEST("->%zd ", realnew);
+        if(indexes[realnew]){ // find least index
+            newval = indexes[realnew];
+            TEST("real: %zd ", newval);
+            break;
+        }
+    }while(assoc[realnew] != realnew);
+    if(newval){
+        TEST(" ==> %zd\n", newval);
+        indexes[i] = newval;
+        continue;
+    }
+    TEST("new index %zd\n", cidx);
+    // enter next label
+    indexes[i] = cidx++;
+}
+/*
+// rebuild indexes
+for(size_t i = 1; i < last_assoc_idx; ++i){
+    printf("%zd\t%zd ",i,assoc[i]);
+    if(i == assoc[i]){
+        if(i == cidx){
+            printf("- keep\n");
+        }else{
+            printf("- change to %zd\n", cidx);
+            size_t _2change = assoc[i];
+            assoc[i] = cidx;
+            for(size_t j = i+1; j < last_assoc_idx; ++j){
+                if(assoc[j] == _2change){
+                    printf("\t%zd\t%zd -> %zd\n", j, assoc[j], cidx);
+                    assoc[j] = cidx;
+                }
+            }
+        }
+        ++cidx;
+    }else printf("\n");
+}*/
+
+--cidx; // cidx now is amount of detected objects
+DBG("amount after rebuild: %zd", cidx);
+
+#ifdef TESTMSGS
+printf("\n\n\nI\tASS[I]\tIDX[I]\n");
+for(size_t i = 1; i < last_assoc_idx; ++i)
+    printf("%zd\t%zd\t%zd\n",i,assoc[i],indexes[i]);
+#endif
+
+int wh = W*H;
+OMP_FOR()
+for(int i = 0; i < wh; ++i) labels[i] = indexes[labels[i]];
+
+if(Nobj) *Nobj = cidx;
+
+FREE(assoc);
+FREE(indexes);

LocCorr/cclabling_ori.h

@@ -0,0 +1,128 @@
+/*
+ * cclabling_1.h - inner part of functions cclabel4 and cclabel8
+ *
+ * Copyright 2015 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+//double t0 = dtime();
+
+    size_t N = 0, // current label
+    Nmax = W*H/4; // max number of labels
+    int w = W - 1, h = H - 1;
+    int y;
+    size_t *assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array
+    size_t last_assoc_idx = 0; // last index filled in assoc array
+    size_t currentnum = 0; // current pixel number
+    inline void remark(size_t old, size_t new){ // remark in assoc[] pixel with value old to assoc[new] or vice versa
+        size_t New = assoc[new], Old = assoc[old];
+        if(Old == New){
+            return;
+        }
+        // now we must check Old: if  Old<New we should swap them
+        if(Old < New){
+            register size_t _tmp_ = Old; Old = New; New = _tmp_; // swap values
+            _tmp_ = old; old = new; new = _tmp_;
+        }
+        // decrement counters for current value (because we make merging)
+        --currentnum;
+        for(size_t i = 1; i < last_assoc_idx; ++i){
+            size_t m = assoc[i];
+            if(m < Old) continue; // lower values
+            if(m == Old){ // change all old markers to new
+                assoc[i] = New;
+            }else{ // decrement all higher values
+                --assoc[i];
+            }
+        }
+    }
+    size_t *ptr = labels;
+    for(y = 0; y < H; y++){
+        bool found = false;
+        size_t curmark = 0; // mark of pixel to the left
+        for(int x = 0; x < W; x++, ptr++){
+            size_t curval = *ptr;
+            if(!curval){ found = false; continue;}
+            size_t *U = (y) ? &ptr[-W] : NULL;
+            size_t upmark = 0; // mark from line above
+            if(!found){ // new pixel, check neighbours above
+                found = true;
+                // now check neighbours in upper string:
+                if(U){
+                    #ifdef LABEL_8
+                    if(x && U[-1]){ // there is something in upper left corner -> use its value
+                        upmark = U[-1];
+                    }else // check point above only if there's nothing in left up
+                    #endif
+                        if(U[0]) upmark = U[0];
+                    #ifdef LABEL_8
+                    if(x < w && U[1]){ // there's something in upper right
+                        if(upmark){ // to the left of it was pixels
+                            remark(U[1], upmark);
+                        }else
+                            upmark = U[1];
+                    }
+                    #endif
+                }
+                if(!upmark){ // there's nothing above - set current pixel to incremented counter
+                    #ifdef LABEL_8  // check, whether pixel is not single
+                    size_t *D = (y < h) ? &ptr[W] : NULL;
+                    if(  !(x && ((D && D[-1]) /*|| ptr[-1]*/))   // no left neighbours
+                      && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours
+                      && !(D && D[0])){ // no neighbour down
+                        *ptr = 0; // erase this hermit!
+                        continue;
+                    }
+                    #else
+                    // no neighbour down & neighbour to the right -> hermit
+                    if((y < h && ptr[W] == 0) && (x < w && ptr[1] == 0)){
+                        *ptr = 0; // erase this hermit!
+                        continue;
+                    }
+                    #endif
+                    upmark = ++N;
+                    assoc[upmark] = ++currentnum; // refresh "assoc"
+                    last_assoc_idx = upmark + 1;
+                }
+                *ptr = upmark;
+                curmark = upmark;
+            }else{ // there was something to the left -> we must chek only U[1]
+                if(U){
+                    if(x < w && U[1]){ // there's something in upper right
+                        remark(U[1], curmark);
+                    }
+                }
+                *ptr = curmark;
+            }
+        }
+    }
+#if 0
+    // Step 2: rename markers
+    // first correct complex assotiations in assoc
+    OMP_FOR()
+    for(y = 0; y < H; y++){
+        size_t *ptr = &labels[y*W];
+        for(int x = 0; x < W; x++, ptr++){
+            size_t p = *ptr;
+            if(!p){continue;}
+            *ptr = assoc[p];
+        }
+    }
+#endif
+    FREE(assoc);
+    if(Nobj) *Nobj = currentnum;
+//printf("%6.4f\t%zd\n", dtime()-t0, currentnum);

LocCorr/cmdlnopts.c

@@ -0,0 +1,107 @@
+/*                                                                                                  geany_encoding=koi8-r
+ * cmdlnopts.c - the only function that parse cmdln args and returns glob parameters
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#include <assert.h>
+#include <stdio.h>
+#include <string.h>
+#include <strings.h>
+#include <math.h>
+#include <usefull_macros.h>
+
+#include "cmdlnopts.h"
+
+/*
+ * here are global parameters initialisation
+ */
+static int help;
+glob_pars *GP = NULL;
+
+// default PID filename:
+#define DEFAULT_PIDFILE "/tmp/loccorr.pid"
+
+//            DEFAULTS
+// default global parameters
+static glob_pars G = {
+    .pidfile = DEFAULT_PIDFILE,
+    .throwpart = 0.5,
+    .medradius = 1.,
+    .ndilations = 2,
+    .nerosions = 2,
+    .minarea = 5,
+    .intensthres = 0.01,
+    .maxexp = 500.,
+    .minexp = 0.001,
+    .xtarget = -1.,
+    .ytarget = -1.
+};
+
+/*
+ * Define command line options by filling structure:
+ *  name        has_arg     flag    val     type        argptr              help
+*/
+static myoption cmdlnopts[] = {
+// common options
+    {"help",    NO_ARGS,    NULL,     0,    arg_int,    APTR(&help),        _("show this help")},
+    {"logfile", NEED_ARG,   NULL,   'l',    arg_string, APTR(&G.logfile),   _("file to save logs (default: none)")},
+    {"pidfile", NEED_ARG,   NULL,   'P',    arg_string, APTR(&G.pidfile),   _("pidfile (default: " DEFAULT_PIDFILE ")")},
+    {"verbose", NO_ARGS,    NULL,   'v',    arg_none,   APTR(&G.verb),      _("increase verbosity level of log file (each -v increased by 1)")},
+    {"input",   NEED_ARG,   NULL,   'i',    arg_string, APTR(&G.inputname), _("file or directory name for monitoring (or grasshopper for capturing)")},
+    {"blackp",  NEED_ARG,   NULL,   'b',    arg_double, APTR(&G.throwpart), _("fraction of black pixels to throw away when make histogram eq")},
+    {"radius",  NEED_ARG,   NULL,   'r',    arg_int,    APTR(&G.medradius), _("radius of median filter (r=1 -> 3x3, r=2 -> 5x5 etc.)")},
+    {"equalize", NO_ARGS,   NULL,   'e',    arg_int,    APTR(&G.equalize),  _("make historam equalization of saved jpeg")},
+    {"ndilat",  NEED_ARG,   NULL,   'D',    arg_int,    APTR(&G.ndilations),_("amount of erosions after thresholding (default: 2)")},
+    {"neros",   NEED_ARG,   NULL,   'E',    arg_int,    APTR(&G.nerosions), _("amount of dilations after erosions (default: 2)")},
+    {"minarea", NEED_ARG,   NULL,   'A',    arg_int,    APTR(&G.minarea),   _("minimal object pixels amount (default: 5)")},
+    {"intthres",NEED_ARG,   NULL,   'T',    arg_double, APTR(&G.intensthres),_("threshold by total object intensity when sorting = |I1-I2|/(I1+I2), default: 0.01")},
+    {"xoff",    NEED_ARG,   NULL,   'x',    arg_int,    APTR(&G.xoff),      _("X offset at grabbed image")},
+    {"yoff",    NEED_ARG,   NULL,   'y',    arg_int,    APTR(&G.yoff),      _("Y offset at grabbed image")},
+    {"width",   NEED_ARG,   NULL,   'w',    arg_int,    APTR(&G.width),     _("grabbed subimage width")},
+    {"height",  NEED_ARG,   NULL,   'h',    arg_int,    APTR(&G.height),    _("grabbed subimage height")},
+    {"maxexp",  NEED_ARG,   NULL,   0,      arg_double, APTR(&G.maxexp),    _("maximal exposition time (ms), default: 500")},
+    {"minexp",  NEED_ARG,   NULL,   0,      arg_double, APTR(&G.minexp),    _("minimal exposition time (ms), default: 0.001")},
+    {"xtarget", NEED_ARG,   NULL,   'X',    arg_double, APTR(&G.xtarget),   _("target point X coordinate")},
+    {"ytarget", NEED_ARG,   NULL,   'Y',    arg_double, APTR(&G.ytarget),   _("target point Y coordinate")},
+    {"logXY",   NEED_ARG,   NULL,   'L',    arg_string, APTR(&G.logXYname), _("file to log XY coordinates of selected star")},
+   end_option
+};
+
+/**
+ * Parse command line options and return dynamically allocated structure
+ *      to global parameters
+ * @param argc - copy of argc from main
+ * @param argv - copy of argv from main
+ * @return allocated structure with global parameters
+ */
+glob_pars *parse_args(int argc, char **argv){
+    size_t hlen = 1024;
+    char helpstring[1024], *hptr = helpstring;
+    snprintf(hptr, hlen, "Usage: %%s [args]\n\n\tWhere args are:\n");
+    // format of help: "Usage: progname [args]\n"
+    change_helpstring(helpstring);
+    // parse arguments
+    parseargs(&argc, &argv, cmdlnopts);
+    if(help) showhelp(-1, cmdlnopts);
+    if(argc > 0){
+        WARNX("Extra parameters!");
+        showhelp(-1, cmdlnopts);
+    }
+    return &G;
+}
+

LocCorr/cmdlnopts.h

@@ -0,0 +1,53 @@
+/*                                                                                                  geany_encoding=koi8-r
+ * cmdlnopts.h - comand line options for parceargs
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef CMDLNOPTS_H__
+#define CMDLNOPTS_H__
+
+/*
+ * here are some typedef's for global data
+ */
+typedef struct{
+    char *pidfile;          // name of PID file
+    char *logfile;          // logging to this file
+    char *inputname;        // input for monitor file or directory name
+    char *logXYname;        // file to log XY coordinates of first point
+    double throwpart;       // fraction of black pixels to throw away when make histogram eq
+    int equalize;           // make historam equalization of saved jpeg
+    int medradius;          // radius of median filter (r=1 -> 3x3, r=2 -> 5x5 etc.)
+    int verb;               // logfile verbosity level
+    int ndilations;         // amount of erosions (default: 2)
+    int nerosions;          // amount of dilations (default: 2)
+    int minarea;            // minimal object pixels amount (default: 5)
+    double intensthres;     // threshold by total object intensity when sorting = |I1-I2|/(I1+I2), default: 0.01
+    double maxexp;          // max exposition time (ms)
+    double minexp;          // min exposition time (ms)
+    int xoff; int yoff;     // offset by X and Y axes
+    int width; int height;  // target width and height of image
+    double xtarget; double ytarget;// target point coordinates
+} glob_pars;
+
+extern glob_pars *GP;  // for GP->pidfile need in `signals`
+
+glob_pars *parse_args(int argc, char **argv);
+
+#endif // CMDLNOPTS_H__

LocCorr/dilation.h

@@ -0,0 +1,73 @@
+/*
+ * dilation.h - inner part of function `dilation`
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+/*
+ *  HERE'S NO ANY "FILE-GUARDS" BECAUSE FILE IS MULTIPLY INCLUDED!
+ *  I use this because don't want to dive into depths of BOOST
+ *
+ *  multi-including with different defines before - is a simplest way to
+ * modify simple code for avoiding extra if-then-else
+ */
+
+
+#if ! defined IM_UP && ! defined IM_DOWN // image without upper and lower borders
+OMP_FOR()
+for(int y = 1; y < h; y++)
+#endif
+{
+    uint8_t *iptr = &image[W0*y];
+    uint8_t *optr = &ret[W0*y];
+    uint8_t p = DIL[*iptr]
+    #ifndef IM_UP
+        | iptr[-W0]
+    #endif
+    #ifndef IM_DOWN
+        | iptr[W0]
+    #endif
+        ;
+    if(iptr[1] & 0x80) p |= 1;
+    *optr = p;
+    optr++; iptr++;
+    for(int x = 1; x < w; x++, iptr++, optr++){
+        p = DIL[*iptr]
+    #ifndef IM_UP
+        | iptr[-W0]
+    #endif
+    #ifndef IM_DOWN
+        | iptr[W0]
+    #endif
+        ;
+        if(iptr[1] & 0x80) p |= 1;
+        if(iptr[-1] & 1) p |= 0x80;
+        *optr = p;
+    }
+    p = DIL[*iptr]
+    #ifndef IM_UP
+        | iptr[-W0]
+    #endif
+    #ifndef IM_DOWN
+        | iptr[W0]
+    #endif
+    ;
+    if(iptr[-1] & 1) p |= 0x80;
+    *optr = p & lastmask; // clear outern pixel
+}
+

LocCorr/draw.c

@@ -0,0 +1,105 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// simplest interface to draw lines & ellipsis
+#include "draw.h"
+#include "fits.h"
+
+#include <usefull_macros.h>
+
+// base colors:
+const uint8_t
+    C_R[] = {255, 0, 0},
+    C_G[] = {0, 255, 0},
+    C_B[] = {0, 0, 255},
+    C_K[] = {0, 0, 0},
+    C_W[] = {255,255,255};
+
+void Img3_free(Img3 **im){
+    if(!im || !*im) return;
+    FREE((*im)->data);
+    FREE(*im);
+}
+
+void Pattern_free(Pattern **p){
+    if(!p || !*p) return;
+    FREE((*p)->data);
+    FREE(*p);
+}
+
+// make a single-channel (opaque) mask for cross; allocated here!!!
+Pattern *Pattern_cross(int h, int w){
+    int s = h*w, hmid = h/2, wmid = w/2;
+    uint8_t *data = MALLOC(uint8_t, s);
+    uint8_t *ptr = &data[wmid];
+    for(int y = 0; y < h; ++y, ptr += w) *ptr = 255;
+    ptr = &data[hmid*w];
+    for(int x = 0; x < w; ++x, ++ptr) *ptr = 255;
+    Pattern *p = MALLOC(Pattern, 1);
+    p->data = data;
+    p->h = h; p->w = w;
+    return p;
+}
+
+/**
+ * @brief draw3_pattern - draw pattern @ 3-channel image
+ * @param img (io)    - image
+ * @param p (i)       - the pattern
+ * @param xc, yc      - coordinates of pattern center @ image
+ * @param colr        - color to draw pattern (when opaque == 255)
+ */
+void Pattern_draw3(Img3 *img, Pattern *p, int xc, int yc, const uint8_t colr[]){
+    int xul = xc - p->w/2, yul = yc - p->h/2;
+    int xdr = xul+p->w-1, ydr = yul+p->h-1;
+    int R = img->w, D = img->h; // right and down border coordinates + 1
+    if(ydr < 0 || xdr < 0 || xul > R-1 || yul > D-1) return; // box outside of image
+
+    int oxlow, oxhigh, oylow, oyhigh; // output limit coordinates
+    int ixlow, iylow; // intput limit coordinates
+    if(xul < 0){
+        oxlow = 0; ixlow = -xul;
+    }else{
+        oxlow = xul; ixlow = 0;
+    }
+    if(yul < 0){
+        oylow = 0; iylow = -yul;
+    }else{
+        oylow = yul; iylow = 0;
+    }
+    if(xdr < R){
+        oxhigh = xdr;
+    }else{
+        oxhigh = R;
+    }
+    if(ydr < D){
+        oyhigh = ydr;
+    }else{
+        oyhigh = D;
+    }
+    OMP_FOR()
+    for(int y = oylow; y < oyhigh; ++y){
+        uint8_t *in = &p->data[(iylow+y-oylow)*p->w + ixlow]; // opaque component
+        uint8_t *out = &img->data[(y*img->w + oxlow)*3]; // 3-colours
+        for(int x = oxlow; x < oxhigh; ++x, ++in, out += 3){
+            float opaque = *in/255.;
+            for(int c = 0; c < 3; ++c){
+                out[c] = (uint8_t)(colr[c] * opaque + out[c]*(1.-opaque));
+            }
+        }
+    }
+}

LocCorr/draw.h

@@ -0,0 +1,56 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+#ifndef DRAW_H__
+#define DRAW_H__
+
+#include <stdint.h>
+
+// 3-channel image for saving into jpg/png
+typedef struct{
+    uint8_t *data;  // image data
+    int w;          // width
+    int h;          // height
+} Img3;
+
+/*
+// box (coordinates from upper left corner)
+typedef struct{
+    int xul;    // coordinates of upper left box corner
+    int yul;
+    int w;      // width and height
+    int h;
+    uint8_t pattern[3]; // pattern for figure filling
+} BBox;
+*/
+
+// opaque pattern for drawing
+typedef struct{
+    uint8_t *data;
+    int w;
+    int h;
+} Pattern;
+
+// base colors
+extern const uint8_t C_R[], C_G[], C_B[], C_K[], C_W[];
+
+void Pattern_free(Pattern **p);
+void Pattern_draw3(Img3 *img, Pattern *p, int xc, int yc, const uint8_t colr[]);
+Pattern *Pattern_cross(int h, int w);
+
+#endif // DRAW_H__

LocCorr/ec_filter.h

@@ -0,0 +1,96 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// inner part of function `filter8`
+
+#if ! defined IM_UP && ! defined IM_DOWN
+OMP_FOR()
+for(int y = 1; y < h; y++)
+#endif
+{
+    uint8_t *iptr = &image[W0*y];
+    uint8_t *optr = &ret[W0*y];
+    // x=0
+    register uint8_t inp = *iptr;
+    register uint8_t p = (inp << 1) | (inp >> 1)
+        #ifndef IM_UP
+            | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1)
+        #endif
+        #ifndef IM_DOWN
+            | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1)
+        #endif
+    ;
+    if((iptr[1]
+        #ifndef IM_UP
+            | iptr[-W0+1]
+        #endif
+        #ifndef IM_DOWN
+           | iptr[W0+1]
+        #endif
+       )& 0x80) p |= 1;
+    *optr++ = inp & p;
+    ++iptr;
+    for(int x = 1; x < w; ++x){
+        inp = *iptr;
+        p = (inp << 1) | (inp >> 1)
+            #ifndef IM_UP
+                | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1)
+            #endif
+            #ifndef IM_DOWN
+                | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1)
+            #endif
+        ;
+        if((iptr[1]
+            #ifndef IM_UP
+                | iptr[-W0+1]
+            #endif
+            #ifndef IM_DOWN
+               | iptr[W0+1]
+            #endif
+           )& 0x80) p |= 1;
+        if((iptr[-1]
+            #ifndef IM_UP
+                | iptr[-W0-1]
+            #endif
+            #ifndef IM_DOWN
+               | iptr[W0-1]
+            #endif
+           )& 1) p |= 0x80;
+        *optr++ = inp & p;
+        ++iptr;
+    }
+    // x = W0-1
+    inp = *iptr;
+    p = (inp << 1) | (inp >> 1)
+        #ifndef IM_UP
+            | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1)
+        #endif
+        #ifndef IM_DOWN
+            | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1)
+        #endif
+    ;
+    if((iptr[-1]
+        #ifndef IM_UP
+            | iptr[-W0-1]
+        #endif
+        #ifndef IM_DOWN
+           | iptr[W0-1]
+        #endif
+       )& 1) p |= 0x80;
+    *optr = inp & p;
+}

LocCorr/fc_filter.h

@@ -0,0 +1,68 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// inner part of function `filter4`
+
+#if ! defined IM_UP && ! defined IM_DOWN
+OMP_FOR()
+for(int y = 1; y < h; y++)
+#endif
+{
+    uint8_t *iptr = &image[W0*y];
+    uint8_t *optr = &ret[W0*y];
+    // x=0
+    register uint8_t inp = *iptr;
+    register uint8_t p = (inp << 1) | (inp >> 1)
+        #ifndef IM_UP
+            | iptr[-W0]
+        #endif
+        #ifndef IM_DOWN
+            | iptr[W0]
+        #endif
+    ;
+    if(iptr[1] & 0x80) p |= 1;
+    *optr++ = inp & p;
+    ++iptr;
+    for(int x = 1; x < w; ++x){
+        inp = *iptr;
+        p = (inp << 1) | (inp >> 1)
+            #ifndef IM_UP
+                | iptr[-W0]
+            #endif
+            #ifndef IM_DOWN
+                | iptr[W0]
+            #endif
+        ;
+        if(iptr[1] & 0x80) p |= 1;
+        if(iptr[-1] & 1) p |= 0x80;
+        *optr++ = inp & p;
+        ++iptr;
+    }
+    // x = W0-1
+    inp = *iptr;
+    p = (inp << 1) | (inp >> 1)
+        #ifndef IM_UP
+            | iptr[-W0]
+        #endif
+        #ifndef IM_DOWN
+            | iptr[W0]
+        #endif
+    ;
+    if(iptr[-1] & 1) p |= 0x80;
+    *optr = inp & p;
+}

LocCorr/fits.c

@@ -0,0 +1,227 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#include <stdint.h>
+#include <string.h>
+#include <stdio.h>
+#include <usefull_macros.h>
+#include <math.h>
+
+#include "fits.h"
+
+static int fitsstatus = 0;
+
+/*
+ * Macros for error processing when working with cfitsio functions
+ */
+#define TRYFITS(f, ...)							\
+do{ fitsstatus = 0;								\
+    f(__VA_ARGS__, &fitsstatus);				\
+    if(fitsstatus){								\
+        fits_report_error(stderr, fitsstatus);	\
+        return FALSE;}							\
+}while(0)
+#define FITSFUN(f, ...)							\
+do{ fitsstatus = 0;								\
+    int ret = f(__VA_ARGS__, &fitsstatus);		\
+    if(ret || fitsstatus)						\
+        fits_report_error(stderr, fitsstatus);	\
+}while(0)
+#define WRITEKEY(...)							\
+do{ fitsstatus = 0;								\
+    fits_write_key(__VA_ARGS__, &fitsstatus);	\
+    if(status) fits_report_error(stderr, status);\
+}while(0)
+
+void Image_free(Image **img){
+    size_t i, sz = (*img)->keynum;
+    char **list = (*img)->keylist;
+    for(i = 0; i < sz; ++i) FREE(list[i]);
+    FREE((*img)->keylist);
+    FREE((*img)->data);
+    FREE(*img);
+}
+
+bool FITS_read(const char *filename, Image **fits){
+    FNAME();
+    bool ret = TRUE;
+    fitsfile *fp;
+    int i, j, hdunum, hdutype, nkeys, keypos;
+    int naxis;
+    long naxes[2];
+    char card[FLEN_CARD];
+    Image *img = MALLOC(Image, 1);
+
+    TRYFITS(fits_open_file, &fp, filename, READONLY);
+    FITSFUN(fits_get_num_hdus, fp, &hdunum);
+    DBG("Got %d HDUs", hdunum);
+    if(hdunum < 1){
+        WARNX(_("Can't read HDU"));
+        ret = FALSE;
+        goto returning;
+    }
+    // get image dimensions
+    TRYFITS(fits_get_img_param, fp, 2, &img->dtype, &naxis, naxes);
+    DBG("Image have %d axis", naxis);
+    if(naxis > 2){
+        WARNX(_("Images with > 2 dimensions are not supported"));
+        ret = FALSE;
+        goto returning;
+    }
+    img->width = naxes[0];
+    img->height = naxes[1];
+    DBG("got image %ldx%ld pix, bitpix=%d", naxes[0], naxes[1], img->dtype);
+    // loop through all HDUs
+    for(i = 1; !(fits_movabs_hdu(fp, i, &hdutype, &fitsstatus)); ++i){
+        TRYFITS(fits_get_hdrpos, fp, &nkeys, &keypos);
+        int oldnkeys = img->keynum;
+        img->keynum += nkeys;
+        if(!(img->keylist = realloc(img->keylist, sizeof(char*) * img->keynum))){
+            ERR(_("Can't realloc"));
+        }
+        char **currec = &(img->keylist[oldnkeys]);
+        DBG("HDU # %d of %d keys", i, nkeys);
+        for(j = 1; j <= nkeys; ++j){
+            FITSFUN(fits_read_record, fp, j, card);
+            if(!fitsstatus){
+                *currec = MALLOC(char, FLEN_CARD);
+                memcpy(*currec, card, FLEN_CARD);
+                DBG("key %d: %s", oldnkeys + j, *currec);
+                ++currec;
+            }
+        }
+    }
+    if(fitsstatus == END_OF_FILE){
+        fitsstatus = 0;
+    }else{
+        fits_report_error(stderr, fitsstatus);
+        ret = FALSE;
+        goto returning;
+    }
+    size_t sz = naxes[0] * naxes[1];
+    img->data = MALLOC(Imtype, sz);
+    int stat = 0;
+    TRYFITS(fits_read_img, fp, FITSDATATYPE, 1, sz, NULL, img->data, &stat);
+    Imtype *d = img->data, min = *d, max = *d;
+    for(size_t x = 0; x < sz; ++x){
+        if(d[x] > max) max = d[x];
+        else if(d[x] < min) min = d[x];
+    }
+    img->maxval = max;
+    img->minval = min;
+    DBG("FITS stat: min=%g, max=%g", min, max);
+    if(stat) WARNX(_("Found %d pixels with undefined value"), stat);
+    DBG("ready");
+
+returning:
+    FITSFUN(fits_close_file, fp);
+    if(!ret){
+        Image_free(&img);
+    }
+    if(fits) *fits = img;
+    return ret;
+}
+
+bool FITS_write(const char *filename, const Image *fits){
+    int w = fits->width, h = fits->height;
+    long naxes[2] = {w, h};
+    size_t sz = w * h, keys = fits->keynum;
+    fitsfile *fp;
+
+    TRYFITS(fits_create_file, &fp, filename);
+    TRYFITS(fits_create_img, fp, fits->dtype, 2, naxes);
+    Imtype *outp = fits->data;
+    bool need2free = FALSE;
+    if(fits->dtype > 0){ // convert floating data into integer
+        Imtype maxval;
+        Imtype minval;
+        switch(fits->dtype){
+            case SBYTE_IMG: // there's a bug in cfitsio, it can't save float->sbyte
+                maxval = (Imtype)INT8_MAX;
+                minval = (Imtype)INT8_MIN;
+            break;
+            case SHORT_IMG:
+                maxval = (Imtype)INT16_MAX;
+                minval = (Imtype)INT16_MIN;
+            break;
+            case USHORT_IMG:
+                maxval = (Imtype)UINT16_MAX;
+                minval = (Imtype)0;
+            break;
+            case LONG_IMG:
+                maxval = (Imtype)INT32_MAX;
+                minval = (Imtype)INT32_MIN;
+            break;
+            case ULONG_IMG:
+                maxval = (Imtype)UINT32_MAX;
+                minval = (Imtype)0;
+            break;
+            case ULONGLONG_IMG:
+                maxval = (Imtype)UINT64_MAX;
+                minval = (Imtype)0;
+            break;
+            case LONGLONG_IMG:
+                maxval = (Imtype)INT64_MAX;
+                minval = (Imtype)INT64_MIN;
+            break;
+            case BYTE_IMG:
+            default:  // byte
+                maxval = (Imtype)UINT8_MAX;
+                minval = (Imtype)0;
+        }
+        DBG("maxval = %g, minval = %g", maxval, minval);
+        int w = fits->width, h = fits->height;
+        Imtype min = fits->minval, max = fits->maxval, W = (maxval - minval)/(max - min);
+        outp = MALLOC(Imtype, w*h);
+        OMP_FOR()
+        for(int y = 0; y < h; ++y){
+            Imtype *o = &outp[y*w], *i = &fits->data[y*w];
+            for(int x = 0; x < w; ++x, ++o, ++i){
+                *o = W*((*i) - min) + minval;
+                //if(*o < minval || *o > maxval) red("o: %g\n", *o);
+                //if(*o < minval) *o = minval;
+                //else if(*o > maxval) *o = maxval;
+            }
+        }
+        need2free = TRUE;
+        DBG("converted");
+    }
+    if(keys){ // there's keys
+        size_t i;
+        char **records = fits->keylist;
+        for(i = 0; i < keys; ++i){
+            char *rec = records[i];
+            if(strncmp(rec, "SIMPLE", 6) == 0 || strncmp(rec, "EXTEND", 6) == 0) // key "file does conform ..."
+                continue;
+            else if(strncmp(rec, "COMMENT", 7) == 0) // comment of obligatory key in FITS head
+                continue;
+            else if(strncmp(rec, "NAXIS", 5) == 0 || strncmp(rec, "BITPIX", 6) == 0) // NAXIS, NAXISxxx, BITPIX
+                continue;
+            FITSFUN(fits_write_record, fp, rec);
+        }
+    }
+    FITSFUN(fits_write_record, fp, "COMMENT  modified by loccorr");
+    fitsstatus = 0;
+    fits_write_img(fp, FITSDATATYPE, 1, sz, outp, &fitsstatus);
+    if(need2free) FREE(outp);
+    if(fitsstatus){
+        fits_report_error(stderr, fitsstatus);
+        return FALSE;
+    }
+    TRYFITS(fits_close_file, fp);
+    return TRUE;
+}

LocCorr/fits.h

@@ -0,0 +1,55 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#ifndef FITS_H__
+#define FITS_H__
+
+#include <fitsio.h>
+#include <omp.h>
+#include <stdio.h>
+#include <stdbool.h>
+
+#define Stringify(x) #x
+#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x))
+
+
+typedef float Imtype; // maybe float or double only
+// this is TFLOAT or TDOUBLE depending on Imtype
+#define FITSDATATYPE    TFLOAT
+/*
+typedef double Imtype;
+#define FITSDATATYPE    TDOUBLE
+*/
+
+typedef struct{
+    int width;			// width
+    int height;			// height
+    int dtype;			// data type for image storage
+    Imtype *data;		// picture data
+    Imtype minval;      // extremal data values
+    Imtype maxval;
+    char **keylist;		// list of options for each key
+    size_t keynum;		// full number of keys (size of *keylist)
+} Image;
+
+void Image_free(Image **ima);
+bool FITS_read(const char *filename, Image **fits);
+bool FITS_write(const char *filename, const Image *fits);
+
+#endif // FITS_H__

LocCorr/grasshopper.c

@@ -0,0 +1,320 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <C/FlyCapture2_C.h>
+#include <C/FlyCapture2Defs_C.h>
+#include <float.h> // FLT_EPSILON
+#include <math.h>
+#include <stdio.h>
+#include <usefull_macros.h>
+
+#include "cmdlnopts.h"
+#include "fits.h"
+#include "grasshopper.h"
+#include "imagefile.h"
+
+static fc2Context context;
+static fc2PGRGuid guid;
+static fc2Error err = FC2_ERROR_OK;
+static float exptime = 10.; // exposition time in milliseconds
+static float gain = 0.;
+
+#define FC2FN(fn, ...) do{err = FC2_ERROR_OK; if(FC2_ERROR_OK != (err=fn(context __VA_OPT__(,) __VA_ARGS__))){ \
+    WARNX(#fn "(): %s", fc2ErrorToDescription(err)); return 0;}}while(0)
+
+/**
+ * @brief setfloat - set absolute property value (float)
+ * @param t        - type of property
+ * @param context  - initialized context
+ * @param f        - new value
+ * @return 1 if all OK
+ */
+static int setfloat(fc2PropertyType t, float f){
+    fc2Property prop = {0};
+    prop.type = t;
+    fc2PropertyInfo i = {0};
+    i.type = t;
+    FC2FN(fc2GetProperty, &prop);
+    FC2FN(fc2GetPropertyInfo, &i);
+    if(!prop.present || !i.present) return 0;
+    if(prop.autoManualMode){
+        if(!i.manualSupported){
+            WARNX("Can't set auto-only property");
+            return 0;
+        }
+        prop.autoManualMode = false;
+    }
+    if(!prop.absControl){
+        if(!i.absValSupported){
+            WARNX("Can't set non-absolute property to absolute value");
+            return 0;
+        }
+        prop.absControl = true;
+    }
+    if(!prop.onOff){
+        if(!i.onOffSupported){
+            WARNX("Can't set property ON");
+            return 0;
+        }
+        prop.onOff = true;
+    }
+    if(prop.onePush && i.onePushSupported) prop.onePush = false;
+    prop.valueA = prop.valueB = 0;
+    prop.absValue = f;
+    FC2FN(fc2SetProperty, &prop);
+    // now check
+    FC2FN(fc2GetProperty, &prop);
+    if(fabsf(prop.absValue - f) > 0.02f){
+        WARNX("Can't set property! Got %g instead of %g.", prop.absValue, f);
+        return 0;
+    }
+    return 1;
+}
+
+int propOnOff(fc2PropertyType t, BOOL onOff){
+    fc2Property prop = {0};
+    prop.type = t;
+    fc2PropertyInfo i = {0};
+    i.type = t;
+    FC2FN(fc2GetPropertyInfo, &i);
+    FC2FN(fc2GetProperty, &prop);
+    if(!prop.present || !i.present) return 0;
+    if(prop.onOff == onOff) return 0;
+    if(!i.onOffSupported){
+        WARNX("Property doesn't support state OFF");
+        return  0;
+    }
+    prop.onOff = onOff;
+    FC2FN(fc2SetProperty, &prop);
+    FC2FN(fc2GetProperty, &prop);
+    if(prop.onOff != onOff){
+        WARNX("Can't change property OnOff state");
+        return 0;
+    }
+    return 1;
+}
+
+#define autoExpOff()        propOnOff(FC2_AUTO_EXPOSURE, false)
+#define whiteBalOff()       propOnOff(FC2_WHITE_BALANCE, false)
+#define gammaOff()          propOnOff(FC2_GAMMA, false)
+#define trigModeOff()       propOnOff(FC2_TRIGGER_MODE, false)
+#define trigDelayOff()      propOnOff(FC2_TRIGGER_DELAY, false)
+#define frameRateOff()      propOnOff(FC2_FRAME_RATE, false)
+// +set: saturation, hue, sharpness
+#define setbrightness(b)    setfloat(FC2_BRIGHTNESS, b)
+#define setexp(e)           setfloat(FC2_SHUTTER, e)
+#define setgain(g)          setfloat(FC2_GAIN, g)
+
+static int connected = 0;
+static void disconnect(){
+    connected = 0;
+    fc2DestroyContext(context);
+}
+
+static int changeformat(){
+    FNAME();
+    BOOL b;
+    fc2Format7Info f7i = {.mode = FC2_MODE_0};
+    FC2FN(fc2GetFormat7Info, &f7i, &b);
+    if(!b) return 0;
+    fc2Format7ImageSettings f7;
+/*    unsigned int packSz; float perc;
+    FC2FN(fc2GetFormat7Configuration, &f7, &packSz, &perc);
+    DBG("packsz=%u, perc=%f, off=%d/%d, HW=%d/%d", packSz, perc, f7.offsetX, f7.offsetY, f7.height, f7.width);
+    */
+    f7.mode = FC2_MODE_0;
+    f7.offsetX = (GP->xoff < (int)f7i.maxWidth && GP->xoff > -1) ? GP->xoff : 0;
+    f7.offsetY = (GP->yoff < (int)f7i.maxHeight && GP->yoff > -1) ? GP->yoff : 0;
+    f7.width = (f7.offsetX+GP->width <= f7i.maxWidth && GP->width > 1) ? (unsigned int)GP->width : f7i.maxWidth - f7.offsetX;
+    f7.height = (f7.offsetY+GP->height <= f7i.maxHeight && GP->height > 1) ? (unsigned int)GP->height : f7i.maxHeight - f7.offsetY;
+    DBG("offx=%d, offy=%d, w=%d, h=%d ", f7.offsetX, f7.offsetY, f7.width, f7.height);
+    f7.pixelFormat = FC2_PIXEL_FORMAT_MONO8;
+    fc2Format7PacketInfo f7p;
+    FC2FN(fc2ValidateFormat7Settings, &f7, &b, &f7p);
+    if(!b) return 0; // invalid
+    FC2FN(fc2SetFormat7Configuration, &f7, f7p.recommendedBytesPerPacket);
+    return 1;
+}
+
+static int connect(){
+    FNAME();
+    if(connected) return 1;
+    unsigned int numCameras = 0;
+    if(FC2_ERROR_OK != (err = fc2CreateContext(&context))){
+        WARNX("fc2CreateContext(): %s", fc2ErrorToDescription(err));
+        return 0;
+    }
+    FC2FN(fc2GetNumOfCameras, &numCameras);
+    if(numCameras == 0){
+        WARNX("No cameras detected!");
+        disconnect();
+        return 0;
+    }
+    DBG("Found %d camera[s]", numCameras);
+    if(numCameras > 1){
+        WARNX("Found %d cameras, will use first", numCameras);
+    }
+    FC2FN(fc2GetCameraFromIndex, 0, &guid);
+    FC2FN(fc2Connect, &guid);
+    // turn off all shit
+    autoExpOff();
+    whiteBalOff();
+    gammaOff();
+    trigModeOff();
+    trigDelayOff();
+    frameRateOff();
+    if(!changeformat()) WARNX("Can't change camera format");
+    connected = 1;
+    return 1;
+}
+
+static int GrabImage(fc2Image *convertedImage){
+    FNAME();
+    fc2Image rawImage;
+    // start capture
+    FC2FN(fc2StartCapture);
+    err = fc2CreateImage(&rawImage);
+    if(err != FC2_ERROR_OK){
+        WARNX("Error in fc2CreateImage: %s", fc2ErrorToDescription(err));
+        return 0;
+    }
+    // Retrieve the image
+    err = fc2RetrieveBuffer(context, &rawImage);
+    if(err != FC2_ERROR_OK){
+        WARNX("Error in fc2RetrieveBuffer: %s", fc2ErrorToDescription(err));
+        return 0;
+    }
+    // Convert image to gray
+    err = fc2ConvertImageTo(FC2_PIXEL_FORMAT_MONO8, &rawImage, convertedImage);
+    if(err != FC2_ERROR_OK){
+        WARNX("Error in fc2ConvertImageTo: %s", fc2ErrorToDescription(err));
+        return 0;
+    }
+    fc2StopCapture(context);
+    fc2DestroyImage(&rawImage);
+    return 1;
+}
+
+static void calcexpgain(float newexp){
+    DBG("recalculate exposition: oldexp=%g, oldgain=%g, newexp=%g", exptime, gain, newexp);
+    if(newexp < exptime){ // first we should make gain lower
+        if(10.*newexp < GP->maxexp){
+            if(gain > 10.){
+                gain = 10.;
+                newexp *= 10.;
+            }else if(gain > 0.){
+                gain = 0.;
+                newexp *= 10.;
+            }
+        }
+    }else{ // if new exposition too large, increase gain
+        if(newexp > GP->maxexp){
+            if(gain < 19.){
+                gain += 10.;
+                newexp /= 10.;
+            }
+        }
+    }
+    if(newexp < GP->minexp) newexp = GP->minexp;
+    else if(newexp > GP->maxexp) newexp = GP->maxexp;
+    exptime = newexp;
+    DBG("New values: exp=%g, gain=%g", exptime, gain);
+}
+
+//convertedImage.pData, convertedImage.cols, convertedImage.rows, convertedImage.stride
+static void recalcexp(fc2Image *img){
+    uint8_t *data = img->pData;
+    int W = img->cols, H = img->rows, S = img->stride;
+    int histogram[256] = {0};
+    DBG("W=%d, H=%d, S=%d", W, H, S);
+    for(int y = 0; y < H; ++y){
+        uint8_t *ptr = &data[y*S];
+        for(int x = 0; x < W; ++x, ++ptr)
+            ++histogram[*ptr];
+    }
+    int idx100, sum100 = 0;
+    for(idx100 = 255; idx100 >= 0; --idx100){
+        sum100 += histogram[idx100];
+        if(sum100 > 100) break;
+    }
+    DBG("Sum100=%d, idx100=%d", sum100, idx100);
+    if(idx100 > 200 && idx100 < 250) return; // good values
+    if(idx100 > 250){ // exposure too long
+        calcexpgain(0.7*exptime);
+    }else{ // exposure too short
+        if(idx100 > 5)
+            calcexpgain(exptime * 210. / (float)idx100);
+        else
+            calcexpgain(exptime * 50.);
+    }
+}
+
+int capture_grasshopper(void (*process)(Image*)){
+    FNAME();
+    static float oldexptime = 0.;
+    static float oldgain = -1.;
+    fc2Image convertedImage;
+    err = fc2CreateImage(&convertedImage);
+    if(err != FC2_ERROR_OK){
+        WARNX("capture_grasshopper(): can't create image, %s", fc2ErrorToDescription(err));
+        return 0;
+    }
+    Image *oIma = NULL;
+    while(1){
+        if(!connect()){ // wait until camera be powered on
+            DBG("Disconnected");
+            sleep(1);
+            continue;
+        }
+        if(fabsf(oldexptime - exptime) > FLT_EPSILON){ // new exsposition value
+            red("Change exptime to %.2fms\n", exptime);
+            if(setexp(exptime)){
+                oldexptime = exptime;
+            }else{
+                WARNX("Can't change exposition time to %gms", exptime);
+                disconnect();
+                continue;
+            }
+        }
+        if(fabs(oldgain - gain) > FLT_EPSILON){ // change gain
+            red("Change gain to %g\n", gain);
+            if(setgain(gain)){
+                oldgain = gain;
+            }else{
+                WARNX("Can't change gain to %g", gain);
+                disconnect();
+                continue;
+            }
+        }
+        if(!GrabImage(&convertedImage)){
+            WARNX("Can't grab image");
+            disconnect();
+            continue;
+        }
+        if(!process) continue;
+        recalcexp(&convertedImage);
+        oIma = u8toImage(convertedImage.pData, convertedImage.cols, convertedImage.rows, convertedImage.stride);
+        if(oIma){
+            process(oIma);
+            FREE(oIma);
+        }
+    }
+    fc2DestroyImage(&convertedImage);
+    fc2DestroyContext(context);
+    return 1;
+}

LocCorr/grasshopper.h

@@ -0,0 +1,28 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+#ifndef GRASSHOPPER_H__
+#define GRASSHOPPER_H__
+
+#include "fits.h" // Image*
+
+#define GRASSHOPPER_CAPT_NAME   "grasshopper"
+
+int capture_grasshopper(void (*process)(Image *));
+
+#endif // GRASSHOPPER_H__

LocCorr/imagefile.c

@@ -0,0 +1,514 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <dirent.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <usefull_macros.h>
+
+#include <stb/stb_image.h>
+#include <stb/stb_image_write.h>
+
+#include "cmdlnopts.h"
+#include "draw.h"
+#include "grasshopper.h"
+#include "imagefile.h"
+#include "median.h"
+
+// weights to convert colour image into gray (sum=1!)
+#define WEIGHT_RED      ()
+#define WEIGHT_GREEN    ()
+#define WEIGHT_BLUE     ()
+
+typedef struct{
+    const char signature[8];
+    uint8_t len;
+    InputType it;
+} imsign;
+
+const imsign signatures[] = {
+    {"BM", 2, T_BMP},
+    {"SIMPLE", 6, T_FITS},
+    {{0x1f, 0x8b, 0x08}, 3, T_GZIP},
+    {"GIF8", 4, T_GIF},
+    {{0xff, 0xd8, 0xff, 0xdb}, 4, T_JPEG},
+    {{0xff, 0xd8, 0xff, 0xe0}, 4, T_JPEG},
+    {{0xff, 0xd8, 0xff, 0xe1}, 4, T_JPEG},
+    {{0x89, 0x50, 0x4e, 0x47}, 4, T_PNG},
+   // {{0x49, 0x49, 0x2a, 0x00}, 4, T_TIFF},
+    {"", 0, T_WRONG}
+};
+
+#ifdef EBUG
+static char *hexdmp(const char sig[8]){
+    static char buf[128];
+    char *bptr = buf;
+    bptr += sprintf(bptr, "[ ");
+    for(int i = 0; i < 7; ++i){
+        bptr += sprintf(bptr, "%02X ", (uint8_t)sig[i]);
+    }
+    bptr += sprintf(bptr, "]");
+    return buf;
+}
+#endif
+
+static InputType imtype(FILE *f){
+    char signature[8];
+    int x = fread(signature, 1, 7, f);
+    DBG("x=%d", x);
+    if(7 != x){
+        WARN("Can't read file signature");
+        return T_WRONG;
+    }
+    signature[7] = 0;
+    const imsign *s = signatures;
+    DBG("Got signature: %s (%s)", hexdmp(signature), signature);
+    while(s->len){
+        DBG("Check %s", s->signature);
+        if(0 == memcmp(s->signature, signature, s->len)){
+            DBG("Found signature %s", s->signature);
+            return s->it;
+        }
+        ++s;
+    }
+    return T_WRONG;
+}
+
+/**
+ * @brief chkinput - check file/directory name
+ * @param name - name of file or directory
+ * @return type of `name`
+ */
+InputType chkinput(const char *name){
+    DBG("input name: %s", name);
+    if(0 == strcmp(name, GRASSHOPPER_CAPT_NAME)) return T_CAPT_GRASSHOPPER;
+    struct stat fd_stat;
+    stat(name, &fd_stat);
+    if(S_ISDIR(fd_stat.st_mode)){
+        DBG("%s is a directory", name);
+        DIR *d = opendir(name);
+        if(!d){
+            WARN("Can't open directory %s", name);
+            return T_WRONG;
+        }
+        closedir(d);
+        return T_DIRECTORY;
+    }
+    FILE *f = fopen(name, "r");
+    if(!f){
+        WARN("Can't open file %s", name);
+        return T_WRONG;
+    }
+    InputType tp = imtype(f);
+    DBG("Image type of %s is %d", name, tp);
+    fclose(f);
+    return tp;
+}
+
+Image *u8toImage(uint8_t *data, int width, int height, int stride){
+    Image *outp = MALLOC(Image, 1);
+    outp->width = width;
+    outp->height = height;
+    outp->dtype = FLOAT_IMG;
+    uint8_t min = *data, max = min;
+    for(int y = 0; y < height; ++y){
+        uint8_t *ptr = &data[y*stride];
+        for(int x = 0; x < width; ++x){
+            uint8_t p = *ptr++;
+            if(p < min) min = p;
+            else if(p > max) max = p;
+        }
+    }
+    outp->minval = (Imtype) min;
+    outp->maxval = (Imtype) max;
+    DBG("\nMAX=%g, MIN=%g\n", outp->maxval, outp->minval);
+    outp->data = MALLOC(Imtype, width*height);
+    // flip image updown for FITS coordinate system
+    OMP_FOR()
+    for(int y = 0; y < height; ++y){
+        Imtype *Out = &outp->data[(height-1-y)*width];
+        uint8_t *In = &data[y*stride];
+        for(int x = 0; x < width; ++x){
+            *Out++ = (Imtype)(*In++);
+        }
+    }
+    return outp;
+}
+
+// load other image file
+static Image *im_load(const char *name){
+    int width, height, channels;
+    uint8_t *img = stbi_load(name, &width, &height, &channels, 1);
+    if(!img){
+        WARNX("Error in loading the image %s\n", name);
+        return NULL;
+    }
+    return u8toImage(img, width, height, width);
+}
+
+/**
+ * @brief Image_read - read image from any supported file type
+ * @param name - path to image
+ * @return image or NULL if failed
+ */
+Image *Image_read(const char *name){
+    InputType tp = chkinput(name);
+    if(tp == T_DIRECTORY || tp == T_WRONG){
+        WARNX("Bad file type to read");
+        return NULL;
+    }
+    Image *outp = NULL;
+    if(tp == T_FITS || tp == T_GZIP){
+        if(!FITS_read(name, &outp)){
+            WARNX("Can't read %s", name);
+            return NULL;
+        }
+    }else outp = im_load(name);
+    return outp;
+}
+
+/**
+ * @brief Image_new - allocate memory for new struct Image & Image->data
+ * @param w, h - image size
+ * @return data allocated here
+ */
+Image *Image_new(int w, int h){
+    Image *outp = MALLOC(Image, 1);
+    outp->width = w;
+    outp->height = h;
+    outp->data = MALLOC(Imtype, w*h);
+    return outp;
+}
+
+/**
+ * @brief Image_sim - allocate memory for new empty Image with similar size & data type
+ * @param i - sample image
+ * @return data allocated here (with empty keylist & zeros in data)
+ */
+Image *Image_sim(const Image *i){
+    if(!i) return NULL;
+    if((i->width * i->height) < 1) return NULL;
+    Image *outp = Image_new(i->width, i->height);
+    outp->dtype = i->dtype;
+    return outp;
+}
+
+/**
+ * @brief linear - linear transform (mirror image upside down!)
+ * @param I - input image
+ * @param nchannels - 1 or 3 colour channels
+ * @return allocated here image for jpeg/png storing
+ */
+uint8_t *linear(const Image *I, int nchannels){ // only 1 and 3 channels supported!
+    if(nchannels != 1 && nchannels != 3) return NULL;
+    int width = I->width, height = I->height;
+    size_t stride = width*nchannels, S = height*stride;
+    uint8_t *outp = MALLOC(uint8_t, S);
+    Imtype min = I->minval, max = I->maxval, W = 255./(max - min);
+    DBG("make linear transform %dx%d, %d channels", I->width, I->height, nchannels);
+    if(nchannels == 3){
+        OMP_FOR()
+        for(int y = 0; y < height; ++y){
+            uint8_t *Out = &outp[(height-1-y)*stride];
+            Imtype *In = &I->data[y*width];
+            for(int x = 0; x < width; ++x){
+                Out[0] = Out[1] = Out[2] = (uint8_t)(W*((*In++) - min));
+                Out += 3;
+            }
+        }
+    }else{
+        OMP_FOR()
+        for(int y = 0; y < height; ++y){
+            uint8_t *Out = &outp[(height-1-y)*width];
+            Imtype *In = &I->data[y*width];
+            for(int x = 0; x < width; ++x){
+                *Out++ = (uint8_t)(W*((*In++) - min));
+            }
+        }
+    }
+    return outp;
+}
+
+/**
+ * @brief equalize - hystogram equalization (mirror image upside down!)
+ * @param I - input image
+ * @param nchannels - 1 or 3 colour channels
+ * @param throwpart - which part of black pixels (from all amount) to throw away
+ * @return allocated here image for jpeg/png storing
+ */
+uint8_t *equalize(const Image *I, int nchannels, double throwpart){
+    if(nchannels != 1 && nchannels != 3) return NULL;
+    int width = I->width, height = I->height;
+    size_t stride = width*nchannels, S = height*stride;
+    uint8_t *lin = linear(I, 1);
+    if(!lin) return NULL;
+    uint8_t *outp = MALLOC(uint8_t, S);
+    int orig_hysto[256] = {0}; // original hystogram (linear)
+    uint8_t eq_levls[256] = {0};   // levels to convert: newpix = eq_levls[oldpix]
+    int s = width*height;
+    for(int i = 0; i < s; ++i){
+        ++orig_hysto[lin[i]];
+    }
+
+    int Nblack = 0, bpart = (int)(throwpart * s);
+    int startidx;
+    // remove first part of black pixels
+    for(startidx = 0; startidx < 256; ++startidx){
+        Nblack += orig_hysto[startidx];
+        if(Nblack >= bpart) break;
+    }
+    ++startidx;
+    /* remove last part of white pixels
+    for(stopidx = 255; stopidx > startidx; --stopidx){
+        Nwhite += orig_hysto[stopidx];
+        if(Nwhite >= wpart) break;
+    }*/
+    DBG("Throw %d (real: %d black) pixels, startidx=%d", bpart, Nblack, startidx);
+/*
+    double part = (double)(s + 1) / 256., N = 0.;
+    for(int i = 0; i < 256; ++i){
+        N += orig_hysto[i];
+        eq_levls[i] = (uint8_t)(N/part);
+    }*/
+    double part = (double)(s + 1. - Nblack) / 256., N = 0.;
+    for(int i = startidx; i < 256; ++i){
+        N += orig_hysto[i];
+        eq_levls[i] = (uint8_t)(N/part);
+    }
+    //for(int i = stopidx; i < 256; ++i) eq_levls[i] = 255;
+#if 0
+    DBG("Original / new histogram");
+    for(int i = 0; i < 256; ++i) printf("%d\t%d\t%d\n", i, orig_hysto[i], eq_levls[i]);
+#endif
+    if(nchannels == 3){
+        OMP_FOR()
+        for(int y = 0; y < height; ++y){
+            uint8_t *Out = &outp[y*stride];
+            uint8_t *In = &lin[y*width];
+            for(int x = 0; x < width; ++x){
+                Out[0] = Out[1] = Out[2] = eq_levls[*In++];
+                Out += 3;
+            }
+        }
+    }else{
+        OMP_FOR()
+        for(int y = 0; y < height; ++y){
+            uint8_t *Out = &outp[y*width];
+            uint8_t *In = &lin[y*width];
+            for(int x = 0; x < width; ++x){
+                *Out++ = eq_levls[*In++];
+            }
+        }
+    }
+    FREE(lin);
+    return outp;
+}
+
+/**
+ * @brief Image_write_jpg - save image as JPG file (flipping upside down)
+ * @param I - image
+ * @param name - filename
+ * @param eq == 0 to write linear, != 0 to write equalized image
+ * @return 0 if failed
+ */
+int Image_write_jpg(const Image *I, const char *name, int eq){
+    if(!I || !I->data) return 0;
+    uint8_t *outp = NULL;
+    if(eq)
+        outp = equalize(I, 1, GP->throwpart);
+    else
+        outp = linear(I, 1);
+    DBG("Try to write %s", name);
+    int r = stbi_write_jpg(name, I->width, I->height, 1, outp, 95);
+    FREE(outp);
+    return r;
+}
+
+// calculate extremal values of image data and store them in it
+void Image_minmax(Image *I){
+    if(!I || !I->data) return;
+    Imtype min = *(I->data), max = min;
+    int wh = I->width * I->height;
+#ifdef EBUG
+    double t0 = dtime();
+#endif
+    /*
+    int N = omp_get_max_threads();
+    Imtype arrmin[N], arrmax[N];
+    for(int i = 0; i < N; ++i){
+        arrmin[i] = min; arrmax[i] = max;
+    }
+    OMP_FOR()
+    for(int i = 0; i < wh; ++i){
+        if(I->data[i] < arrmin[omp_get_thread_num()]) arrmin[omp_get_thread_num()] = I->data[i];
+        else if(I->data[i] > arrmax[omp_get_thread_num()]) arrmax[omp_get_thread_num()] = I->data[i];
+    }
+
+    for(int i = 0; i < N; ++i){
+        if(min > arrmin[i]) min = arrmin[i];
+        if(max < arrmax[i]) max = arrmax[i];
+    }*/
+    #pragma omp parallel shared(min, max)
+    {
+        int min_p = min, max_p = min;
+        #pragma omp for nowait
+        for(int i = 0; i < wh; ++i){
+            if(I->data[i] < min_p) min_p = I->data[i];
+            else if(I->data[i] > max_p) max_p = I->data[i];
+        }
+        #pragma omp critical
+        {
+            if(min > min_p) min = min_p;
+            if(max < max_p) max = max_p;
+        }
+    }
+    I->maxval = max;
+    I->minval = min;
+    DBG("Image_minmax(): Min=%g, Max=%g, time: %gms", min, max, (dtime()-t0)*1e3);
+}
+
+/*
+ * =================== CONVERT IMAGE TYPES ===================>
+ */
+
+// convert binarized image into floating
+/**
+ * @brief bin2Im - convert binarized image into floating
+ * @param image - binarized image
+ * @param W, H  - its size (in pixels!)
+ * @return Image structure
+ */
+Image *bin2Im(uint8_t *image, int W, int H){
+    Image *ret = Image_new(W, H);
+    int stride = (W + 7) / 8, s1 = (stride*8 == W) ? stride : stride - 1;
+    OMP_FOR()
+    for(int y = 0; y < H; y++){
+        Imtype *optr = &ret->data[y*W];
+        uint8_t *iptr = &image[y*stride];
+        for(int x = 0; x < s1; x++){
+            register uint8_t inp = *iptr++;
+            for(int i = 0; i < 8; ++i){
+                *optr++ = (inp & 0x80) ? 1. : 0;
+                inp <<= 1;
+            }
+        }
+        int rest = W - s1*8;
+        register uint8_t inp = *iptr;
+        for(int i = 0; i < rest; ++i){
+            *optr++ = (inp & 0x80) ? 1. : 0;
+            inp <<= 1;
+        }
+    }
+    ret->dtype = BYTE_IMG;
+    ret->minval = 0;
+    ret->maxval = 1;
+    return ret;
+}
+
+/**
+ * Convert floatpoint image into pseudo-packed (1 char == 8 pixels), all values > 0 will be 1, else - 0
+ * @param im (i)     - image to convert
+ * @param stride (o) - new width of image
+ * @param bk         - background level (all values < bk will be 0, other will be 1)
+ * @return allocated memory area with "packed" image
+ */
+uint8_t *Im2bin(const Image *im, Imtype bk){
+    if(!im) return NULL;
+    int W = im->width, H = im->height;
+    if(W < 2 || H < 2) return NULL;
+    int y, W0 = (W + 7) / 8, s1 = (W/8 == W0) ? W0 : W0 - 1;
+    uint8_t *ret = MALLOC(uint8_t, W0 * H);
+    OMP_FOR()
+    for(y = 0; y < H; y++){
+        Imtype *iptr = &im->data[y*W];
+        uint8_t *optr = &ret[y*W0];
+        register uint8_t o;
+        for(int x = 0; x < s1; ++x){
+            o = 0;
+            for(int i = 0; i < 8; ++i){
+                o <<= 1;
+                if(*iptr++ > bk) o |= 1;
+            }
+            *optr++ = o;
+        }
+        int rest = 7 - (W - s1*8);
+        if(rest < 7){
+            o = 0;
+            for(int x = 7; x > rest; --x){
+                if(*iptr++ > 0.) o |= 1<<x;
+            }
+            *optr = o;
+            //*optr = o | 1<<rest; // mark outern right edge (for good CC labeling)
+        }
+    }
+    return ret;
+}
+
+// convert size_t labels into floating
+Image *ST2Im(size_t *image, int W, int H){
+    Image *ret = Image_new(W, H);
+    OMP_FOR()
+    for(int y = 0; y < H; ++y){
+        Imtype *optr = &ret->data[y*W];
+        size_t *iptr = &image[y*W];
+        for(int x = 0; x < W; ++x){
+            *optr++ = (Imtype)*iptr++;
+        }
+    }
+    ret->dtype = FLOAT_IMG;
+    Image_minmax(ret);
+    return ret;
+}
+
+/**
+ * Convert "packed" image into size_t array for conncomp procedure
+ * @param image (i) - input image
+ * @param W, H      - size of image in pixels
+ * @return allocated memory area with copy of an image
+ */
+size_t *bin2ST(uint8_t *image, int W, int H){
+    size_t *ret = MALLOC(size_t, W * H);
+    int W0 = (W + 7) / 8, s1 = W0 - 1;
+    OMP_FOR()
+    for(int y = 0; y < H; y++){
+        size_t *optr = &ret[y*W];
+        uint8_t *iptr = &image[y*W0];
+        for(int x = 0; x < s1; ++x){
+            register uint8_t inp = *iptr++;
+            for(int i = 0; i < 8; ++i){
+                *optr++ = (inp & 0x80) ? 1 : 0;
+                inp <<= 1;
+            }
+        }
+        int rest = W - s1*8;
+        register uint8_t inp = *iptr;
+        for(int i = 0; i < rest; ++i){
+            *optr++ = (inp & 0x80) ? 1 : 0;
+            inp <<= 1;
+        }
+    }
+    return ret;
+}
+
+
+/*
+ * <=================== CONVERT IMAGE TYPES ===================
+ */

LocCorr/imagefile.h

@@ -0,0 +1,55 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+#ifndef IMAGEFILE_H__
+#define IMAGEFILE_H__
+
+#include <stdint.h>
+
+#include "fits.h"
+
+// input file/directory type
+typedef enum{
+    T_WRONG,
+    T_DIRECTORY,
+    T_BMP,
+    T_FITS,
+    T_GZIP,
+    T_GIF,
+    T_JPEG,
+    T_PNG,
+    T_CAPT_GRASSHOPPER // capture grasshopper
+} InputType;
+
+void Image_minmax(Image *I);
+uint8_t *linear(const Image *I, int nchannels);
+uint8_t *equalize(const Image *I, int nchannels, double throwpart);
+InputType chkinput(const char *name);
+Image *Image_read(const char *name);
+Image *Image_new(int w, int h);
+Image *Image_sim(const Image *i);
+void Image_free(Image **I);
+int Image_write_jpg(const Image *I, const char *name, int equalize);
+
+Image *u8toImage(uint8_t *data, int width, int height, int stride);
+Image *bin2Im(uint8_t *image, int W, int H);
+uint8_t *Im2bin(const Image *im, Imtype bk);
+size_t *bin2ST(uint8_t *image, int W, int H);
+Image *ST2Im(size_t *image, int W, int H);
+
+#endif // IMAGEFILE_H__

LocCorr/improc.c

@@ -0,0 +1,303 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+#include <stb/stb_image_write.h>
+#include <stdbool.h>
+#include <string.h>
+#include <time.h>
+#include <usefull_macros.h>
+
+#include "binmorph.h"
+#include "cmdlnopts.h"
+#include "draw.h"
+#include "grasshopper.h"
+#include "fits.h"
+#include "improc.h"
+#include "inotify.h"
+#include "median.h"
+
+// how many frames will be averaged to count image deviation
+#define AVERAGING_ARRAY_SIZE        (5)
+// tolerance of deviations by X and Y axis
+#define XY_TOLERANCE                (1.)
+
+static FILE *fXYlog = NULL;
+
+static double Xtarget = -1., Ytarget = -1.; // target coordinates
+static double tstart = 0.; // time of logging start
+
+typedef struct{
+    uint32_t area;      // object area in pixels
+    double Isum;        // total object's intensity over background
+    double WdivH;       // width of object's box divided by height
+    double xc;          // centroid coordinates
+    double yc;
+    double xsigma;      // STD by horizontal and vertical axes
+    double ysigma;
+} object;
+
+static bool save_fits(Image *I, const char *name){
+    uint8_t *outp = NULL;
+    if(GP->equalize)
+        outp = equalize(I, 1, GP->throwpart);
+    else
+        outp = linear(I, 1);
+    char fname[PATH_MAX];
+    char *flname = strdup(name);
+    char *pt = strrchr(flname, '.');
+    if(pt) *pt = 0;
+    snprintf(fname, PATH_MAX, "%s.jpg", flname);
+    stbi_write_jpg(fname, I->width, I->height, 1, outp, 95);
+    FREE(outp);
+    unlink(name);
+    return FITS_write(name, I);
+}
+
+static void savebin(uint8_t *b, int W, int H, const char *name){
+    Image *I = bin2Im(b, W, H);
+    if(I){
+        save_fits(I, name);
+        Image_free(&I);
+    }
+}
+
+// function for Qsort
+static int compObjs(const void *a, const void *b){
+    const object *oa = (const object*)a;
+    const object *ob = (const object*)b;
+    double idiff = (oa->Isum - ob->Isum)/(oa->Isum + ob->Isum);
+    if(fabs(idiff) > GP->intensthres) return (idiff > 0) ? -1:1;
+    double r2a = oa->xc * oa->xc + oa->yc * oa->yc;
+    double r2b = ob->xc * ob->xc + ob->yc * ob->yc;
+    return (r2a < r2b) ? -1 : 1;
+}
+
+static void getDeviation(object *curobj){
+    if(Xtarget < 0. || Ytarget < 0.){ // init from user parameters
+        Xtarget = GP->xtarget; Ytarget = GP->ytarget;
+    }
+    static double Xc[AVERAGING_ARRAY_SIZE], Yc[AVERAGING_ARRAY_SIZE];
+    static int counter = 0;
+    Xc[counter] = curobj->xc; Yc[counter] = curobj->yc;
+    if(fXYlog){ // make log record
+        fprintf(fXYlog, "%.2f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t",
+                dtime() - tstart, curobj->xc, curobj->yc,
+                curobj->xsigma, curobj->ysigma, curobj->WdivH);
+    }
+    if(++counter != AVERAGING_ARRAY_SIZE){
+        if(fXYlog) fprintf(fXYlog, "\n");
+        return;
+    }
+    // it's time to calculate average deviations
+    counter = 0;
+    double xdev = 0., ydev = 0., xsum2 = 0., ysum2 = 0.;
+    for(int i = 0; i < AVERAGING_ARRAY_SIZE; ++i){
+        double dx = Xc[i] - Xtarget, dy = Yc[i] - Ytarget;
+        xdev += dx; ydev += dy;
+        xsum2 += dx*dx; ysum2 += dy*dy;
+    }
+    xdev /= AVERAGING_ARRAY_SIZE; ydev /= AVERAGING_ARRAY_SIZE;
+    double Sx = sqrt(xsum2/AVERAGING_ARRAY_SIZE - xdev*xdev);
+    double Sy = sqrt(ysum2/AVERAGING_ARRAY_SIZE - ydev*ydev);
+    DBG("xtag=%g, ytag=%g", Xtarget, Ytarget);
+    green("\n\n\n Deviations: dX=%g (+-%g), dY=%g (+-%g)\n", xdev, Sx, ydev, Sy);
+    if(Xtarget < 0. || Ytarget < 0. || fabs(xdev) < XY_TOLERANCE || fabs(ydev) < XY_TOLERANCE){
+        DBG("Target coordinates not defined or correction too small");
+        if(fXYlog) fprintf(fXYlog, "\n");
+        return;
+    }
+    if(fXYlog){
+        fprintf(fXYlog, "%.1f\t%.1f\t%.1f\t%.1f\n", xdev, ydev, Sx, Sy);
+    }
+}
+
+void process_file(Image *I){
+#ifdef EBUG
+    double t0 = dtime(), tlast = t0;
+#define DELTA(p) do{double t = dtime(); DBG("---> %s @ %gms (delta: %gms)", p, (t-t0)*1e3, (t-tlast)*1e3); tlast = t;}while(0)
+#else
+#define DELTA(x)
+#endif
+    // I - original image
+    // M - median filtering
+    // mean - local mean
+    // std  - local STD
+    /**** read original image ****/
+    DELTA("Imread");
+    if(!I){
+        WARNX("No image");
+        return;
+    }
+    int W = I->width, H = I->height;
+    if(!I->dtype) I->dtype = FLOAT_IMG;
+    save_fits(I, "fitsout.fits");
+    DELTA("Save original");
+/*
+    uint8_t *outp = NULL;
+    if(GP->equalize)
+        outp = equalize(I, 1, GP->throwpart);
+    else
+        outp = linear(I, 1);
+    stbi_write_jpg("jpegout.jpg", I->width, I->height, 3, outp, 95);
+    FREE(outp);
+*/
+    Imtype bk;
+    if(calc_background(I, &bk)){
+        DBG("backgr = %g", bk);
+        DELTA("Got background");
+        uint8_t *ibin = Im2bin(I, bk);
+        DELTA("Made binary");
+        if(ibin){
+            savebin(ibin, W, H, "binary.fits");
+            DELTA("save binary.fits");
+            uint8_t *er = erosionN(ibin, W, H, GP->nerosions);
+            FREE(ibin);
+            DELTA("Erosion");
+            savebin(er, W, H, "erosion.fits");
+            DELTA("Save erosion");
+            uint8_t *opn = dilationN(er, W, H, GP->ndilations);
+            FREE(er);
+            DELTA("Opening");
+            savebin(opn, W, H, "opening.fits");
+            DELTA("Save opening");
+            ConnComps *cc;
+            size_t *S = cclabel4(opn, W, H, &cc);
+            FREE(opn);
+            if(cc->Nobj > 1){
+                object *Objects = MALLOC(object, cc->Nobj-1);
+                int objctr = 0;
+                for(size_t i = 1; i < cc->Nobj; ++i){
+                    Box *b = &cc->boxes[i];
+                    double wh = ((double)b->xmax - b->xmin)/(b->ymax - b->ymin);
+                    //DBG("Obj# %zd: wh=%g, area=%d", i, wh, b->area);
+                    // TODO: change magick numbers to parameters
+                    if(wh < 0.2 || wh > 5.) continue;
+                    if((int)b->area < GP->minarea) continue;
+                    double xc = 0., yc = 0.;
+                    double x2c = 0., y2c = 0., Isum = 0.;
+                    for(size_t y = b->ymin; y <= b->ymax; ++y){
+                        size_t idx = y*W + b->xmin;
+                        size_t *maskptr = &S[idx];
+                        Imtype *Iptr = &I->data[idx];
+                        for(size_t x = b->xmin; x <= b->xmax; ++x, ++maskptr, ++Iptr){
+                            if(*maskptr != i) continue;
+                            double intens = (double) (*Iptr - bk);
+                            if(intens < 0.) continue;
+                            double xw = x * intens, yw = y * intens;
+                            xc += xw;
+                            yc += yw;
+                            x2c += xw * x;
+                            y2c += yw * y;
+                            Isum += intens;
+                        }
+                    }
+                    xc /= Isum; yc /= Isum;
+                    x2c = x2c/Isum - xc*xc;
+                    y2c = y2c/Isum - yc*yc;
+                    Objects[objctr++] = (object){
+                        .area = b->area, .Isum = Isum,
+                        .WdivH = wh, .xc = xc, .yc = yc,
+                        .xsigma = sqrt(x2c), .ysigma = sqrt(y2c)
+                    };
+                }
+                DELTA("Labeling");
+                printf("T%zd, N=%d\n", time(NULL), objctr);
+                qsort(Objects, objctr, sizeof(object), compObjs);
+                object *o = Objects;
+                green("%6s\t%6s\t%6s\t%6s\t%6s\t%6s\t%6s\t%6s\n",
+                      "N", "Area", "Mv", "W/H", "Xc", "Yc", "Sx", "Sy");
+                for(int i = 0; i < objctr; ++i, ++o){
+                    // 1.0857 = 2.5/ln(10)
+                    printf("%6d\t%6d\t%6.1f\t%6.1f\t%6.1f\t%6.1f\t%6.1f\t%6.1f\n",
+                           i, o->area, 20.-1.0857*log(o->Isum), o->WdivH, o->xc, o->yc, o->xsigma, o->ysigma);
+                }
+                getDeviation(Objects);
+                if(objctr){ // draw crosses @ objects' centers
+                    uint8_t *outp = NULL;
+                    if(GP->equalize)
+                        outp = equalize(I, 3, GP->throwpart);
+                    else
+                        outp = linear(I, 3);
+                    Pattern *cross = Pattern_cross(33, 33);
+                    int H = I->height;
+                    Img3 i3 = {.data = outp, .w = I->width, .h = H};
+                    Pattern_draw3(&i3, cross, Objects[0].xc, H-Objects[0].yc, C_G);
+                    for(int i = 1; i < objctr; ++i)
+                        Pattern_draw3(&i3, cross, Objects[i].xc, H-Objects[i].yc, C_R);
+                    Pattern_free(&cross);
+                    stbi_write_jpg("outpWcrosses.jpg", I->width, I->height, 3, outp, 95);
+                    FREE(outp);
+                }
+                FREE(cc);
+                FREE(Objects);
+                Image *c = ST2Im(S, W, H);
+                DELTA("conv size_t -> Ima");
+                save_fits(c, "size_t.fits");
+                Image_free(&c);
+                DELTA("Save size_t");
+                Image *obj = Image_sim(I);
+                OMP_FOR()
+                for(int y = 0; y < H; ++y){
+                    size_t idx = y*W;
+                    Imtype *optr = &obj->data[idx];
+                    Imtype *iptr = &I->data[idx];
+                    size_t *mask = &S[idx];
+                    for(int x = 0; x < W; ++x, ++mask, ++iptr, ++optr){
+                        if(*mask) *optr = *iptr - bk;
+                    }
+                }
+                Image_minmax(obj);
+                save_fits(obj, "object.fits");
+                Image_free(&obj);
+            }
+            FREE(S);
+        }
+    }
+    DELTA("End");
+}
+
+int process_input(InputType tp, char *name){
+    if(tp == T_DIRECTORY) return watch_directory(name, process_file);
+    else if(tp == T_CAPT_GRASSHOPPER) return capture_grasshopper(process_file);
+    return watch_file(name, process_file);
+}
+
+/**
+ * @brief openXYlog - open file to log XY values
+ * @param name - filename
+ */
+void openXYlog(const char *name){
+    closeXYlog();
+    fXYlog = fopen(name, "a");
+    if(!fXYlog){
+        char *e = strerror(errno);
+        WARNX("Can't create file %s: %s", name, e);
+        LOGERR("Can't create file %s: %s", name, e);
+        return;
+    }
+    time_t t = time(NULL);
+    fprintf(fXYlog, "# Start at: %s", ctime(&t));
+    fprintf(fXYlog, "# time,s\tXc\tYc\t\tSx\tSy\tW/H\tcorrX\tcorrY\tSXcorr\tSYcorr\n");
+    tstart = dtime();
+}
+void closeXYlog(){
+    if(!fXYlog) return;
+    fclose(fXYlog);
+    fXYlog = NULL;
+}

LocCorr/improc.h

@@ -0,0 +1,29 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+#ifndef IMPROC_H__
+#define IMPROC_H__
+
+#include "imagefile.h"
+
+void process_file(Image *I);
+int  process_input(InputType tp, char *name);
+void openXYlog(const char *name);
+void closeXYlog();
+
+#endif // IMPROC_H__

LocCorr/inotify.c

@@ -0,0 +1,129 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <errno.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h> // strlen
+#include <sys/inotify.h>
+#include <sys/select.h>
+#include <unistd.h>
+#include <usefull_macros.h>
+
+#include "imagefile.h"
+
+static char filenm[FILENAME_MAX];
+
+// handle inotify event IN_CLOSE_WRITE
+/**
+ * @brief handle_events - handle inotify event IN_CLOSE_WRITE
+ * @param name  - directory or file name
+ * @param fd    - inotify init descriptor
+ * @return 1 if file was changed and ready to be processed, 0 if not changed, -1 if removed
+ */
+static int changed(const char *name, int fd, uint32_t mask){
+    struct inotify_event buf[10] = {0};
+    ssize_t len = read(fd, buf, sizeof buf);
+    if(len == -1 && errno != EAGAIN){
+        WARN("inotify read()");
+        return -1;
+    }
+    if(len < 1) return 0; // not ready
+    uint32_t bm = buf[0].mask;
+    buf[0].mask = 0;
+    for(int i = 0; i < 10; ++i)
+        DBG("CHANGED: %s (%d)\n", buf[i].name, bm);
+    if(bm & mask){
+        if(name){
+            snprintf(filenm, FILENAME_MAX, "%s/%s", name, buf[0].name); // full path
+        }else{
+            snprintf(filenm, FILENAME_MAX, "%s", buf[0].name); // file name
+        }
+        return 1;
+    }
+    DBG("IN_IGNORED");
+    return -1; // IN_IGNORED (file removed)
+}
+
+static int initinot(const char *name, int *filed, uint32_t mask){
+    static int fd = -1, wd = -1;
+    if(wd > -1){
+        inotify_rm_watch(fd, wd);
+        wd = -1; fd = -1;
+    }
+    if(fd == -1) fd = inotify_init();//1(IN_NONBLOCK);
+    if(fd == -1){
+       WARN("inotify_init()");
+       return 0;
+    }
+    if(-1 == (wd = inotify_add_watch(fd, name, mask))){
+        WARN("inotify_add_watch()");
+        return 0;
+    }
+    if(filed) *filed = fd;
+    return wd;
+}
+
+static int watch_any(const char *name, void (*process)(Image*), uint32_t mask){
+    int fd = -1, wd = -1;
+    while(1){
+        if(wd < 1 || fd < 1){
+            wd = initinot(name, &fd, mask);
+            if(wd < 1){
+                usleep(1000);
+                continue;
+            }
+        }
+        int ch;
+        if(mask == IN_CLOSE_WRITE)
+            ch = changed(NULL, fd, mask); // only file name
+        else
+            ch = changed(name, fd, mask); // full path
+        if(ch == -1){ sleep(1); wd = -1; continue; } // removed
+        if(ch == 1){ // changed
+            if(process){
+                Image *I = Image_read(filenm);
+                process(I);
+            }
+        }
+    }
+    close(fd);
+    return 0;
+}
+
+
+
+int watch_file(const char *name, void (*process)(Image*)){
+    FNAME();
+    if(!name){
+        WARNX("Need filename");
+        return 1;
+    }
+    return watch_any(name, process, IN_CLOSE_WRITE);
+}
+
+int watch_directory(char *name, void (*process)(Image*)){
+    FNAME();
+    if(!name){
+        WARNX("Need directory name");
+        return 1;
+    }
+    int l = strlen(name) - 1;
+    if(name[l] == '/') name[l] = 0;
+    return watch_any(name, process, IN_CLOSE_WRITE|IN_ISDIR);
+}

LocCorr/inotify.h

@@ -0,0 +1,28 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+#ifndef INOTIFY_H__
+#define INOTIFY_H__
+
+#include "fits.h" // Image*
+
+int watch_file(const char *name, void (*process)(Image*));
+int watch_directory(char *name, void (*process)(Image*));
+
+#endif // INOTIFY_H__

LocCorr/main.c

@@ -0,0 +1,117 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <signal.h>         // signal
+#include <string.h>         // strdup
+#include <usefull_macros.h>
+
+#include "cmdlnopts.h"
+#include "improc.h"
+
+/**
+ * We REDEFINE the default WEAK function of signal processing
+ */
+void signals(int sig){
+    DBG("exit %d", sig);
+    if(sig){
+        signal(sig, SIG_IGN);
+        DBG("Get signal %d, quit.\n", sig);
+    }
+    closeXYlog();
+    if(GP && GP->pidfile) // remove unnesessary PID file
+        unlink(GP->pidfile);
+    LOGERR("Exit with status %d", sig);
+    exit(sig);
+}
+
+void iffound_default(pid_t pid){
+    ERRX("Another copy of this process found, pid=%d. Exit.", pid);
+}
+
+static InputType chk_inp(const char *name){
+    if(!name) ERRX("Point file or directory name to monitor");
+    InputType tp = chkinput(GP->inputname);
+    if(T_WRONG == tp) return T_WRONG;
+    green("\n%s is a ", name);
+    switch(tp){
+        case T_DIRECTORY:
+            printf("directory");
+        break;
+        case T_JPEG:
+            printf("jpeg");
+        break;
+        case T_PNG:
+            printf("png");
+        break;
+        case T_GIF:
+            printf("gif");
+        break;
+        case T_FITS:
+            printf("fits");
+        break;
+        case T_BMP:
+            printf("bmp");
+        break;
+        case T_GZIP:
+            printf("maybe fits.gz?");
+        break;
+        case T_CAPT_GRASSHOPPER:
+            printf("capture grasshopper camera");
+        break;
+        default:
+            printf("Unsupported type\n");
+            return T_WRONG;
+    }
+    printf("\n");
+    return tp;
+}
+
+int main(int argc, char *argv[]){
+    initial_setup();
+    char *self = strdup(argv[0]);
+    GP = parse_args(argc, argv);
+    if(GP->throwpart < 0. || GP->throwpart > 0.99){
+        ERRX("Fraction of black pixels should be in [0., 0.99]");
+    }
+    InputType tp = chk_inp(GP->inputname);
+    if(tp == T_WRONG) ERRX("Enter correct image file or directory name");
+    check4running(self, GP->pidfile);
+    DBG("%s started, snippets library version is %s\n", self, sl_libversion());
+    free(self);
+    signal(SIGTERM, signals); // kill (-15) - quit
+    signal(SIGHUP, SIG_IGN);  // hup - ignore
+    signal(SIGINT, signals);  // ctrl+C - quit
+    signal(SIGQUIT, signals); // ctrl+\ - quit
+    signal(SIGTSTP, SIG_IGN); // ignore ctrl+Z
+    if(GP->logfile){
+        sl_loglevel lvl = LOGLEVEL_ERR; // default log level - errors
+        int v = GP->verb;
+        while(v--){ // increase loglevel for each "-v"
+            if(++lvl == LOGLEVEL_ANY) break;
+        }
+        OPENLOG(GP->logfile, lvl, 1);
+        DBG("Opened log file @ level %d", lvl);
+    }
+    if(GP->logXYname) openXYlog(GP->logXYname);
+    LOGMSG("Start application...");
+    LOGDBG("xtag=%g, ytag=%g", GP->xtarget, GP->ytarget);
+    int p = process_input(tp, GP->inputname);
+    // never reached
+    signals(p); // clean everything
+    return p;
+}

LocCorr/median.c

@@ -0,0 +1,519 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// FOR MEDIATOR:
+// Copyright (c) 2011 ashelly.myopenid.com under <http://www.opensource.org/licenses/mit-license>
+
+
+// TODO: resolve problem with borders
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <usefull_macros.h>
+
+#include "imagefile.h"
+#include "median.h"
+
+
+#define ELEM_SWAP(a, b) {register Imtype t = a; a = b; b = t;}
+#define PIX_SORT(a, b)  {if (p[a] > p[b]) ELEM_SWAP(p[a], p[b]);}
+
+static Imtype opt_med2(Imtype *p){
+    return (p[0] + p[1]) * 0.5;
+}
+static Imtype opt_med3(Imtype *p){
+    PIX_SORT(0, 1); PIX_SORT(1, 2); PIX_SORT(0, 1);
+    return(p[1]) ;
+}
+static Imtype opt_med4(Imtype *p){
+    PIX_SORT(0, 2); PIX_SORT(1, 3);
+    PIX_SORT(0, 1); PIX_SORT(2, 3);
+    return(p[1] + p[2]) * 0.5;
+}
+static Imtype opt_med5(Imtype *p){
+    PIX_SORT(0, 1); PIX_SORT(3, 4); PIX_SORT(0, 3);
+    PIX_SORT(1, 4); PIX_SORT(1, 2); PIX_SORT(2, 3) ;
+    PIX_SORT(1, 2);
+    return(p[2]) ;
+}
+// even values are from "FAST, EFFICIENT MEDIAN FILTERS WITH EVEN LENGTH WINDOWS", J.P. HAVLICEK, K.A. SAKADY, G.R.KATZ
+static Imtype opt_med6(Imtype *p){
+    PIX_SORT(1, 2); PIX_SORT(3, 4);
+    PIX_SORT(0, 1); PIX_SORT(2, 3); PIX_SORT(4, 5);
+    PIX_SORT(1, 2); PIX_SORT(3, 4);
+    PIX_SORT(0, 1); PIX_SORT(2, 3); PIX_SORT(4, 5);
+    PIX_SORT(1, 2); PIX_SORT(3, 4);
+    return ( p[2] + p[3] ) * 0.5;
+}
+static Imtype opt_med7(Imtype *p){
+    PIX_SORT(0, 5); PIX_SORT(0, 3); PIX_SORT(1, 6);
+    PIX_SORT(2, 4); PIX_SORT(0, 1); PIX_SORT(3, 5);
+    PIX_SORT(2, 6); PIX_SORT(2, 3); PIX_SORT(3, 6);
+    PIX_SORT(4, 5); PIX_SORT(1, 4); PIX_SORT(1, 3);
+    PIX_SORT(3, 4); return (p[3]);
+}
+// optimal Batcher's sort for 8 elements (http://myopen.googlecode.com/svn/trunk/gtkclient_tdt/include/fast_median.h)
+static Imtype opt_med8(Imtype *p){
+    PIX_SORT(0, 4); PIX_SORT(1, 5); PIX_SORT(2, 6);
+    PIX_SORT(3, 7); PIX_SORT(0, 2); PIX_SORT(1, 3);
+    PIX_SORT(4, 6); PIX_SORT(5, 7); PIX_SORT(2, 4);
+    PIX_SORT(3, 5); PIX_SORT(0, 1); PIX_SORT(2, 3);
+    PIX_SORT(4, 5); PIX_SORT(6, 7); PIX_SORT(1, 4);
+    PIX_SORT(3, 6);
+    return(p[3] + p[4]) * 0.5;
+}
+static Imtype opt_med9(Imtype *p){
+    PIX_SORT(1, 2); PIX_SORT(4, 5); PIX_SORT(7, 8);
+    PIX_SORT(0, 1); PIX_SORT(3, 4); PIX_SORT(6, 7);
+    PIX_SORT(1, 2); PIX_SORT(4, 5); PIX_SORT(7, 8);
+    PIX_SORT(0, 3); PIX_SORT(5, 8); PIX_SORT(4, 7);
+    PIX_SORT(3, 6); PIX_SORT(1, 4); PIX_SORT(2, 5);
+    PIX_SORT(4, 7); PIX_SORT(4, 2); PIX_SORT(6, 4);
+    PIX_SORT(4, 2); return(p[4]);
+}
+static Imtype opt_med16(Imtype *p){
+    PIX_SORT(0, 8); PIX_SORT(1, 9); PIX_SORT(2, 10); PIX_SORT(3, 11);
+    PIX_SORT(4, 12); PIX_SORT(5, 13); PIX_SORT(6, 14); PIX_SORT(7, 15);
+    PIX_SORT(0, 4); PIX_SORT(1, 5); PIX_SORT(2, 6); PIX_SORT(3, 7);
+    PIX_SORT(8, 12); PIX_SORT(9, 13); PIX_SORT(10, 14); PIX_SORT(11, 15);
+    PIX_SORT(4, 8); PIX_SORT(5, 9); PIX_SORT(6, 10); PIX_SORT(7, 11);
+    PIX_SORT(0, 2); PIX_SORT(1, 3); PIX_SORT(4, 6); PIX_SORT(5, 7);
+    PIX_SORT(8, 10); PIX_SORT(9, 11); PIX_SORT(12, 14); PIX_SORT(13, 15);
+    PIX_SORT(2, 8); PIX_SORT(3, 9); PIX_SORT(6, 12); PIX_SORT(7, 13);
+    PIX_SORT(2, 4); PIX_SORT(3, 5); PIX_SORT(6, 8); PIX_SORT(7, 9);
+    PIX_SORT(10, 12); PIX_SORT(11, 13); PIX_SORT(0, 1); PIX_SORT(2, 3);
+    PIX_SORT(4, 5); PIX_SORT(6, 7); PIX_SORT(8, 9); PIX_SORT(10, 11);
+    PIX_SORT(12, 13); PIX_SORT(14, 15); PIX_SORT(1, 8); PIX_SORT(3, 10);
+    PIX_SORT(5, 12); PIX_SORT(7, 14); PIX_SORT(5, 8); PIX_SORT(7, 10);
+    return (p[7] + p[8]) * 0.5;
+}
+static Imtype opt_med25(Imtype *p){
+    PIX_SORT(0, 1)  ; PIX_SORT(3, 4)  ; PIX_SORT(2, 4) ;
+    PIX_SORT(2, 3)  ; PIX_SORT(6, 7)  ; PIX_SORT(5, 7) ;
+    PIX_SORT(5, 6)  ; PIX_SORT(9, 10) ; PIX_SORT(8, 10) ;
+    PIX_SORT(8, 9)  ; PIX_SORT(12, 13); PIX_SORT(11, 13) ;
+    PIX_SORT(11, 12); PIX_SORT(15, 16); PIX_SORT(14, 16) ;
+    PIX_SORT(14, 15); PIX_SORT(18, 19); PIX_SORT(17, 19) ;
+    PIX_SORT(17, 18); PIX_SORT(21, 22); PIX_SORT(20, 22) ;
+    PIX_SORT(20, 21); PIX_SORT(23, 24); PIX_SORT(2, 5) ;
+    PIX_SORT(3, 6)  ; PIX_SORT(0, 6)  ; PIX_SORT(0, 3) ;
+    PIX_SORT(4, 7)  ; PIX_SORT(1, 7)  ; PIX_SORT(1, 4) ;
+    PIX_SORT(11, 14); PIX_SORT(8, 14) ; PIX_SORT(8, 11) ;
+    PIX_SORT(12, 15); PIX_SORT(9, 15) ; PIX_SORT(9, 12) ;
+    PIX_SORT(13, 16); PIX_SORT(10, 16); PIX_SORT(10, 13) ;
+    PIX_SORT(20, 23); PIX_SORT(17, 23); PIX_SORT(17, 20) ;
+    PIX_SORT(21, 24); PIX_SORT(18, 24); PIX_SORT(18, 21) ;
+    PIX_SORT(19, 22); PIX_SORT(8, 17) ; PIX_SORT(9, 18) ;
+    PIX_SORT(0, 18) ; PIX_SORT(0, 9)  ; PIX_SORT(10, 19) ;
+    PIX_SORT(1, 19) ; PIX_SORT(1, 10) ; PIX_SORT(11, 20) ;
+    PIX_SORT(2, 20) ; PIX_SORT(2, 11) ; PIX_SORT(12, 21) ;
+    PIX_SORT(3, 21) ; PIX_SORT(3, 12) ; PIX_SORT(13, 22) ;
+    PIX_SORT(4, 22) ; PIX_SORT(4, 13) ; PIX_SORT(14, 23) ;
+    PIX_SORT(5, 23) ; PIX_SORT(5, 14) ; PIX_SORT(15, 24) ;
+    PIX_SORT(6, 24) ; PIX_SORT(6, 15) ; PIX_SORT(7, 16) ;
+    PIX_SORT(7, 19) ; PIX_SORT(13, 21); PIX_SORT(15, 23) ;
+    PIX_SORT(7, 13) ; PIX_SORT(7, 15) ; PIX_SORT(1, 9) ;
+    PIX_SORT(3, 11) ; PIX_SORT(5, 17) ; PIX_SORT(11, 17) ;
+    PIX_SORT(9, 17) ; PIX_SORT(4, 10) ; PIX_SORT(6, 12) ;
+    PIX_SORT(7, 14) ; PIX_SORT(4, 6)  ; PIX_SORT(4, 7) ;
+    PIX_SORT(12, 14); PIX_SORT(10, 14); PIX_SORT(6, 7) ;
+    PIX_SORT(10, 12); PIX_SORT(6, 10) ; PIX_SORT(6, 17) ;
+    PIX_SORT(12, 17); PIX_SORT(7, 17) ; PIX_SORT(7, 10) ;
+    PIX_SORT(12, 18); PIX_SORT(7, 12) ; PIX_SORT(10, 18) ;
+    PIX_SORT(12, 20); PIX_SORT(10, 20); PIX_SORT(10, 12) ;
+    return (p[12]);
+}
+#undef PIX_SORT
+#define PIX_SORT(a, b)  {if (a > b) ELEM_SWAP(a, b);}
+/**
+ * quick select - algo for approximate median calculation for array idata of size n
+ */
+static Imtype quick_select(Imtype *idata, int n){
+    int low, high;
+    int median;
+    int middle, ll, hh;
+    Imtype *arr = MALLOC(Imtype, n);
+    memcpy(arr, idata, n*sizeof(Imtype));
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for(;;){
+        if(high <= low) // One element only
+            break;
+        if(high == low + 1){ // Two elements only
+            PIX_SORT(arr[low], arr[high]) ;
+            break;
+        }
+        // Find median of low, middle and high Imtypes; swap into position low
+        middle = (low + high) / 2;
+        PIX_SORT(arr[middle], arr[high]) ;
+        PIX_SORT(arr[low], arr[high]) ;
+        PIX_SORT(arr[middle], arr[low]) ;
+        // Swap low Imtype (now in position middle) into position (low+1)
+        ELEM_SWAP(arr[middle], arr[low+1]) ;
+        // Nibble from each end towards middle, swapping Imtypes when stuck
+        ll = low + 1;
+        hh = high;
+        for(;;){
+            do ll++; while (arr[low] > arr[ll]);
+            do hh--; while (arr[hh] > arr[low]);
+            if(hh < ll) break;
+            ELEM_SWAP(arr[ll], arr[hh]) ;
+        }
+        // Swap middle Imtype (in position low) back into correct position
+        ELEM_SWAP(arr[low], arr[hh]) ;
+        // Re-set active partition
+        if (hh <= median) low = ll;
+        if (hh >= median) high = hh - 1;
+    }
+    Imtype ret = arr[median];
+    FREE(arr);
+    return ret;
+}
+#undef PIX_SORT
+#undef ELEM_SWAP
+
+/**
+ * calculate median of array idata with size n
+ */
+Imtype calc_median(Imtype *idata, int n){
+    if(!idata || n < 1){
+        WARNX("calc_median(): wrong dta"); return 0.;
+    }
+    typedef Imtype (*medfunc)(Imtype *p);
+    medfunc fn = NULL;
+    const medfunc fnarr[] = {opt_med2, opt_med3, opt_med4, opt_med5, opt_med6,
+            opt_med7, opt_med8, opt_med9};
+    if(n == 1) return *idata;
+    if(n < 10) fn = fnarr[n - 1];
+    else if(n == 16) fn = opt_med16;
+    else if(n == 25) fn = opt_med25;
+    if(fn){
+        return fn(idata);
+    }else{
+        return quick_select(idata, n);
+    }
+}
+
+#define ImtypeLess(a,b) ((a)<(b))
+#define ImtypeMean(a,b) (((a)+(b))/2)
+
+typedef struct Mediator_t{
+    Imtype* data;   // circular queue of values
+    int* pos;       // index into `heap` for each value
+    int* heap;      // max/median/min heap holding indexes into `data`.
+    int N;          // allocated size.
+    int idx;        // position in circular queue
+    int ct;         // count of Imtypes in queue
+} Mediator;
+
+/*--- Helper Functions ---*/
+
+#define minCt(m) (((m)->ct-1)/2) //count of Imtypes in minheap
+#define maxCt(m) (((m)->ct)/2) //count of Imtypes in maxheap
+
+//returns 1 if heap[i] < heap[j]
+inline int mmless(Mediator* m, int i, int j){
+    return ImtypeLess(m->data[m->heap[i]],m->data[m->heap[j]]);
+}
+
+//swaps Imtypes i&j in heap, maintains indexes
+inline int mmexchange(Mediator* m, int i, int j){
+    int t = m->heap[i];
+    m->heap[i] = m->heap[j];
+    m->heap[j] = t;
+    m->pos[m->heap[i]] = i;
+    m->pos[m->heap[j]] = j;
+    return 1;
+}
+
+//swaps Imtypes i&j if i<j; returns true if swapped
+inline int mmCmpExch(Mediator* m, int i, int j){
+    return (mmless(m,i,j) && mmexchange(m,i,j));
+}
+
+//maintains minheap property for all Imtypes below i/2.
+void minSortDown(Mediator* m, int i){
+    for(; i <= minCt(m); i*=2){
+        if(i>1 && i < minCt(m) && mmless(m, i+1, i)) ++i;
+        if(!mmCmpExch(m,i,i/2)) break;
+    }
+}
+
+//maintains maxheap property for all Imtypes below i/2. (negative indexes)
+void maxSortDown(Mediator* m, int i){
+    for(; i >= -maxCt(m); i*=2){
+        if(i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) --i;
+    if(!mmCmpExch(m,i/2,i)) break;
+    }
+}
+
+//maintains minheap property for all Imtypes above i, including median
+//returns true if median changed
+int minSortUp(Mediator* m, int i){
+    while (i > 0 && mmCmpExch(m, i, i/2)) i /= 2;
+    return (i == 0);
+}
+
+//maintains maxheap property for all Imtypes above i, including median
+//returns true if median changed
+int maxSortUp(Mediator* m, int i){
+    while (i < 0 && mmCmpExch(m, i/2, i)) i /= 2;
+    return (i == 0);
+}
+
+/*--- Public Interface ---*/
+
+//creates new Mediator: to calculate `nImtypes` running median.
+//mallocs single block of memory, caller must free.
+Mediator* MediatorNew(int nImtypes){
+    int size = sizeof(Mediator) + nImtypes*(sizeof(Imtype)+sizeof(int)*2);
+    Mediator* m = malloc(size);
+    m->data = (Imtype*)(m + 1);
+    m->pos = (int*) (m->data + nImtypes);
+    m->heap = m->pos + nImtypes + (nImtypes / 2); //points to middle of storage.
+    m->N = nImtypes;
+    m->ct = m->idx = 0;
+    while (nImtypes--){ //set up initial heap fill pattern: median,max,min,max,...
+        m->pos[nImtypes] = ((nImtypes+1)/2) * ((nImtypes&1)? -1 : 1);
+        m->heap[m->pos[nImtypes]] = nImtypes;
+    }
+    return m;
+}
+
+//Inserts Imtype, maintains median in O(lg nImtypes)
+void MediatorInsert(Mediator* m, Imtype v){
+    int isNew=(m->ct<m->N);
+    int p = m->pos[m->idx];
+    Imtype old = m->data[m->idx];
+    m->data[m->idx]=v;
+    m->idx = (m->idx+1) % m->N;
+    m->ct+=isNew;
+    if(p>0){ //new Imtype is in minHeap
+        if (!isNew && ImtypeLess(old,v)) minSortDown(m,p*2);
+        else if (minSortUp(m,p)) maxSortDown(m,-1);
+    }else if (p<0){ //new Imtype is in maxheap
+        if (!isNew && ImtypeLess(v,old)) maxSortDown(m,p*2);
+        else if (maxSortUp(m,p)) minSortDown(m, 1);
+    }else{ //new Imtype is at median
+        if (maxCt(m)) maxSortDown(m,-1);
+        if (minCt(m)) minSortDown(m, 1);
+    }
+}
+
+//returns median Imtype (or average of 2 when Imtype count is even)
+Imtype MediatorMedian(Mediator* m){
+    Imtype v = m->data[m->heap[0]];
+    if ((m->ct&1) == 0) v = ImtypeMean(v, m->data[m->heap[-1]]);
+    return v;
+}
+
+// median + min/max
+Imtype MediatorStat(Mediator* m, Imtype *minval, Imtype *maxval){
+    Imtype v= m->data[m->heap[0]];
+    if ((m->ct&1) == 0) v = ImtypeMean(v, m->data[m->heap[-1]]);
+    Imtype min = v, max = v;
+    int i;
+    for(i = -maxCt(m); i < 0; ++i){
+        int v = m->data[m->heap[i]];
+        if(v < min) min = v;
+    }
+    *minval = min;
+    for(i = 1; i <= minCt(m); ++i){
+        int v = m->data[m->heap[i]];
+        if(v > max) max = v;
+    }
+    *maxval = max;
+    return v;
+}
+
+/**
+ * filter image by median (seed*2 + 1) x (seed*2 + 1)
+ */
+Image *get_median(const Image *img, int seed){
+    if(seed < 1) return NULL;
+    size_t w = img->width, h = img->height;
+    Image *out = Image_sim(img);
+    Imtype *med = out->data, *inputima = img->data;
+
+    size_t blksz = seed * 2 + 1, fullsz = blksz * blksz;
+#ifdef EBUG
+    double t0 = dtime();
+#endif
+    OMP_FOR(shared(inputima, med))
+    for(size_t x = seed; x < w - seed; ++x){
+        size_t xx, yy, xm = x + seed + 1, y, ymax = blksz - 1, xmin = x - seed;
+        Mediator* m = MediatorNew(fullsz);
+        // initial fill
+        for(yy = 0; yy < ymax; ++yy)
+            for(xx = xmin; xx < xm; ++xx)
+                MediatorInsert(m, inputima[xx + yy*w]);
+        ymax = 2*seed*w;
+        xmin += ymax;
+        xm += ymax;
+        ymax = h - seed;
+        size_t medidx = x + seed * w;
+        for(y = seed; y < ymax; ++y, xmin += w, xm += w, medidx += w){
+            for(xx = xmin; xx < xm; ++xx)
+                MediatorInsert(m, inputima[xx]);
+            med[medidx] = MediatorMedian(m);
+        }
+        FREE(m);
+    }
+    Image_minmax(out);
+    DBG("time for median filtering %zdx%zd of image %zdx%zd: %gs", blksz, blksz, w, h,
+        dtime() - t0);
+    return out;
+}
+
+/**
+ * @brief get_stat - calculate floating statistics in (seed*2+1)^2
+ * @param in (i) - input image
+ * @param seed - radius of box
+ * @param mean (o) - mean by box (excluding borders)
+ * @param std  (o) - STD by box (excluding borders)
+ * @retur 0 if error
+ */
+int get_stat(const Image *in, int seed, Image **mean, Image **std){
+    if(!in) return 0;
+    if(seed < 1 || seed > (in->width - 1)/2 || seed > (in->height - 1)/2) return 0;
+#ifdef EBUG
+    double t0 = dtime();
+#endif
+    Image *M = NULL, *S = NULL;
+    if(mean) M = Image_sim(in);
+    if(std) S = Image_sim(in);
+    int ymax = in->height - seed, xmax = in->width - seed;
+    int hsz = (seed*2 + 1), sz = hsz * hsz, w = in->width;
+    OMP_FOR()
+    for(int y = seed; y < ymax; ++y){ // dumb calculations
+        int startidx = y*w + seed;
+        Imtype *om = (M) ? &M->data[startidx] : NULL;
+        Imtype *os = (S) ? &S->data[startidx] : NULL;
+        for(int x = seed; x < xmax; ++x){
+            Imtype sum = 0, sum2 = 0;
+            int yb = y + seed + 1, xm = x - seed;
+            for(int yy = y - seed; yy < yb; ++yy){
+                Imtype *ptr = &in->data[yy * w + xm];
+                for(int xx = 0; xx < hsz; ++xx){
+                    Imtype d = *ptr++;
+                    sum += d;
+                    sum2 += d*d;
+                }
+            }
+            //DBG("sum=%g, sum2=%g, sz=%d", sum, sum2, sz);
+            sum /= sz;
+            if(om){
+                *om++ = sum;
+                //DBG("mean (%d, %d): %g", x, y, sum);
+            }
+            if(os) *os++ = sqrt(sum2/sz - sum*sum);
+        }
+    }
+    if(mean){
+        Image_minmax(M);
+        *mean = M;
+    }
+    if(std){
+        Image_minmax(S);
+        *std = S;
+    }
+    DBG("time for mean/sigma computation: %gs", dtime() - t0);
+    return 1;
+}
+
+/**
+ * @brief calc_background - Simple background calculation by histogram
+ * @param img (i) - input image (here will be modified its top2proc field)
+ * @param bk (o)  - background value
+ * @return 0 if error
+ */
+int calc_background(Image *img, Imtype *bk){
+    if(img->maxval - img->minval < DBL_EPSILON){
+        WARNX("Zero image!");
+        return 0;
+    }
+    int w = img->width, h = img->height, wh = w*h;
+    Imtype min = img->minval, ampl = img->maxval - min;
+    int histogram[256] = {0};
+    //DBG("min: %g, max: %g, ampl: %g", min, img->maxval, ampl);
+    #pragma omp parallel
+    {
+        int histogram_private[256] = {0};
+        #pragma omp for nowait
+        for(int i = 0; i < wh; ++i){
+            int newval = (int)((((img->data[i]) - min)/ampl)*255. + 0.5);
+            ++histogram_private[newval];
+        }
+        #pragma omp critical
+        {
+            for(int i=0; i<256; ++i) histogram[i] += histogram_private[i];
+        }
+    }
+    int modeidx = 0, modeval = 0;
+    for(int i = 0; i < 256; ++i)
+        if(modeval < histogram[i]){
+            modeval = histogram[i];
+            modeidx = i;
+        }
+    //DBG("Mode=%g @ idx%d (N=%d)", ((Imtype)modeidx / 255.)*ampl, modeidx, modeval);
+/*
+    int diff[256] = {0};
+    for(int i = 1; i < 255; ++i) diff[i] = (histogram[i+1]-histogram[i-1])/2;
+    if(modeidx == 0) modeidx = 1;
+    if(modeidx > 253) return NULL; // very bad image: overilluminated
+    int borderidx = modeidx;
+    green("1\n");
+    for(int i = modeidx; i < 255; ++i){ // search bend-point by first derivate
+        printf("%d: %d, %d\n", i, diff[i], diff[i+1]);
+        if(diff[i] >= 0 && diff[i+1] >=0){
+            borderidx = i; break;
+        }
+    }
+*/
+    int diff2[256] = {0};
+    for(int i = 2; i < 254; ++i) diff2[i] = (histogram[i+2]+histogram[i-2]-2*histogram[i])/4;
+    if(modeidx < 2) modeidx = 2;
+    if(modeidx > 253) return 0; // very bad image: overilluminated
+    int borderidx = modeidx;
+   // green("2\n");
+    for(int i = modeidx; i < 254; ++i){ // search bend-point by second derivate
+     //   printf("%d: %d, %d\n", i, diff2[i], diff2[i+1]);
+        if(diff2[i] <= 0 && diff2[i+1] <=0){
+            borderidx = i; break;
+        }
+    }
+    //DBG("borderidx=%d -> %d", borderidx, (borderidx+modeidx)/2);
+    //borderidx = (borderidx + modeidx) / 2;
+    Imtype borderval = ((Imtype)borderidx / 255.)*ampl + min;
+    if(bk) *bk = borderval;
+    //green("HISTO:\n");
+    //for(int i = 0; i < 256; ++i) printf("%d:\t%d\t%d\n", i, histogram[i], diff2[i]);
+    // calculate values of upper 2% border
+#if 0
+    Image *out = Image_sim(img);
+    //DBG("found border: %g @ %d", borderval, borderidx);
+    OMP_FOR()
+    for(int i = 0; i < wh; ++i){
+        register Imtype val = img->data[i];
+        if(val > borderval) out->data[i] = val - borderval;
+    }
+    Image_minmax(out);
+#endif
+    return 1;
+}

LocCorr/median.h

@@ -0,0 +1,32 @@
+/*
+ * median.h
+ *
+ * Copyright 2015 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#pragma once
+#ifndef __MEDIAN_H__
+#define __MEDIAN_H__
+
+#include "fits.h"
+
+Imtype calc_median(Imtype *idata, int n);
+Image *get_median(const Image *img, int seed);
+int  get_stat(const Image *in, int seed, Image **mean, Image **std);
+int calc_background(Image *img, Imtype *bk);
+
+#endif // __MEDIAN_H__

LocCorr/stbimpl.c

@@ -0,0 +1,29 @@
+/*
+ * This file is part of the loccorr project.
+ * Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// implementation of STB library (to accelerate compilation)
+#define STBI_NO_PSD
+#define STBI_NO_TGA
+#define STBI_NO_HDR
+#define STBI_NO_PIC
+#define STBI_NO_PNM
+
+#define STB_IMAGE_IMPLEMENTATION
+#include <stb/stb_image.h>
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+#include <stb/stb_image_write.h>