diff --git a/FITSmanip.h b/FITSmanip.h index c3cf9e3..4f53353 100644 --- a/FITSmanip.h +++ b/FITSmanip.h @@ -17,9 +17,19 @@ */ #include +#include #include #include +#define Stringify(x) #x +#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x)) +#ifndef MAX +#define MAX(x,y) ((x) > (y) ? (x) : (y)) +#endif +#ifndef MIN +#define MIN(x,y) ((x) < (y) ? (x) : (y)) +#endif + /************************************************************************************** * fits.c * **************************************************************************************/ @@ -176,6 +186,7 @@ void table_print_all(FITS *fits); void image_free(FITSimage **ima); FITSimage *image_read(FITS *fits); +FITSimage *image_rebuild(FITSimage *img, double *dimg); int image_datatype_size(int bitpix, int *dtype); void *image_data_malloc(long totpix, int pxbytes); FITSimage *image_new(int naxis, long *naxes, int bitpix); diff --git a/examples/imstat.c b/examples/imstat.c index 4a2e4ed..1bc31d4 100644 --- a/examples/imstat.c +++ b/examples/imstat.c @@ -24,11 +24,15 @@ #include #include #include +#include typedef struct{ char *outfile; // output file name (for single operation) char **infiles; // input files list int Ninfiles; // input files number + char *add; // add some value to all pixels + double mult; // multiply all pixels by some value + int rmneg; // remove negative values (assign them to 0) } glob_pars; /* @@ -36,6 +40,7 @@ typedef struct{ */ int help; glob_pars G = { + .mult = 1., }; /* @@ -46,6 +51,9 @@ myoption cmdlnopts[] = { // common options {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")}, {"outfile", NEED_ARG, NULL, 'o', arg_string, APTR(&G.outfile), _("output file name (collect all input files)")}, + {"add", NEED_ARG, NULL, 'a', arg_string, APTR(&G.add), _("add some value (double, or 'mean', 'std', 'min', 'max')")}, + {"multiply",NEED_ARG, NULL, 'm', arg_double, APTR(&G.mult), _("multiply by some value (double, operation run after adding)")}, + {"rmneg", NO_ARGS, NULL, 'z', arg_none, APTR(&G.rmneg), _("remove negative values (assign them to 0)")}, end_option }; @@ -105,6 +113,39 @@ void printstat(imgstat *stat){ printf("MEAN=%g\nSTD=%g\nMIN=%g\nMAX=%g\n", stat->mean, stat->std, stat->min, stat->max); } +bool addsomething(FITSimage *img, double *dimg, imgstat *stat){ + if(!G.add || !img || !stat) return FALSE; + // parser: + char *eptr; + double val = 1., addval = 0.; + val = strtod(G.add, &eptr); + if(eptr == G.add && *G.add == '-'){++eptr; val = -1.;} + if(*eptr){ // something more + if(strncasecmp(eptr, "mean", 4) == 0) addval = stat->mean; + else if(strncasecmp(eptr, "std", 3) == 0) addval = stat->std; + else if(strncasecmp(eptr, "min", 3) == 0) addval = stat->min; + else if(strncasecmp(eptr, "max", 3) == 0) addval = stat->max; + } + if(fabs(val) > DBL_EPSILON) addval *= val; + if(fabs(addval) > DBL_EPSILON) val = addval; + DBG("Add value %g", val); + if(fabs(val) < 2.*DBL_EPSILON) return FALSE; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i) + dimg[i] += val; + return TRUE; +} + +bool multbysomething(FITSimage *img, double *dimg){ + if(!img || !dimg) return FALSE; + if(fabs(G.mult) < 2*DBL_EPSILON) return FALSE; + DBG("multiply by %g", G.mult); + OMP_FOR() + for(long i = 0; i < img->totpix; ++i) + dimg[i] *= G.mult; + return TRUE; +} + bool process_fitsfile(char *inname, FITS *output){ DBG("File %s", inname); bool mod = FALSE; @@ -129,16 +170,39 @@ bool process_fitsfile(char *inname, FITS *output){ } if(!f->curHDU){ WARNX("Didn't find image HDU in %s", inname); - }else{ // 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); - // calculate image statistics - imgstat *stat = get_imgstat(dImg, img->totpix); - printstat(stat); - //mod = TRUE; - FREE(dImg); + FITS_free(&f); + return FALSE; } + // 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); + // calculate image statistics + imgstat *stat = get_imgstat(dImg, img->totpix); + printstat(stat); + DBG("i[1000] = %d, o[1000]=%g", ((uint16_t*)img->data)[1000], dImg[1000]); + if(G.add){ + if(addsomething(img, dImg, stat)) mod = TRUE; + } + DBG("ADD: i[1000]=%g", dImg[1000]); + if(fabs(G.mult - 1.) > DBL_EPSILON){ + if(multbysomething(img, dImg)) mod = TRUE; + } + DBG("MUL: i[1000]=%g", dImg[1000]); + if(G.rmneg){ + DBG("REMOVE negative values"); + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + if(dImg[i] < 0.) dImg[i] = 0.; + mod = TRUE; + } + } + DBG("NEG: i[1000]=%g", dImg[1000]); + if(mod){ // modified -> change file type + image_rebuild(img, dImg); + DBG("i[1000]=%g", dImg[1000]); + } + FREE(dImg); if(mod || output){ // file modified (or output file pointed), write differences if(!output){ green("Rewrite file %s.\n", f->filename); @@ -166,6 +230,7 @@ int main(int argc, char *argv[]){ ofits = MALLOC(FITS, 1); ofits->filename = G.outfile; } + initomp(); for(int i = 0; i < G.Ninfiles; ++i){ if(process_fitsfile(G.infiles[i], ofits)) mod = 1; } diff --git a/fits.c b/fits.c index 45e261c..f49d0c0 100644 --- a/fits.c +++ b/fits.c @@ -940,7 +940,7 @@ int image_datatype_size(int bitpix, int *dtype){ if(dtype){ switch(s){ case 1: // BYTE_IMG - *dtype = TBYTE; + *dtype = TBYTE; // uchar break; case 2: // SHORT_IMG *dtype = TUSHORT; @@ -1009,6 +1009,104 @@ FITSimage *image_new(int naxis, long *naxes, int bitpix){ return out; } +// function for qsort +static int cmpdbl(const void *d1, const void *d2){ + register double D1 = *(double*)d1, D2 = *(double*)d2; + if(fabs(D1 - D2) < DBL_EPSILON) return 0; + if(D1 > D2) return 1; + else return -1; +} + +// functions to convert double to different datatypes +static void convu8(FITSimage *img, double *dimg){ + uint8_t *dptr = (uint8_t*) img->data; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + dptr[i] = (uint8_t) dimg[i]; + } +} +static void convu16(FITSimage *img, double *dimg){ + uint16_t *dptr = (uint16_t*) img->data; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + dptr[i] = (uint16_t) dimg[i]; + } +} +static void convu32(FITSimage *img, double *dimg){ + uint32_t *dptr = (uint32_t*) img->data; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + dptr[i] = (uint32_t) dimg[i]; + } +} +static void convu64(FITSimage *img, double *dimg){ + uint64_t *dptr = (uint64_t*) img->data; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + dptr[i] = (uint64_t) dimg[i]; + } +} +static void convf(FITSimage *img, double *dimg){ + float *dptr = (float*) img->data; + OMP_FOR() + for(long i = 0; i < img->totpix; ++i){ + dptr[i] = (float) dimg[i]; + } +} + +/** + * @brief image_rebuild substitute content of image with array dimg, change its output type + * @param img - input image + * @param dimg - data to change + * @return rebuilt image; array dimg no longer used an can be FREEd + */ +FITSimage *image_rebuild(FITSimage *img, double *dimg){ + if(!img || !dimg) return NULL; + // first we should calculate statistics of new image + double *sr = MALLOC(double, img->totpix); + memcpy(sr, dimg, sizeof(double)*img->totpix); + qsort(sr, img->totpix, sizeof(double), cmpdbl); + double mindiff = DBL_MAX; + bool isint = TRUE; + for(long i = 1; i < img->totpix; ++i){ + double d = fabs(sr[i] - sr[i-1]); + if(d > DBL_EPSILON && d < mindiff) mindiff = d; + if(fabs(d - floor(d)) > DBL_EPSILON) isint = FALSE; + } + // now we know min, max and minimal difference between elements + double min = sr[0], max = sr[img->totpix-1]; + FREE(sr); + DBG("min: %g, max: %g, mindiff: %g", min, max, mindiff); + int bitpix = -64; // double by default + void (*convdata)(FITSimage*, double*) = NULL; + if(isint)do{ // check which integer type will suits better + DBG("INTEGER?"); + if(min < 0){ isint = FALSE; break;} // TODO: correct with BZERO + if(max < UINT8_MAX){bitpix = 8; convdata = convu8; break; } + if(max < UINT16_MAX){bitpix = 16; convdata = convu16; break; } + if(max < UINT32_MAX){bitpix = 32; convdata = convu32; break; } + if(max < UINT64_MAX){bitpix = 64; convdata = convu64; break; } + }while(0); + if(!isint){ // float or double + DBG("mindiff: %g(%g), min: %g(%g), max: %g(%g)", mindiff,FLT_EPSILON, + min, FLT_MIN, max, FLT_MAX); + if(mindiff > FLT_EPSILON && (min > -FLT_MAX) && max < FLT_MAX){ + bitpix = -32; + convdata = convf; + } + } + DBG("NOW: bitpix = %d", bitpix); + img->bitpix = bitpix; + void *data = calloc(img->totpix, abs(bitpix/8)); + if(!data){ WARN("calloc()"); return NULL; } + FREE(img->data); + img->data = data; + if(convdata) convdata(img, dimg); + else memcpy(img->data, dimg, sizeof(double)*img->totpix); + img->pxsz = image_datatype_size(bitpix, &img->dtype); + return img; +} + /** * build IMAGE image from data array indata * @@ -1157,8 +1255,8 @@ double *image2double(FITSimage *img){ } uint8_t *din = img->data; OMP_FOR() - for(size_t i = 0; i < tot; i += img->pxsz){ - ret[i] = fconv(&din[i]); + for(size_t i = 0; i < tot; ++i){ + ret[i] = fconv(&din[i*img->pxsz]); } return ret; } diff --git a/local.h b/local.h index 8205473..823493e 100644 --- a/local.h +++ b/local.h @@ -15,6 +15,7 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ +#include // xx_EPSILON etc. #if defined GETTEXT #include @@ -33,12 +34,4 @@ #define DBL_MAX (1.7976931348623157e+308) #endif -#define Stringify(x) #x -#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x)) -#ifndef MAX -#define MAX(x,y) ((x) > (y) ? (x) : (y)) -#endif -#ifndef MIN -#define MIN(x,y) ((x) < (y) ? (x) : (y)) -#endif