start adding image manipulation functions

FITSmanip.c

@@ -50,3 +50,128 @@ void FITS_reporterr(int *errcode){
         fprintf(stderr, COLOR_OLD);
     *errcode = 0;
 }
+
+/**
+ * @brief mktransform - make image intensity transformation
+ * @param dimg (io) - double image
+ * @param st   (i)  - image statistics
+ * @param transf    - type of transformation
+ * @return NULL if failed
+ *      Be carefull: image should be equalized before some types of transform
+ */
+doubleimage *mktransform(doubleimage *im, imgstat *st, intens_transform transf){
+    if(!im || !im->data || !st) return NULL;
+    double max = st->max, min = st->min;
+    if((max-min) < 2.*DBL_EPSILON){
+        WARNX(_("Data range is too small"));
+        return NULL;
+    }
+    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);
+    double logtrans(double in){ // logaryphmic
+        return log(1. + in);
+    }
+    double exptrans(double in){ // exponential
+        return exp(in);
+    }
+    double powtrans(double in){ // x^2
+        return in*in;
+    }
+    double sqrtrans(double in){ // square root
+        return sqrt(in);
+    }
+    switch(transf){
+        case TRANSF_EXP:
+            transfn = exptrans;
+        break;
+        case TRANSF_LOG:
+            transfn = logtrans;
+        break;
+        case TRANSF_POW:
+            transfn = powtrans;
+        break;
+        case TRANSF_SQR:
+            transfn = sqrtrans;
+        break;
+        default: return NULL;
+    }
+    size_t totpix = im->totpix;
+    OMP_FOR()
+    for(size_t i = 0; i < totpix; ++i){
+        double d = dimg[i] - min;
+        dimg[i] = transfn(d);
+    }
+    return im;
+}
+
+/**
+ * @brief palette_gray - simplest gray conversion
+ * @param gray  - nornmalized double value
+ * @param rgb - red, green and blue components
+ */
+static void palette_gray(double gray, uint8_t *rgb){
+    rgb[0] = rgb[1] = rgb[2] = (uint8_t)(255.*gray);
+}
+
+/**
+ * @brief palette_BR - palette from blue to red
+ * @param gray  - nornmalized double value
+ * @param rgb - red, green and blue components
+ */
+static void palette_BR(double gray, uint8_t *rgb){
+	int i = (int)(gray * 4.);
+	double x = gray - (double)i * .25;
+	uint8_t r = 0, g = 0, b = 0;
+	switch(i){
+		case 0:
+			g = (uint8_t)(255. * x);
+			b = 255;
+		break;
+		case 1:
+			g = 255;
+			b = (uint8_t)(255. * (1. - x));
+		break;
+		case 2:
+			r = (uint8_t)(255. * x);
+			g = 255;
+		break;
+		case 3:
+			r = 255;
+			g = (uint8_t)(255. * (1. - x));
+		break;
+		default:
+			r = 255;
+	}
+	rgb[0] = r;
+	rgb[1] = g;
+	rgb[2] = b;
+}
+
+typedef void (*palette)(double, uint8_t[3]); // pointer to palette function
+static palette palette_F[PALETTE_COUNT] = {
+    [PALETTE_GRAY] = palette_gray,
+    [PALETTE_BR] = palette_BR
+};
+
+/**
+ * @brief convert2palette - convert normalized double image into colour using some palette
+ * @param im (i) - image to convert
+ * @param cmap   - palette (colormap) used
+ * @return allocated here array with color image
+ */
+uint8_t *convert2palette(doubleimage *im, image_palette cmap){
+    if(!im || !im->data || cmap < 0 || cmap >= PALETTE_COUNT) return NULL;
+    palette impalette = palette_F[cmap];
+    size_t totpix = im->totpix;
+    if(totpix == 0) return NULL;
+    double *inarr = im->data;
+    uint8_t *colored = MALLOC(uint8_t, totpix * 3);
+    initomp();
+    OMP_FOR()
+    for(size_t i = 0; i < totpix; ++i){
+        impalette(inarr[i], &colored[i*3]);
+    }
+    return colored;
+}

FITSmanip.h

@@ -102,6 +102,40 @@ typedef struct{
 	void *data; 	    // picture data
 } FITSimage;
 
+// 2-dimensional image data as double value
+typedef struct{
+    size_t height;
+    size_t width;
+    size_t totpix;
+    double *data;
+} doubleimage;
+
+// simplest statistics
+typedef struct{
+    double mean;
+    double std;
+    double min;
+    double max;
+} imgstat;
+
+// type of intensity transformation
+typedef enum{
+    TRANSF_WRONG = 0,
+    TRANSF_LINEAR,
+    TRANSF_LOG,
+    TRANSF_EXP,
+    TRANSF_POW,
+    TRANSF_SQR,
+    TRANSF_HISTEQ
+} intens_transform;
+
+// different color maps
+typedef enum{
+    PALETTE_GRAY = 0,
+    PALETTE_BR,
+    PALETTE_COUNT
+} image_palette;
+
 typedef union{
     FITSimage *image;
     FITStable *table;
@@ -163,7 +197,10 @@ 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);
-double *image2double(FITSimage *img);
+void dblima_free(doubleimage **im);
+doubleimage *image2double(FITSimage *img);
+imgstat *get_imgstat(const doubleimage *dimg, imgstat *est);
+doubleimage *normalize_dbl(doubleimage *dimg, imgstat *st);
 //FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata);
 
 /**************************************************************************************
@@ -176,13 +213,15 @@ FITS *FITS_open(char *filename);
 bool FITS_write(char *filename, FITS *fits);
 bool FITS_rewrite(FITS *fits);
 char* make_filename(char *buff, size_t buflen, char *prefix, char *suffix);
-bool file_is_absent(char *name);
+bool file_absent(char *name);
 
 /**************************************************************************************
  *                                   FITSmanip.c                                      *
  **************************************************************************************/
 void FITS_reporterr(int *errcode);
 void initomp();
+doubleimage *mktransform(doubleimage *im, imgstat *st, intens_transform transf);
+uint8_t *convert2palette(doubleimage *im, image_palette cmap);
 
 /*
 // pointer to image conversion function

examples/CMakeLists.txt

@@ -7,3 +7,5 @@ link_libraries(FITSmanip cfitsio m)
 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)

examples/Readme.md

@@ -4,6 +4,22 @@ Examples
 ## common.h
 Common files for all
 
+# gd.c
+
+Usage: gd [args]
+
+        Where args are:
+
+  -T, --transform=arg   type of intensity transformation (log, sqr, exp, pow)
+  -h, --help            show this help
+  -i, --inname=arg      name of input file
+  -n, --hdunumber=arg   open image from given HDU number
+  -o, --outpname=arg    output file name (jpeg)
+  -p, --palette=arg     convert as given palette
+  -r, --rewrite         rewrite output file
+  -t, --textline=arg    add text line to output image (at bottom)
+
+
 ## imstat.c
 
 Usage: imstat [args] input files
@@ -19,7 +35,6 @@ Get statistics and modify images from first image HDU of each input file
 
 ## keylist.c
 
-
 Usage: keylist [args] infile.fits
 
         Where args are:

examples/gd.c

@@ -0,0 +1,245 @@
+/*
+ * 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"
+#include <gd.h>
+
+/*
+ * Read FITS image, convert it to double and save as JPEG
+ * with given pallette.
+ * WARNING! Supports only 2-dimensional images
+ */
+
+typedef struct{
+    char *fitsname;     // input file name
+    char *outfile;      // output file name
+    char *text;         // add this text to image
+    char *transform;    // type of intensity transform
+    char *palette;      // palette to convert FITS image
+    int   nhdu;         // HDU number to read image from
+    int   rewrite;      // rewrite output file
+} glob_pars;
+
+/*
+ * here are global parameters initialisation
+ */
+static int help;
+static glob_pars G = {
+    .nhdu = 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")},
+    {"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")},
+    {"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)")},
+    {"rewrite", NO_ARGS,    NULL,   'r',    arg_none,   APTR(&G.rewrite),   _("rewrite output file")},
+    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;
+}
+
+/**
+ * @brief write_jpeg - save JPEG file
+ * @param fname - file name
+ * @param data  - RGB array with data (order: R,G,B)
+ * @param str   - string to add @ left top corner of image
+ * @param img   - image structure (for size parameters)
+ * @return true if all OK
+ */
+static bool write_jpeg(const char *fname, const uint8_t *data, const char *str, FITSimage *img){
+    if(!img) return FALSE;
+    int W = img->naxes[0], H = img->naxes[1];
+    gdImagePtr im = gdImageCreateTrueColor(W, H);
+    if(!im) return FALSE;
+    DBG("Create image %dx%d", W,H);
+    //OMP_FOR()
+    for(int y = H - 1; y >= 0; --y){ // flip up-down
+        for(int x = 0; x < (int)W; ++x){
+            im->tpixels[y][x] = (data[0] << 16) | (data[1] << 8) | data[2];
+            data += 3;
+        }
+    }
+    DBG("Converted");
+    FILE *fp = fopen(fname, "w");
+    if(!fp){
+        WARN(_("Can't save jpg image %s\n"), fname);
+        gdImageDestroy(im);
+        return FALSE;
+    }
+    if(str){
+        gdFTUseFontConfig(1);
+        char *font = (char*)"monotype";
+        char *ret = gdImageStringFT(im, NULL, 0xffffff, font, 10, 0., 2, 12, (char*)str);
+        if(ret) WARNX(_("Error: %s"), ret);
+    }
+    gdImageJpeg(im, fp, 90);
+    fclose(fp);
+    gdImageDestroy(im);
+    return TRUE;
+}
+
+/**
+ * @brief conv2rgb - convert grayscale normalized image to grayscale RGB
+ * @param inarr (i) - input array
+ * @param totpix    - its size
+ * @return allocated here array with data
+ *
+static uint8_t *conv2rgb(doubleimage *in){
+    if(!in) return NULL;
+    size_t totpix = in->totpix;
+    if(totpix == 0) return NULL;
+    double *inarr = in->data;
+    uint8_t *colored = MALLOC(uint8_t, totpix * 3);
+    OMP_FOR()
+    for(size_t i = 0; i < totpix; ++i){
+        uint8_t *pcl = &colored[i*3];
+        pcl[0] = pcl[1] = pcl[2] = (uint8_t)(inarr[i] * 255.);
+    }
+    return colored;
+}*/
+
+/**
+ * @brief gettransf - convert string with transformation type into intens_transform
+ * @param transf - type of transformation
+ * @return TRANSF_WRONG if arg is wrong or appropriate transformation type
+ * this function tests only first char[s] of type, so instead of, e.g. "log", you can
+ * write "lo", "looo", "Logogo", etc.
+ */
+static intens_transform gettransf(const char *transf){
+    if(!transf) return TRANSF_WRONG;
+    switch(transf[0]){
+        case 'e': // exp
+            return TRANSF_EXP;
+        break;
+        case 'l': // linear, log
+        case 'L':
+            switch(transf[1]){
+                case 'i':
+                case 'I':
+                    return TRANSF_LINEAR;
+                break;
+                case 'o':
+                case 'O':
+                    return TRANSF_LOG;
+                break;
+                default:
+                    return TRANSF_WRONG;
+            }
+        break;
+        case 'p':
+            return TRANSF_POW;
+        break;
+        case 's':
+            return TRANSF_SQR;
+        break;
+        default:
+            return TRANSF_WRONG;
+    }
+}
+
+/**
+ * @brief palette_transform - transform string with colormap name into its number
+ * @param p - colormap name
+ * @return PALETTE_COUNT if wrong or appropriate palette
+ */
+static image_palette palette_transform(char *p){
+    if(!p) return PALETTE_COUNT;
+    switch(p[0]){
+        case 'B':
+        case 'b':
+            return PALETTE_BR;
+        break;
+        default:
+            return PALETTE_COUNT;
+    }
+}
+
+int main(int argc, char *argv[]){
+    intens_transform tr = TRANSF_LINEAR;
+    initial_setup();
+    parse_args(argc, argv);
+    if(!G.fitsname) ERRX(_("No input filename given!"));
+    if(!G.outfile)  ERRX(_("Point the name of output file!"));
+    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);
+    DBG("Open file %s", G.fitsname);
+    FITS *f = FITS_read(G.fitsname);
+    DBG("HERE");
+    green("got file %s, HDUs: %d, working HDU #%d\n", G.fitsname, f->NHDUs, G.nhdu);
+    if(f->NHDUs < G.nhdu) ERRX(_("File %s consists %d HDUs!"), G.fitsname, f->NHDUs);
+    f->curHDU = &f->HDUs[G.nhdu];
+    if(f->curHDU->hdutype != IMAGE_HDU) ERRX(_("HDU %d is not image!"), G.nhdu);
+    FITSimage *img = f->curHDU->contents.image;
+    if(img->naxis != 2) ERRX(_("Support only 2-dimensional images"));
+    DBG("convert image from HDU #%d into double", G.nhdu);
+    doubleimage *dblimg = image2double(img);
+    if(!dblimg) ERRX(_("Can't convert image from HDU %s"), G.nhdu);
+    DBG("Done");
+    imgstat *st = get_imgstat(dblimg, NULL);
+    DBG("Image statistics: MIN=%g, MAX=%g, AVR=%g, STD=%g", st->min, st->max, st->mean, st->std);
+    if(!normalize_dbl(dblimg, st)) ERRX(_("Can't normalize image!"));
+#ifdef EBUG
+    st = get_imgstat(dblimg, NULL);
+#endif
+    DBG("NOW: MIN=%g, MAX=%g, AVR=%g, STD=%g", st->min, st->max, st->mean, st->std);
+    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);
+    image_palette colormap = PALETTE_GRAY;
+    if(G.palette){ // convert normalized image due to choosen palette
+        colormap = palette_transform(G.palette);
+        if(colormap == PALETTE_COUNT) ERRX(_("Wrong colormap name"));
+    }
+    //uint8_t *colored = conv2rgb(dblimg);
+    uint8_t *colored = convert2palette(dblimg, colormap);
+    DBG("Save jpeg to %s", G.outfile);
+    if(!write_jpeg(G.outfile, colored, G.text, img)) ERRX(_("Can't save modified file %s"), G.outfile);
+    green("File %s saved\n", G.outfile);
+    return 0;
+}

examples/imstat.c

@@ -57,13 +57,6 @@ static myoption cmdlnopts[] = {
     end_option
 };
 
-typedef struct{
-    double mean;
-    double std;
-    double min;
-    double max;
-} imgstat;
-
 /**
  * Parse command line options and return dynamically allocated structure
  *      to global parameters
@@ -89,25 +82,6 @@ static glob_pars *parse_args(int argc, char **argv){
     return &G;
 }
 
-static imgstat *get_imgstat(double *dimg, long totpix){
-    static imgstat st;
-    if(!dimg || !totpix) return &st; // return some trash if wrong data
-    st.min = dimg[0];
-    st.max = dimg[0];
-    double sum = dimg[0], sum2 = dimg[0];
-    for(long i = 1; i < totpix; ++i){
-        double val = dimg[i];
-        if(st.min > val) st.min = val;
-        if(st.max < val) st.max = val;
-        sum += val;
-        sum2 += val*val;
-    }
-    DBG("tot:%ld, sum=%g, sum2=%g, min=%g, max=%g", totpix, sum, sum2, st.min, st.max);
-    st.mean = sum / totpix;
-    st.std = sqrt(sum2/totpix - st.mean*st.mean);
-    return &st;
-}
-
 static void printstat(imgstat *stat){
     green("Statistics:\n");
     printf("MEAN=%g\nSTD=%g\nMIN=%g\nMAX=%g\n", stat->mean, stat->std, stat->min, stat->max);
@@ -176,10 +150,11 @@ static bool process_fitsfile(char *inname, FITS *output){
     // OK, we have an image and can do something with it
     green("\tGet image from this HDU.\n");
     FITSimage *img = f->curHDU->contents.image;
-    double *dImg = image2double(img);
+    doubleimage *dblim = image2double(img);
     // calculate image statistics
-    imgstat *stat = get_imgstat(dImg, img->totpix);
+    imgstat *stat = get_imgstat(dblim, NULL);
     printstat(stat);
+    double *dImg = dblim->data;
     DBG("i[1000] = %d, o[1000]=%g", ((uint16_t*)img->data)[1000], dImg[1000]);
     if(G.add){
         if(addsomething(img, dImg, stat)) mod = TRUE;

fitsimages.c

@@ -170,6 +170,7 @@ FITSimage *image_rebuild(FITSimage *img, double *dimg){
     // first we should calculate statistics of new image
     double *sr = MALLOC(double, img->totpix);
     memcpy(sr, dimg, sizeof(double)*img->totpix);
+    initomp();
     qsort(sr, img->totpix, sizeof(double), cmpdbl);
     double mindiff = DBL_MAX;
     bool isint = TRUE;
@@ -318,21 +319,31 @@ FITSimage *image_read(FITS *fits){
     return img;
 }
 
+void dblima_free(doubleimage **im){
+    FREE((*im)->data);
+    FREE(*im);
+}
+
 /**
  * @brief image2double convert image values to double
  * @param img - input image
- * @return array of double with size imt->totpix
+ * @return array of double with size img->totpix
  */
-double *image2double(FITSimage *img){
+doubleimage *image2double(FITSimage *img){
     size_t tot = img->totpix;
     double *ret = MALLOC(double, tot);
+    doubleimage *dblim = MALLOC(doubleimage, 1);
+    dblim->data = ret;
+    dblim->width = img->naxes[0];
+    dblim->height = img->naxes[1];
+    dblim->totpix = tot;
+    DBG("image: %ldx%ld=%ld", dblim->width, dblim->height, tot);
     double (*fconv)(uint8_t *x);
     double ubyteconv(uint8_t *data){return (double)*data;}
     double ushortconv(uint8_t *data){return (double)*((uint16_t*)data);}
     double ulongconv(uint8_t *data){return (double)*((uint32_t*)data);}
     double ulonglongconv(uint8_t *data){return (double)*((uint64_t*)data);}
     double floatconv(uint8_t *data){return (double)*((float*)data);}
-    initomp();
     switch(img->dtype){
         case TBYTE:
             fconv = ubyteconv;
@@ -351,17 +362,77 @@ double *image2double(FITSimage *img){
         break;
         case TDOUBLE:
             memcpy(ret, img->data, sizeof(double)*img->totpix);
-            return ret;
+            return dblim;
         break;
         default:
             WARNX(_("Undefined image type, cant convert to double"));
             FREE(ret);
+            FREE(dblim);
             return NULL;
     }
     uint8_t *din = img->data;
+    initomp();
     OMP_FOR()
     for(size_t i = 0; i < tot; ++i){
         ret[i] = fconv(&din[i*img->pxsz]);
     }
-    return ret;
+    return dblim;
+}
+
+/**
+ * @brief get_imgstat - calculate simplest statistics: mean/std/min/max
+ * @param dimg   - double array
+ * @param totpix - total amount of pixels
+ * @param est    - structure for output data (for thread-safe operations)
+ * @return structure with statistics data
+ */
+imgstat *get_imgstat(const doubleimage *im, imgstat *est){
+    static imgstat st;
+    if(!im || !im->totpix) return &st; // return some trash if wrong data
+    double *dimg = im->data;
+    size_t totpix = im->totpix;
+    st.min = dimg[0];
+    st.max = dimg[0];
+    double sum = dimg[0], sum2 = dimg[0];
+    for(size_t i = 1; i < totpix; ++i){
+        double val = dimg[i];
+        if(st.min > val) st.min = val;
+        if(st.max < val) st.max = val;
+        sum += val;
+        sum2 += val*val;
+    }
+    DBG("tot:%ld, sum=%g, sum2=%g, min=%g, max=%g", totpix, sum, sum2, st.min, st.max);
+    st.mean = sum / totpix;
+    st.std = sqrt(sum2/totpix - st.mean*st.mean);
+    if(est){
+        memcpy(est, &st, sizeof(imgstat));
+        return est;
+    }
+    return &st;
+}
+
+/**
+ * @brief normalize_dbl - convert double image array to normalized (0..1)
+ * @param dimg (io) - array with image pixels
+ * @param st   (i)  - image statistics (maybe NULL, then calculates here)
+ * @return pointer to dimg
+ */
+doubleimage *normalize_dbl(doubleimage *im, imgstat *st){
+    if(!im || !im->data || !st) return NULL;
+    double *dimg = im->data;
+    size_t totpix = im->totpix;
+    if(totpix < 1) return NULL;
+    imgstat imst;
+    if(!st) st = get_imgstat(im, &imst);
+    double rng = st->max - st->min;
+    if(rng < 2*DBL_EPSILON){
+        WARNX(_("Data range is too small"));
+        return NULL;
+    }
+    initomp();
+    OMP_FOR()
+    for(size_t i = 0; i < totpix; ++i){
+        dimg[i] = (dimg[i] - st->min) / rng;
+    }
+    return im;
 }