First commit with changes to enable reading of tile-compressed images using CFitsIO (3f06d390) · Commits · csst-pipeline / msc / sextractor

src/back.c

+48 −15

Original line number	Diff line number	Diff line
		@@ -52,8 +52,9 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)

		{
		backstruct backmesh,wbackmesh, bm,wbm;
		tabstruct tab, wtab;
		PIXTYPE buf,wbuf, buft,wbuft;
		OFF_T fcurpos,wfcurpos, wfcurpos2,fcurpos2, bufshift, jumpsize;
		OFF_T2 fcurpos,wfcurpos, wfcurpos2,fcurpos2, bufshift, jumpsize;
		size_t bufsize, bufsize2,
		size,meshsize;
		int i,j,k,m,n, step, nlines,
		@@ -65,7 +66,11 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		/* If the weight-map is not an external one, no stats are needed for it */
		if (wfield && wfield->flags&(INTERP_FIELD\|BACKRMS_FIELD))
		wfield= NULL;

		tab = field->tab;
		if (wfield)
		wtab = wfield->tab;
		else
		wtab = NULL; /* to avoid gcc -Wall warnings */
		w = field->width;
		bw = field->backw;
		bh = field->backh;
		@@ -83,12 +88,17 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)

		wfcurpos = wfcurpos2 = 0; /* to avoid gcc -Wall warnings */
		QFTELL(field->file, fcurpos, field->filename);
		tab->currentElement = 1; // CFITSIO

		if (wfield)
		{
		QFTELL(wfield->file, wfcurpos, wfield->filename);
		wtab->currentElement = 1; // CFITSIO
		}

		/* Allocate a correct amount of memory to store pixels */

		bufsize = (OFF_T)w*bh;
		bufsize = (OFF_T2)w*bh;
		meshsize = (size_t)bufsize;
		nlines = 0;
		if (bufsize > (size_t)BACK_BUFSIZE)
		@@ -96,8 +106,8 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		nlines = BACK_BUFSIZE/w;
		step = (field->backh-1)/nlines+1;
		bufsize = (size_t)(nlines = field->backh/step)*w;
		bufshift = (step/2)*(OFF_T)w;
		jumpsize = (step-1)*(OFF_T)w;
		bufshift = (step/2)*(OFF_T2)w;
		jumpsize = (step-1)*(OFF_T2)w;
		}
		else
		bufshift = jumpsize = 0; /* to avoid gcc -Wall warnings */
		@@ -180,8 +190,8 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		nlines = BACK_BUFSIZE/w;
		step = (n-1)/nlines+1;
		bufsize = (nlines = n/step)*(size_t)w;
		bufshift = (step/2)*(OFF_T)w;
		jumpsize = (step-1)*(OFF_T)w;
		bufshift = (step/2)*(OFF_T2)w;
		jumpsize = (step-1)*(OFF_T2)w;
		free(buf);
		QMALLOC(buf, PIXTYPE, bufsize); /* pixel buffer */
		if (wfield)
		@@ -192,35 +202,47 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		}

		/---- Read and skip, read and skip, etc... /
		QFSEEK(field->file, bufshift*(OFF_T)field->bytepix, SEEK_CUR,
		QFSEEK(field->file, bufshift*(OFF_T2)field->bytepix, SEEK_CUR,
		field->filename);
		tab->currentElement += bufshift; // CFITSIO

		buft = buf;
		for (i=nlines; i--; buft += w)
		{
		read_body(field->tab, buft, w);
		if (i)
		QFSEEK(field->file, jumpsize*(OFF_T)field->bytepix, SEEK_CUR,
		if (i) {
		QFSEEK(field->file, jumpsize*(OFF_T2)field->bytepix, SEEK_CUR,
		field->filename);

		tab->currentElement += jumpsize; // CFITSIO
		}
		}

		if (wfield)
		{
		/------ Read and skip, read and skip, etc... now on the weight-map /
		QFSEEK(wfield->file, bufshift*(OFF_T)wfield->bytepix, SEEK_CUR,
		QFSEEK(wfield->file, bufshift*(OFF_T2)wfield->bytepix, SEEK_CUR,
		wfield->filename);
		wtab->currentElement += bufshift; // CFITSIO

		wbuft = wbuf;
		for (i=nlines; i--; wbuft += w)
		{
		read_body(wfield->tab, wbuft, w);
		weight_to_var(wfield, wbuft, w);
		if (i)
		QFSEEK(wfield->file, jumpsize*(OFF_T)wfield->bytepix, SEEK_CUR,
		if (i){
		QFSEEK(wfield->file, jumpsize*(OFF_T2)wfield->bytepix, SEEK_CUR,
		wfield->filename);
		wtab->currentElement += jumpsize; // CFITSIO

		}
		}
		}
		backstat(backmesh, wbackmesh, buf, wbuf, bufsize, nx, w, bw,
		wfield?wfield->weight_thresh:0.0);
		QFSEEK(field->file, fcurpos2, SEEK_SET, field->filename);
		tab->currentElement = (fcurpos2 == 0) ? 1 : fcurpos2; // CFITSIO

		bm = backmesh;
		for (m=nx; m--; bm++)
		if (bm->mean <= -BIG)
		@@ -230,6 +252,8 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		if (wfield)
		{
		QFSEEK(wfield->file, wfcurpos2, SEEK_SET, wfield->filename);
		wtab->currentElement = (wfcurpos2 == 0) ? 1 : wfcurpos2; // CFITSIO

		wbm = wbackmesh;
		for (m=nx; m--; wbm++)
		if (wbm->mean <= -BIG)
		@@ -242,10 +266,14 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)
		{
		if (bufsize2>size)
		bufsize2 = size;

		read_body(field->tab, buf, bufsize2);

		if (wfield)
		{

		read_body(wfield->tab, wbuf, bufsize2);

		weight_to_var(wfield, wbuf, bufsize2);
		}
		backhisto(backmesh, wbackmesh, buf, wbuf, bufsize2, nx, w, bw,
		@@ -284,8 +312,13 @@ void makeback(picstruct field, picstruct wfield, int wscale_flag)

		/* Go back to the original position */
		QFSEEK(field->file, fcurpos, SEEK_SET, field->filename);
		if (wfield)
		field->tab->currentElement = fcurpos; // CFITSIO

		if (wfield) {
		QFSEEK(wfield->file, wfcurpos, SEEK_SET, wfield->filename);
		wfield->tab->currentElement = (wfcurpos == 0) ? 1 : wfcurpos; // CFITSIO

		}

		/* Median-filter and check suitability of the background map */
		NFPRINTF(OUTPUT, "Filtering background map(s)");

src/catout.c

+1 −1

Original line number	Diff line number	Diff line
		@@ -1050,7 +1050,7 @@ void reendcat()
		{
		keystruct *key;
		tabstruct *tab;
		OFF_T pos;
		OFF_T2 pos;
		char *head;

		switch(prefs.cat_type)

src/field.c

+2 −1

Original line number	Diff line number	Diff line
		@@ -70,6 +70,7 @@ picstruct newfield(char filename, int flags, int ext)
		field->cat = cat;
		nok = 0;
		tab = cat->tab;
		if(tab->isTileCompressed) nok++;
		if (tab->naxis >= 2
		&& strncmp(tab->xtension, "BINTABLE", 8)
		&& strncmp(tab->xtension, "ASCTABLE", 8))
		@@ -78,6 +79,7 @@ picstruct newfield(char filename, int flags, int ext)
		for (ntab=cat->ntab; ext2-- && ntab--;)
		{
		tab=tab->nexttab;
		if(tab->isTileCompressed) nok++;
		if (tab->naxis >= 2
		&& strncmp(tab->xtension, "BINTABLE", 8)
		&& strncmp(tab->xtension, "ASCTABLE", 8))
		@@ -85,7 +87,6 @@ picstruct newfield(char filename, int flags, int ext)
		}
		if (!nok)
		error(EXIT_FAILURE, "Not a valid FITS image in ",filename);

		field->tab = tab;
		if (ntab<0)
		error(EXIT_FAILURE, "Not enough valid FITS image extensions in ",filename);

src/fits/fitsbody.c

+78 −4

Original line number	Diff line number	Diff line
		@@ -85,7 +85,10 @@ PIXTYPE alloc_body(tabstruct tab, void (func)(PIXTYPE ptr, int npix))
		tab->extname);

		/* Decide if the data will go in physical memory or on swap-space */
		npix = tab->tabsize/tab->bytepix;


		//npix = tab->tabsize/tab->bytepix;
		npix = tab->naxisn[0] * tab->naxisn[1];
		size = npix*sizeof(PIXTYPE);
		if (size < body_ramleft)
		{
		@@ -93,6 +96,7 @@ PIXTYPE alloc_body(tabstruct tab, void (func)(PIXTYPE ptr, int npix))
		if ((tab->bodybuf = malloc(size)))
		{
		QFSEEK(tab->cat->file, tab->bodypos, SEEK_SET, tab->cat->filename);
		tab->currentElement = 1; // CFITSIO
		read_body(tab, (PIXTYPE *)tab->bodybuf, npix);
		/---- Apply pixel processing /
		if (func)
		@@ -120,6 +124,7 @@ PIXTYPE alloc_body(tabstruct tab, void (func)(PIXTYPE ptr, int npix))
		if (!spoonful)
		spoonful = DATA_BUFSIZE;
		QFSEEK(tab->cat->file, tab->bodypos, SEEK_SET, tab->cat->filename);
		tab->currentElement = 1; // CFITSIO
		read_body(tab, buffer, spoonful/sizeof(PIXTYPE));
		/-- Apply pixel processing /
		if (func)
		@@ -323,14 +328,20 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		int curval, dval, blankflag, ival, iblank;

		size_t i, bowl, spoonful, npix;
		PIXTYPE bs,bz;
		//PIXTYPE bs,bz;
		double bs,bz;

		/* a NULL cat structure indicates that no data can be read */
		if (!(cat = tab->cat))
		return;

		bs = (PIXTYPE)tab->bscale;
		bz = (PIXTYPE)tab->bzero;
		// this cast from double to float loses precision
		//bs = (PIXTYPE)tab->bscale;
		//bz = (PIXTYPE)tab->bzero;

		bs = tab->bscale;
		bz = tab->bzero;

		blankflag = tab->blankflag;

		switch(tab->compress_type)
		@@ -344,7 +355,41 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		if (spoonful>size)
		spoonful = size;
		bufdata = (char *)bufdata0;

		// CFITSIO
		if (tab->isTileCompressed) {

		int status, hdutype;

		// first of all, move to correct HDU
		status = 0; fits_movabs_hdu(tab->infptr, tab->hdunum, &hdutype, &status);
		if (status != 0) {

		printf("Error moving to HDU %d\n", tab->hdunum);
		fits_report_error(stderr, status);
		}

		// pixels count from 1
		if (tab->currentElement == 0) tab->currentElement = 1;

		// now read section of image
		status = 0; fits_read_img(tab->infptr, TFLOAT, tab->currentElement, spoonful, NULL, bufdata, NULL, &status);

		//printf("CFITSIO OK reading start=%d end=%d absolute end=%d\n", tab->currentElement, (tab->currentElement + spoonful) , (tab->naxisn[0]*tab->naxisn[1]));

		// report reading error
		if (status != 0) {

		printf("CFITSIO ERROR reading start=%d end=%d absolute end=%d\n", tab->currentElement, (tab->currentElement + spoonful) , (tab->naxisn[0]*tab->naxisn[1]));
		fits_report_error(stderr, status);
		}

		// update file 'pointer'
		tab->currentElement += spoonful;
		}
		else
		QFREAD(bufdata, spoonful*tab->bytepix, cat->file, cat->filename);

		switch(tab->bitpix)
		{
		case BP_BYTE:
		@@ -381,6 +426,7 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		break;

		case BP_SHORT:
		if (!tab->isTileCompressed)
		if (bswapflag)
		swapbytes(bufdata, 2, spoonful);
		if (blankflag)
		@@ -416,6 +462,7 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		break;

		case BP_LONG:
		if (!tab->isTileCompressed)
		if (bswapflag)
		swapbytes(bufdata, 4, spoonful);
		if (blankflag)
		@@ -452,6 +499,7 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)

		#ifdef HAVE_LONG_LONG_INT
		case BP_LONGLONG:
		if (!tab->isTileCompressed)
		if (bswapflag)
		swapbytes(bufdata, 8, spoonful);
		if (blankflag)
		@@ -487,6 +535,7 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		break;
		#endif
		case BP_FLOAT:
		if (!tab->isTileCompressed)
		if (bswapflag)
		swapbytes(bufdata, 4, spoonful);
		#pragma ivdep
		@@ -497,6 +546,7 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		case BP_DOUBLE:
		if (bswapflag)
		{
		if (!tab->isTileCompressed)
		swapbytes(bufdata, 8, spoonful);
		#pragma ivdep
		for (i=spoonful; i--; bufdata += sizeof(double))
		@@ -615,8 +665,12 @@ void read_body(tabstruct tab, PIXTYPE ptr, size_t size)
		default:
		error(EXIT_FAILURE,"Internal Error: unknown compression mode in ",
		"read_body()");


		}

		//printf("SSSS %f %f %f %f %f\n", bufdata[0], bufdata[10], bufdata[100], bufdata[1000], bufdata[spoonful-1]);

		return;
		}

		@@ -922,6 +976,10 @@ void write_body(tabstruct tab, PIXTYPE ptr, size_t size)
		#pragma ivdep
		for (i=spoonful; i--;)
		(bufdata++) = ((ptr++)-bz)/bs;


		// CFITSIO - only perform byte-swap if we are NOT writing a tile-compressed format using cfitsio
		if (0 && tab->infptr == NULL) // TODO
		if (bswapflag)
		swapbytes(cbufdata0, 4, spoonful);
		}
		@@ -943,6 +1001,22 @@ void write_body(tabstruct tab, PIXTYPE ptr, size_t size)
		"read_body()");
		break;
		}

		// CFITSIO - if cfitsio output file has been set up, then proceed to write using cfitsio
		if (0 && tab->infptr != NULL) { // TODO

		int status = 0; fits_write_img(tab->infptr, TFLOAT, tab->currentElement, spoonful, cbufdata0, &status);

		if (status != 0) {

		printf("CFITSIO ERROR writing start=%d end=%d absolute end=%d\n", tab->currentElement, (tab->currentElement + spoonful) , (tab->naxisn[0]*tab->naxisn[1]));
		fits_report_error(stderr, status);
		}

		tab->currentElement = tab->currentElement + spoonful;
		}
		// otherwise, continue with usual AstrOmatic fits writing routine
		else
		QFWRITE(cbufdata0, spoonful*tab->bytepix, cat->file, cat->filename);
		}
		break;

src/fits/fitscat.c

+30 −2

Original line number	Diff line number	Diff line
		@@ -323,14 +323,39 @@ int map_cat(catstruct *cat)
		QCALLOC(tab, tabstruct, 1);
		tab->cat = cat;
		QFTELL(cat->file, tab->headpos, cat->filename);

		// CFITSIO
		fitsfile *infptr;
		int status, hdutype, hdunum;
		status = 0; fits_open_file(&infptr, cat->filename, READONLY, &status);
		if (status != 0) {
		fits_report_error(stderr, status);
		printf("ERROR could not open FITS file with cfitsio: %s\n", cat->filename);
		}
		hdunum = 1;

		for (ntab=0; !get_head(tab); ntab++)
		{
		readbasic_head(tab);
		readbintabparam_head(tab);
		QFTELL(cat->file, tab->bodypos, cat->filename);
		tab->nseg = tab->seg = 1;
		if (tab->tabsize)
		QFSEEK(cat->file, PADTOTAL(tab->tabsize), SEEK_CUR, cat->filename);

		// CFITSIO
		tab->hdunum = hdunum;
		tab->infptr = infptr;
		status = 0; fits_movabs_hdu(tab->infptr, tab->hdunum, &hdutype, &status);
		if (status != 0) printf("ERROR could not move to hdu %d in file %s\n", tab->hdunum, cat->filename);
		tab->tabsize = infptr->Fptr->rowlength;

		if (tab->tabsize) {

		// IMPORTANT: moving to start of next header using fseek and cfitsio position rather than table size, as done previously
		fseek(cat->file, infptr->Fptr->headstart[hdunum], SEEK_SET);

		// this is how it was done previously
		//QFSEEK(cat->file, PADTOTAL(tab->tabsize), SEEK_CUR, cat->filename);
		}
		if (prevtab)
		{
		tab->prevtab = prevtab;
		@@ -342,6 +367,9 @@ int map_cat(catstruct *cat)
		QCALLOC(tab, tabstruct, 1);
		tab->cat = cat;
		QFTELL(cat->file, tab->headpos, cat->filename);

		// CFITSIO
		hdunum++;
		}

		cat->ntab = ntab;