eddy/fitsmaniplib
1
0
Fork 0
mirror of https://github.com/eddyem/fitsmaniplib.git synced 2026-02-03 05:45:09 +03:00
Code Issues Packages Projects Releases Wiki Activity

modified FITS structure, start debugging image reading

Browse Source
This commit is contained in:
eddyem 2019-02-18 23:29:47 +03:00
parent 620f2c0612
commit 82b214db95
3 changed files with 206 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
FITSmanip.h
View File

@ -37,6 +37,7 @@ cfitsio.h BITPIX code values for FITS image types:
#define FLOAT_IMG -32 #define FLOAT_IMG -32
#define DOUBLE_IMG -64 #define DOUBLE_IMG -64
*/ */
/*
// FilterType (not only convolution!) // FilterType (not only convolution!)
typedef enum{ typedef enum{
FILTER_NONE = 0 // simple start FILTER_NONE = 0 // simple start
@ -58,7 +59,11 @@ typedef struct{
double *data; double *data;
size_t size; size_t size;
}Itmarray; }Itmarray;
*/
/**
Keylist: all keys from given HDU
*/
typedef struct klist_{ typedef struct klist_{
int keyclass; // key class [look int CFITS_API ffgkcl(char *tcard) ] int keyclass; // key class [look int CFITS_API ffgkcl(char *tcard) ]
char *record; // record itself char *record; // record itself
@ -66,6 +71,9 @@ typedef struct klist_{
struct klist_ *last; // previous record struct klist_ *last; // previous record
} KeyList; } KeyList;
/**
Table column
*/
typedef struct{ typedef struct{
void *contents; // contents of table void *contents; // contents of table
int coltype; // type of columns int coltype; // type of columns
@ -74,20 +82,26 @@ typedef struct{
char colname[FLEN_KEYWORD]; // column name (arg ttype of fits_create_tbl) char colname[FLEN_KEYWORD]; // column name (arg ttype of fits_create_tbl)
char format[FLEN_FORMAT]; // format codes (tform) char format[FLEN_FORMAT]; // format codes (tform)
char unit[FLEN_CARD]; // units (tunit) char unit[FLEN_CARD]; // units (tunit)
}table_column; } table_column;
/**
FITS table
*/
typedef struct{ typedef struct{
int ncols; // amount of columns int ncols; // amount of columns
long nrows; // max amount of rows long nrows; // max amount of rows
char tabname[FLEN_CARD]; // table name char tabname[FLEN_CARD]; // table name
table_column *columns; // array of structures 'table_column' table_column *columns; // array of structures 'table_column'
}FITStable; } FITStable;
/* /*
typedef struct{ typedef struct{
size_t amount; // amount of tables in file size_t amount; // amount of tables in file
FITStable **tables; // array of pointers to tables FITStable **tables; // array of pointers to tables
} FITStables; } FITStables;
*/ */
/**
FITS image
*/
typedef struct{ typedef struct{
int width; // width int width; // width
int height; // height int height; // height
@ -95,16 +109,29 @@ typedef struct{
void *data; // picture data void *data; // picture data
} FITSimage; } FITSimage;
typedef union{
FITSimage *image;
FITStable *table;
} FITSimtbl;
/**
One of HDU
*/
typedef struct{
int hdutype; // type of current HDU: image/binary table/ACSII table/bad data
FITSimtbl contents;// data contents of HDU
KeyList *keylist; // keylist of given HDU
} FITSHDU;
typedef struct{ typedef struct{
fitsfile *fp; // cfitsio file structure fitsfile *fp; // cfitsio file structure
char *filename; // filename char *filename; // filename
int Nimages; // amount of images in file int NHDUs; // HDU amount
FITSimage **images; // image array FITSHDU *HDUs; // HDUs array itself
int Ntables; // amount of tables in file FITSHDU *curHDU; // pointer to current HDU
FITStable **tables; // table array
KeyList *keylist; // list of options for each key
} FITS; } FITS;
/*
typedef struct _Filter{ typedef struct _Filter{
char *name; // filter name char *name; // filter name
FType FilterType; // filter type FType FilterType; // filter type
@ -123,6 +150,7 @@ typedef enum{
,MATH_MEAN // calculate mean for all files ,MATH_MEAN // calculate mean for all files
,MATH_DIFF // difference of first and rest files ,MATH_DIFF // difference of first and rest files
} MathOper; } MathOper;
*/
void keylist_free(KeyList **list); void keylist_free(KeyList **list);
KeyList *keylist_add_record(KeyList **list, char *rec); KeyList *keylist_add_record(KeyList **list, char *rec);
@ -144,10 +172,13 @@ void table_print(FITStable *tbl);
void table_print_all(FITS *fits); void table_print_all(FITS *fits);
void image_free(FITSimage **ima); void image_free(FITSimage **ima);
FITSimage *image_new(size_t h, size_t w, int dtype); FITSimage *image_read(FITS *fits);
FITSimage *image_mksimilar(FITSimage *in, int dtype); #define image_datatype_size(d) (abs(d)/8)
void *image_data_malloc(size_t w, size_t h, int dtype);
FITSimage *image_new(size_t w, size_t h, int dtype);
FITSimage *image_mksimilar(FITSimage *in);
FITSimage *image_copy(FITSimage *in); FITSimage *image_copy(FITSimage *in);
FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata); //FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata);
void FITS_free(FITS **fits); void FITS_free(FITS **fits);
FITS *FITS_read(char *filename); FITS *FITS_read(char *filename);
@ -163,6 +194,7 @@ bool file_is_absent(char *name);
/*
// pointer to image conversion function // pointer to image conversion function
typedef FITS* (*imfuncptr)(FITS *in, Filter *f, Itmarray *i); typedef FITS* (*imfuncptr)(FITS *in, Filter *f, Itmarray *i);
*/

35
examples/keylist.c
View File

@ -20,9 +20,10 @@
#include <string.h> #include <string.h>
typedef struct{ typedef struct{
char *fitsname; char *fitsname; // input file name
int list; int list; // print keyword list
char **addrec; int contents; // print short description of file contents
char **addrec; // add some records to keywords in first HDU
} glob_pars; } glob_pars;
/* /*
@ -40,6 +41,7 @@ glob_pars G; /* = {
myoption cmdlnopts[] = { myoption cmdlnopts[] = {
// common options // common options
{"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")}, {"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), _("show this help")},
{"contents",NO_ARGS, NULL, 'c', arg_none, APTR(&G.contents), _("show short file contents")},
{"list", NO_ARGS, NULL, 'l', arg_none, APTR(&G.list), _("list all keywords")}, {"list", NO_ARGS, NULL, 'l', arg_none, APTR(&G.list), _("list all keywords")},
{"addrec", MULT_PAR, NULL, 'a', arg_string, APTR(&G.addrec), _("add record to file (you can add more than one record in once, point more -a)")}, {"addrec", MULT_PAR, NULL, 'a', arg_string, APTR(&G.addrec), _("add record to file (you can add more than one record in once, point more -a)")},
end_option end_option
@ -74,8 +76,33 @@ int main(int argc, char *argv[]){
green("Open file %s\n", G.fitsname); green("Open file %s\n", G.fitsname);
FITS *f = FITS_read(G.fitsname); FITS *f = FITS_read(G.fitsname);
if(!f) ERRX(_("Can't open file %s"), G.fitsname); if(!f) ERRX(_("Can't open file %s"), G.fitsname);
int N = f->NHDUs;
if(G.list){ if(G.list){
green("List of keywords:\n"); green("\n\nList of keywords:\n");
for(int i = 1; i <= N; ++i){
green("\nHDU #%d\n", i);
keylist_print(f->HDUs[i].keylist);
}
}
if(G.contents){
green("\n\nFile consists of %d HDUs:\n", N);
for(int i = 0; i <= N; ++i){
printf("\tHDU #%d - ", i);
switch(f->HDUs[i].hdutype){
case IMAGE_HDU:
printf("Image");
break;
case ASCII_TBL:
printf("ASCII table");
break;
case BINARY_TBL:
printf("Binary table");
break;
default:
printf("Some shit");
}
printf("\n");
}
} }
if(G.addrec){ if(G.addrec){
char **ptr = G.addrec; char **ptr = G.addrec;

239
fits.c
View File

@ -22,6 +22,7 @@
#include "FITSmanip.h" #include "FITSmanip.h"
#include <errno.h> #include <errno.h>
#include <stdlib.h>
#include <string.h> #include <string.h>
#include <sys/stat.h> #include <sys/stat.h>
#include <usefull_macros.h> #include <usefull_macros.h>
@ -249,10 +250,11 @@ void keylist_print(KeyList *list){
* @return keylist read * @return keylist read
*/ */
KeyList *keylist_read(FITS *fits){ KeyList *keylist_read(FITS *fits){
if(!fits || !fits->fp) return NULL; if(!fits || !fits->fp || !fits->curHDU) return NULL;
int fst, nkeys = -1, keypos = -1; int fst = 0, nkeys = -1, keypos = -1;
KeyList *list = fits->keylist; KeyList *list = fits->curHDU->keylist;
fits_get_hdrpos(fits->fp, &nkeys, &keypos, &fst); fits_get_hdrpos(fits->fp, &nkeys, &keypos, &fst);
DBG("nkeys=%d, keypos=%d, status=%d", nkeys, keypos, fst);
if(nkeys < 1){ if(nkeys < 1){
WARNX(_("No keywords in given HDU")); WARNX(_("No keywords in given HDU"));
return NULL; return NULL;
@ -349,6 +351,7 @@ FITStable *table_read(FITS *fits){
DBG("Table named %s with %ld rows and %d columns", extname, nrows, ncols); DBG("Table named %s with %ld rows and %d columns", extname, nrows, ncols);
FITStable *tbl = table_new(extname); FITStable *tbl = table_new(extname);
if(!tbl) return NULL; if(!tbl) return NULL;
fits->curHDU->contents.table = tbl;
for(i = 1; i <= ncols; ++i){ for(i = 1; i <= ncols; ++i){
int typecode; int typecode;
long repeat, width; long repeat, width;
@ -390,10 +393,6 @@ FITStable *table_read(FITS *fits){
//table_addcolumn(tbl, &col); //table_addcolumn(tbl, &col);
FREE(array); FREE(array);
} }
int N = fits->Ntables + 1;
if(!(fits->tables = realloc(fits->tables, sizeof(FITStable**)*N)))
ERR("realloc()");
fits->tables[fits->Ntables++] = tbl;
return tbl; return tbl;
} }
@ -600,55 +599,56 @@ void table_print(FITStable *tbl){
* @param fits - pointer to given file structure * @param fits - pointer to given file structure
*/ */
void table_print_all(FITS *fits){ void table_print_all(FITS *fits){
size_t i, N = fits->Ntables; size_t i, N = fits->NHDUs+1;
if(N == 0) return; if(N == 0) return;
for(i = 0; i < N; ++i) for(i = 1; i < N; ++i){
table_print(fits->tables[i]); if(fits->HDUs[i].hdutype == BINARY_TBL || fits->HDUs[i].hdutype == ASCII_TBL)
table_print(fits->HDUs[i].contents.table);
}
} }
/** /**
* @brief table_write - write tables to FITS file * @brief table_write - write tables to FITS file from current HDU
* @param fits (i) - pointer to FITS file structure * @param fits (i) - pointer to FITS file structure
*/ */
bool table_write(FITS *file){ bool table_write(FITS *file){
FNAME(); FNAME();
size_t N = file->Ntables, i;
if(N == 0) return FALSE;
fitsfile *fp = file->fp; fitsfile *fp = file->fp;
for(i = 0; i < N; ++i){ int fst = 0;
int hdutype = file->curHDU->hdutype;
if(hdutype != BINARY_TBL || hdutype != ASCII_TBL)
return FALSE;
FITStable *tbl = file->curHDU->contents.table;
size_t c, cols = tbl->ncols;
char **columns = MALLOC(char*, cols);
char **formats = MALLOC(char*, cols);
char **units = MALLOC(char*, cols);
table_column *col = tbl->columns;
for(c = 0; c < cols; ++c, ++col){
columns[c] = col->colname;
formats[c] = col->format;
units[c] = col->unit;
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,
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){
WARNX(_("Can't write table %s!"), tbl->tabname);
return FALSE;
}
//col = tbl->columns;
for(c = 0; c < cols; ++c, ++col){
DBG("write column %zd", c);
int fst = 0; int fst = 0;
FITStable *tbl = file->tables[i]; fits_write_col(fp, col->coltype, c+1, 1, 1, col->repeat, col->contents, &fst);
size_t c, cols = tbl->ncols; if(fst){
char **columns = MALLOC(char*, cols); fits_report_error(stderr, fst);
char **formats = MALLOC(char*, cols); WARNX(_("Can't write column %s!"), col->colname);
char **units = MALLOC(char*, cols);
table_column *col = tbl->columns;
for(c = 0; c < cols; ++c, ++col){
columns[c] = col->colname;
formats[c] = col->format;
units[c] = col->unit;
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, BINARY_TBL, 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){
WARNX(_("Can't write table %s!"), tbl->tabname);
return FALSE; return FALSE;
} }
//col = tbl->columns;
for(c = 0; c < cols; ++c, ++col){
DBG("write column %zd", c);
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);
WARNX(_("Can't write column %s!"), col->colname);
return FALSE;
}
}
} }
return TRUE; return TRUE;
} }
@ -661,14 +661,24 @@ void FITS_free(FITS **fits){
if(!fits || !*fits) return; if(!fits || !*fits) return;
FITS *f = *fits; FITS *f = *fits;
FREE(f->filename); FREE(f->filename);
int fst; int fst = 0;
fits_close_file(f->fp, &fst); fits_close_file(f->fp, &fst);
keylist_free(&f->keylist); int n, N = f->NHDUs;
int n; for(n = 1; n < N; ++n){
for(n = 0; n < f->Ntables; ++n) FITSHDU *hdu = &f->HDUs[n];
table_free(&f->tables[n]); keylist_free(&hdu->keylist);
for(n = 0; n < f->Nimages; ++n) switch(hdu->hdutype){
image_free(&f->images[n]); case IMAGE_HDU:
image_free(&hdu->contents.image);
break;
case BINARY_TBL:
case ASCII_TBL:
table_free(&hdu->contents.table);
break;
default: // do nothing
break;
}
}
FREE(*fits); FREE(*fits);
} }
@ -689,7 +699,7 @@ void FITS_free(FITS **fits){
*/ */
FITS *FITS_open(char *filename){ FITS *FITS_open(char *filename){
FITS *fits = MALLOC(FITS, 1); FITS *fits = MALLOC(FITS, 1);
int fst; int fst = 0;
// use fits_open_diskfile instead of fits_open_file to prevent using of extended name syntax // use fits_open_diskfile instead of fits_open_file to prevent using of extended name syntax
fits_open_diskfile(&fits->fp, filename, READWRITE, &fst); // READWRITE allows to modify files on-the-fly fits_open_diskfile(&fits->fp, filename, READWRITE, &fst); // READWRITE allows to modify files on-the-fly
if(fst){ if(fst){
@ -718,15 +728,21 @@ FITS *FITS_read(char *filename){
WARNX(_("Can't read HDU")); WARNX(_("Can't read HDU"));
goto returning; goto returning;
} }
fits->NHDUs = hdunum;
fits->HDUs = MALLOC(FITSHDU, hdunum+1);
int hdutype; int hdutype;
for(int i = 1; !(fits_movabs_hdu(fits->fp, i, &hdutype, &fst)); ++i){ for(int i = 1; i <= hdunum && !(fits_movabs_hdu(fits->fp, i, &hdutype, &fst)); ++i){
DBG("try to read keys"); FITSHDU *curHDU = &fits->HDUs[i];
keylist_read(fits); fits->curHDU = curHDU;
DBG("try to read keys from HDU #%d (type: %d)", i, hdutype);
curHDU->hdutype = hdutype;
curHDU->keylist = keylist_read(fits);
// types: IMAGE_HDU , ASCII_TBL, BINARY_TBL // types: IMAGE_HDU , ASCII_TBL, BINARY_TBL
DBG("HDU[%d] type %d", i, hdutype); DBG("HDU[%d] type %d", i, hdutype);
switch(hdutype){ switch(hdutype){
case IMAGE_HDU: case IMAGE_HDU:
DBG("Image"); DBG("Image");
curHDU->contents.image = image_read(fits);
break; break;
case BINARY_TBL: case BINARY_TBL:
DBG("Binary table"); DBG("Binary table");
@ -742,38 +758,6 @@ FITS *FITS_read(char *filename){
if(fst == END_OF_FILE){ if(fst == END_OF_FILE){
fst = 0; fst = 0;
}else goto returning; }else goto returning;
#if 0
// TODO: open not only images (liststruc.c from cexamples)!
// TODO: open not only 2-dimensional files!
// get image dimensions
int i, j, hdunum = 0, hdutype, nkeys, keypos, naxis;
long naxes[2];
char card[FLEN_CARD];
fits_get_img_param(fp, 2, &img->dtype, &naxis, naxes, &fst);
if(fst){fits_report_error(stderr, fst); goto returning;}
if(naxis > 2){
WARNX(_("Images with > 2 dimensions are not supported"));
fst = 1;
goto returning;
}
img->width = naxes[0];
img->height = naxes[1];
DBG("got image %ldx%ld pix, bitpix=%d", naxes[0], naxes[1], img->dtype);
// loop through all HDUs
if(fits_movabs_hdu(fp, imghdu, &hdutype, &fst)){
WARNX(_("Can't open image HDU #%d"), imghdu);
fst = 1;
goto returning;
}
size_t sz = naxes[0] * naxes[1];
img->data = MALLOC(double, sz);
int stat = 0;
fits_read_img(fp, TDOUBLE, 1, sz, NULL, img->data, &stat, &fst);
if(fst){fits_report_error(stderr, fst);}
if(stat) WARNX(_("Found %d pixels with undefined value"), stat);
DBG("ready");
#endif
returning: returning:
if(fst){ if(fst){
fits_report_error(stderr, fst); fits_report_error(stderr, fst);
@ -835,14 +819,28 @@ void image_free(FITSimage **img){
} }
/** /**
* create an empty image without headers, assign data type to "dtype" * @brief image_malloc - allocate memory for given data type
* @param w - image width
* @param h - image height
* @param dtype - data type
* @return allocated memory
*/ */
FITSimage *image_new(size_t h, size_t w, int dtype){ void *image_data_malloc(size_t w, size_t h, int dtype){
size_t bufsiz = w*h; void *data = calloc(w*h, image_datatype_size(dtype));
if(!data) ERR(_("calloc()"));
return data;
}
/**
* @brief image_new - create an empty image without headers, assign data type to "dtype"
* @param w - image width
* @param h - image height
* @param dtype - image data type
* @return
*/
FITSimage *image_new(size_t w, size_t h, int dtype){
FITSimage *out = MALLOC(FITSimage, 1); FITSimage *out = MALLOC(FITSimage, 1);
// TODO: make allocation when reading! out->data = image_data_malloc(w, h, dtype);
// TODO: allocate input data type ?
out->data = MALLOC(double, bufsiz);
out->width = w; out->width = w;
out->height = h; out->height = h;
out->dtype = dtype; out->dtype = dtype;
@ -851,13 +849,13 @@ FITSimage *image_new(size_t h, size_t w, int dtype){
/** /**
* build IMAGE image from data array indata * build IMAGE image from data array indata
*/ *
FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata){ FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata){
size_t stride = 0; size_t stride = 0;
double (*fconv)(uint8_t *data) = NULL; double (*fconv)(uint8_t *data) = NULL;
double ubyteconv(uint8_t *data){return (double)*data;} double ubyteconv(uint8_t *data){return (double)*data;}
double ushortconv(uint8_t *data){return (double)*(int16_t*)data;} double ushortconv(uint8_t *data){return (double)*(short*)data;}
double ulongconv(uint8_t *data){return (double)*(uint32_t*)data;} double ulongconv(uint8_t *data){return (double)*(unsigned long*)data;}
double ulonglongconv(uint8_t *data){return (double)*(uint64_t*)data;} double ulonglongconv(uint8_t *data){return (double)*(uint64_t*)data;}
double floatconv(uint8_t *data){return (double)*(float*)data;} double floatconv(uint8_t *data){return (double)*(float*)data;}
FITSimage *out = image_new(h, w, dtype); FITSimage *out = image_new(h, w, dtype);
@ -866,6 +864,7 @@ FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata){
stride = 1; stride = 1;
fconv = ubyteconv; fconv = ubyteconv;
break; break;
//case USHORT_IMG?
case SHORT_IMG: case SHORT_IMG:
stride = 2; stride = 2;
fconv = ushortconv; fconv = ushortconv;
@ -901,29 +900,55 @@ FITSimage *image_build(size_t h, size_t w, int dtype, uint8_t *indata){
*dout++ = fconv(din); *dout++ = fconv(din);
} }
return out; return out;
} }*/
/** /**
* create an empty copy of image "in" without headers, assign data type to "dtype" * create an empty copy of image "in" without headers, assign data type to "dtype"
*/ */
FITSimage *image_mksimilar(FITSimage *img, int dtype){ FITSimage *image_mksimilar(FITSimage *img){
size_t w = img->width, h = img->height, bufsiz = w*h; if(!img || img->height < 1 || img->width < 1) return NULL;
FITSimage *out = MALLOC(FITSimage, 1); return image_new(img->width, img->height, img->dtype);
// TODO: allocate buffer as in original!
out->data = MALLOC(double, bufsiz);
out->width = w;
out->height = h;
out->dtype = dtype;
return out;
} }
/** /**
* make full copy of image 'in' * make full copy of image 'in'
*/ */
FITSimage *image_copy(FITSimage *in){ FITSimage *image_copy(FITSimage *in){
FITSimage *out = image_mksimilar(in, in->dtype); FITSimage *out = image_mksimilar(in);
if(!out) return NULL;
// TODO: size of data as in original! // TODO: size of data as in original!
memcpy(out->data, in->data, sizeof(double)*in->width*in->height); memcpy(out->data, in->data, sizeof(double)*in->width*in->height);
return out; return out;
} }
/**
* @brief image_read - read image from current HDU
* @param fits - fits structure pointer
* @return - pointer to allocated image structure or NULL if failed
*/
FITSimage *image_read(FITS *fits){
// TODO: open not only 2-dimensional files!
// get image dimensions
int naxis, fst = 0, dtype;
long naxes[2] = {0,0};
fits_get_img_param(fits->fp, 2, &dtype, &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;
}
DBG("got image %ldx%ld pix, bitpix=%d", naxes[0], naxes[1], dtype);
FITSimage *img = image_new(naxes[0], naxes[1], dtype);
int stat = 0;
fits_read_img(fits->fp, dtype, 1, naxes[0]*naxes[1], NULL, img->data, &stat, &fst);
if(fst){
fits_report_error(stderr, fst);
image_free(&img);
return NULL;
}
if(stat) WARNX(_("Found %d pixels with undefined value"), stat);
DBG("ready");
return img;
}