diff --git a/FITSmanip.c b/FITSmanip.c index 9d0aa7e..0d51cf2 100644 --- a/FITSmanip.c +++ b/FITSmanip.c @@ -20,10 +20,33 @@ #include #include #include +#include -static inline void initomp(){ +/* + * Here are some common functions used in library + */ + +/** + * @brief initomp init optimal CPU number for OPEN MP + */ +void initomp(){ + static bool got = FALSE; + if(got) return; int cpunumber = sysconf(_SC_NPROCESSORS_ONLN); if(omp_get_max_threads() != cpunumber) omp_set_num_threads(cpunumber); + got = TRUE; } +/** + * @brief FITS_reporterr - show error in red + * @param errcode + */ +void FITS_reporterr(int *errcode){ + if(isatty(STDERR_FILENO)) + fprintf(stderr, COLOR_RED); + fits_report_error(stderr, *errcode); + if(isatty(STDERR_FILENO)) + fprintf(stderr, COLOR_OLD); + *errcode = 0; +} diff --git a/FITSmanip.h b/FITSmanip.h index b322b18..c3cf9e3 100644 --- a/FITSmanip.h +++ b/FITSmanip.h @@ -103,8 +103,9 @@ typedef struct{ FITS image */ typedef struct{ - int width; // width - int height; // height + int naxis; // amount of image dimensions + long *naxes; // dimensions + long totpix; // total pixels amount int bitpix; // original bitpix int dtype; // type of stored data int pxsz; // number of bytes for one pixel data @@ -176,12 +177,14 @@ void table_print_all(FITS *fits); void image_free(FITSimage **ima); FITSimage *image_read(FITS *fits); int image_datatype_size(int bitpix, int *dtype); -void *image_data_malloc(size_t w, size_t h, int pxbytes); -FITSimage *image_new(size_t w, size_t h, int bitpix); +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); //FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata); +FITSHDU *FITS_addHDU(FITS *fits); void FITS_free(FITS **fits); FITS *FITS_read(char *filename); FITS *FITS_open(char *filename); @@ -195,7 +198,11 @@ bool FITS_rewrite(FITS *fits); char* make_filename(char *buff, size_t buflen, char *prefix, char *suffix); bool file_is_absent(char *name); - +/************************************************************************************** + * FITSmanip.c * + **************************************************************************************/ +void FITS_reporterr(int *errcode); +void initomp(); /* // pointer to image conversion function diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 40eda2c..8baeace 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,7 +1,8 @@ cmake_minimum_required(VERSION 3.9) project(examples) include_directories(../) -link_libraries(usefull_macros FITSmanip cfitsio) +link_libraries(usefull_macros FITSmanip cfitsio m) #add_executable(fitsstat fitsstat.c) add_executable(keylist keylist.c) +add_executable(imstat imstat.c) diff --git a/examples/imstat.c b/examples/imstat.c new file mode 100644 index 0000000..4a2e4ed --- /dev/null +++ b/examples/imstat.c @@ -0,0 +1,178 @@ +/* + * 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 . + */ + +/* + * This example allows to calculate simple statistics over fits file, + * add constant number or multiply by constant all images in file + */ +#include "common.h" +#include +#include +#include + +typedef struct{ + char *outfile; // output file name (for single operation) + char **infiles; // input files list + int Ninfiles; // input files number +} glob_pars; + +/* + * here are global parameters initialisation + */ +int help; +glob_pars G = { +}; + +/* + * Define command line options by filling structure: + * name has_arg flag val type argptr help +*/ +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)")}, + 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 + * @param argc - copy of argc from main + * @param argv - copy of argv from main + * @return allocated structure with global parameters + */ +glob_pars *parse_args(int argc, char **argv){ + int i; + char *helpstring = "Usage: %%s [args] input files\n" + "Get statistics and modify images from first image HDU of each input file\n" + "\tWhere args are:\n"; + change_helpstring(helpstring); + // parse arguments + parseargs(&argc, &argv, cmdlnopts); + if(help) showhelp(-1, cmdlnopts); + if(argc > 0){ + G.Ninfiles = argc; + G.infiles = MALLOC(char*, argc); + for (i = 0; i < argc; i++) + G.infiles[i] = strdup(argv[i]); + } + return &G; +} + +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; +} + +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); +} + +bool process_fitsfile(char *inname, FITS *output){ + DBG("File %s", inname); + bool mod = FALSE; + FITS *f = FITS_read(inname); + if(!f){WARNX("Can't read %s", inname); return FALSE;} + // search first image HDU + f->curHDU = NULL; + for(int i = 1; i <= f->NHDUs; ++i){ + green("File %s, %dth HDU, type: %d\n", inname, i, f->HDUs[i].hdutype); + if(f->HDUs[i].hdutype != IMAGE_HDU) continue; + FITSimage *img = f->HDUs[i].contents.image; + if(!img){ + WARNX("empty image"); + continue; + } + if(img->totpix == 0){ + WARNX("totpix=0"); + continue; // no data - only header + } + f->curHDU = &f->HDUs[i]; + break; + } + 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); + } + if(mod || output){ // file modified (or output file pointed), write differences + if(!output){ + green("Rewrite file %s.\n", f->filename); + if(!FITS_rewrite(f)) WARNX("Can't rewrite %s", inname); + }else{ // copy first image HDU and its header to output file + green("Add image to %s.\n", output->filename); + FITSHDU *hdu = FITS_addHDU(output); + if(!hdu) return NULL; + hdu->hdutype = IMAGE_HDU; + hdu->keylist = keylist_copy(f->curHDU->keylist); + hdu->contents.image = image_copy(f->curHDU->contents.image); + } + } + FITS_free(&f); + return mod; +} + +int main(int argc, char *argv[]){ + FITS *ofits = NULL; + bool mod = FALSE; + initial_setup(); + parse_args(argc, argv); + if(!G.Ninfiles) ERRX(_("No input filename[s] given!")); + if(G.outfile){ + ofits = MALLOC(FITS, 1); + ofits->filename = G.outfile; + } + for(int i = 0; i < G.Ninfiles; ++i){ + if(process_fitsfile(G.infiles[i], ofits)) mod = 1; + } + if(ofits && mod){ + green("\nWrite all modified images to output file %s\n", ofits->filename); + FITS_write(ofits->filename, ofits); + } + return 0; +} + diff --git a/examples/keylist.c b/examples/keylist.c index a008464..0793657 100644 --- a/examples/keylist.c +++ b/examples/keylist.c @@ -55,6 +55,7 @@ myoption cmdlnopts[] = { {"addrec", MULT_PAR, NULL, 'a', arg_string, APTR(&G.addrec), _("add record to first HDU (you can add more than one record in once, point more -a)")}, {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&G.outfile), _("save result to file (else save to same file)")}, {"modify", MULT_PAR, NULL, 'm', arg_string, APTR(&G.modify), _("modify values values of given keys (each param should be \"key = new_value\")")}, + {"infile", NEED_ARG, NULL, 'i', arg_string, APTR(&G.fitsname), _("input file name (you can also point it without any keys)")}, end_option }; @@ -72,7 +73,7 @@ glob_pars *parse_args(int argc, char **argv){ // parse arguments parseargs(&argc, &argv, cmdlnopts); if(help) showhelp(-1, cmdlnopts); - if(argc)(G.fitsname = strdup(argv[0])); + if(argc && !G.fitsname)(G.fitsname = strdup(argv[0])); if(argc > 1){ for (i = 1; i < argc; i++) printf("Ignore extra argument: %s\n", argv[i]); @@ -86,6 +87,16 @@ void ch(int s){ signal(s, ch); } +void print_imgHDU(FITSimage *image){ + printf("Image: naxis=%d, totpix=%ld, ", image->naxis, image->totpix); + printf("naxes=("); + for(int i = 0; i < image->naxis; ++i) + printf("%s%ld", i?", ":"", image->naxes[i]); + if(image->naxis) printf("), "); + else printf("none), "); + printf("bitpix=%d, dtype=%d", image->bitpix, image->dtype); +} + int main(int argc, char *argv[]){ initial_setup(); parse_args(argc, argv); @@ -107,7 +118,7 @@ int main(int argc, char *argv[]){ printf("\tHDU #%d - ", i); switch(f->HDUs[i].hdutype){ case IMAGE_HDU: - printf("Image"); + print_imgHDU(f->HDUs[i].contents.image); break; case ASCII_TBL: printf("ASCII table"); @@ -121,11 +132,12 @@ int main(int argc, char *argv[]){ printf("\n"); } } + int differs = 0; if(G.addrec){ char **ptr = G.addrec; while(*ptr){ printf("record: %s\n", *ptr); - keylist_add_record(&(f->HDUs[1].keylist), *ptr++, 1); + if(keylist_add_record(&(f->HDUs[1].keylist), *ptr++, 1)) differs = 1; } } if(G.modify){ @@ -140,28 +152,30 @@ int main(int argc, char *argv[]){ *val++ = 0; // now `val` is value + comment; ptr is key if(!keylist_modify_key(f->HDUs[1].keylist, *ptr, val)){ WARNX("key %s not found", *ptr); - } + }else differs = 1; ++ptr; } } - // test for signals handler: ctrl+c, ctrl+z - signal(SIGINT, ch); - signal(SIGTSTP, ch); - // protect critical zone blocking all possible signals: - sigset_t mask, oldmask; - sigfillset(&mask); - sigprocmask(SIG_SETMASK, &mask, &oldmask); - if(G.outfile){ // save result to new file - FITS_write(G.outfile, f); - }else{ - FITS_rewrite(f); + if(differs){ // file differs, need to save new file + // test for signals handler: ctrl+c, ctrl+z + signal(SIGINT, ch); + signal(SIGTSTP, ch); + // protect critical zone blocking all possible signals: + sigset_t mask, oldmask; + sigfillset(&mask); + sigprocmask(SIG_SETMASK, &mask, &oldmask); + if(G.outfile){ // save result to new file + FITS_write(G.outfile, f); + }else{ + FITS_rewrite(f); + } + DBG("Written! Sleep for 2 seconds in ctitical section"); + sleep(2); + // return ignoring + DBG("Unblock signals, sleep for 2 seconds"); + sigprocmask(SIG_SETMASK, &oldmask, NULL); + sleep(2); } - DBG("Written! Sleep for 2 seconds in ctitical section"); - sleep(2); - // return ignoring - DBG("Unblock signals, sleep for 2 seconds"); - sigprocmask(SIG_SETMASK, &oldmask, NULL); - sleep(2); /* * Do something here */ diff --git a/fits.c b/fits.c index 83636e1..45e261c 100644 --- a/fits.c +++ b/fits.c @@ -78,7 +78,7 @@ KeyList *keylist_add_record(KeyList **list, char *rec, int check){ char card[FLEN_CARD]; fits_parse_template(rec, card, &tp, &st); if(st){ - fits_report_error(stderr, st); + FITS_reporterr(&st); return NULL; } DBG("\n WAS: %s\nBECOME: %s\ntp=%d", rec, card, tp); @@ -141,7 +141,7 @@ KeyList *keylist_modify_key(KeyList *list, char *key, char *newval){ int tp, st = 0; fits_parse_template(test, buf, &tp, &st); if(st){ - fits_report_error(stderr, st); + FITS_reporterr(&st); return NULL; } DBG("new record:\n%s", buf); @@ -283,14 +283,14 @@ KeyList *keylist_read(FITS *fits){ return NULL; } if(fst){ - fits_report_error(stderr, fst); + FITS_reporterr(&fst); return NULL; } DBG("Find %d keys, keypos=%d", nkeys, keypos); for(int j = 1; j <= nkeys; ++j){ char card[FLEN_CARD]; fits_read_record(fits->fp, j, card, &fst); - if(fst) fits_report_error(stderr, fst); + if(fst) FITS_reporterr(&fst); else{ KeyList *kl = keylist_add_record(&list, card, 0); if(!kl){ @@ -361,16 +361,16 @@ FITStable *table_copy(FITStable *intab){ * @return */ FITStable *table_read(FITS *fits){ - int ncols, i, fst = 0, ret; + int ncols, i, fst = 0; long nrows; char extname[FLEN_VALUE]; fitsfile *fp = fits->fp; fits_get_num_rows(fp, &nrows, &fst); - if(fst){fits_report_error(stderr, fst); return NULL;} + if(fst){FITS_reporterr(&fst); return NULL;} fits_get_num_cols(fp, &ncols, &fst); - if(fst){fits_report_error(stderr, fst); return NULL;} + if(fst){FITS_reporterr(&fst); return NULL;} fits_read_key(fp, TSTRING, "EXTNAME", extname, NULL, &fst); - if(fst){fits_report_error(stderr, fst); return NULL;} + if(fst){FITS_reporterr(&fst); return NULL;} DBG("Table named %s with %ld rows and %d columns", extname, nrows, ncols); FITStable *tbl = table_new(extname); if(!tbl) return NULL; @@ -378,9 +378,9 @@ FITStable *table_read(FITS *fits){ for(i = 1; i <= ncols; ++i){ int typecode; long repeat, width; - ret = fits_get_coltype(fp, i, &typecode, &repeat, &width, &fst); - if(fst){fits_report_error(stderr, fst); ret = fst; fst = 0;} - if(ret){ + fits_get_coltype(fp, i, &typecode, &repeat, &width, &fst); + if(fst){ + FITS_reporterr(&fst); WARNX(_("Can't read column %d!"), i); continue; } @@ -392,9 +392,9 @@ FITStable *table_read(FITS *fits){ int64_t nullval = 0; int j; for(j = 0; j < repeat; ++j){ - ret = fits_read_col(fp, typecode, i, j=1, 1, 1, (void*)nullval, array, &anynul, &fst); - if(fst){fits_report_error(stderr, fst); ret = fst; fst = 0;} - if(ret){ + fits_read_col(fp, typecode, i, j=1, 1, 1, (void*)nullval, array, &anynul, &fst); + if(fst){ + FITS_reporterr(&fst); WARNX(_("Can't read column %d row %d!"), i, j); continue; } @@ -654,11 +654,11 @@ bool table_write(FITS *file){ DBG("col: %s, form: %s, unit: %s", columns[c], formats[c], units[c]); } //fits_movabs_hdu(fptr, 2, &hdutype, &status) - int ret = fits_create_tbl(fp, hdutype, tbl->nrows, cols, + fits_create_tbl(fp, hdutype, tbl->nrows, cols, columns, formats, units, tbl->tabname, &fst); - if(fst){fits_report_error(stderr, fst); ret = fst; fst = 0;} FREE(columns); FREE(formats); FREE(units); - if(ret){ + if(fst){ + FITS_reporterr(&fst); WARNX(_("Can't write table %s!"), tbl->tabname); return FALSE; } @@ -668,7 +668,7 @@ bool table_write(FITS *file){ int fst = 0; fits_write_col(fp, col->coltype, c+1, 1, 1, col->repeat, col->contents, &fst); if(fst){ - fits_report_error(stderr, fst); + FITS_reporterr(&fst); WARNX(_("Can't write column %s!"), col->colname); return FALSE; } @@ -680,6 +680,29 @@ bool table_write(FITS *file){ * FITS files * **************************************************************************************/ +/** + * @brief FITS_addHDU - add new HDU to FITS file structure + * @param fits (io) - fits file to use + * @return pointer to new HDU or NULL in case of error + */ +FITSHDU *FITS_addHDU(FITS *fits){ + int hdunum = fits->NHDUs + 1; + // add 1 to `hdunum` because HDU numbering starts @1 + FITSHDU *newhdu = realloc(fits->HDUs, sizeof(FITSHDU)*(1+hdunum)); + if(!newhdu){ + WARN("FITS_addHDU, realloc() failed"); + return NULL; + } + fits->HDUs = newhdu; + fits->curHDU = &fits->HDUs[hdunum]; + fits->NHDUs = hdunum; + return fits->curHDU; +} + +/** + * @brief FITS_free - delete FITS structure from memory + * @param fits - address of FITS pointer + */ void FITS_free(FITS **fits){ if(!fits || !*fits) return; FITS *f = *fits; @@ -689,7 +712,9 @@ void FITS_free(FITS **fits){ int n, N = f->NHDUs; for(n = 1; n < N; ++n){ FITSHDU *hdu = &f->HDUs[n]; + if(!hdu) continue; keylist_free(&hdu->keylist); + if(!hdu->contents.image) continue; switch(hdu->hdutype){ case IMAGE_HDU: image_free(&hdu->contents.image); @@ -706,13 +731,7 @@ void FITS_free(FITS **fits){ } /** - - TODO: READWRITE allows to modify files on-the-fly, need to use it! TODO: what's about HCOMPRESS? - - * read FITS file and fill 'IMAGE' structure (with headers and tables) - * can't work with image stack - opens the first image met - * works only with binary tables */ /** @@ -726,7 +745,7 @@ FITS *FITS_open(char *filename){ // use fits_open_diskfile instead of fits_open_file to prevent using of extended name syntax fits_open_diskfile(&fits->fp, filename, READONLY, &fst); if(fst){ - fits_report_error(stderr, fst); + FITS_reporterr(&fst); FITS_free(&fits); return NULL; } @@ -780,11 +799,10 @@ FITS *FITS_read(char *filename){ } if(fst == END_OF_FILE){ fst = 0; - }else goto returning; + } returning: if(fst){ - fits_report_error(stderr, fst); - fst = 0; + FITS_reporterr(&fst); FITS_free(&fits); } return fits; @@ -798,7 +816,7 @@ static bool keylist_write(KeyList *kl, fitsfile *fp){ if(kl->keyclass > TYP_CMPRS_KEY){ // this record should be written fits_write_record(fp, kl->record, &st); DBG("Write %s, st = %d", kl->record, st); - if(st){fits_report_error(stderr, st); st = 0; ret = FALSE;} + if(st){FITS_reporterr(&st); ret = FALSE;} } kl = kl->next; } @@ -817,13 +835,12 @@ bool FITS_write(char *filename, FITS *fits){ int fst = 0; fits_create_file(&fp, filename, &fst); DBG("create file %s", filename); - if(fst){fits_report_error(stderr, fst); return FALSE;} + if(fst){FITS_reporterr(&fst); return FALSE;} int N = fits->NHDUs; for(int i = 1; i <= N; ++i){ FITSHDU *hdu = &fits->HDUs[i]; if(!hdu) continue; FITSimage *img; - long naxes[2] = {0,0}; KeyList *records = hdu->keylist; DBG("HDU #%d (type %d)", i, hdu->hdutype); switch(hdu->hdutype){ @@ -831,28 +848,26 @@ bool FITS_write(char *filename, FITS *fits){ img = hdu->contents.image; if(!img && records){ // something wrong - just write keylist DBG("create empty image with records"); - fits_create_img(fp, SHORT_IMG, 0, naxes, &fst); - if(fst){fits_report_error(stderr, fst); fst = 0; continue;} + fits_create_img(fp, SHORT_IMG, 0, NULL, &fst); + if(fst){FITS_reporterr(&fst); continue;} keylist_write(records, fp); DBG("OK"); continue; } - naxes[0] = img->width, naxes[1] = img->height; - DBG("create, bitpix: %d, naxes = {%zd, %zd}", img->bitpix, naxes[0], naxes[1]); - //fits_create_img(fp, img->bitpix, 2, naxes, &fst); - fits_create_img(fp, SHORT_IMG, 2, naxes, &fst); - if(fst){fits_report_error(stderr, fst); fst = 0; continue;} + DBG("create, bitpix: %d, naxis = %d, totpix = %ld", img->bitpix, img->naxis, img->totpix); + fits_create_img(fp, img->bitpix, img->naxis, img->naxes, &fst); + //fits_create_img(fp, SHORT_IMG, img->naxis, img->naxes, &fst); + if(fst){FITS_reporterr(&fst); continue;} keylist_write(records, fp); DBG("OK, now write image"); - // TODO: change to original data type according to bitpix or more whide according to min/max - DBG("bitpix: %d, dtype: %d", SHORT_IMG, img->dtype); //int bscale = 1, bzero = 32768, status = 0; //fits_set_bscale(fp, bscale, bzero, &status); if(img && img->data){ - fits_write_img(fp, TUSHORT, 1, img->width * img->height, img->data, &fst); + DBG("bitpix: %d, dtype: %d", img->bitpix, img->dtype); + fits_write_img(fp, img->dtype, 1, img->totpix, img->data, &fst); //fits_write_img(fp, img->dtype, 1, img->width * img->height, img->data, &fst); DBG("status: %d", fst); - if(fst){fits_report_error(stderr, fst); fst = 0; continue;} + if(fst){FITS_reporterr(&fst); continue;} } break; case BINARY_TBL: @@ -863,7 +878,7 @@ bool FITS_write(char *filename, FITS *fits){ } fits_close_file(fp, &fst); - if(fst){fits_report_error(stderr, fst);} + if(fst){FITS_reporterr(&fst);} return TRUE; } @@ -910,6 +925,7 @@ bool FITS_rewrite(FITS *fits){ void image_free(FITSimage **img){ FREE((*img)->data); + FREE((*img)->naxes); FREE(*img); } @@ -923,63 +939,70 @@ int image_datatype_size(int bitpix, int *dtype){ int s = bitpix/8; if(dtype){ switch(s){ - case 1: - case 2: - case 4: - *dtype = TINT; - s = 4; + case 1: // BYTE_IMG + *dtype = TBYTE; break; - case 8: - *dtype = TLONG; + case 2: // SHORT_IMG + *dtype = TUSHORT; break; - case -4: - case -8: + case 4: // LONG_IMG + *dtype = TUINT; + break; + case 8: // LONGLONG_IMG + *dtype = TULONG; + break; + case -4: // FLOAT_IMG + *dtype = TFLOAT; + break; + case -8: // DOUBLE_IMG *dtype = TDOUBLE; - s = 8; break; default: return 0; // wrong bitpix } } DBG("bitpix: %d, dtype=%d, imgs=%d", bitpix, *dtype, s); - return s; + return abs(s); } /** - * @brief image_malloc - allocate memory for given bitpix - * @param w - image width - * @param h - image height - * @param pxbytes- amount of bytes to store data from one pixel + * @brief image_data_malloc - allocate memory for given bitpix + * @param totpix - total pixels amount + * @param pxbytes - number of bytes for each pixel * @return allocated memory */ -void *image_data_malloc(size_t w, size_t h, int pxbytes){ - if(!pxbytes || !w || !h) return NULL; - void *data = calloc(w*h, pxbytes); - DBG("Allocate %zd members of size %d", w*h, pxbytes); - if(!data) ERR(_("calloc()")); +void *image_data_malloc(long totpix, int pxbytes){ + if(!pxbytes || !totpix) return NULL; + void *data = calloc(totpix, pxbytes); + DBG("Allocate %zd members of size %d", totpix, pxbytes); + if(!data) ERR("calloc()"); return data; } /** * @brief image_new - create an empty image without headers, assign BITPIX to "bitpix" - * @param w - image width - * @param h - image height - * @param dtype - image data type - * @return + * @param naxis - number of dimensions + * @param naxes - sizes by each dimension + * @param bitpix - BITPIX for given image + * @return allocated structure or NULL */ -FITSimage *image_new(size_t w, size_t h, int bitpix){ +FITSimage *image_new(int naxis, long *naxes, int bitpix){ FITSimage *out = MALLOC(FITSimage, 1); int dtype, pxsz = image_datatype_size(bitpix, &dtype); - if(w && h){ - out->data = image_data_malloc(w, h, pxsz); + long totpix = 0; + if(naxis){ // not empty image + totpix = 1; + for(int i = 0; i < naxis; ++i) if(naxes[i]) totpix *= naxes[i]; + out->data = image_data_malloc(totpix, pxsz); if(!out->data){ - WARNX(_("Bad w, h or pxsz")); FREE(out); return NULL; } } - out->width = w; - out->height = h; + out->totpix = totpix; + out->naxes = MALLOC(long, naxis); + memcpy(out->naxes, naxes, sizeof(long)*naxis); + out->naxis = naxis; out->pxsz = pxsz; out->bitpix = bitpix; out->dtype = dtype; @@ -1045,8 +1068,7 @@ FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata){ * create an empty copy of image "in" without headers, assign data type to "dtype" */ FITSimage *image_mksimilar(FITSimage *img){ - if(!img || img->height < 1 || img->width < 1) return NULL; - return image_new(img->width, img->height, img->bitpix); + return image_new(img->naxis, img->naxes, img->bitpix); } /** @@ -1055,8 +1077,7 @@ FITSimage *image_mksimilar(FITSimage *img){ FITSimage *image_copy(FITSimage *in){ FITSimage *out = image_mksimilar(in); if(!out) return NULL; - // TODO: size of data as in original! - memcpy(out->data, in->data, (in->pxsz)*(in->width)*(in->height)); + memcpy(out->data, in->data, (in->pxsz)*(in->totpix)); return out; } @@ -1069,31 +1090,23 @@ FITSimage *image_read(FITS *fits){ // TODO: open not only 2-dimensional files! // get image dimensions int naxis, fst = 0, bitpix; - long naxes[2] = {0,0}; - fits_get_img_param(fits->fp, 2, &bitpix, &naxis, naxes, &fst); - if(fst){fits_report_error(stderr, fst); return NULL;} - if(naxis > 2){ - WARNX(_("Images with > 2 dimensions are not supported")); - return NULL; - }/* - if(naxis < 2){ - WARNX(_("Not an image: NAXIS = %d"), naxis); - return NULL; - }*/ - DBG("got image %ldx%ld pix, bitpix=%d", naxes[0], naxes[1], bitpix); - - FITSimage *img = image_new(naxes[0], naxes[1], bitpix); - + fits_get_img_dim(fits->fp, &naxis, &fst); + if(fst){FITS_reporterr(&fst); return NULL;} + long *naxes = MALLOC(long, naxis); + fits_get_img_param(fits->fp, naxis, &bitpix, &naxis, naxes, &fst); + if(fst){FITS_reporterr(&fst); return NULL;} + FITSimage *img = image_new(naxis, naxes, bitpix); + FREE(naxes); int stat = 0; if(!img) return NULL; if(!img->data) return img; // empty "image" - no data inside - DBG("try to read, dt=%d, sz=%ld", img->dtype, naxes[0]*naxes[1]); + DBG("try to read, dt=%d, sz=%ld", img->dtype, img->totpix); //int bscale = 1, bzero = 32768, status = 0; //fits_set_bscale(fits->fp, bscale, bzero, &status); - //fits_read_img(fits->fp, img->dtype, 1, naxes[0]*naxes[1], NULL, img->data, &stat, &fst); - fits_read_img(fits->fp, TUSHORT, 1, naxes[0]*naxes[1], NULL, img->data, &stat, &fst); + fits_read_img(fits->fp, img->dtype, 1, img->totpix, NULL, img->data, &stat, &fst); + //fits_read_img(fits->fp, TUSHORT, 1, img->totpix, NULL, img->data, &stat, &fst); if(fst){ - fits_report_error(stderr, fst); + FITS_reporterr(&fst); image_free(&img); return NULL; } @@ -1101,3 +1114,51 @@ FITSimage *image_read(FITS *fits){ DBG("ready"); return img; } + +/** + * @brief image2double convert image values to double + * @param img - input image + * @return array of double with size imt->totpix + */ +double *image2double(FITSimage *img){ + size_t tot = img->totpix; + double *ret = MALLOC(double, 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; + break; + case TUSHORT: + fconv = ushortconv; + break; + case TUINT: + fconv = ulongconv; + break; + case TULONG: + fconv = ulonglongconv; + break; + case TFLOAT: + fconv = floatconv; + break; + case TDOUBLE: + memcpy(ret, img->data, sizeof(double)*img->totpix); + return ret; + break; + default: + WARNX(_("Undefined image type, cant convert to double")); + FREE(ret); + return NULL; + } + uint8_t *din = img->data; + OMP_FOR() + for(size_t i = 0; i < tot; i += img->pxsz){ + ret[i] = fconv(&din[i]); + } + return ret; +} diff --git a/locale/ru/LC_MESSAGES/FITSmanip.mo b/locale/ru/LC_MESSAGES/FITSmanip.mo index 1025262..35fc8e5 100644 Binary files a/locale/ru/LC_MESSAGES/FITSmanip.mo and b/locale/ru/LC_MESSAGES/FITSmanip.mo differ diff --git a/locale/ru/messages.po b/locale/ru/messages.po index e48776d..fe32c89 100644 --- a/locale/ru/messages.po +++ b/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-02-21 20:26+0300\n" +"POT-Creation-Date: 2019-02-22 22:46+0300\n" "PO-Revision-Date: YEAR-MO-DA HO:MI+ZONE\n" "Last-Translator: FULL NAME \n" "Language-Team: LANGUAGE \n" @@ -67,28 +67,24 @@ msgstr "" msgid "Can't write column %s!" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:751 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:770 msgid "Can't read HDU" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:778 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:797 msgid "Unknown HDU type" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:900 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:915 #, c-format msgid "Can't get real path for %s, use cfitsio to rewrite" msgstr "" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:976 -msgid "Bad w, h or pxsz" -msgstr "" - -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1076 -msgid "Images with > 2 dimensions are not supported" -msgstr "" - -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1100 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1113 #, c-format msgid "Found %d pixels with undefined value" msgstr "" + +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1154 +msgid "Undefined image type, cant convert to double" +msgstr "" diff --git a/locale/ru/ru.po b/locale/ru/ru.po index 494774a..3d364f4 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-02-21 20:26+0300\n" + "POT-Creation-Date: 2019-02-22 22:46+0300\n" "PO-Revision-Date: YEAR-MO-DA HO:MI+ZONE\n" "Last-Translator: FULL NAME \n" "Language-Team: LANGUAGE \n" @@ -16,10 +16,6 @@ msgstr "Project-Id-Version: PACKAGE VERSION\n" "Content-Type: text/plain; charset=koi8-r\n" "Content-Transfer-Encoding: 8bit\n" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:976 -msgid "Bad w, h or pxsz" -msgstr "Неверные значения w, h или pxsz" - #. / "Не могу добавить запись в список" #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:298 msgid "Can't add record to list" @@ -30,13 +26,13 @@ msgstr " msgid "Can't copy data" msgstr "Не могу скопировать данные" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:900 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:915 #, c-format msgid "Can't get real path for %s, use cfitsio to rewrite" msgstr "Не могу определить путь (realpath) к %s, использую cfitsio для " "перезаписи" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:751 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:770 msgid "Can't read HDU" msgstr "Не могу прочесть HDU" @@ -68,20 +64,20 @@ msgstr " msgid "Can't write table %s!" msgstr "Не могу записать таблицу %s!" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1100 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1113 #, c-format msgid "Found %d pixels with undefined value" msgstr "Найдено %d пикселей с неопределенными значениями" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1076 -msgid "Images with > 2 dimensions are not supported" -msgstr "Изображения с более чем двумя измерениями не поддерживаются" - #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:282 msgid "No keywords in given HDU" msgstr "В данном HDU ключи отсутствуют" -#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:778 +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:1154 +msgid "Undefined image type, cant convert to double" +msgstr "" + +#: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:797 msgid "Unknown HDU type" msgstr "Неизвестный тип HDU" @@ -92,3 +88,9 @@ msgstr " #: /home/eddy/C-files/FITSmaniplib/sharedlib_template/fits.c:351 msgid "strdup() failed!" msgstr "Не удалось сделать strdup()!" + +#~ msgid "Bad w, h or pxsz" +#~ msgstr "Неверные значения w, h или pxsz" + +#~ msgid "Images with > 2 dimensions are not supported" +#~ msgstr "Изображения с более чем двумя измерениями не поддерживаются"