diff --git a/LocCorr/CMakeLists.txt b/LocCorr/CMakeLists.txt new file mode 100644 index 0000000..d978d0e --- /dev/null +++ b/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") diff --git a/LocCorr/FindFLYCAP.cmake b/LocCorr/FindFLYCAP.cmake new file mode 100644 index 0000000..6eafdc1 --- /dev/null +++ b/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 +# +# 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) diff --git a/LocCorr/Readme.md b/LocCorr/Readme.md new file mode 100644 index 0000000..166f6bb --- /dev/null +++ b/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. diff --git a/LocCorr/binmorph.c b/LocCorr/binmorph.c new file mode 100644 index 0000000..3541415 --- /dev/null +++ b/LocCorr/binmorph.c @@ -0,0 +1,571 @@ +/* + * binmorph.c - functions for morphological operations on binary image + * + * Copyright 2015 Edward V. Emelianoff + * + * 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 +#include +#include // memcpy +#include +#include +#include +#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< + */ +#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 ===================> + */ diff --git a/LocCorr/binmorph.h b/LocCorr/binmorph.h new file mode 100644 index 0000000..933fae0 --- /dev/null +++ b/LocCorr/binmorph.h @@ -0,0 +1,85 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#pragma once +#ifndef BINMORPH_H__ +#define BINMORPH_H__ + +#include +#include +#include +#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__ diff --git a/LocCorr/cclabling.h b/LocCorr/cclabling.h new file mode 100644 index 0000000..e01f3c7 --- /dev/null +++ b/LocCorr/cclabling.h @@ -0,0 +1,136 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +//#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); diff --git a/LocCorr/cclabling_ori.h b/LocCorr/cclabling_ori.h new file mode 100644 index 0000000..76fe3bb --- /dev/null +++ b/LocCorr/cclabling_ori.h @@ -0,0 +1,128 @@ +/* + * cclabling_1.h - inner part of functions cclabel4 and cclabel8 + * + * Copyright 2015 Edward V. Emelianoff + * + * 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 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); diff --git a/LocCorr/cmdlnopts.c b/LocCorr/cmdlnopts.c new file mode 100644 index 0000000..8af09e5 --- /dev/null +++ b/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 + * + * 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 +#include +#include +#include +#include +#include + +#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; +} + diff --git a/LocCorr/cmdlnopts.h b/LocCorr/cmdlnopts.h new file mode 100644 index 0000000..5c78b8a --- /dev/null +++ b/LocCorr/cmdlnopts.h @@ -0,0 +1,53 @@ +/* geany_encoding=koi8-r + * cmdlnopts.h - comand line options for parceargs + * + * Copyright 2013 Edward V. Emelianoff + * + * 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__ diff --git a/LocCorr/dilation.h b/LocCorr/dilation.h new file mode 100644 index 0000000..7fdba20 --- /dev/null +++ b/LocCorr/dilation.h @@ -0,0 +1,73 @@ +/* + * dilation.h - inner part of function `dilation` + * + * Copyright 2013 Edward V. Emelianoff + * + * 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 +} + diff --git a/LocCorr/draw.c b/LocCorr/draw.c new file mode 100644 index 0000000..ce95a20 --- /dev/null +++ b/LocCorr/draw.c @@ -0,0 +1,105 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +// simplest interface to draw lines & ellipsis +#include "draw.h" +#include "fits.h" + +#include + +// 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)); + } + } + } +} diff --git a/LocCorr/draw.h b/LocCorr/draw.h new file mode 100644 index 0000000..cbea4e7 --- /dev/null +++ b/LocCorr/draw.h @@ -0,0 +1,56 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ +#pragma once +#ifndef DRAW_H__ +#define DRAW_H__ + +#include + +// 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__ diff --git a/LocCorr/ec_filter.h b/LocCorr/ec_filter.h new file mode 100644 index 0000000..c95adec --- /dev/null +++ b/LocCorr/ec_filter.h @@ -0,0 +1,96 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +// 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; +} diff --git a/LocCorr/fc_filter.h b/LocCorr/fc_filter.h new file mode 100644 index 0000000..5d30d3c --- /dev/null +++ b/LocCorr/fc_filter.h @@ -0,0 +1,68 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +// 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; +} diff --git a/LocCorr/fits.c b/LocCorr/fits.c new file mode 100644 index 0000000..3e06ec7 --- /dev/null +++ b/LocCorr/fits.c @@ -0,0 +1,227 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ +#include +#include +#include +#include +#include + +#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; +} diff --git a/LocCorr/fits.h b/LocCorr/fits.h new file mode 100644 index 0000000..865c49d --- /dev/null +++ b/LocCorr/fits.h @@ -0,0 +1,55 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#pragma once +#ifndef FITS_H__ +#define FITS_H__ + +#include +#include +#include +#include + +#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__ diff --git a/LocCorr/grasshopper.c b/LocCorr/grasshopper.c new file mode 100644 index 0000000..c4354ac --- /dev/null +++ b/LocCorr/grasshopper.c @@ -0,0 +1,320 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#include +#include +#include // FLT_EPSILON +#include +#include +#include + +#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; +} diff --git a/LocCorr/grasshopper.h b/LocCorr/grasshopper.h new file mode 100644 index 0000000..de10feb --- /dev/null +++ b/LocCorr/grasshopper.h @@ -0,0 +1,28 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ +#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__ diff --git a/LocCorr/imagefile.c b/LocCorr/imagefile.c new file mode 100644 index 0000000..d2dde61 --- /dev/null +++ b/LocCorr/imagefile.c @@ -0,0 +1,514 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#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<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 =================== + */ diff --git a/LocCorr/imagefile.h b/LocCorr/imagefile.h new file mode 100644 index 0000000..07c3696 --- /dev/null +++ b/LocCorr/imagefile.h @@ -0,0 +1,55 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ +#pragma once +#ifndef IMAGEFILE_H__ +#define IMAGEFILE_H__ + +#include + +#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__ diff --git a/LocCorr/improc.c b/LocCorr/improc.c new file mode 100644 index 0000000..0805e20 --- /dev/null +++ b/LocCorr/improc.c @@ -0,0 +1,303 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#include +#include +#include +#include +#include +#include + +#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; +} diff --git a/LocCorr/improc.h b/LocCorr/improc.h new file mode 100644 index 0000000..f106056 --- /dev/null +++ b/LocCorr/improc.h @@ -0,0 +1,29 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ +#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__ diff --git a/LocCorr/inotify.c b/LocCorr/inotify.c new file mode 100644 index 0000000..542cb13 --- /dev/null +++ b/LocCorr/inotify.c @@ -0,0 +1,129 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#include +#include +#include +#include // strlen +#include +#include +#include +#include + +#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); +} diff --git a/LocCorr/inotify.h b/LocCorr/inotify.h new file mode 100644 index 0000000..19caecc --- /dev/null +++ b/LocCorr/inotify.h @@ -0,0 +1,28 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#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__ diff --git a/LocCorr/main.c b/LocCorr/main.c new file mode 100644 index 0000000..0a5d340 --- /dev/null +++ b/LocCorr/main.c @@ -0,0 +1,117 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +#include // signal +#include // strdup +#include + +#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; +} diff --git a/LocCorr/median.c b/LocCorr/median.c new file mode 100644 index 0000000..1dbf11e --- /dev/null +++ b/LocCorr/median.c @@ -0,0 +1,519 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +// FOR MEDIATOR: +// Copyright (c) 2011 ashelly.myopenid.com under + + +// TODO: resolve problem with borders + +#include +#include +#include +#include +#include + +#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 i1 && 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->ctN); + 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; +} diff --git a/LocCorr/median.h b/LocCorr/median.h new file mode 100644 index 0000000..92fe1a5 --- /dev/null +++ b/LocCorr/median.h @@ -0,0 +1,32 @@ +/* + * median.h + * + * Copyright 2015 Edward V. Emelianov + * + * 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__ diff --git a/LocCorr/stbimpl.c b/LocCorr/stbimpl.c new file mode 100644 index 0000000..4c76251 --- /dev/null +++ b/LocCorr/stbimpl.c @@ -0,0 +1,29 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * 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 . + */ + +// 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 +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include