eddy/astrovideoguide_v3
1
0
Fork 0
mirror of https://github.com/eddyem/astrovideoguide_v3.git synced 2025-12-06 02:35:11 +03:00
Code Issues Packages Projects Releases Wiki Activity

initial commit

Browse Source
This commit is contained in:
Edward Emelianov 2021-03-29 18:27:24 +03:00
parent 2711e8586c
commit ae3f52db70
22 changed files with 3342 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
.gitignore vendored Normal file
View File

@ -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

74
CMakeLists.txt Normal file
View File

@ -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")

571
binmorph.c Normal file
View File

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

85
binmorph.h Normal file
View File

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

136
cclabling.h Normal file
View File

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

94
cmdlnopts.c Normal file
View File

@ -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 <eddy@sao.ru>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
*/
#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <strings.h>
#include <math.h>
#include <usefull_macros.h>
#include "cmdlnopts.h"
/*
* here are global parameters initialisation
*/
static int help;
glob_pars *GP = NULL;
// default PID filename:
#define DEFAULT_PIDFILE "/tmp/loccorr.pid"
// DEFAULTS
// default global parameters
static glob_pars G = {
.pidfile = DEFAULT_PIDFILE,
.throwpart = 0.5,
.medradius = 1.,
.ndilations = 2,
.nerosions = 2,
.minarea = 5,
.intensthres = 0.01
};
/*
* 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;
}

47
cmdlnopts.h Normal file
View File

@ -0,0 +1,47 @@
/* geany_encoding=koi8-r
* cmdlnopts.h - comand line options for parceargs
*
* Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
*/
#pragma once
#ifndef CMDLNOPTS_H__
#define CMDLNOPTS_H__
/*
* here are some typedef's for global data
*/
typedef struct{
char *pidfile; // name of PID file
char *logfile; // logging to this file
char *inputname; // input for monitor file or directory name
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__

73
dilation.h Normal file
View File

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

105
draw.c Normal file
View File

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

56
draw.h Normal file
View File

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

96
ec_filter.h Normal file
View File

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

68
fc_filter.h Normal file
View File

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

227
fits.c Normal file
View File

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

55
fits.h Normal file
View File

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

504
imagefile.c Normal file
View File

@ -0,0 +1,504 @@
/*
* This file is part of the loccorr project.
* Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <dirent.h>
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <usefull_macros.h>
#include <stb/stb_image.h>
#include <stb/stb_image_write.h>
#include "cmdlnopts.h"
#include "draw.h"
#include "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<<x;
}
*optr = o;
//*optr = o | 1<<rest; // mark outern right edge (for good CC labeling)
}
}
return ret;
}
// convert size_t labels into floating
Image *ST2Im(size_t *image, int W, int H){
Image *ret = Image_new(W, H);
OMP_FOR()
for(int y = 0; y < H; ++y){
Imtype *optr = &ret->data[y*W];
size_t *iptr = &image[y*W];
for(int x = 0; x < W; ++x){
*optr++ = (Imtype)*iptr++;
}
}
ret->dtype = FLOAT_IMG;
Image_minmax(ret);
return ret;
}
/**
* Convert "packed" image into size_t array for conncomp procedure
* @param image (i) - input image
* @param W, H - size of image in pixels
* @return allocated memory area with copy of an image
*/
size_t *bin2ST(uint8_t *image, int W, int H){
size_t *ret = MALLOC(size_t, W * H);
int W0 = (W + 7) / 8, s1 = W0 - 1;
OMP_FOR()
for(int y = 0; y < H; y++){
size_t *optr = &ret[y*W];
uint8_t *iptr = &image[y*W0];
for(int x = 0; x < s1; ++x){
register uint8_t inp = *iptr++;
for(int i = 0; i < 8; ++i){
*optr++ = (inp & 0x80) ? 1 : 0;
inp <<= 1;
}
}
int rest = W - s1*8;
register uint8_t inp = *iptr;
for(int i = 0; i < rest; ++i){
*optr++ = (inp & 0x80) ? 1 : 0;
inp <<= 1;
}
}
return ret;
}
/*
* <=================== CONVERT IMAGE TYPES ===================
*/

53
imagefile.h Normal file
View File

@ -0,0 +1,53 @@
/*
* This file is part of the loccorr project.
* Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#ifndef IMAGEFILE_H__
#define IMAGEFILE_H__
#include <stdint.h>
#include "fits.h"
// input file/directory type
typedef enum{
T_WRONG,
T_DIRECTORY,
T_BMP,
T_FITS,
T_GZIP,
T_GIF,
T_JPEG,
T_PNG,
} 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__

125
inotify.c Normal file
View File

@ -0,0 +1,125 @@
/*
* This file is part of the loccorr project.
* Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h> // strlen
#include <sys/inotify.h>
#include <sys/select.h>
#include <unistd.h>
#include <usefull_macros.h>
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);
}

26
inotify.h Normal file
View File

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

345
main.c Normal file
View File

@ -0,0 +1,345 @@
/*
* This file is part of the loccorr project.
* Copyright 2021 Edward V. Emelianov <edward.emelianoff@gmail.com>.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <signal.h> // signal
#include <stb/stb_image_write.h>
#include <stdio.h> // printf
#include <stdlib.h> // exit, free
#include <string.h> // strdup
#include <time.h>
#include <unistd.h> // sleep
#include <usefull_macros.h>
#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;
}

518
median.c Normal file
View File

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

32
median.h Normal file
View File

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

29
stbimpl.c Normal file
View File

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