From ae3f52db7041ef950d4ce3bc2d151b780baf6875 Mon Sep 17 00:00:00 2001 From: Edward Emelianov Date: Mon, 29 Mar 2021 18:27:24 +0300 Subject: [PATCH] initial commit --- .gitignore | 23 ++ CMakeLists.txt | 74 +++++++ binmorph.c | 571 +++++++++++++++++++++++++++++++++++++++++++++++++ binmorph.h | 85 ++++++++ cclabling.h | 136 ++++++++++++ cmdlnopts.c | 94 ++++++++ cmdlnopts.h | 47 ++++ dilation.h | 73 +++++++ draw.c | 105 +++++++++ draw.h | 56 +++++ ec_filter.h | 96 +++++++++ fc_filter.h | 68 ++++++ fits.c | 227 ++++++++++++++++++++ fits.h | 55 +++++ imagefile.c | 504 +++++++++++++++++++++++++++++++++++++++++++ imagefile.h | 53 +++++ inotify.c | 125 +++++++++++ inotify.h | 26 +++ main.c | 345 ++++++++++++++++++++++++++++++ median.c | 518 ++++++++++++++++++++++++++++++++++++++++++++ median.h | 32 +++ stbimpl.c | 29 +++ 22 files changed, 3342 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 binmorph.c create mode 100644 binmorph.h create mode 100644 cclabling.h create mode 100644 cmdlnopts.c create mode 100644 cmdlnopts.h create mode 100644 dilation.h create mode 100644 draw.c create mode 100644 draw.h create mode 100644 ec_filter.h create mode 100644 fc_filter.h create mode 100644 fits.c create mode 100644 fits.h create mode 100644 imagefile.c create mode 100644 imagefile.h create mode 100644 inotify.c create mode 100644 inotify.h create mode 100644 main.c create mode 100644 median.c create mode 100644 median.h create mode 100644 stbimpl.c diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..60663b5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,23 @@ +# Prerequisites +*.d + +# Object files +*.o + +# Libraries +*.lib +*.a +*.la +*.lo + +# Shared objects (inc. Windows DLLs) +*.so +*.so.* + +# qt-creator +*.config +*.cflags +*.cxxflags +*.creator* +*.files +*.includes diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..80a9da1 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,74 @@ +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: +find_package(PkgConfig 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}) +# -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} -lm) + +# Installation of the program +INSTALL(TARGETS ${PROJ} DESTINATION "bin") diff --git a/binmorph.c b/binmorph.c new file mode 100644 index 0000000..3541415 --- /dev/null +++ b/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/binmorph.h b/binmorph.h new file mode 100644 index 0000000..933fae0 --- /dev/null +++ b/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/cclabling.h b/cclabling.h new file mode 100644 index 0000000..e01f3c7 --- /dev/null +++ b/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/cmdlnopts.c b/cmdlnopts.c new file mode 100644 index 0000000..55a5237 --- /dev/null +++ b/cmdlnopts.c @@ -0,0 +1,94 @@ +/* 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 +}; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +static myoption cmdlnopts[] = { +// common options + {"help", NO_ARGS, NULL, 'h', 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")}, + {"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",NEED_ARG, 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")}, + 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/cmdlnopts.h b/cmdlnopts.h new file mode 100644 index 0000000..b300568 --- /dev/null +++ b/cmdlnopts.h @@ -0,0 +1,47 @@ +/* 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 + 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 +} 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/dilation.h b/dilation.h new file mode 100644 index 0000000..7fdba20 --- /dev/null +++ b/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/draw.c b/draw.c new file mode 100644 index 0000000..ce95a20 --- /dev/null +++ b/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/draw.h b/draw.h new file mode 100644 index 0000000..cbea4e7 --- /dev/null +++ b/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/ec_filter.h b/ec_filter.h new file mode 100644 index 0000000..c95adec --- /dev/null +++ b/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/fc_filter.h b/fc_filter.h new file mode 100644 index 0000000..5d30d3c --- /dev/null +++ b/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/fits.c b/fits.c new file mode 100644 index 0000000..3e06ec7 --- /dev/null +++ b/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/fits.h b/fits.h new file mode 100644 index 0000000..865c49d --- /dev/null +++ b/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/imagefile.c b/imagefile.c new file mode 100644 index 0000000..90422fd --- /dev/null +++ b/imagefile.c @@ -0,0 +1,504 @@ +/* + * 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 "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){ + 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; +} + +// 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; + } + Image *outp = MALLOC(Image, 1); + int wh = width * height; + outp->width = width; + outp->height = height; + outp->dtype = USHORT_IMG; // let it be 16-bit @ save + uint8_t min = *img, max = min; + for(int i = 0; i < wh; ++i){ + if(img[i] < min) min = img[i]; + else if(img[i] > max) max = img[i]; + } + outp->minval = (Imtype) min; + outp->maxval = (Imtype) max; + outp->data = MALLOC(Imtype, wh); + // 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 = &img[y*width]; + for(int x = 0; x < width; ++x){ + *Out++ = (Imtype)(*In++); + } + } + FREE(img); + return outp; +} + +/** + * @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 + * @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 + * @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/imagefile.h b/imagefile.h new file mode 100644 index 0000000..75887a2 --- /dev/null +++ b/imagefile.h @@ -0,0 +1,53 @@ +/* + * 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, +} 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 *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/inotify.c b/inotify.c new file mode 100644 index 0000000..c02f8a7 --- /dev/null +++ b/inotify.c @@ -0,0 +1,125 @@ +/* + * 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 + +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)(const char*), 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) process(filenm); + } + } + close(fd); + return 1; +} + + + +int watch_file(const char *name, void (*process)(const char*)){ + FNAME(); + if(!name){ + WARNX("Need filename"); + return 0; + } + return watch_any(name, process, IN_CLOSE_WRITE); +} + +int watch_directory(char *name, void (*process)(const char*)){ + FNAME(); + if(!name){ + WARNX("Need directory name"); + return 0; + } + int l = strlen(name) - 1; + if(name[l] == '/') name[l] = 0; + return watch_any(name, process, IN_CLOSE_WRITE|IN_ISDIR); + //return watch_any(name, process, IN_CREATE); +} diff --git a/inotify.h b/inotify.h new file mode 100644 index 0000000..577db95 --- /dev/null +++ b/inotify.h @@ -0,0 +1,26 @@ +/* + * 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__ + +int watch_file(const char *name, void (*process)(const char *n)); +int watch_directory(char *name, void (*process)(const char*)); + +#endif // INOTIFY_H__ diff --git a/main.c b/main.c new file mode 100644 index 0000000..af084ef --- /dev/null +++ b/main.c @@ -0,0 +1,345 @@ +/* + * 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 // signal +#include +#include // printf +#include // exit, free +#include // strdup +#include +#include // sleep +#include + +#include "binmorph.h" +#include "cmdlnopts.h" +#include "draw.h" +#include "inotify.h" +#include "fits.h" +#include "imagefile.h" +#include "median.h" + +/** + * We REDEFINE the default WEAK function of signal processing + */ +void signals(int sig){ + if(sig){ + signal(sig, SIG_IGN); + DBG("Get signal %d, quit.\n", sig); + } + LOGERR("Exit with status %d", sig); + if(GP && GP->pidfile) // remove unnesessary PID file + unlink(GP->pidfile); + 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; + default: + printf("Unsupported type\n"); + return T_WRONG; + } + printf("\n"); + return tp; +} + +static bool save_fits(Image *I, const char *name){ + 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); + } +} +//++npoints, b->area, Isum, wh, xc, yc, x2c, y2c +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; double yc;// centroid coordinates + double xsigma; // STD by horizontal and vertical axes + double ysigma; +} object; + +// function for Qsort +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 process_file(const char *name){ +#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 ****/ + Image *I = Image_read(name); + DELTA("Imread"); + if(!I){ + WARNX("Can't read"); + return; + } + int W = I->width, H = I->height; + I->dtype = FLOAT_IMG; + save_fits(I, "fitsout.fits"); + DELTA("Save original"); + /* + uint8_t *outp = NULL; + if(GP->equalize) + outp = equalize(I, 3, GP->throwpart); + else + outp = linear(I, 3); + // draw test crosses + Pattern *cross = Pattern_cross(301, 301); + Img3 i3 = {.data = outp, .w = I->width, .h = I->height}; + Pattern_draw3(&i3, cross, I->width-100, I->height-100, C_R); + Pattern_draw3(&i3, cross, I->width/2, I->height/2, C_G); + Pattern_draw3(&i3, cross, 100, 100, C_W); + Pattern_free(&cross); + DBG("Try to write %s", name); + stbi_write_jpg("jpegout.jpg", I->width, I->height, 3, outp, 95); + FREE(outp);*/ +#ifdef GETMEDIAN + /**** get median image ****/ + Image *M = get_median(I, GP->medradius); + if(M){ + DELTA("Got median"); + /*outp = linear(M, 3); + stbi_write_jpg("median.jpg", I->width, I->height, 3, outp, 95); + FREE(outp);*/ + save_fits(M, "median.fits"); + DELTA("Save median"); +#endif +/* + Image *mean = NULL, *std = NULL; + if(get_stat(I, GP->medradius, &mean, &std)){ + DBG("Save std & mean"); + save_fits(mean, "mean.fits"); + save_fits(std, "std.fits"); + int wh = I->width*I->height; + Image *diff = Image_sim(I); + OMP_FOR() + for(int i = 0; i < wh; ++i){ + Imtype pixval = fabs(I->data[i] - M->data[i]); + //register Imtype mode = fabs(2.5*M->data[i] - 1.5*mean->data[i]); + //register Imtype pixval = fabs(mode - I->data[i]); + if(pixval > std->data[i]) diff->data[i] = 100.; + } + save_fits(diff, "val_med.fits"); + Image_free(&diff); + }else WARNX("Can't calculate statistics"); + Image_free(&mean); + Image_free(&std); +*/ + Imtype bk; + if(calc_background(I, &bk)){ + DBG("backgr = %g", bk); + DELTA("Got background"); + //uint8_t *ibin = Im2bin(I, 1960.); + 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); + 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); + //double averW = 0., averH = 0.; + //green("Found %zd objects\n", cc->Nobj-1); + //printf("%6s\t%6s\t%6s\t%6s\t%6s\t%6s\t%6s\n", "N", "area", "w/h", "Xc", "Yc", "Xw", "Yw"); + object *Objects = MALLOC(object, cc->Nobj-1); + int objctr = 0; + for(size_t i = 1; i < cc->Nobj; ++i){ + Box *b = &cc->boxes[i]; + //DBG("BOX %zd: area=%d, x:(%d-%d), y:(%d:%d)", i, b->area, b->xmin, b->xmax, b->ymin, b->ymax); + double wh = ((double)b->xmax - b->xmin)/(b->ymax - b->ymin); + if(wh < 0.77 || wh > 1.3) 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){ + //DBG("(%zd, %zd): %zd / %g", x, y, *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) + }; + ///object *o = &Objects[objctr-1]; + ///printf("%6d\t%6d\t%14.1f\t%6.1f\t%6.1f\t%6.1f\t%6.1f\t%6.1f\n", i, o->area, 40.-2.5*log(o->Isum), o->WdivH, o->xc, o->yc, o->xsigma, o->ysigma); + //averH += y2c; averW += x2c; + } + DELTA("Labeling"); + printf("%zd %d\n", time(NULL), objctr); + qsort(Objects, objctr, sizeof(object), compObjs); + object *o = Objects; + 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); + } + FREE(cc); + FREE(Objects); + Image *c = ST2Im(S, W, H); + FREE(S); + DELTA("conv size_t -> Ima"); + save_fits(c, "size_t.fits"); + Image_free(&c); + DELTA("Save size_t"); +#if 0 + uint8_t *f4 = filter8(ibin, W, H); + DELTA("calc f8"); + savebin(f4, W, H, "f8.fits"); + DELTA("save f8.fits"); + uint8_t *e1 = erosion(ibin, W, H); + DELTA("Get e"); + savebin(e1, W, H, "er.fits"); + DELTA("Save er.fits"); + uint8_t *e2 = dilation(ibin, W, H); + DELTA("Get di"); + savebin(e2, W, H, "dil.fits"); + DELTA("Save dil.fits"); + FREE(e1); + FREE(e2); +#endif +#if 0 + e1 = openingN(ibin, W, H, 1); + savebin(e1, W, H, "op.fits"); + e2 = closingN(ibin, W, H, 1); + savebin(e2, W, H, "cl.fits"); + FREE(e1); + FREE(e2); +#endif + FREE(ibin); + FREE(opn); + } + } +#ifdef GETMEDIAN + Image_free(&M); + } +#endif + DELTA("End"); + Image_free(&I); +} + +static int process_input(InputType tp, char *name){ + if(tp == T_DIRECTORY) return watch_directory(name, process_file); + return watch_file(name, process_file); +} + +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); + } + LOGMSG("Start application..."); + int p = process_input(tp, GP->inputname); + // never reached + signals(p); // clean everything + return p; +} diff --git a/median.c b/median.c new file mode 100644 index 0000000..e88de48 --- /dev/null +++ b/median.c @@ -0,0 +1,518 @@ +/* + * 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 + * @param bk (o) - background value + * @return 0 if error + */ +int calc_background(const 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; +#if 0 + //green("HISTO:\n"); + //for(int i = 0; i < 256; ++i) printf("%d:\t%d\t%d\t%d\n", i, histogram[i], diff[i], diff2[i]); + 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/median.h b/median.h new file mode 100644 index 0000000..c0a72c8 --- /dev/null +++ b/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(const Image *img, Imtype *bk); + +#endif // __MEDIAN_H__ diff --git a/stbimpl.c b/stbimpl.c new file mode 100644 index 0000000..4c76251 --- /dev/null +++ b/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