From 88ae442d9392f24c44ac7656012bbc622b34db41 Mon Sep 17 00:00:00 2001 From: eddyem Date: Tue, 2 Apr 2019 14:01:09 +0300 Subject: [PATCH] Add simple median filtering, TODO: create adaptive median --- FITSmanip.c | 2 - FITSmanip.h | 62 +-- examples/CMakeLists.txt | 10 +- examples/Readme.md | 9 +- examples/gd.c | 70 ++- examples/med.c | 107 +++++ fitsimages.c | 17 +- locale/ru/LC_MESSAGES/FITSmanip.mo | Bin 2189 -> 2332 bytes locale/ru/messages.po | 16 +- locale/ru/ru.po | 16 +- median.c | 698 +++++++++++++++++++++++++++++ 11 files changed, 919 insertions(+), 88 deletions(-) create mode 100644 examples/med.c create mode 100644 median.c diff --git a/FITSmanip.c b/FITSmanip.c index dbbfa7e..2eae981 100644 --- a/FITSmanip.c +++ b/FITSmanip.c @@ -74,7 +74,6 @@ double sqrtrans(double in){ // square root typedef double (*transfunct)(double in); static transfunct tfunctions[TRANSF_COUNT] = { [TRANSF_EXP] = exptrans, - [TRANSF_HISTEQ] = NULL, // another type of transform! [TRANSF_LOG] = logtrans, [TRANSF_LINEAR] = lintrans, [TRANSF_POW] = powtrans, @@ -98,7 +97,6 @@ doubleimage *mktransform(doubleimage *im, imgstat *st, intens_transform transf){ } double *dimg = im->data; if(transf == TRANSF_LINEAR) return im; // identity -if(transf == TRANSF_HISTEQ) return NULL; // histogram equalization; TODO: add this option too! double (*transfn)(double in) = tfunctions[transf]; if(!transfn) ERRX(_("Given transform type not supported yet")); size_t totpix = im->totpix; diff --git a/FITSmanip.h b/FITSmanip.h index 8d467b4..bfa690a 100644 --- a/FITSmanip.h +++ b/FITSmanip.h @@ -94,7 +94,7 @@ typedef struct{ */ typedef struct{ int naxis; // amount of image dimensions - long *naxes; // dimensions + long *naxes; // dimensions (x, y, z, etc) long totpix; // total pixels amount int bitpix; // original bitpix int dtype; // type of stored data @@ -126,7 +126,6 @@ typedef enum{ TRANSF_EXP, TRANSF_POW, TRANSF_SQR, - TRANSF_HISTEQ, TRANSF_COUNT // amount of transforms } intens_transform; @@ -209,7 +208,8 @@ void *image_data_malloc(long totpix, int pxbytes); FITSimage *image_new(int naxis, long *naxes, int bitpix); FITSimage *image_mksimilar(FITSimage *in); FITSimage *image_copy(FITSimage *in); -void dblima_free(doubleimage **im); +doubleimage *doubleimage_new(size_t w, size_t h); +void doubleimage_free(doubleimage **im); doubleimage *image2double(FITSimage *img); imgstat *get_imgstat(const doubleimage *dimg, imgstat *est); doubleimage *normalize_dbl(doubleimage *dimg, imgstat *st); @@ -243,54 +243,12 @@ histogram *dbl2histogram(doubleimage *im, size_t nvalues); doubleimage *dbl_histcutoff(doubleimage *im, size_t nlevls, double fracbtm, double fractop); doubleimage *dbl_histeq(doubleimage *im, size_t nlevls); -/* -// pointer to image conversion function -typedef FITS* (*imfuncptr)(FITS *in, Filter *f, Itmarray *i); -*/ - -/* -// FilterType (not only convolution!) -typedef enum{ - FILTER_NONE = 0 // simple start - ,MEDIAN // median filter - ,ADPT_MEDIAN // simple adaptive median - ,LAPGAUSS // laplasian of gaussian - ,GAUSS // gaussian - ,SOBELH // Sobel horizontal - ,SOBELV // -//- vertical - ,SIMPLEGRAD // simple gradient (by Sobel) - ,PREWITTH // Prewitt (horizontal) - simple derivative - ,PREWITTV // -//- (vertical) - ,SCHARRH // Scharr (modified Sobel) - ,SCHARRV - ,STEP // "posterisation" -} FType; - -typedef struct{ - double *data; - size_t size; -}Itmarray; -*/ - -/* -typedef struct _Filter{ - char *name; // filter name - FType FilterType; // filter type - int w; // filter width - int h; // height - double sx; // x half-width - double sy; // y half-width (sx, sy - for Gaussian-type filters) - FITS* (*imfunc)(FITS *in, struct _Filter *f, Itmarray *i); // image function for given conversion type -} Filter; - -// mathematical operations when there's no '-i' parameter (for >1 FITS-files) -typedef enum{ - MATH_NONE = 0 - ,MATH_SUM // make sum of all files - ,MATH_MEDIAN // calculate median by all files - ,MATH_MEAN // calculate mean for all files - ,MATH_DIFF // difference of first and rest files -} MathOper; -*/ +/************************************************************************************** + * median.c * + **************************************************************************************/ +doubleimage *get_median(const doubleimage *img, size_t radius); +//doubleimage *get_adaptive_median(const doubleimage *img, size_t radius); +double quick_select(const double *idata, int n); +double calc_median(const double *idata, int n); #endif // FITSMANIP_H__ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2a5b436..0aa9e24 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -3,9 +3,9 @@ project(examples) include_directories(../) link_libraries(FITSmanip cfitsio m) -#add_executable(fitsstat fitsstat.c) -add_executable(keylist keylist.c) -add_executable(imstat imstat.c) -add_executable(listtable listtable.c) add_executable(gd gd.c) -target_link_libraries(gd -lgd) \ No newline at end of file +target_link_libraries(gd -lgd) +add_executable(imstat imstat.c) +add_executable(keylist keylist.c) +add_executable(listtable listtable.c) +add_executable(median med.c) diff --git a/examples/Readme.md b/examples/Readme.md index 5e2c7e3..a021ce0 100644 --- a/examples/Readme.md +++ b/examples/Readme.md @@ -10,16 +10,19 @@ Usage: gd [args] Where args are: - -T, --transform=arg type of intensity transformation (log, sqr, exp, pow) + -E, --histeq histogram equalisation + -H, --hcuthigh=arg histogram cut-off high limit + -L, --hcutlow=arg histogram cut-off low limit + -T, --transform=arg type of intensity transformation (exp, lin, log, pow, sqrt) -h, --help show this help -i, --inname=arg name of input file + -l, --histlvl=arg amount of levels for histogram calculation -n, --hdunumber=arg open image from given HDU number -o, --outpname=arg output file name (jpeg) - -p, --palette=arg convert as given palette + -p, --palette=arg convert as given palette (br, cold, gray, hot, jet) -r, --rewrite rewrite output file -t, --textline=arg add text line to output image (at bottom) - ## imstat.c Usage: imstat [args] input files diff --git a/examples/gd.c b/examples/gd.c index c85be5a..6147dc3 100644 --- a/examples/gd.c +++ b/examples/gd.c @@ -35,6 +35,9 @@ typedef struct{ int nhdu; // HDU number to read image from int rewrite; // rewrite output file int nlvl; // amount of histogram levels + int histeq; // histogram equalisation + double histcutlow; // low limit of histogram cut-off + double histcuthigh; // top limit -//- } glob_pars; /* @@ -52,15 +55,18 @@ static glob_pars G = { */ static myoption cmdlnopts[] = { // common options - {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")}, + {"help", NO_ARGS, NULL, 'h', arg_none, APTR(&help), _("show this help")}, {"inname", NEED_ARG, NULL, 'i', arg_string, APTR(&G.fitsname), _("name of input file")}, {"outpname",NEED_ARG, NULL, 'o', arg_string, APTR(&G.outfile), _("output file name (jpeg)")}, {"textline",NEED_ARG, NULL, 't', arg_string, APTR(&G.text), _("add text line to output image (at bottom)")}, - {"palette", NEED_ARG, NULL, 'p', arg_string, APTR(&G.palette), _("convert as given palette (linear, br, hot)")}, + {"palette", NEED_ARG, NULL, 'p', arg_string, APTR(&G.palette), _("convert as given palette (br, cold, gray, hot, jet)")}, {"hdunumber",NEED_ARG, NULL, 'n', arg_int, APTR(&G.nhdu), _("open image from given HDU number")}, - {"transform",NEED_ARG, NULL, 'T', arg_string, APTR(&G.transform), _("type of intensity transformation (log, sqr, exp, pow)")}, + {"transform",NEED_ARG, NULL, 'T', arg_string, APTR(&G.transform), _("type of intensity transformation (exp, lin, log, pow, sqrt)")}, {"rewrite", NO_ARGS, NULL, 'r', arg_none, APTR(&G.rewrite), _("rewrite output file")}, {"histlvl", NEED_ARG, NULL, 'l', arg_int, APTR(&G.nlvl), _("amount of levels for histogram calculation")}, + {"hcutlow", NEED_ARG, NULL, 'L', arg_double, APTR(&G.histcutlow),_("histogram cut-off low limit")}, + {"hcuthigh",NEED_ARG, NULL, 'H', arg_double, APTR(&G.histcuthigh),_("histogram cut-off high limit")}, + {"histeq", NO_ARGS, NULL, 'E', arg_none, APTR(&G.histeq), _("histogram equalisation")}, end_option }; @@ -159,9 +165,14 @@ static intens_transform gettransf(const char *transf){ case 's': return TRANSF_SQR; break; - default: - return TRANSF_WRONG; } + fprintf(stderr, "Possible arguments of " COLOR_RED "\"Transformation\"" COLOR_OLD ":\n"); + fprintf(stderr, "exp - exponential transform\n"); + fprintf(stderr, "linear (default) - linear transform (do nothing)\n"); + fprintf(stderr, "log - logariphmic transform\n"); + fprintf(stderr, "pow - x^2\n"); + fprintf(stderr, "sqrt - sqrt(x)\n"); + return TRANSF_WRONG; } /** @@ -184,17 +195,27 @@ static image_palette palette_transform(char *p){ case 'G': return PALETTE_GRAY; break; - case 'h': + case 'h': // hot, help case 'H': - return PALETTE_HOT; + switch(p[1]){ + case 'o': + case 'O': + return PALETTE_HOT; + break; + } break; case 'j': case 'J': return PALETTE_JET; break; - default: - return PALETTE_WRONG; } + fprintf(stderr, "Possible arguments of " COLOR_RED "\"palette\"" COLOR_OLD ":\n"); + fprintf(stderr, "br - blue->red->yellow->white\n"); + fprintf(stderr, "cold - black->blue->cyan->white\n"); + fprintf(stderr, "gray (default) - simple gray\n"); + fprintf(stderr, "hot - black->red->yellow->white\n"); + fprintf(stderr, "jet - black->white->blue\n"); + return PALETTE_WRONG; } void print_histo(histogram *H){ @@ -212,17 +233,19 @@ void print_histo(histogram *H){ int main(int argc, char *argv[]){ initial_setup(); parse_args(argc, argv); - if(!G.fitsname) ERRX(_("No input filename given!")); - if(!G.outfile) ERRX(_("Point the name of output file!")); - intens_transform tr = TRANSF_LINEAR; - if(G.transform) tr = gettransf(G.transform); - if(tr == TRANSF_WRONG) ERRX(_("Wrong transform: %s"), G.transform); - if(!file_absent(G.outfile) && !G.rewrite) ERRX(_("File %s exists"), G.outfile); image_palette colormap = PALETTE_GRAY; if(G.palette){ // convert normalized image due to choosen palette colormap = palette_transform(G.palette); - if(colormap == PALETTE_WRONG) ERRX(_("Wrong colormap name")); + if(colormap == PALETTE_WRONG) ERRX(_("Wrong colormap: %s"), G.palette); } + intens_transform tr = TRANSF_LINEAR; + if(G.transform){ + tr = gettransf(G.transform); + if(tr == TRANSF_WRONG) ERRX(_("Wrong transform: %s"), G.transform); + } + if(!G.fitsname) ERRX(_("No input filename given!")); + if(!G.outfile) ERRX(_("Point the name of output file!")); + if(!file_absent(G.outfile) && !G.rewrite) ERRX(_("File %s exists"), G.outfile); DBG("Open file %s", G.fitsname); FITS *f = FITS_read(G.fitsname); if(!f) ERRX(_("Failed to open")); @@ -244,12 +267,25 @@ int main(int argc, char *argv[]){ st = get_imgstat(dblimg, NULL); #endif DBG("NOW: MIN=%g, MAX=%g, AVR=%g, STD=%g", st->min, st->max, st->mean, st->std); + green("Histogram before transformations:\n"); + histogram *h = dbl2histogram(dblimg, G.nlvl); + print_histo(h); + histogram_free(&h); + if(G.histeq){ // equalize histogram + if(!dbl_histeq(dblimg, G.nlvl)) + ERRX(_("Can't do histogram equalization")); + } + if(G.histcutlow > DBL_EPSILON || G.histcuthigh > DBL_EPSILON){ + if(!dbl_histcutoff(dblimg, G.nlvl, G.histcutlow, G.histcuthigh)) + ERRX(_("Can't make histogram cut-off")); + } if(!mktransform(dblimg, st, tr)) ERRX(_("Can't do given transform")); #ifdef EBUG st = get_imgstat(dblimg, NULL); #endif DBG("After transformation: MIN=%g, MAX=%g, AVR=%g, STD=%g", st->min, st->max, st->mean, st->std); - histogram *h = dbl2histogram(dblimg, G.nlvl); + green("Histogram after transformations:\n"); + h = dbl2histogram(dblimg, G.nlvl); print_histo(h); histogram_free(&h); uint8_t *colored = convert2palette(dblimg, colormap); diff --git a/examples/med.c b/examples/med.c new file mode 100644 index 0000000..22ea6ef --- /dev/null +++ b/examples/med.c @@ -0,0 +1,107 @@ +/* + * This file is part of the FITSmaniplib project. + * Copyright 2019 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 "common.h" + +/* + * Median filtering of image + */ + +typedef struct{ + char *fitsname; // input file name + char *outfile; // output file name + int rewrite; // rewrite existing file + int medr; // radius of filter +} glob_pars; + +/* + * here are global parameters initialisation + */ +static int help; +static glob_pars G = { + .medr = 1, +}; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +static myoption cmdlnopts[] = { +// common options + {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")}, + {"fitsname",NEED_ARG, NULL, 'i', arg_string, APTR(&G.fitsname), _("name of input file")}, + {"outpname",NEED_ARG, NULL, 'o', arg_string, APTR(&G.outfile), _("output file name (jpeg)")}, + {"rewrite", NO_ARGS, NULL, 'r', arg_none, APTR(&G.rewrite), _("rewrite output file")}, + {"radius", NEED_ARG, NULL, 'R', arg_int, APTR(&G.medr), _("radius of median (0 for cross 3x3)")}, + end_option +}; + +/** + * Parse command line options and return dynamically allocated structure + * to global parameters + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +static glob_pars *parse_args(int argc, char **argv){ + int i; + char *helpstring = "Usage: %%s [args]\n\n\tWhere args are:\n"; + change_helpstring(helpstring); + // parse arguments + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(argc > 0){ + for (i = 0; i < argc; i++) + printf("Ignore extra argument: %s\n", argv[i]); + } + return &G; +} + +int main(int argc, char *argv[]){ + initial_setup(); + parse_args(argc, argv); + if(!G.fitsname) ERRX(_("No input filename given!")); + if(!G.outfile) ERRX(_("No output filename given!")); + if(G.medr < 0) ERRX(_("Median radius should be >= 0")); + if(!file_absent(G.outfile) && !G.rewrite) ERRX(_("File %s exists"), G.outfile); + FITS *f = FITS_read(G.fitsname); + if(!f) ERRX(_("Failed to open %s"), G.fitsname); + f->curHDU = NULL; + int i; + for(i = 1; i <= f->NHDUs; ++i){ + if(f->HDUs[i].hdutype == IMAGE_HDU){f->curHDU = &f->HDUs[i]; break;} + } + if(!f->curHDU) ERRX(_("No image HDUs in %s"), G.fitsname); + green("First HDU with image: #%d\n", i); + FITSimage *img = f->curHDU->contents.image; + if(img->naxis != 2) ERRX(_("Support only 2-dimensional images")); + doubleimage *dblimg = image2double(img); + if(!dblimg) ERRX(_("Can't convert image to double")); + doubleimage *filtered = get_median(dblimg, G.medr); + if(!filtered) ERRX(_("WTF?")); + if(!image_rebuild(img, filtered->data)) ERRX(_("Can't rebuild image")); + doubleimage_free(&dblimg); + f->filename = G.outfile; + bool w = FALSE; + if(file_absent(G.outfile)) w = FITS_write(G.outfile, f); + else w = FITS_rewrite(f); + if(!w) ERRX(_("Can't write %s"), f->filename); + FITS_free(&f); + return 0; +} + diff --git a/fitsimages.c b/fitsimages.c index 663d3b0..66c8237 100644 --- a/fitsimages.c +++ b/fitsimages.c @@ -317,7 +317,7 @@ FITSimage *image_read(FITS *fits){ return img; } -void dblima_free(doubleimage **im){ +void doubleimage_free(doubleimage **im){ FREE((*im)->data); FREE(*im); } @@ -434,3 +434,18 @@ doubleimage *normalize_dbl(doubleimage *im, imgstat *st){ } return im; } + +/** + * @brief new_doubleimage - create image of double numbers + * @param w - width + * @param h - height + * @return empty image + */ +doubleimage *doubleimage_new(size_t w, size_t h){ + doubleimage *out = MALLOC(doubleimage, 1); + out->height = h; + out->width = w; + out->totpix = w*h; + out->data = MALLOC(double, out->totpix); + return out; +} diff --git a/locale/ru/LC_MESSAGES/FITSmanip.mo b/locale/ru/LC_MESSAGES/FITSmanip.mo index 9aee15b4036813ec3947054ebd520e27095767ee..ebc70126cbc2df031364353434b39cd99d611d3c 100644 GIT binary patch delta 647 zcmXZZOG{fp6u|M5G@AIR=4E`O#7lfYu>~Ur)FMS(xzkk_brZuSxG)JZ;zqCwS0z~m zS9MV+^b2Go2*MQuMrtByFh)=)-E`Ye(Ed-n2hRP?%$<8@=1kFd;y--yyN-nzCSQ?F za-1wUcn}_^NE~}`0^j2ZCUFoCa1;xu|M-q`c#Z@37l%+@h=d@7+JB48Ka^D(lMH;s zD1O3e{D!*FHO}I1G%)BAd4qBE;+E~M?Gdi9`D0w$R=`$9W!s{c1yFaZTFlJi!H_D z<*%-SBmA@eaA`H~jaK`rep#zmjOs)EqG8=tYSnVId)*nlZCI`9kNaA?V%%B{>!#9d Jv`ZOJY5P;#?^OAUe=Zr>;%3-1)T6jVrDvBQ<*jrjCSSY0M2MB3Wh?wfI5kUk? z!Qv`GQVI5AVyMX~ehAlZoo?K_(2saIsv8Tke^OZ_=AiWW2&Jy0\n" "Language-Team: LANGUAGE \n" @@ -17,16 +17,16 @@ msgstr "" "Content-Type: text/plain; charset=koi8-r\n" "Content-Transfer-Encoding: 8bit\n" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:96 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:95 #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fitsimages.c:427 msgid "Data range is too small" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:103 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:101 msgid "Given transform type not supported yet" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:267 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:265 msgid "Given colormap doesn't support yet" msgstr "" @@ -113,3 +113,11 @@ msgstr "" #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/histogram.c:122 msgid "Can't find bottom index" msgstr "" + +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/median.c:211 +msgid "Wrong parameters" +msgstr "" + +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/median.c:486 +msgid "Can't create output image" +msgstr "" diff --git a/locale/ru/ru.po b/locale/ru/ru.po index 9a3aa10..37791a3 100644 --- a/locale/ru/ru.po +++ b/locale/ru/ru.po @@ -7,7 +7,7 @@ msgid "" msgstr "Project-Id-Version: PACKAGE VERSION\n" "Report-Msgid-Bugs-To: \n" - "POT-Creation-Date: 2019-04-01 11:58+0300\n" + "POT-Creation-Date: 2019-04-02 11:55+0300\n" "PO-Revision-Date: YEAR-MO-DA HO:MI+ZONE\n" "Last-Translator: FULL NAME \n" "Language-Team: LANGUAGE \n" @@ -34,6 +34,10 @@ msgstr " msgid "Can't copy data" msgstr "Не могу скопировать данные" +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/median.c:486 +msgid "Can't create output image" +msgstr "Не могу создать выходное изображение" + #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/histogram.c:122 msgid "Can't find bottom index" msgstr "Не могу найти нижний индекс" @@ -73,7 +77,7 @@ msgstr " msgid "Can't write table %s!" msgstr "Не могу записать таблицу %s!" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:96 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:95 #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fitsimages.c:427 msgid "Data range is too small" msgstr "Диапазон данных слишком мал" @@ -83,11 +87,11 @@ msgstr " msgid "Found %d pixels with undefined value" msgstr "Найдено %d пикселей с неопределенными значениями" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:267 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:265 msgid "Given colormap doesn't support yet" msgstr "Данная палитра пока не поддерживается" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:103 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/FITSmanip.c:101 msgid "Given transform type not supported yet" msgstr "Данный тип преобразований пока не поддерживается" @@ -112,6 +116,10 @@ msgstr " msgid "Unknown HDU type" msgstr "Неизвестный тип HDU" +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/median.c:211 +msgid "Wrong parameters" +msgstr "Неправильные параметры" + #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fitstables.c:119 msgid "strdup() failed!" msgstr "Не удалось сделать strdup()!" diff --git a/median.c b/median.c new file mode 100644 index 0000000..5542a67 --- /dev/null +++ b/median.c @@ -0,0 +1,698 @@ +/* + * This file is part of the FITSmaniplib project. + * Copyright 2019 Edward V. Emelianov , . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +// FOR MEDIATOR: +// Copyright (c) 2011 ashelly.myopenid.com under +// FOR opt_medXX: +// Copyright (c) 1998 Nicolas Devillard. Public domain. +// FOR qickselect: +// "Numerical recipes in C", Second Edition, +// Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5 +// Code by Nicolas Devillard - 1998. Public domain. + +// TODO: resolve problem with borders + +#include "FITSmanip.h" +#include "local.h" + +// largest radius for adaptive median filter +#define LARGEST_ADPMED_RADIUS (3) + +#define ELEM_SWAP(a, b) {register double t = a; a = b; b = t;} +#define PIX_SORT(a, b) {if (p[a] > p[b]) ELEM_SWAP(p[a], p[b]);} + +/* + * simplest short functions for median calculation + */ +// even values are from "FAST, EFFICIENT MEDIAN FILTERS WITH EVEN LENGTH WINDOWS", +// J.P. HAVLICEK, K.A. SAKADY, G.R.KATZ +static double opt_med2(double *p){ + return (p[0] + p[1]) * 0.5; +} +static double opt_med3(double *p){ + PIX_SORT(0, 1); PIX_SORT(1, 2); PIX_SORT(0, 1); + return (p[1]); +} +static double opt_med4(double *p){ + PIX_SORT(0, 2); PIX_SORT(1, 3); + PIX_SORT(0, 1); PIX_SORT(2, 3); + return (p[1] + p[2]) / 2.; +} +static double opt_med5(double *p){ + PIX_SORT(0, 1); PIX_SORT(3, 4); PIX_SORT(0, 3); + PIX_SORT(1, 4); PIX_SORT(1, 2); PIX_SORT(2, 3) ; + PIX_SORT(1, 2); + return (p[2]); +} +static double opt_med6(double *p){ + PIX_SORT(1, 2); PIX_SORT(3, 4); + PIX_SORT(0, 1); PIX_SORT(2, 3); PIX_SORT(4, 5); + PIX_SORT(1, 2); PIX_SORT(3, 4); + PIX_SORT(0, 1); PIX_SORT(2, 3); PIX_SORT(4, 5); + PIX_SORT(1, 2); PIX_SORT(3, 4); + return (p[2] + p[3]) / 2.; +} +static double opt_med7(double *p){ + PIX_SORT(0, 5); PIX_SORT(0, 3); PIX_SORT(1, 6); + PIX_SORT(2, 4); PIX_SORT(0, 1); PIX_SORT(3, 5); + PIX_SORT(2, 6); PIX_SORT(2, 3); PIX_SORT(3, 6); + PIX_SORT(4, 5); PIX_SORT(1, 4); PIX_SORT(1, 3); + PIX_SORT(3, 4); + return (p[3]); +} +// optimal Batcher's sort for 8 elements (http://myopen.googlecode.com/svn/trunk/gtkclient_tdt/include/fast_median.h) +static double opt_med8(double *p){ + PIX_SORT(0, 4); PIX_SORT(1, 5); PIX_SORT(2, 6); + PIX_SORT(3, 7); PIX_SORT(0, 2); PIX_SORT(1, 3); + PIX_SORT(4, 6); PIX_SORT(5, 7); PIX_SORT(2, 4); + PIX_SORT(3, 5); PIX_SORT(0, 1); PIX_SORT(2, 3); + PIX_SORT(4, 5); PIX_SORT(6, 7); PIX_SORT(1, 4); + PIX_SORT(3, 6); + return (p[3] + p[4]) / 2.; +} +static double opt_med9(double *p){ + PIX_SORT(1, 2); PIX_SORT(4, 5); PIX_SORT(7, 8); + PIX_SORT(0, 1); PIX_SORT(3, 4); PIX_SORT(6, 7); + PIX_SORT(1, 2); PIX_SORT(4, 5); PIX_SORT(7, 8); + PIX_SORT(0, 3); PIX_SORT(5, 8); PIX_SORT(4, 7); + PIX_SORT(3, 6); PIX_SORT(1, 4); PIX_SORT(2, 5); + PIX_SORT(4, 7); PIX_SORT(4, 2); PIX_SORT(6, 4); + PIX_SORT(4, 2); + return (p[4]); +} +static double opt_med16(double *p){ + PIX_SORT(0, 8); PIX_SORT(1, 9); PIX_SORT(2, 10); PIX_SORT(3, 11); + PIX_SORT(4, 12); PIX_SORT(5, 13); PIX_SORT(6, 14); PIX_SORT(7, 15); + PIX_SORT(0, 4); PIX_SORT(1, 5); PIX_SORT(2, 6); PIX_SORT(3, 7); + PIX_SORT(8, 12); PIX_SORT(9, 13); PIX_SORT(10, 14); PIX_SORT(11, 15); + PIX_SORT(4, 8); PIX_SORT(5, 9); PIX_SORT(6, 10); PIX_SORT(7, 11); + PIX_SORT(0, 2); PIX_SORT(1, 3); PIX_SORT(4, 6); PIX_SORT(5, 7); + PIX_SORT(8, 10); PIX_SORT(9, 11); PIX_SORT(12, 14); PIX_SORT(13, 15); + PIX_SORT(2, 8); PIX_SORT(3, 9); PIX_SORT(6, 12); PIX_SORT(7, 13); + PIX_SORT(2, 4); PIX_SORT(3, 5); PIX_SORT(6, 8); PIX_SORT(7, 9); + PIX_SORT(10, 12); PIX_SORT(11, 13); PIX_SORT(0, 1); PIX_SORT(2, 3); + PIX_SORT(4, 5); PIX_SORT(6, 7); PIX_SORT(8, 9); PIX_SORT(10, 11); + PIX_SORT(12, 13); PIX_SORT(14, 15); PIX_SORT(1, 8); PIX_SORT(3, 10); + PIX_SORT(5, 12); PIX_SORT(7, 14); PIX_SORT(5, 8); PIX_SORT(7, 10); + return (p[7] + p[8]) / 2.; +} +static double opt_med25(double *p){ + PIX_SORT(0, 1) ; PIX_SORT(3, 4) ; PIX_SORT(2, 4) ; + PIX_SORT(2, 3) ; PIX_SORT(6, 7) ; PIX_SORT(5, 7) ; + PIX_SORT(5, 6) ; PIX_SORT(9, 10) ; PIX_SORT(8, 10) ; + PIX_SORT(8, 9) ; PIX_SORT(12, 13); PIX_SORT(11, 13) ; + PIX_SORT(11, 12); PIX_SORT(15, 16); PIX_SORT(14, 16) ; + PIX_SORT(14, 15); PIX_SORT(18, 19); PIX_SORT(17, 19) ; + PIX_SORT(17, 18); PIX_SORT(21, 22); PIX_SORT(20, 22) ; + PIX_SORT(20, 21); PIX_SORT(23, 24); PIX_SORT(2, 5) ; + PIX_SORT(3, 6) ; PIX_SORT(0, 6) ; PIX_SORT(0, 3) ; + PIX_SORT(4, 7) ; PIX_SORT(1, 7) ; PIX_SORT(1, 4) ; + PIX_SORT(11, 14); PIX_SORT(8, 14) ; PIX_SORT(8, 11) ; + PIX_SORT(12, 15); PIX_SORT(9, 15) ; PIX_SORT(9, 12) ; + PIX_SORT(13, 16); PIX_SORT(10, 16); PIX_SORT(10, 13) ; + PIX_SORT(20, 23); PIX_SORT(17, 23); PIX_SORT(17, 20) ; + PIX_SORT(21, 24); PIX_SORT(18, 24); PIX_SORT(18, 21) ; + PIX_SORT(19, 22); PIX_SORT(8, 17) ; PIX_SORT(9, 18) ; + PIX_SORT(0, 18) ; PIX_SORT(0, 9) ; PIX_SORT(10, 19) ; + PIX_SORT(1, 19) ; PIX_SORT(1, 10) ; PIX_SORT(11, 20) ; + PIX_SORT(2, 20) ; PIX_SORT(2, 11) ; PIX_SORT(12, 21) ; + PIX_SORT(3, 21) ; PIX_SORT(3, 12) ; PIX_SORT(13, 22) ; + PIX_SORT(4, 22) ; PIX_SORT(4, 13) ; PIX_SORT(14, 23) ; + PIX_SORT(5, 23) ; PIX_SORT(5, 14) ; PIX_SORT(15, 24) ; + PIX_SORT(6, 24) ; PIX_SORT(6, 15) ; PIX_SORT(7, 16) ; + PIX_SORT(7, 19) ; PIX_SORT(13, 21); PIX_SORT(15, 23) ; + PIX_SORT(7, 13) ; PIX_SORT(7, 15) ; PIX_SORT(1, 9) ; + PIX_SORT(3, 11) ; PIX_SORT(5, 17) ; PIX_SORT(11, 17) ; + PIX_SORT(9, 17) ; PIX_SORT(4, 10) ; PIX_SORT(6, 12) ; + PIX_SORT(7, 14) ; PIX_SORT(4, 6) ; PIX_SORT(4, 7) ; + PIX_SORT(12, 14); PIX_SORT(10, 14); PIX_SORT(6, 7) ; + PIX_SORT(10, 12); PIX_SORT(6, 10) ; PIX_SORT(6, 17) ; + PIX_SORT(12, 17); PIX_SORT(7, 17) ; PIX_SORT(7, 10) ; + PIX_SORT(12, 18); PIX_SORT(7, 12) ; PIX_SORT(10, 18) ; + PIX_SORT(12, 20); PIX_SORT(10, 20); PIX_SORT(10, 12) ; + return (p[12]); +} +#undef PIX_SORT + +#define PIX_SORT(a, b) {if (a > b) ELEM_SWAP(a, b);} +/** + * @brief quick_select - algorithm for approximate median calculation for array idata of size n + * @param idata (i) - input data array + * @param n - size of `idata` + * @return median value + */ +double quick_select(const double *idata, int n){ + int low, high; + int median; + int middle, ll, hh; + double *arr = MALLOC(double, n); + memcpy(arr, idata, n*sizeof(double)); + low = 0 ; high = n-1 ; median = (low + high) / 2; + for(;;){ + if(high <= low) // One element only + break; + if(high == low + 1){ // Two elements only + PIX_SORT(arr[low], arr[high]) ; + break; + } + // Find median of low, middle and high doubles; swap into position low + middle = (low + high) / 2; + PIX_SORT(arr[middle], arr[high]) ; + PIX_SORT(arr[low], arr[high]) ; + PIX_SORT(arr[middle], arr[low]) ; + // Swap low double (now in position middle) into position (low+1) + ELEM_SWAP(arr[middle], arr[low+1]) ; + // Nibble from each end towards middle, swapping doubles when stuck + ll = low + 1; + hh = high; + for(;;){ + do ll++; while (arr[low] > arr[ll]); + do hh--; while (arr[hh] > arr[low]); + if(hh < ll) break; + ELEM_SWAP(arr[ll], arr[hh]) ; + } + // Swap middle double (in position low) back into correct position + ELEM_SWAP(arr[low], arr[hh]) ; + // Re-set active partition + if (hh <= median) low = ll; + if (hh >= median) high = hh - 1; + } + double ret = arr[median]; + FREE(arr); + return ret; +} +#undef PIX_SORT +#undef ELEM_SWAP + +/** + * @brief calc_median - calculate median of array idata with size n + * the specific type of algorythm is choosen according to `n` + * @param idata (i) - input data array + * @param n - size of array `idata` + * @return median value + */ +double calc_median(const double *idata, int n){ + if(!idata || n < 1){ + WARNX(_("Wrong parameters")); + return 0.; + } + typedef double (*medfunc)(double *p); + medfunc fn = NULL; + const medfunc fnarr[] = {opt_med2, opt_med3, opt_med4, opt_med5, opt_med6, + opt_med7, opt_med8, opt_med9}; + if(n == 1) return *idata; + if(n < 10) fn = fnarr[n - 2]; + else if(n == 16) fn = opt_med16; + else if(n == 25) fn = opt_med25; + if(fn){ + // copy data to new buffer - `idata` should leave unchanged + double *dataarr = MALLOC(double, n); + memcpy(dataarr, idata, sizeof(double)*n); + double medval = fn(dataarr); + FREE(dataarr); + return medval; + }else{ + return quick_select(idata, n); + } +} + +#define doubleLess(a,b) ((a)<(b)) +#define doubleMean(a,b) (((a)+(b))/2) + +typedef struct Mediator_t{ + double* data; // circular queue of values + int* pos; // index into `heap` for each value + int* heap; // max/median/min heap holding indexes into `data`. + int N; // allocated size. + int idx; // position in circular queue + int ct; // count of doubles in queue +} Mediator; + +/*--- Helper Functions ---*/ + +#define minCt(m) (((m)->ct-1)/2) //count of doubles in minheap +#define maxCt(m) (((m)->ct)/2) //count of doubles in maxheap + +//returns 1 if heap[i] < heap[j] +static inline int mmless(Mediator* m, int i, int j){ + return doubleLess(m->data[m->heap[i]],m->data[m->heap[j]]); +} + +//swaps doubles i&j in heap, maintains indexes +static inline int mmexchange(Mediator* m, int i, int j){ + int t = m->heap[i]; + m->heap[i] = m->heap[j]; + m->heap[j] = t; + m->pos[m->heap[i]] = i; + m->pos[m->heap[j]] = j; + return 1; +} + +//swaps doubles i&j if i1 && i < minCt(m) && mmless(m, i+1, i)) ++i; + if(!mmCmpExch(m,i,i/2)) break; + } +} + +//maintains maxheap property for all doubles below i/2. (negative indexes) +static void maxSortDown(Mediator* m, int i){ + for(; i >= -maxCt(m); i*=2){ + if(i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) --i; + if(!mmCmpExch(m,i/2,i)) break; + } +} + +//maintains minheap property for all doubles above i, including median +//returns true if median changed +static int minSortUp(Mediator* m, int i){ + while (i > 0 && mmCmpExch(m, i, i/2)) i /= 2; + return (i == 0); +} + +//maintains maxheap property for all doubles above i, including median +//returns true if median changed +static int maxSortUp(Mediator* m, int i){ + while (i < 0 && mmCmpExch(m, i/2, i)) i /= 2; + return (i == 0); +} + +/*--- Public Interface ---*/ + +//creates new Mediator: to calculate `ndoubles` running median. +//mallocs single block of memory, caller must free. +static Mediator* MediatorNew(int ndoubles){ + int size = sizeof(Mediator) + ndoubles*(sizeof(double)+sizeof(int)*2); + Mediator* m = malloc(size); + m->data = (double*)(m + 1); + m->pos = (int*) (m->data + ndoubles); + m->heap = m->pos + ndoubles + (ndoubles / 2); //points to middle of storage. + m->N = ndoubles; + m->ct = m->idx = 0; + while (ndoubles--){ //set up initial heap fill pattern: median,max,min,max,... + m->pos[ndoubles] = ((ndoubles+1)/2) * ((ndoubles&1)? -1 : 1); + m->heap[m->pos[ndoubles]] = ndoubles; + } + return m; +} + +//Inserts double, maintains median in O(lg ndoubles) +static void MediatorInsert(Mediator* m, double v){ + int isNew=(m->ctN); + int p = m->pos[m->idx]; + double old = m->data[m->idx]; + m->data[m->idx]=v; + m->idx = (m->idx+1) % m->N; + m->ct+=isNew; + if(p>0){ //new double is in minHeap + if (!isNew && doubleLess(old,v)) minSortDown(m,p*2); + else if (minSortUp(m,p)) maxSortDown(m,-1); + }else if (p<0){ //new double is in maxheap + if (!isNew && doubleLess(v,old)) maxSortDown(m,p*2); + else if (maxSortUp(m,p)) minSortDown(m, 1); + }else{ //new double is at median + if (maxCt(m)) maxSortDown(m,-1); + if (minCt(m)) minSortDown(m, 1); + } +} + +//returns median double (or average of 2 when double count is even) +static double MediatorMedian(Mediator* m){ + double v = m->data[m->heap[0]]; + if ((m->ct&1) == 0) v = doubleMean(v, m->data[m->heap[-1]]); + return v; +} +/* +// median + min/max +static double MediatorStat(Mediator* m, double *minval, double *maxval){ + double v= m->data[m->heap[0]]; + if ((m->ct&1) == 0) v = doubleMean(v,m->data[m->heap[-1]]); + double min = v, max = v; + int i; + for(i = -maxCt(m); i < 0; ++i){ + int v = m->data[m->heap[i]]; + if(v < min) min = v; + } + *minval = min; + for(i = 1; i <= minCt(m); ++i){ + int v = m->data[m->heap[i]]; + if(v > max) max = v; + } + *maxval = max; + return v; +}*/ + +// TODO: add adaptive filtering +/** + * @brief get_adp_median_cross - adaptive median filter by cross 3x3 + * We have 5 datapoints and 4 inserts @ each step, so it's better to use opt_med5 instead of Mediator + * @param img (i) - input image + * @param out (o) - output image (allocated outside) + * @param adp - TRUE for adaptive filtering and FALSE for regular + */ +static void get_adp_median_cross(const doubleimage *img, doubleimage *out, _U_ bool adp){ + size_t w = img->width, h = img->height; + double *med = out->data, *inputima = img->data, *iptr; +#ifdef EBUG + double t0 = dtime(); +#endif + OMP_FOR() + for(size_t x = 1; x < w - 1; ++x){ + double buffer[5]; + size_t curpix = x + w, // index of current pixel image arrays + y, ymax = h - 1; + for(y = 1; y < ymax; ++y, curpix += w){ + double md, *I = &inputima[curpix]; //, Ival = *I; + memcpy(buffer, I - 1, 3*sizeof(double)); + buffer[3] = I[-w]; buffer[4] = I[w]; + md = opt_med5(buffer); + /* + if(adp){ + double s, l; + s = DBL_EPSILON + MIN(buffer[0], buffer[1]); + l = MAX(buffer[3], buffer[4]) - DBL_EPSILON; + if(s < md && md < l){ + if(s < Ival && Ival < l) med[curpix] = Ival; + else med[curpix] = md; + }else{ + med[curpix] = adp_med_5by5(img, x, y); + } + }else */ + med[curpix] = md; + } + } + // process borders & corners (without adaptive) + double buf[5]; + // left top + buf[0] = inputima[0]; buf[1] = inputima[0]; + buf[2] = inputima[1]; buf[3] = inputima[w]; + buf[4] = inputima[w + 1]; + med[0] = opt_med5(buf); + // right top + iptr = &inputima[w - 1]; + buf[0] = iptr[0]; buf[1] = iptr[0]; + buf[2] = iptr[-1]; buf[3] = iptr[w - 1]; + buf[4] = iptr[w]; + med[w - 1] = opt_med5(buf); + // left bottom + iptr = &inputima[(h - 1) * w]; + buf[0] = iptr[0]; buf[1] = iptr[0]; + buf[2] = iptr[-w]; buf[3] = iptr[1 - w]; + buf[4] = iptr[1]; + med[(h - 1) * w] = opt_med5(buf); + // right bottom + iptr = &inputima[h * w - 1]; + buf[0] = iptr[0]; buf[1] = iptr[0]; + buf[2] = iptr[-w-1]; buf[3] = iptr[-w]; + buf[4] = iptr[-1]; + med[h * w - 1] = opt_med5(buf); + // process borders without corners + // top + OMP_FOR(shared(med)) + for(size_t x = 1; x < w - 1; ++x){ + double *iptr = &inputima[x]; + buf[0] = buf[1] = *iptr; + buf[2] = iptr[-1]; buf[3] = iptr[2]; + buf[4] = iptr[w]; + med[x] = opt_med5(buf); + } + // bottom + size_t curidx = (h-2)*w; + OMP_FOR(shared(curidx, med)) + for(size_t x = 1; x < w - 1; --x){ + double *iptr = &inputima[curidx + x]; + buf[0] = buf[1] = *iptr; + buf[2] = iptr[-w]; buf[3] = iptr[-1]; + buf[4] = iptr[1]; + med[curidx + x] = opt_med5(buf); + } + // left + OMP_FOR(shared(med)) + for(size_t y = 1; y < h - 1; ++y){ + size_t cur = y * w; + double *iptr = &inputima[cur]; + buf[0] = buf[1] = *iptr; + buf[2] = iptr[-w]; buf[3] = iptr[1]; + buf[4] = iptr[w]; + med[cur] = opt_med5(buf); + } + // right + curidx = w - 1; + OMP_FOR(shared(curidx, med)) + for(size_t y = 1; y < h - 1; ++y){ + size_t cur = curidx + y * w; + double *iptr = &inputima[cur]; + buf[0] = buf[1] = *iptr; + buf[2] = iptr[-w]; buf[3] = iptr[-1]; + buf[4] = iptr[w]; + med[cur] = opt_med5(buf); + } + DBG("time for median filtering by cross 3x3 of image %zdx%zd: %gs", w, h, + dtime() - t0); +} + +// TODO: add borders and corners +/** + * @brief get_median - filter image by median (radius*2 + 1) x (radius*2 + 1) + * @param img (i) - input image + * @param radius - zone radius (0 for cross 3x3) + * @return image filtered by median (allocated here) + */ +doubleimage *get_median(const doubleimage *img, size_t radius){ + size_t w = img->width, h = img->height; + doubleimage *out = doubleimage_new(img->width, img->height); + if(!out){ + WARNX(_("Can't create output image")); + return NULL; + } + memcpy(out->data, img->data, sizeof(double)*img->totpix); + double *med = out->data, *inputima = img->data; + if(radius == 0){ + get_adp_median_cross(img, out, 0); + return out; + } + + size_t blksz = radius * 2 + 1, fullsz = blksz * blksz; +#ifdef EBUG + double t0 = dtime(); +#endif + OMP_FOR(shared(inputima, med)) + for(size_t x = radius; x < w - radius; ++x){ + size_t xx, yy, xm = x + radius + 1, y, ymax = blksz - 1, xmin = x - radius; + Mediator* m = MediatorNew(fullsz); + // initial fill + for(yy = 0; yy < ymax; ++yy) + for(xx = xmin; xx < xm; ++xx) + MediatorInsert(m, inputima[xx + yy*w]); + ymax = 2*radius*w; + xmin += ymax; + xm += ymax; + ymax = h - radius; + size_t medidx = x + radius * w; + for(y = radius; y < ymax; ++y, xmin += w, xm += w, medidx += w){ + for(xx = xmin; xx < xm; ++xx) + MediatorInsert(m, inputima[xx]); + med[medidx] = MediatorMedian(m); + } + FREE(m); + } + DBG("time for median filtering %zdx%zd of image %zdx%zd: %gs", blksz, blksz, w, h, + dtime() - t0); + return out; +} + +#if 0 + +/** + * procedure for finding median value in window 5x5 + * PROBLEM: bounds + */ +static double adp_med_5by5(const IMAGE *img, size_t x, size_t y){ + size_t blocklen, w = img->width, h = img->height, yy, _2w = 2 * w; + double arr[25], *arrptr = arr, *dataptr, *currpix; + int position = ((x < 1) ? 1 : 0) // left columns + + ((x > w - 2) ? 2 : 0) // right columns + + ((y < 1) ? 4 : 0) // top rows + + ((y > w - 2) ? 8 : 0); // bottom rows + /* Now by value of "position" we know where is the point: + *************************** + * 5 * 4 * 6 * + *************************** + * * * * + * * * * + * 1 * 0 * 2 * + * * * * + * * * * + *************************** + * 9 * 8 *10 * + ***************************/ + currpix = &img->data[x + y * w]; // pointer to current pixel + dataptr = currpix - _2w - 2; // pointer to left upper corner of 5x5 square + inline void copy5times(double val){ + for(int i = 0; i < 5; ++i) *arrptr++ = val; + } + inline void copy9times(double val){ + for(int i = 0; i < 9; ++i) *arrptr++ = val; + } + void copycolumn(double *startpix){ + for(int i = 0; i < 5; ++i, startpix += w) *arrptr++ = *startpix; + } + inline void copyvertblock(size_t len){ + for(int i = 0; i < 5; ++i, dataptr += w, arrptr += len) + memcpy(arrptr, dataptr, len * sizeof(double)); + } + inline void copyhorblock(size_t len){ + for(size_t i = 0; i < len; ++i, dataptr += w, arrptr += 5) + memcpy(arrptr, dataptr, 5 * sizeof(double)); + } + inline void copyblock(){ + for(size_t i = 0; i < 4; ++i, dataptr += w, arrptr += 4) + memcpy(arrptr, dataptr, 4 * sizeof(double)); + } + switch(position){ + case 1: // left + copy5times(*currpix); // make 5 copies of current pixel + if(x == 0){ // copy 1st column too + dataptr += 2; + copycolumn(dataptr); + blocklen = 3; + }else{ // 2nd column - no copy need + ++dataptr; + blocklen = 4; + } + copyvertblock(blocklen); + break; + case 2: // right + copy5times(*currpix); + if(x == w - 1){ // copy last column too + copycolumn(dataptr + 2); + blocklen = 3; + }else{ // 2nd column - no copy need + blocklen = 4; + } + copyvertblock(blocklen); + break; + case 4: // top + copy5times(*currpix); + if(y == 0){ + dataptr += _2w; + memcpy(arrptr, dataptr, 5 * sizeof(double)); + blocklen = 3; + }else{ + dataptr += w; + blocklen = 4; + } + copyhorblock(blocklen); + break; + case 8: // bottom + copy5times(*currpix); + if(y == h - 1){ + memcpy(arrptr, dataptr + _2w, 5 * sizeof(double)); + blocklen = 3; + }else{ + blocklen = 4; + } + copyhorblock(blocklen); + break; + case 5: // top left corner: in all corners we just copy 4x4 square & 9 times this pixel + copy9times(*currpix); + dataptr = img->data; + copyblock(); + break; + case 6: // top right corner + copy9times(*currpix); + dataptr = &img->data[w - 4]; + copyblock(); + break; + case 9: // bottom left cornet + copy9times(*currpix); + dataptr = &img->data[(y - 4) * w]; + copyblock(); + break; + case 10: // bottom right cornet + copy9times(*currpix); + dataptr = &img->data[(y - 3) * w - 4]; + copyblock(); + break; + default: // 0 + for(yy = 0; yy < 5; ++yy, dataptr += w, arrptr += 5) + memcpy(arrptr, dataptr, 5*sizeof(double)); + } + return opt_med25(arr); +} + + +/** + * filter image by median (radius*2 + 1) x (radius*2 + 1) + */ +doubleimage *get_adaptive_median(const doubleimage *img, size_t radius);{ + int radius = f->w; + size_t w = img->width, h = img->height, siz = w*h, bufsiz = siz*sizeof(double); + IMAGE *out = similarFITS(img, img->dtype); + double *med = out->data, *inputima = img->data; + memcpy(med, inputima, bufsiz); + if(radius == 0){ + get_adp_median_cross(img, out, 1); + return out; + } + size_t blksz = radius * 2 + 1, fullsz = blksz * blksz; +#ifdef EBUG + double t0 = dtime(); +#endif + OMP_FOR(shared(inputima, med)) + for(size_t x = radius; x < w - radius; ++x){ + size_t xx, yy, xm = x + radius + 1, y, ymax = blksz - 1, xmin = x - radius; + Mediator* m = MediatorNew(fullsz); + // initial fill + for(yy = 0; yy < ymax; ++yy) + for(xx = xmin; xx < xm; ++xx) + MediatorInsert(m, inputima[xx + yy*w]); + ymax = 2*radius*w; + xmin += ymax; + xm += ymax; + ymax = h - radius; + size_t curpos = x + radius * w; + for(y = radius; y < ymax; ++y, xmin += w, xm += w, curpos += w){ + for(xx = xmin; xx < xm; ++xx) + MediatorInsert(m, inputima[xx]); + double s, l, md, I = inputima[curpos]; + md = MediatorStat(m, &s, &l); + s += ITM_EPSILON, l -= ITM_EPSILON; + if(s < md && md < l){ + if(s < I && I < l) med[curpos] = I; + else med[curpos] = md; + }else{ + if(radius > LARGEST_ADPMED_RADIUS) + med[curpos] = I; + else + med[curpos] = adp_med_5by5(img, x, y); + } + } + FREE(m); + } + DBG("time for adadptive median filtering %zdx%zd of image %zdx%zd: %gs", blksz, blksz, w, h, + dtime() - t0); + return out; +} +#endif