Add simple median filtering, TODO: create adaptive median

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;

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
+ *                                     median.c                                       *
-typedef FITS* (*imfuncptr)(FITS *in, Filter *f, Itmarray *i);
+ **************************************************************************************/
-*/
+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);
-// FilterType (not only convolution!)
+double calc_median(const double *idata, int n);
-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;
-*/
 
 #endif // FITSMANIP_H__

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)
+add_executable(imstat imstat.c)
+add_executable(keylist keylist.c)
+add_executable(listtable listtable.c)
+add_executable(median med.c)

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

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':
+            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);

examples/med.c

@@ -0,0 +1,107 @@
+/*
+ * This file is part of the FITSmaniplib project.
+ * Copyright 2019  Edward V. Emelianov <edward.emelianoff@gmail.com>, <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 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 "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;
+}
+

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;
+}

locale/ru/messages.po

@@ -8,7 +8,7 @@ msgid ""
 msgstr ""
 "Project-Id-Version: PACKAGE VERSION\n"
 "Report-Msgid-Bugs-To: \n"
-"POT-Creation-Date: 2019-04-01 13:42+0300\n"
+"POT-Creation-Date: 2019-04-02 13:36+0300\n"
 "PO-Revision-Date: YEAR-MO-DA HO:MI+ZONE\n"
 "Last-Translator: FULL NAME <EMAIL@ADDRESS>\n"
 "Language-Team: LANGUAGE <LL@li.org>\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 ""

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 <EMAIL@ADDRESS>\n"
         "Language-Team: LANGUAGE <LL@li.org>\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()!"

median.c

@@ -0,0 +1,698 @@
+/*
+ * This file is part of the FITSmaniplib project.
+ * Copyright 2019  Edward V. Emelianov <edward.emelianoff@gmail.com>, <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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// FOR MEDIATOR:
+// Copyright (c) 2011 ashelly.myopenid.com under <http://www.opensource.org/licenses/mit-license>
+// 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 i<j; returns true if swapped
+static inline int mmCmpExch(Mediator* m, int i, int j){
+	return (mmless(m,i,j) && mmexchange(m,i,j));
+}
+
+//maintains minheap property for all doubles below i/2.
+static void minSortDown(Mediator* m, int i){
+	for(; i <= minCt(m); i*=2){
+		if(i>1 && i < minCt(m) && mmless(m, i+1, i)) ++i;
+		if(!mmCmpExch(m,i,i/2)) break;
+	}
+}
+
+//maintains maxheap property for all 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->ct<m->N);
+	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