From f3570498e175aa410c2615f046bb46b91f5e0ae6 Mon Sep 17 00:00:00 2001 From: Edward Emelianov Date: Thu, 13 Apr 2023 18:03:05 +0300 Subject: [PATCH] add some binary morphological operations --- binmorph.c | 492 ++++++++++++++++++++++++++++++++++++++++ cclabling.h | 136 +++++++++++ cclabling_ori.h | 128 +++++++++++ converttypes.c | 231 +++++++++++++++++++ dilation.inc | 79 +++++++ ec_filter.inc | 102 +++++++++ examples/CMakeLists.txt | 3 + examples/Readme.md | 30 ++- examples/equalize.c | 2 + examples/gauss.c | 64 ++++++ examples/generate.c | 23 +- examples/objdet.c | 95 ++++++++ examples/poisson.c | 55 +++++ fc_filter.inc | 75 ++++++ imagefile.c | 224 +----------------- improclib.files | 9 + improclib.h | 71 +++++- random.c | 159 +++++++++++++ 18 files changed, 1743 insertions(+), 235 deletions(-) create mode 100644 binmorph.c create mode 100644 cclabling.h create mode 100644 cclabling_ori.h create mode 100644 converttypes.c create mode 100644 dilation.inc create mode 100644 ec_filter.inc create mode 100644 examples/gauss.c create mode 100644 examples/objdet.c create mode 100644 examples/poisson.c create mode 100644 fc_filter.inc create mode 100644 random.c diff --git a/binmorph.c b/binmorph.c new file mode 100644 index 0000000..39ba9a8 --- /dev/null +++ b/binmorph.c @@ -0,0 +1,492 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include // memcpy +#include +#include + +#include "openmp.h" +#include "improclib.h" + +//#define TESTMSGS +#ifdef TESTMSGS +#define TEST(...) printf(__VA_ARGS__) +#else +#define TEST(...) +#endif + +// 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 *ilfilter4(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; + { + // top of image, y = 0 + #define IM_UP + #include "fc_filter.inc" + #undef IM_UP + } + { + // mid of image, y = 1..h-1 + #include "fc_filter.inc" + } + { + // image bottom, y = h + #define IM_DOWN + #include "fc_filter.inc" + #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 *ilfilter8(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; + { + #define IM_UP + #include "ec_filter.inc" + #undef IM_UP + } + { + #include "ec_filter.inc" + } + { + #define IM_DOWN + #include "ec_filter.inc" + #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 *ildilation(const 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< + */ +/** + * 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 *ilimand(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 *ilsubstim(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; +} +/* + * <=================== LOGICAL OPERATIONS =================== + */ + +/* + * =================== CONNECTED COMPONENTS LABELING ===================> + */ + +// check table and rename all "oldval" into "newval" +static inline void remark(size_t newval, size_t oldval, size_t *assoc){ + 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); + } +} + +/** + * 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 (numeration starts from 1!!!), so first box is CC->box[1], amount of boxes is CC->Nobj-1 !!! + * @return an array of labeled components + */ +size_t *ilCClabel4(uint8_t *Img, int W, int H, ilConnComps **CC){ + size_t *assoc; + if(W < MINWIDTH || H < MINHEIGHT) return NULL; + uint8_t *f = ilfilter4(Img, W, H); // remove all non 4-connected pixels + //DBG("convert to size_t"); + size_t *labels = ilbin2sizet(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, assoc); + 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 + ilBox *boxes = MALLOC(ilBox, cidx); + OMP_FOR() + for(size_t i = 1; i < cidx; ++i){ // init borders + boxes[i].xmin = W; + boxes[i].ymin = H; + } +#pragma omp parallel shared(boxes) + { + ilBox *l_boxes = MALLOC(ilBox, 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; + ilBox *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){ + ilBox *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(ilConnComps, 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 + +/* + * <=================== CONNECTED COMPONENTS LABELING =================== + */ + + +/* + * <=================== template ===================> + */ diff --git a/cclabling.h b/cclabling.h new file mode 100644 index 0000000..e01f3c7 --- /dev/null +++ b/cclabling.h @@ -0,0 +1,136 @@ +/* + * This file is part of the loccorr project. + * Copyright 2021 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +//#define TESTMSGS + +#ifdef TESTMSGS +#define TEST(...) printf(__VA_ARGS__) +#else +#define TEST(...) +#endif + +size_t Nmax = W*H/4; // max number of 4-connected labels +size_t *assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array +size_t last_assoc_idx = 1; // last index filled in assoc array +// was: 35 +// check table and rename all "oldval" into "newval" +inline void remark(size_t newval, size_t oldval){ + TEST("\tnew = %zd, old=%zd; ", newval, oldval); + // find the least values + do{newval = assoc[newval];}while(assoc[newval] != newval); + do{oldval = assoc[oldval];}while(assoc[oldval] != oldval); + TEST("\trealnew = %zd, realold=%zd ", newval, oldval); + // now change larger value to smaller + if(newval > oldval){ + assoc[newval] = oldval; + TEST("change %zd to %zd\n", newval, oldval); + }else{ + assoc[oldval] = newval; + TEST("change %zd to %zd\n", oldval, newval); + } +} + +for(int y = 0; y < H; ++y){ + bool found = false; + size_t *ptr = &labels[y*W]; + size_t curmark = 0; // mark of pixel to the left + for(int x = 0; x < W; ++x, ++ptr){ + if(!*ptr){found = false; continue;} // empty pixel + size_t U = (y) ? ptr[-W] : 0; // upper mark + if(found){ // there's a pixel to the left + if(U && U != curmark){ // meet old mark -> remark one of them in assoc[] + TEST("(%d, %d): remark %zd --> %zd\n", x, y, U, curmark); + remark(U, curmark); + curmark = U; // change curmark to upper mark (to reduce further checks) + } + }else{ // new mark -> change curmark + found = true; + if(U) curmark = U; // current mark will copy upper value + else{ // current mark is new value + curmark = last_assoc_idx++; + assoc[curmark] = curmark; + TEST("(%d, %d): new mark=%zd\n", x, y, curmark); + } + } + *ptr = curmark; + } +} +size_t *indexes = MALLOC(size_t, last_assoc_idx); // new indexes +size_t cidx = 1; +TEST("\n\n\nRebuild indexes\n\n"); +for(size_t i = 1; i < last_assoc_idx; ++i){ + TEST("%zd\t%zd ",i,assoc[i]); + // search new index + register size_t realnew = i, newval = 0; + do{ + realnew = assoc[realnew]; + TEST("->%zd ", realnew); + if(indexes[realnew]){ // find least index + newval = indexes[realnew]; + TEST("real: %zd ", newval); + break; + } + }while(assoc[realnew] != realnew); + if(newval){ + TEST(" ==> %zd\n", newval); + indexes[i] = newval; + continue; + } + TEST("new index %zd\n", cidx); + // enter next label + indexes[i] = cidx++; +} +/* +// rebuild indexes +for(size_t i = 1; i < last_assoc_idx; ++i){ + printf("%zd\t%zd ",i,assoc[i]); + if(i == assoc[i]){ + if(i == cidx){ + printf("- keep\n"); + }else{ + printf("- change to %zd\n", cidx); + size_t _2change = assoc[i]; + assoc[i] = cidx; + for(size_t j = i+1; j < last_assoc_idx; ++j){ + if(assoc[j] == _2change){ + printf("\t%zd\t%zd -> %zd\n", j, assoc[j], cidx); + assoc[j] = cidx; + } + } + } + ++cidx; + }else printf("\n"); +}*/ + +--cidx; // cidx now is amount of detected objects +DBG("amount after rebuild: %zd", cidx); + +#ifdef TESTMSGS +printf("\n\n\nI\tASS[I]\tIDX[I]\n"); +for(size_t i = 1; i < last_assoc_idx; ++i) + printf("%zd\t%zd\t%zd\n",i,assoc[i],indexes[i]); +#endif + +int wh = W*H; +OMP_FOR() +for(int i = 0; i < wh; ++i) labels[i] = indexes[labels[i]]; + +if(Nobj) *Nobj = cidx; + +FREE(assoc); +FREE(indexes); diff --git a/cclabling_ori.h b/cclabling_ori.h new file mode 100644 index 0000000..76fe3bb --- /dev/null +++ b/cclabling_ori.h @@ -0,0 +1,128 @@ +/* + * cclabling_1.h - inner part of functions cclabel4 and cclabel8 + * + * Copyright 2015 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +//double t0 = dtime(); + + size_t N = 0, // current label + Nmax = W*H/4; // max number of labels + int w = W - 1, h = H - 1; + int y; + size_t *assoc = MALLOC(size_t, Nmax); // allocate memory for "remark" array + size_t last_assoc_idx = 0; // last index filled in assoc array + size_t currentnum = 0; // current pixel number + inline void remark(size_t old, size_t new){ // remark in assoc[] pixel with value old to assoc[new] or vice versa + size_t New = assoc[new], Old = assoc[old]; + if(Old == New){ + return; + } + // now we must check Old: if Old use its value + upmark = U[-1]; + }else // check point above only if there's nothing in left up + #endif + if(U[0]) upmark = U[0]; + #ifdef LABEL_8 + if(x < w && U[1]){ // there's something in upper right + if(upmark){ // to the left of it was pixels + remark(U[1], upmark); + }else + upmark = U[1]; + } + #endif + } + if(!upmark){ // there's nothing above - set current pixel to incremented counter + #ifdef LABEL_8 // check, whether pixel is not single + size_t *D = (y < h) ? &ptr[W] : NULL; + if( !(x && ((D && D[-1]) /*|| ptr[-1]*/)) // no left neighbours + && !(x < w && ((D && D[1]) || ptr[1])) // no right neighbours + && !(D && D[0])){ // no neighbour down + *ptr = 0; // erase this hermit! + continue; + } + #else + // no neighbour down & neighbour to the right -> hermit + if((y < h && ptr[W] == 0) && (x < w && ptr[1] == 0)){ + *ptr = 0; // erase this hermit! + continue; + } + #endif + upmark = ++N; + assoc[upmark] = ++currentnum; // refresh "assoc" + last_assoc_idx = upmark + 1; + } + *ptr = upmark; + curmark = upmark; + }else{ // there was something to the left -> we must chek only U[1] + if(U){ + if(x < w && U[1]){ // there's something in upper right + remark(U[1], curmark); + } + } + *ptr = curmark; + } + } + } +#if 0 + // Step 2: rename markers + // first correct complex assotiations in assoc + OMP_FOR() + for(y = 0; y < H; y++){ + size_t *ptr = &labels[y*W]; + for(int x = 0; x < W; x++, ptr++){ + size_t p = *ptr; + if(!p){continue;} + *ptr = assoc[p]; + } + } +#endif + FREE(assoc); + if(Nobj) *Nobj = currentnum; +//printf("%6.4f\t%zd\n", dtime()-t0, currentnum); diff --git a/converttypes.c b/converttypes.c new file mode 100644 index 0000000..7f316a7 --- /dev/null +++ b/converttypes.c @@ -0,0 +1,231 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include + +#include "improclib.h" +#include "openmp.h" + +/** + * @brief bin2Im - convert binarized image into uint8t + * @param image - binarized image + * @param W, H - its size (in pixels!) + * @return Image structure + */ +ilImage *ilbin2Image(const uint8_t *image, int W, int H){ + ilImage *ret = ilImage_new(W, H, IMTYPE_U8); + int stride = (W + 7) / 8, s1 = (stride*8 == W) ? stride : stride - 1; + uint8_t *data = (uint8_t*) ret->data; + int rest = W - s1*8; + OMP_FOR() + for(int y = 0; y < H; y++){ + uint8_t *optr = &data[y*W]; + const 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) ? 255 : 0; + inp <<= 1; + } + } + if(rest){ + register uint8_t inp = *iptr; + for(int i = 0; i < rest; ++i){ + *optr++ = (inp & 0x80) ? 255 : 0; + inp <<= 1; + } + } + } + ret->minval = 0; + ret->maxval = 255; + return ret; +} + +/** + * Convert image into pseudo-packed (1 char == 8 pixels), all values > bk 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 *ilImage2bin(const ilImage *im, double bk){ + if(!im) return NULL; + if(im->type != IMTYPE_U8){ + WARNX("ilImage2bin(): supported only 8-bit images"); + 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); + uint8_t *data = (uint8_t*) im->data; + int rest = W - s1*8; + //OMP_FOR() + for(y = 0; y < H; ++y){ + uint8_t *iptr = &data[y*W]; + uint8_t *optr = &ret[y*W0]; + for(int x = 0; x < s1; ++x){ + register uint8_t o = 0; + for(int i = 0; i < 8; ++i){ + o <<= 1; + if(*iptr++ > bk) o |= 1; + } + *optr++ = o; + } + if(rest){ + register uint8_t o = 0; + for(int x = 0; x < rest; ++x){ + o <<= 1; + if(*iptr++ > bk) o |= 1; + } + *optr = o << (8 - rest); + } + } + return ret; +} + +// transformation functions for Im2u8 +#define TRANSMACRO(datatype) \ + int height = I->height, width = I->width, stride = width * nchannels; \ + uint8_t *outp = MALLOC(uint8_t, height * stride); \ + double min = I->minval, max = I->maxval, W = 255./(max - min); \ + datatype *Idata = (datatype*) I->data; \ + if(nchannels == 3){ \ + OMP_FOR() \ + for(int y = 0; y < height; ++y){ \ + uint8_t *Out = &outp[y*stride]; \ + datatype *In = &Idata[y*width]; \ + for(int x = 0; x < width; ++x){ \ + Out[0] = Out[1] = Out[2] = (uint8_t)(W*((double)(*In++) - min)); \ + Out += 3; \ + } \ + } \ + }else{ \ + OMP_FOR() \ + for(int y = 0; y < height; ++y){ \ + uint8_t *Out = &outp[y*stride]; \ + datatype *In = &Idata[y*width]; \ + for(int x = 0; x < width; ++x){ \ + *Out++ = (uint8_t)(W*((double)(*In++) - min)); \ + } \ + } \ + } + +static uint8_t *Iu8(const ilImage *I, int nchannels){ + TRANSMACRO(uint8_t); + return outp; +} +static uint8_t *Iu16(const ilImage *I, int nchannels){ + TRANSMACRO(uint16_t); + return outp; +} +static uint8_t *Iu32(const ilImage *I, int nchannels){ + TRANSMACRO(uint32_t); + return outp; +} +static uint8_t *If(const ilImage *I, int nchannels){ + TRANSMACRO(float); + return outp; +} +static uint8_t *Id(const ilImage *I, int nchannels){ + TRANSMACRO(double); + return outp; +} + +/** + * @brief ilImage2u8 - linear transform for preparing file to save as JPEG or other type + * @param I - input image + * @param nchannels - 1 or 3 colour channels + * @return allocated here image for jpeg/png storing + */ +uint8_t *ilImage2u8(ilImage *I, int nchannels){ // only 1 and 3 channels supported! + if(!I || !I->data || (nchannels != 1 && nchannels != 3)) return NULL; + ilImage_minmax(I); + //DBG("make linear transform %dx%d, %d channels", I->width, I->height, nchannels); + switch(I->type){ + case IMTYPE_U8: + return Iu8(I, nchannels); + break; + case IMTYPE_U16: + return Iu16(I, nchannels); + break; + case IMTYPE_U32: + return Iu32(I, nchannels); + break; + case IMTYPE_F: + return If(I, nchannels); + break; + case IMTYPE_D: + return Id(I, nchannels); + break; + default: + WARNX("Im2u8: unsupported image type %d", I->type); + } + return NULL; +} + +#if 0 +UNUSED function! Need to be refactored +// convert size_t labels into Image +Image *ST2Im(const 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]; + const size_t *iptr = &image[y*W]; + for(int x = 0; x < W; ++x){ + *optr++ = (Imtype)*iptr++; + } + } + Image_minmax(ret); + return ret; +} +#endif + +/** + * 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 *ilbin2sizet(const 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]; + const 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; + if(rest){ + register uint8_t inp = *iptr; + for(int i = 0; i < rest; ++i){ + *optr++ = (inp & 0x80) ? 1 : 0; + inp <<= 1; + } + } + } + return ret; +} diff --git a/dilation.inc b/dilation.inc new file mode 100644 index 0000000..8dcca97 --- /dev/null +++ b/dilation.inc @@ -0,0 +1,79 @@ +/* + * dilation.h - inner part of function `dilation` + * + * Copyright 2013 Edward V. Emelianoff + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, + * MA 02110-1301, USA. + */ + +/* + * HERE'S NO ANY "FILE-GUARDS" BECAUSE FILE IS MULTIPLY INCLUDED! + * I use this because don't want to dive into depths of BOOST + * + * multi-including with different defines before - is a simplest way to + * modify simple code for avoiding extra if-then-else + */ + + +#if ! defined IM_UP && ! defined IM_DOWN // image without upper and lower borders +OMP_FOR() +for(int y = 1; y < h; y++) +#endif +{ +#ifdef IM_UP + int y = 0; +#endif +#ifdef IM_DOWN + int y = h; +#endif + const uint8_t *iptr = &image[W0*y]; + uint8_t *optr = &ret[W0*y]; + uint8_t p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + *optr = p; + optr++; iptr++; + for(int x = 1; x < w; x++, iptr++, optr++){ + p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + if(iptr[-1] & 1) p |= 0x80; + *optr = p; + } + p = DIL[*iptr] + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[-1] & 1) p |= 0x80; + *optr = p & lastmask; // clear outern pixel +} + diff --git a/ec_filter.inc b/ec_filter.inc new file mode 100644 index 0000000..b5216a5 --- /dev/null +++ b/ec_filter.inc @@ -0,0 +1,102 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +// inner part of function `filter8` + +#if ! defined IM_UP && ! defined IM_DOWN +OMP_FOR() +for(int y = 1; y < h; y++) +#endif +{ +#ifdef IM_UP + int y = 0; +#endif +#ifdef IM_DOWN + int y = h; +#endif + uint8_t *iptr = &image[W0*y]; + uint8_t *optr = &ret[W0*y]; + // x=0 + register uint8_t inp = *iptr; + register uint8_t p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1) + #endif + #ifndef IM_DOWN + | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1) + #endif + ; + if((iptr[1] + #ifndef IM_UP + | iptr[-W0+1] + #endif + #ifndef IM_DOWN + | iptr[W0+1] + #endif + )& 0x80) p |= 1; + *optr++ = inp & p; + ++iptr; + for(int x = 1; x < w; ++x){ + inp = *iptr; + p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1) + #endif + #ifndef IM_DOWN + | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1) + #endif + ; + if((iptr[1] + #ifndef IM_UP + | iptr[-W0+1] + #endif + #ifndef IM_DOWN + | iptr[W0+1] + #endif + )& 0x80) p |= 1; + if((iptr[-1] + #ifndef IM_UP + | iptr[-W0-1] + #endif + #ifndef IM_DOWN + | iptr[W0-1] + #endif + )& 1) p |= 0x80; + *optr++ = inp & p; + ++iptr; + } + // x = W0-1 + inp = *iptr; + p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] | (iptr[-W0]<<1) | (iptr[-W0]>>1) + #endif + #ifndef IM_DOWN + | iptr[W0] | (iptr[W0]<<1) | (iptr[W0]>>1) + #endif + ; + if((iptr[-1] + #ifndef IM_UP + | iptr[-W0-1] + #endif + #ifndef IM_DOWN + | iptr[W0-1] + #endif + )& 1) p |= 0x80; + *optr = inp & p; +} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index edb8897..b9a1990 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -9,3 +9,6 @@ link_libraries(usefull_macros improc m) add_executable(equalize equalize.c) add_executable(generate generate.c) add_executable(genu16 genu16.c) +add_executable(gauss gauss.c) +add_executable(poisson poisson.c) +add_executable(objdet objdet.c) diff --git a/examples/Readme.md b/examples/Readme.md index 8edbab8..359b8d6 100644 --- a/examples/Readme.md +++ b/examples/Readme.md @@ -16,11 +16,37 @@ Usage: genu16 [args] x1,y1[,w1] x2,y2[,w2] ... xn,yn[,w3] - draw 'stars' at coor -b, --beta=arg beta Moffat parameter of 'star' images (default: 1) -h, --height=arg resulting image height (default: 1024) -i, --input=arg input file with coordinates and amplitudes (comma separated) - -o, --output=arg output file name (default: output.png) + -l, --lambda=arg lambda of Poisson noice (default: 10) + -o, --output=arg output file name (default: output.jpg) -s, --halfwidth=arg FWHM of 'star' images (default: 3.5) -w, --width=arg resulting image width (default: 1024) ## genu16 -The same as 'generate', but works with 16-bit image and save it as 1-channel png (`ampi` now is weight). +The same as 'generate', but works with 16-bit image and save it as 1-channel png (`ampi` now is weight). No noice. + +## gauss +Simulator of falling photons with given gaussian distribution. + +Usage: gauss [args], where args are: + + -?, --help show this help + -X, --xstd=arg STD of 'photons' distribution by X (default: 10) + -Y, --ystd=arg STD of 'photons' distribution by Y (default: 10) + -h, --height=arg resulting image height (default: 1024) + -n, --niter=arg iterations ("falling photons") number (default: 1000000) + -o, --output=arg output file name (default: output.png) + -w, --width=arg resulting image width (default: 1024) + -x, --xcenter=arg X coordinate of 'image' center (default: 512) + -y, --ycenter=arg Y coordinate of 'image' center (default: 512) + + +## poisson +Add poisson noice to every pixel of image. + + -?, --help show this help + -h, --height=arg resulting image height (default: 1024) + -l, --lambda=arg mean (and dispersion) of distribution (default: 15.) + -o, --output=arg output file name (default: output.png) + -w, --width=arg resulting image width (default: 1024) diff --git a/examples/equalize.c b/examples/equalize.c index 64b9026..becb738 100644 --- a/examples/equalize.c +++ b/examples/equalize.c @@ -31,7 +31,9 @@ int main(int argc, char **argv){ return 2; } int w = I->width, h = I->height; + double t0 = dtime(); uint8_t *eq = ilequalize8(I, 3, 0.1); + green("Equalize: %g ms\n", (dtime() - t0)*1e3); ilImage_free(&I); if(!eq) return 3; ilImg3 *I3 = MALLOC(ilImg3, 1); diff --git a/examples/gauss.c b/examples/gauss.c new file mode 100644 index 0000000..5048102 --- /dev/null +++ b/examples/gauss.c @@ -0,0 +1,64 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "improclib.h" +#include +#include + +static int help = 0, w = 1024, h = 1024, Niter = 1000000; +static double xsigma = 10., ysigma = 10., x0 = 512., y0 = 512.; +static char *outp = "output.png"; + +static myoption cmdlnopts[] = { + {"help", NO_ARGS, NULL, '?', arg_int, APTR(&help), "show this help"}, + {"width", NEED_ARG, NULL, 'w', arg_int, APTR(&w), "resulting image width (default: 1024)"}, + {"height", NEED_ARG, NULL, 'h', arg_int, APTR(&h), "resulting image height (default: 1024)"}, + {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&outp), "output file name (default: output.png)"}, + {"xstd", NEED_ARG, NULL, 'X', arg_double, APTR(&xsigma), "STD of 'photons' distribution by X (default: 10)"}, + {"ystd", NEED_ARG, NULL, 'Y', arg_double, APTR(&ysigma), "STD of 'photons' distribution by Y (default: 10)"}, + {"xcenter", NEED_ARG, NULL, 'x', arg_double, APTR(&x0), "X coordinate of 'image' center (default: 512)"}, + {"ycenter", NEED_ARG, NULL, 'y', arg_double, APTR(&y0), "Y coordinate of 'image' center (default: 512)"}, + {"niter", NEED_ARG, NULL, 'n', arg_int, APTR(&Niter), "iterations (\"falling photons\") number (default: 1000000)"}, + end_option +}; + +int main(int argc, char **argv){ + initial_setup(); + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(w < 1 || h < 1) ERRX("Wrong image size"); + if(xsigma < DBL_EPSILON || ysigma < DBL_EPSILON) ERRX("STD should be >0"); + if(Niter < 1) ERRX("Iteration number should be a large positive number"); + ilImage *I = ilImage_new(w, h, IMTYPE_U8); + if(!I) ERRX("Can't create image %dx%d pixels", w, h); + int hits = 0; + for(int i = 0; i < Niter; ++i){ + //int x = (int)ilNormal(x0, sigma), y = (int)ilNormal(y0, sigma); + double x, y; + ilNormalPair(&x, &y, x0, y0, xsigma, ysigma); + if(x < 0 || x >= I->width || y < 0 || y >= I->height) continue; + uint8_t *pix = I->data + (int)x + ((int)y)*I->width; + if(*pix < 255) ++*pix; + ++hits; + } + int ret = ilwrite_png(outp, I->width, I->height, 1, I->data); + ilImage_free(&I); + if(!ret) return 1; + printf("File %s ready; %d hits of %d\n", outp, hits, Niter); + return 0; +} diff --git a/examples/generate.c b/examples/generate.c index 6e689af..810468c 100644 --- a/examples/generate.c +++ b/examples/generate.c @@ -22,7 +22,7 @@ #include static int w = 1024, h = 1024, help = 0; -static double fwhm = 3.5, beta = 1.; +static double fwhm = 3.5, beta = 1., lambda = 10.; static char *outp = "output.jpg", *inp = NULL; static ilPattern *star = NULL, *cross = NULL; @@ -34,6 +34,7 @@ static myoption cmdlnopts[] = { {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&outp), "output file name (default: output.jpg)"}, {"halfwidth",NEED_ARG, NULL, 's', arg_double, APTR(&fwhm), "FWHM of 'star' images (default: 3.5)"}, {"beta", NEED_ARG, NULL, 'b', arg_double, APTR(&beta), "beta Moffat parameter of 'star' images (default: 1)"}, + {"lambda", NEED_ARG, NULL, 'l', arg_double, APTR(&lambda), "lambda of Poisson noice (default: 10)"}, {"input", NEED_ARG, NULL, 'i', arg_string, APTR(&inp), "input file with coordinates and amplitudes (comma separated)"}, end_option }; @@ -63,10 +64,15 @@ static void addstar(ilImg3 *I, const char *str){ printf("Add 'star' at %d,%d (ampl=%d)\n", x,y,a); uint8_t c[3] = {a,a,a}; ilImg3_drawpattern(I, star, x, y, c); +} +static void addcross(ilImg3 *I, const char *str){ + int x, y, a; + if(!getpars(str, &x, &y, &a)) return; + printf("Add 'cross' at %d,%d (ampl=%d)\n", x,y,a); ilImg3_drawpattern(I, cross, x, y, ilColor_red); } -static void addfromfile(ilImg3 *I){ +static void addfromfile(ilImg3 *I, void (*fn)(ilImg3*, const char*)){ FILE *f = fopen(inp, "r"); if(!f){ WARN("Can't open %s", inp); @@ -74,7 +80,7 @@ static void addfromfile(ilImg3 *I){ } char *line = NULL; size_t n = 0; - while(getline(&line, &n, f) > 0) addstar(I, line); + while(getline(&line, &n, f) > 0) fn(I, line); fclose(f); } @@ -92,8 +98,15 @@ int main(int argc, char **argv){ star = ilPattern_star(par, par, fwhm, beta); cross = ilPattern_xcross(25, 25); for(int i = 0; i < argc; ++i) addstar(I, argv[i]); - if(inp) addfromfile(I); + if(inp) addfromfile(I, addstar); ilPattern_free(&star); + double t0 = dtime(); + ilImg3_addPoisson(I, lambda); + green("Poisson noice took %gms\n", (dtime()-t0) * 1e3); + if(!ilImg3_jpg(outp, I, 95)) WARNX("Can't save %s", outp); + for(int i = 0; i < argc; ++i) addcross(I, argv[i]); + if(inp) addfromfile(I, addcross); + ilPattern_free(&cross); uint8_t color[] = {255, 0, 100}; //uint8_t color[] = {0, 0, 0}; ilImg3_putstring(I, "Test string", 450, 520, color); @@ -106,7 +119,7 @@ int main(int argc, char **argv){ ilImg3_jpg("outpsubimage.jpg", s, 95); ilImg3_free(&s); }else WARNX("Bad subimage parameters"); - int ret = ilImg3_jpg(outp, I, 95); + int ret = ilImg3_jpg("crosses.jpg", I, 95); //int ret = ilImg3_png(outp, I); ilImg3_free(&I); if(!ret) return 4; diff --git a/examples/objdet.c b/examples/objdet.c new file mode 100644 index 0000000..7c4ead5 --- /dev/null +++ b/examples/objdet.c @@ -0,0 +1,95 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +// detect objects on an image + +#include "improclib.h" +#include +#include +#include + +static int help = 0, ndilat = 0, neros = 0; +double bg = -1.; +static char *infile = NULL, *outbg = NULL, *outbin = NULL; + +static myoption cmdlnopts[] = { + {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), "show this help"}, + {"input", NEED_ARG, NULL, 'i', arg_string, APTR(&infile), "input file name"}, + {"obg", NEED_ARG, NULL, 0, arg_string, APTR(&outbg), "input minus bg jpeg filename"}, + {"background",NEED_ARG, NULL, 'b', arg_double, APTR(&bg), "background level (default: auto)"}, + {"obin", NEED_ARG, NULL, 0, arg_string, APTR(&outbin), "--obg after binarizing"}, + {"ndilat", NEED_ARG, NULL, 'd', arg_int, APTR(&ndilat), "amount of dilations after erosions"}, + {"neros", NEED_ARG, NULL, 'e', arg_int, APTR(&neros), "amount of image erosions"}, + end_option +}; + +int main(int argc, char **argv){ + initial_setup(); + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(!infile) ERRX("Point name of input file"); + ilImage *I = ilImage_read(infile); + if(!I) ERR("Can't read %s", infile); + if(bg < 0. && !ilImage_background(I, &bg)) ERRX("Can't calculate background"); + uint8_t ibg = (int)(bg + 0.5); + printf("Background level: %d\n", ibg); + int w = I->width, h = I->height, wh = w*h; + ilImage *Ibg = ilImage_sim(I); + memcpy(Ibg->data, I->data, wh); + uint8_t *idata = (uint8_t*) Ibg->data; + for(int i = 0; i < wh; ++i) idata[i] = (idata[i] > ibg) ? idata[i] - ibg : 0; + if(outbg) ilwrite_jpg(outbg, Ibg->width, Ibg->height, 1, idata, 95); + double t0 = dtime(); + uint8_t *Ibin = ilImage2bin(I, bg); + if(!Ibin) ERRX("Can't binarize image"); + green("Binarization: %gms\n", 1e3*(dtime()-t0)); + if(neros > 0){ + t0 = dtime(); + uint8_t *eros = ilerosionN(Ibin, w, h, neros); + FREE(Ibin); + Ibin = eros; + green("%d erosions: %gms\n", neros, 1e3*(dtime()-t0)); + } + if(ndilat > 0){ + t0 = dtime(); + uint8_t *dilat = ildilationN(Ibin, w, h, ndilat); + FREE(Ibin); + Ibin = dilat; + green("%d dilations: %gms\n", ndilat, 1e3*(dtime()-t0)); + } + if(outbin){ + ilImage *tmp = ilbin2Image(Ibin, w, h); + ilwrite_jpg(outbin, tmp->width, tmp->height, 1, (uint8_t*)tmp->data, 95); + ilImage_free(&tmp); + } + ilConnComps *comps; + t0 = dtime(); + size_t *labels = ilCClabel4(Ibin, w, h, &comps); + green("Labeling: %gms\n", 1e3*(dtime()-t0)); + if(labels && comps->Nobj > 1){ + printf("Detected %zd components\n", comps->Nobj-1); + ilBox *box = comps->boxes + 1; + for(size_t i = 1; i < comps->Nobj; ++i, ++box){ + printf("\t%4zd: s=%d, LU=(%d, %d), RD=(%d, %d)\n", i, box->area, box->xmin, box->ymin, box->xmax, box->ymax); + } + FREE(labels); + FREE(comps); + } + ; + return 0; +} diff --git a/examples/poisson.c b/examples/poisson.c new file mode 100644 index 0000000..1e02629 --- /dev/null +++ b/examples/poisson.c @@ -0,0 +1,55 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "improclib.h" +#include +#include + +static int help = 0, w = 1024, h = 1024; +static double lambda = 15.; +static char *outp = "output.png"; + +static myoption cmdlnopts[] = { + {"help", NO_ARGS, NULL, '?', arg_int, APTR(&help), "show this help"}, + {"width", NEED_ARG, NULL, 'w', arg_int, APTR(&w), "resulting image width (default: 1024)"}, + {"height", NEED_ARG, NULL, 'h', arg_int, APTR(&h), "resulting image height (default: 1024)"}, + {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&outp), "output file name (default: output.png)"}, + {"lambda", NEED_ARG, NULL, 'l', arg_double, APTR(&lambda), "mean (and dispersion) of distribution (default: 15.)"}, + end_option +}; + +int main(int argc, char **argv){ + initial_setup(); + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(w < 1 || h < 1) ERRX("Wrong image size"); + if(lambda < 1.) ERRX("LAMBDA should be >=1"); + ilImage *I = ilImage_new(w, h, IMTYPE_U8); + if(!I) ERRX("Can't create image %dx%d pixels", w, h); + int npix = I->height * I->width; + uint8_t *d = I->data; + for(int i = 0; i < npix; ++i, ++d){ + int ampl = ilPoisson(lambda); + *d = ampl < 255 ? ampl : 255; + } + int ret = ilwrite_png(outp, I->width, I->height, 1, I->data); + ilImage_free(&I); + if(!ret) return 1; + printf("File %s ready\n", outp); + return 0; +} diff --git a/fc_filter.inc b/fc_filter.inc new file mode 100644 index 0000000..555d580 --- /dev/null +++ b/fc_filter.inc @@ -0,0 +1,75 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + + +// inner part of function `filter4` + +#if ! defined IM_UP && ! defined IM_DOWN +OMP_FOR() +for(int y = 1; y < h; y++) +#endif +{ +#ifdef IM_UP + int y = 0; +#endif +#ifdef IM_DOWN + int y = h; +#endif + uint8_t *iptr = &image[W0*y]; + uint8_t *optr = &ret[W0*y]; + // x=0 + register uint8_t inp = *iptr; + register uint8_t p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + *optr++ = inp & p; + ++iptr; + for(int x = 1; x < w; ++x){ + inp = *iptr; + p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[1] & 0x80) p |= 1; + if(iptr[-1] & 1) p |= 0x80; + *optr++ = inp & p; + ++iptr; + } + // x = W0-1 + inp = *iptr; + p = (inp << 1) | (inp >> 1) + #ifndef IM_UP + | iptr[-W0] + #endif + #ifndef IM_DOWN + | iptr[W0] + #endif + ; + if(iptr[-1] & 1) p |= 0x80; + *optr = inp & p; +} diff --git a/imagefile.c b/imagefile.c index eb3f1c3..1779578 100644 --- a/imagefile.c +++ b/imagefile.c @@ -143,7 +143,7 @@ ilInputType ilchkinput(const char *name){ * @param stride - image width with alignment * @return Image structure (fully allocated, you can FREE(data) after it) */ -ilImage *ilu8toImage(const uint8_t *data, int width, int height){ +ilImage *ilu82Image(const uint8_t *data, int width, int height){ FNAME(); ilImage *outp = ilImage_new(width, height, IMTYPE_U8); memcpy(outp->data, data, width*height); @@ -268,6 +268,11 @@ size_t *ilhistogram8(const ilImage *I){ for(int i = 0; i < 256; ++i) histogram[i] += histogram_private[i]; } } +#if 0 + for(int i = 0; i < wh; ++i){ + ++histogram[data[i]]; + } +#endif return histogram; } @@ -363,86 +368,6 @@ int ilImage_background(ilImage *img, double *bkg){ return TRUE; } -// transformation functions for Im2u8 -#define TRANSMACRO(datatype) \ - int height = I->height, width = I->width, stride = width * nchannels; \ - uint8_t *outp = MALLOC(uint8_t, height * stride); \ - double min = I->minval, max = I->maxval, W = 255./(max - min); \ - datatype *Idata = (datatype*) I->data; \ - if(nchannels == 3){ \ - OMP_FOR() \ - for(int y = 0; y < height; ++y){ \ - uint8_t *Out = &outp[y*stride]; \ - datatype *In = &Idata[y*width]; \ - for(int x = 0; x < width; ++x){ \ - Out[0] = Out[1] = Out[2] = (uint8_t)(W*((double)(*In++) - min)); \ - Out += 3; \ - } \ - } \ - }else{ \ - OMP_FOR() \ - for(int y = 0; y < height; ++y){ \ - uint8_t *Out = &outp[y*stride]; \ - datatype *In = &Idata[y*width]; \ - for(int x = 0; x < width; ++x){ \ - *Out++ = (uint8_t)(W*((double)(*In++) - min)); \ - } \ - } \ - } - -static uint8_t *Iu8(const ilImage *I, int nchannels){ - TRANSMACRO(uint8_t); - return outp; -} -static uint8_t *Iu16(const ilImage *I, int nchannels){ - TRANSMACRO(uint16_t); - return outp; -} -static uint8_t *Iu32(const ilImage *I, int nchannels){ - TRANSMACRO(uint32_t); - return outp; -} -static uint8_t *If(const ilImage *I, int nchannels){ - TRANSMACRO(float); - return outp; -} -static uint8_t *Id(const ilImage *I, int nchannels){ - TRANSMACRO(double); - return outp; -} - -/** - * @brief ilImage2u8 - linear transform for preparing file to save as JPEG or other type - * @param I - input image - * @param nchannels - 1 or 3 colour channels - * @return allocated here image for jpeg/png storing - */ -uint8_t *ilImage2u8(ilImage *I, int nchannels){ // only 1 and 3 channels supported! - if(!I || !I->data || (nchannels != 1 && nchannels != 3)) return NULL; - ilImage_minmax(I); - //DBG("make linear transform %dx%d, %d channels", I->width, I->height, nchannels); - switch(I->type){ - case IMTYPE_U8: - return Iu8(I, nchannels); - break; - case IMTYPE_U16: - return Iu16(I, nchannels); - break; - case IMTYPE_U32: - return Iu32(I, nchannels); - break; - case IMTYPE_F: - return If(I, nchannels); - break; - case IMTYPE_D: - return Id(I, nchannels); - break; - default: - WARNX("Im2u8: unsupported image type %d", I->type); - } - return NULL; -} - /** * @brief equalize8 - 8bit image hystogram equalization * @param I - input image @@ -700,143 +625,6 @@ void ilImage_minmax(ilImage *I){ DBG("Image_minmax(): Min=%g, Max=%g, time: %gms", I->minval, I->maxval, (dtime()-t0)*1e3); } -/* - * =================== CONVERT IMAGE TYPES ===================> - */ - -/** - * @brief bin2Im - convert binarized image into uint8t - * @param image - binarized image - * @param W, H - its size (in pixels!) - * @return Image structure - */ -ilImage *ilbin2Image(const uint8_t *image, int W, int H){ - ilImage *ret = ilImage_new(W, H, IMTYPE_U8); - int stride = (W + 7) / 8, s1 = (stride*8 == W) ? stride : stride - 1; - uint8_t *data = (uint8_t*) ret->data; - OMP_FOR() - for(int y = 0; y < H; y++){ - uint8_t *optr = &data[y*W]; - const 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) ? 255 : 0; - inp <<= 1; - } - } - int rest = W - s1*8; - if(rest){ - register uint8_t inp = *iptr; - for(int i = 0; i < rest; ++i){ - *optr++ = (inp & 0x80) ? 255 : 0; - inp <<= 1; - } - } - } - ret->minval = 0; - ret->maxval = 255; - return ret; -} - -/** - * Convert image into pseudo-packed (1 char == 8 pixels), all values > bk 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 *ilImage2bin(const ilImage *im, double bk){ - if(!im) return NULL; - if(im->type != IMTYPE_U8){ - WARNX("ilImage2bin(): supported only 8-bit images"); - 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); - uint8_t *data = (uint8_t*) im->data; - OMP_FOR() - for(y = 0; y < H; ++y){ - uint8_t *iptr = &data[y*W]; - uint8_t *optr = &ret[y*W0]; - for(int x = 0; x < s1; ++x){ - register uint8_t o = 0; - for(int i = 0; i < 8; ++i){ - o <<= 1; - if(*iptr++ > bk) o |= 1; - } - *optr++ = o; - } - int rest = W - s1*8; - if(rest){ - register uint8_t o = 0; - for(int x = 0; x < rest; ++x){ - o <<= 1; - if(*iptr++ > bk) o |= 1; - } - *optr = o << (8 - rest); - } - } - return ret; -} - -#if 0 -UNUSED function! Need to be refactored -// convert size_t labels into Image -Image *ST2Im(const 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]; - const size_t *iptr = &image[y*W]; - for(int x = 0; x < W; ++x){ - *optr++ = (Imtype)*iptr++; - } - } - Image_minmax(ret); - return ret; -} -#endif - -/** - * 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 *ilbin2sizet(const 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]; - const 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; - if(rest){ - register uint8_t inp = *iptr; - for(int i = 0; i < rest; ++i){ - *optr++ = (inp & 0x80) ? 1 : 0; - inp <<= 1; - } - } - } - return ret; -} - - -/* - * <=================== CONVERT IMAGE TYPES =================== - */ - /* * =================== SAVE IMAGES ===========================> */ diff --git a/improclib.files b/improclib.files index 6ae92b8..cad183b 100644 --- a/improclib.files +++ b/improclib.files @@ -1,11 +1,20 @@ +binmorph.c +converttypes.c +dilation.inc draw.c +ec_filter.inc examples/equalize.c +examples/gauss.c examples/generate.c examples/genu16.c +examples/objdet.c +examples/poisson.c +fc_filter.inc imagefile.c improclib.h letters.c openmp.h +random.c stb/stb_image.h stb/stb_image_write.h stbimpl.c diff --git a/improclib.h b/improclib.h index 357525e..77ad034 100644 --- a/improclib.h +++ b/improclib.h @@ -24,8 +24,6 @@ /*================================================================================* * basic types * *================================================================================*/ - - // 3-channel image for saving into jpg/png typedef struct{ uint8_t *data; // image data @@ -72,6 +70,15 @@ typedef enum{ T_AMOUNT } ilInputType; +/*================================================================================* + * converttypes.c * + *================================================================================*/ +ilImage *ilu82Image(const uint8_t *data, int width, int height); +uint8_t *ilImage2u8(ilImage *I, int nchannels); +ilImage *ilbin2Image(const uint8_t *image, int W, int H); +uint8_t *ilImage2bin(const ilImage *im, double bk); +size_t *ilbin2sizet(const uint8_t *image, int W, int H); + /*================================================================================* * draw.c * *================================================================================*/ @@ -104,10 +111,8 @@ ilImg3 *ilImg3_subimage(const ilImg3 *I, int x0, int y0, int x1, int y1); /*================================================================================* * imagefile.c * *================================================================================*/ - int ilgetpixbytes(ilimtype_t type); void ilImage_minmax(ilImage *I); -uint8_t *ilImage2u8(ilImage *I, int nchannels); uint8_t *ilequalize8(ilImage *I, int nchannels, double throwpart); uint8_t *ilequalize16(ilImage *I, int nchannels, double throwpart); @@ -121,11 +126,6 @@ size_t *ilhistogram8(const ilImage *I); size_t *ilhistogram16(const ilImage *I); int ilImage_background(ilImage *img, double *bkg); -ilImage *ilu8toImage(const uint8_t *data, int width, int height); -ilImage *ilbin2Image(const uint8_t *image, int W, int H); -uint8_t *ilImage2bin(const ilImage *im, double bk); -size_t *ilbin2sizet(const uint8_t *image, int W, int H); - ilImg3 *ilImg3_read(const char *name); int ilImg3_jpg(const char *name, ilImg3 *I3, int quality); int ilImg3_png(const char *name, ilImg3 *I3); @@ -140,8 +140,59 @@ int ilImg3_putstring(ilImg3 *I, const char *str, int x, int y, const uint8_t col /*================================================================================* - * * + * random.c * *================================================================================*/ +double ilNormalBase(); +double ilNormal(double mean, double std); +void ilNormalPair(double *x, double *y, double xmean, double ymean, double xstd, double ystd); + +int ilPoissonSetStep(double s); +int ilPoisson(double lambda); +void ilImage_addPoisson(ilImage *I, double lambda); +void ilImg3_addPoisson(ilImg3 *I, double lambda); + +/*================================================================================* + * binmorph.c * + *================================================================================*/ +// 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 +} ilBox; + +typedef struct{ + size_t Nobj; + ilBox *boxes; +} ilConnComps; + +// morphological operations: +uint8_t *ildilation(const uint8_t *image, int W, int H); +uint8_t *ildilationN(const uint8_t *image, int W, int H, int N); +uint8_t *ilerosion(const uint8_t *image, int W, int H); +uint8_t *ilerosionN(const uint8_t *image, int W, int H, int N); +uint8_t *ilopeningN(uint8_t *image, int W, int H, int N); +uint8_t *ilclosingN(uint8_t *image, int W, int H, int N); +uint8_t *iltopHat(uint8_t *image, int W, int H, int N); +uint8_t *ilbotHat(uint8_t *image, int W, int H, int N); + +// logical operations +uint8_t *ilimand(uint8_t *im1, uint8_t *im2, int W, int H); +uint8_t *ilsubstim(uint8_t *im1, uint8_t *im2, int W, int H); + +// clear non 4-connected pixels +uint8_t *ilfilter4(uint8_t *image, int W, int H); +// clear single pixels +uint8_t *ilfilter8(uint8_t *image, int W, int H); + +size_t *ilCClabel4(uint8_t *Img, int W, int H, ilConnComps **CC); +//size_t *ilcclabel8(uint8_t *Img, int W, int H, size_t *Nobj); /*================================================================================* * * diff --git a/random.c b/random.c new file mode 100644 index 0000000..be85d12 --- /dev/null +++ b/random.c @@ -0,0 +1,159 @@ +/* + * This file is part of the improclib project. + * Copyright 2023 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "improclib.h" +#include "openmp.h" +#include +#include + +// Box&Muller method for gaussian RNG (mean=0, std=1) +double ilNormalBase(){ + double U = drand48(), V = drand48(); + double S = sqrt(-2*log(U)); + double X = S * cos(2*M_PI*V); + double Y = S * sin(2*M_PI*V); + return X * Y; +} + +double ilNormal(double mean, double std){ + return ilNormalBase() * std + mean; +} + +/** + * @brief ilNormalPair - gaussian distributed pair of coordinates + * @param x,y - coordinates + * @param xmean, ymean - mean of coordinated + * @param xstd, ystd - sigma + */ +void ilNormalPair(double *x, double *y, double xmean, double ymean, double xstd, double ystd){ + if(!x || !y) return; + double U = drand48(), V = drand48(); + double S = sqrt(-2*log(U)); + *x = xmean + xstd * S * cos(2*M_PI*V); + *y = ymean + ystd * S * sin(2*M_PI*V); +} + +static double step = 500., es = -1.; +/** + * @brief ilPoissonSetStep - change step value for ilPoisson algo + * @param s - new step (>1) + * @return TRUE if OK + */ +int ilPoissonSetStep(double s){ + if(s < 1.){ + WARNX("ilPoissonSetStep(): step should be > 1."); + return FALSE; + } + step = s; + es = exp(step); + return TRUE; +} + +/** + * @brief ilPoisson - integer number with Poisson distribution + * @param lambda - mean (and dispersion) of distribution + * @return number + */ +int ilPoisson(double lambda){ + if(es < 0.) es = exp(step); + double L = lambda, p = 1.; + int k = 0; + do{ + ++k; + p *= drand48(); + while(p < 1. && L > 0.){ + if(L > step){ + p *= es; + L -= step; + }else{ + p *= exp(L); + L = 0; + } + } + }while(p > 1.); + return k-1; +} + +#define ADDPU(type, max) \ + type *d = (type*)I->data; \ + int wh = I->width * I->height; \ + for(int i = 0; i < wh; ++i, ++d){ \ + type res = *d + ilPoisson(lambda); \ + *d = (res >= *d) ? res : max; \ + } +#define ADDPF(type) \ + type *d = (type*)I->data; \ + int wh = I->width * I->height; \ + for(int i = 0; i < wh; ++i, ++d){ \ + *d += ilPoisson(lambda); \ + } +static void add8p(ilImage *I, double lambda){ + ADDPU(uint8_t, 0xff); +} +static void add16p(ilImage *I, double lambda){ + ADDPU(uint16_t, 0xffff); +} +static void add32p(ilImage *I, double lambda){ + ADDPU(uint32_t, 0xffffffff); +} +static void addfp(ilImage *I, double lambda){ + ADDPF(float); +} +static void adddp(ilImage *I, double lambda){ + ADDPF(double); +} +/** + * @brief ilImage_addPoisson - add poisson noice to each image pixel + * @param I (inout) - image + * @param lambda - lambda of noice + */ +void ilImage_addPoisson(ilImage *I, double lambda){ + switch(I->type){ + case IMTYPE_U8: + return add8p(I, lambda); + break; + case IMTYPE_U16: + return add16p(I, lambda); + break; + case IMTYPE_U32: + return add32p(I, lambda); + break; + case IMTYPE_F: + return addfp(I, lambda); + break; + case IMTYPE_D: + return adddp(I, lambda); + break; + default: + ERRX("ilImage_addPoisson(): invalid data type"); + } +} + +// the same as ilImage_addPoisson but for coloured image (add same noice to all three pixel colour components) +void ilImg3_addPoisson(ilImg3 *I, double lambda){ + int wh = I->width * I->height * 3; + uint8_t *id = (uint8_t*)I->data; + //OMP_FOR() - only will be more slowly + for(int i = 0; i < wh; i += 3){ + uint8_t n = ilPoisson(lambda), *d = &id[i]; + for(int j = 0; j < 3; ++j){ + uint8_t newval = d[j] + n; + d[j] = (newval >= d[j]) ? newval : 255; + } + } +}