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

add some binary morphological operations

Browse Source
This commit is contained in:
Edward Emelianov 2023-04-13 18:03:05 +03:00
parent 03eecf2595
commit f3570498e1
18 changed files with 1743 additions and 235 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

492
binmorph.c Normal file
View File

@ -0,0 +1,492 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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 <stdio.h>
#include <stdlib.h>
#include <string.h> // memcpy
#include <math.h>
#include <sys/time.h>
#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<<rest);
if(!DIL) morph_init();
uint8_t *ret = MALLOC(uint8_t, W0*H);
{
// top of image, y = 0
#define IM_UP
#include "dilation.inc"
#undef IM_UP
}
{
// mid of image, y = 1..h-1
#include "dilation.inc"
}
{
// image bottom, y = h
#define IM_DOWN
#include "dilation.inc"
#undef IM_DOWN
}
return ret;
}
static void mkerosion(const uint8_t *in, uint8_t *out, int W, int H){
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);
// clear first and last rows:
bzero(out, W0);
bzero(out + W0*h, W0);
//DBG("rest=%d, mask:0x%x", rest, lastmask);
OMP_FOR()
for(int y = 1; y < h; y++){ // reset first & last rows of image
const uint8_t *iptr = &in[W0*y];
uint8_t *optr = &out[W0*y];
register uint8_t p;
// x == 0, reset left bit
p = ER[*iptr] & iptr[-W0] & iptr[W0];
if(!(*(++iptr) & 0x80)) p &= 0xfe;
*optr++ = p & 0x7f;
for(int x = 1; x < w; x++, iptr++, optr++){
p = iptr[-W0];
p = ER[*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[-W0] & iptr[W0];
if(!(iptr[-1] & 1)) p &= 0x7f;
*optr = p & lastmask; // reset right bit
}
}
/**
* 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 *ilerosion(const uint8_t *image, int W, int H){
if(W < MINWIDTH || H < MINHEIGHT) return NULL;
int W0 = (W + 7) / 8; // width in bytes
uint8_t *ret = MALLOC(uint8_t, W0*H);
mkerosion(image, ret, W, H);
return ret;
}
// Make erosion N times
uint8_t *ilerosionN(const uint8_t *image, int W, int H, int N){
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
int W0 = (W + 7) / 8, sz = W0*H;
uint8_t *in = MALLOC(uint8_t, sz), *out = MALLOC(uint8_t, sz);
memcpy(out, image, sz);
for(int i = 0; i < N; ++i){
register uint8_t *tmp = in;
in = out; out = tmp;
mkerosion(in, out, W, H);
}
FREE(in);
return out;
}
// Make dilation N times
uint8_t *ildilationN(const uint8_t *image, int W, int H, int N){
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
uint8_t *cur = (uint8_t*)image, *next = NULL;
for(int i = 0; i < N; ++i){
next = ildilation(cur, W, H);
if(cur != image) FREE(cur);
cur = next;
}
return next;
}
// Ntimes opening
uint8_t *ilopeningN(uint8_t *image, int W, int H, int N){
//FNAME();
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
uint8_t *er = ilerosionN(image, W, H, N);
uint8_t *op = ildilationN(er, W, H, N);
FREE(er);
return op;
}
// Ntimes closing
uint8_t *ilclosingN(uint8_t *image, int W, int H, int N){
//FNAME();
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
uint8_t *di = ildilationN(image, W, H, N);
uint8_t *cl = ilerosionN(di, W, H, N);
FREE(di);
return cl;
}
// top hat operation: image - opening(image)
uint8_t *iltopHat(uint8_t *image, int W, int H, int N){
//FNAME();
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
uint8_t *op = ilopeningN(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 *ilbotHat(uint8_t *image, int W, int H, int N){
//FNAME();
if(W < MINWIDTH || H < MINHEIGHT || N < 1) return NULL;
uint8_t *op = ilclosingN(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 ===================>
*/
/**
* 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 ===================>
*/

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);

128
cclabling_ori.h Normal file
View File

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

231
converttypes.c Normal file
View File

@ -0,0 +1,231 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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 <usefull_macros.h>
#include <stdio.h>
#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;
}

79
dilation.inc Normal file
View File

@ -0,0 +1,79 @@
/*
* 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
{
#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
}

102
ec_filter.inc Normal file
View File

@ -0,0 +1,102 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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
{
#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;
}

3
examples/CMakeLists.txt
View File

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

30
examples/Readme.md
View File

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

2
examples/equalize.c
View File

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

64
examples/gauss.c Normal file
View File

@ -0,0 +1,64 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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 "improclib.h"
#include <usefull_macros.h>
#include <stdio.h>
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;
}

23
examples/generate.c
View File

@ -22,7 +22,7 @@
#include <string.h>
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;

95
examples/objdet.c Normal file
View File

@ -0,0 +1,95 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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/>.
*/
// detect objects on an image
#include "improclib.h"
#include <usefull_macros.h>
#include <stdio.h>
#include <string.h>
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;
}

55
examples/poisson.c Normal file
View File

@ -0,0 +1,55 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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 "improclib.h"
#include <usefull_macros.h>
#include <stdio.h>
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;
}

75
fc_filter.inc Normal file
View File

@ -0,0 +1,75 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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
{
#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;
}

224
imagefile.c
View File

@ -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 ===========================>
*/

9
improclib.files
View File

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

71
improclib.h
View File

@ -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);
/*================================================================================*
* *

159
random.c Normal file
View File

@ -0,0 +1,159 @@
/*
* This file is part of the improclib project.
* Copyright 2023 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 "improclib.h"
#include "openmp.h"
#include <math.h>
#include <usefull_macros.h>
// 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;
}
}
}