Merge pull request #10 from teake/teake/merge-roy

Rice compression support

configure.ac

@@ -53,6 +53,7 @@ sinclude(acx_atlas.m4)
 sinclude(acx_openblas.m4)
 sinclude(acx_fftw.m4)
 sinclude(acx_mkl.m4)
+sinclude(acx_cfitsio.m4)
 sinclude(acx_prog_cc_optim.m4)
 sinclude(acx_pthread.m4)
 sinclude(acx_urbi_resolve_dir.m4)
@@ -180,6 +181,20 @@ AC_ARG_WITH(openblas-incdir,
 	[AS_HELP_STRING([--with-openblas-incdir=<OpenBLAS header dir>],
 	[Provide an alternative path to the OpenBLAS header directory])])
 
+# Provide special option for CFITSIO
+AC_MSG_CHECKING([whether CFITSIO support is enabled])
+AC_ARG_ENABLE(cfitsio,
+	[AS_HELP_STRING([--enable-cfitsio],
+	[Enable support for compressed FITS files through the CFITSIO library (default = no)])],
+	AC_MSG_RESULT([yes]),
+	AC_MSG_RESULT([no]))
+AC_ARG_WITH(cfitsio-libdir,
+	[AS_HELP_STRING([--with-cfitsio-libdir=<CFITSIO library path>],
+	[Provide an alternative path to the CFITSIO library])])
+AC_ARG_WITH(cfitsio-incdir,
+	[AS_HELP_STRING([--with-cfitsio-incdir=<CFITSIO header dir>],
+	[Provide an alternative path to the CFITSIO header directory])])
+
 # Provide a special option for the default XSLT URL
 with_xsl_url="file://"$(URBI_RESOLVE_DIR([$datadir]))"/$PACKAGE_TARNAME/$PACKAGE_TARNAME.xsl"
 AC_ARG_WITH(xsl_url,
@@ -305,6 +320,20 @@ if test "$enable_model_fitting" != "no"; then
   fi
 fi
 
+AS_IF([test "x$enable_cfitsio" = "xyes"], [
+	ACX_CFITSIO($with_cfitsio_libdir, $with_cfitsio_incdir,, no,
+    [
+      AM_CFLAGS="$AM_CFLAGS $CFITSIO_CFLAGS "
+      AM_LDFLAGS="$AM_LDFLAGS $CFITSIO_LDFLAGS "
+      LIBS="$CFITSIO_LIBS $LIBS"
+      if test "$CFITSIO_WARN" != ""; then
+        AC_MSG_WARN([$CFITSIO_WARN])
+      fi
+    ],
+    AC_MSG_ERROR([$CFITSIO_ERROR Exiting.])
+    )
+])
+
 AM_CONDITIONAL(USE_MODEL, [test "$enable_model_fitting" != "no"])
 
 # Compile with profiling option

m4/acx_cfitsio.m4

+dnl
+dnl				acx_cfitsio.m4
+dnl
+dnl Figure out if the CFITSIO library and header files are installed.
+dnl
+dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dnl
+dnl	This file part of:	AstrOmatic software
+dnl
+dnl	Copyright:		(C) 2019 IAP/CNRS/UPMC
+dnl
+dnl	License:		GNU General Public License
+dnl
+dnl	AstrOmatic software is free software: you can redistribute it and/or
+dnl	modify it under the terms of the GNU General Public License as
+dnl	published by the Free Software Foundation, either version 3 of the
+dnl	License, or (at your option) any later version.
+dnl	AstrOmatic software is distributed in the hope that it will be useful,
+dnl	but WITHOUT ANY WARRANTY; without even the implied warranty of
+dnl	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+dnl	GNU General Public License for more details.
+dnl	You should have received a copy of the GNU General Public License
+dnl	along with AstrOmatic software.
+dnl	If not, see <http://www.gnu.org/licenses/>.
+dnl
+dnl	Last modified:		16/07/2019
+dnl
+dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dnl
+dnl @synopsis ACX_CFITSIO([CFITSIO_LIBSDIR, CFITSIO_INCDIR, CFITSIO_PFLAG,
+dnl                  ILP64_FLAG, [ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]]])
+dnl
+dnl You may wish to use these variables in your default LIBS:
+dnl
+dnl        LIBS="$CFITSIO_LIBS $LIBS"
+dnl
+dnl ACTION-IF-FOUND is a list of shell commands to run if CFITSIO
+dnl is found (HAVE_CFITSIO is defined first), and ACTION-IF-NOT-FOUND
+dnl is a list of commands to run it if it is not found.
+
+AC_DEFUN([ACX_CFITSIO], [
+AC_REQUIRE([AC_CANONICAL_HOST])
+
+dnl --------------------
+dnl Search include files
+dnl --------------------
+
+CFITSIO_ERROR=""
+if test x$2 = x; then
+  [acx_cfitsio_incdir="/"]
+  AC_CHECK_HEADERS(
+    [${acx_cfitsio_incdir}fitsio.h],,
+    [
+      [acx_cfitsio_incdir=""]
+      AC_CHECK_HEADER(
+        [fitsio.h],,
+        [CFITSIO_ERROR="CFITSIO header files not found!"]
+      )
+    ]
+  )
+else
+  acx_cfitsio_incdir="$2/"
+  AC_CHECK_HEADER(
+    [${acx_cfitsio_incdir}fitsio.h],,
+    [
+      [acx_cfitsio_incdir="$2/include/"]
+      AC_CHECK_HEADERS(
+        [${acx_cfitsio_incdir}fitsio.h],,
+        [CFITSIO_ERROR="CFITSIO header files not found in "$2"!"]
+    )]
+  )
+fi
+
+if test "x$CFITSIO_ERROR" = "x"; then
+  AC_DEFINE_UNQUOTED(FITSIO_H, "${acx_cfitsio_incdir}fitsio.h", [CFITSIO header filename.])
+
+dnl ----------------------------
+dnl Search CFITSIO library file
+dnl ----------------------------
+
+  OLIBS="$LIBS"
+  LIBS=""
+  if test x$4 = xyes; then
+    acx_cfitsio_suffix="64"
+    CFITSIO_CFLAGS="-DCFITSIO_USE64BITINT -DLAPACK_ILP64"
+  else
+    acx_cfitsio_suffix=""
+    CFITSIO_CFLAGS=""
+  fi
+  if test x$1 = x; then
+    acx_cfitsio_libopt=""
+  else
+    acx_cfitsio_libopt="-L$1"
+  fi
+  AC_SEARCH_LIBS(
+      [ffmahd], ["cfitsio"$acx_cfitsio_suffix],,
+      [CFITSIO_ERROR="CFITSIO"$acx_cfitsio_suffix" library file not found!"],
+      $acx_cfitsio_libopt
+    )
+  LIBS="$OLIBS"
+fi
+
+dnl -------------------------------------------------------------------------
+dnl Finally execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND
+dnl -------------------------------------------------------------------------
+
+if test "x$CFITSIO_ERROR" = "x"; then
+  AC_DEFINE(HAVE_CFITSIO,1, [Define if you have the CFITSIO library and header files.])
+  CFITSIO_LIBS="$acx_cfitsio_libopt $ac_cv_search_ffmahd"
+  AC_SUBST(CFITSIO_CFLAGS)
+  AC_SUBST(CFITSIO_LDFLAGS, "")
+  AC_SUBST(CFITSIO_LIBS)
+  AC_SUBST(CFITSIO_WARN)
+  $5
+else
+  AC_SUBST(CFITSIO_ERROR)
+  $6
+fi
+
+])dnl ACX_CFITSIO
+

src/analyse.c

@@ -7,7 +7,7 @@
 *
 *	This file part of:	SExtractor
 *
-*	Copyright:		(C) 1993-2016 IAP/CNRS/UPMC
+*	Copyright:		(C) 1993-2016 Emmanuel Bertin -- IAP/CNRS/UPMC
 *
 *	License:		GNU General Public License
 *
@@ -280,8 +280,7 @@ void  examineiso(picstruct *field, picstruct *dfield, objstruct *obj,
     emy2 /= flux2;	/* variance of ym */
     emxy /= flux2;	/* covariance */
 
-/*-- Handle fully correlated profile
-s (which cause a singularity...) */
+/*-- Handle fully correlated profiles (which cause a singularity...) */
     esum *= 0.08333/flux2;
     if (obj->singuflag && (emx2*emy2-emxy*emxy) < esum*esum)
       {
@@ -718,10 +717,10 @@ void	endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
         dgeo_copy(dgeofield, outobj2.vignet_dgeox , outobj2.vignet_dgeoy,
 		prefs.vignet_dgeoxsize[0],prefs.vignet_dgeoxsize[1], ix,iy);
       else {
-        if (FLAG(obj2.vignet_dgeox)) 
+        if (FLAG(obj2.vignet_dgeox))
           dgeo_copy(dgeofield, outobj2.vignet_dgeox, NULL,
 		prefs.vignet_dgeoxsize[0],prefs.vignet_dgeoxsize[1], ix,iy);
-        if (FLAG(obj2.vignet_dgeoy)) 
+        if (FLAG(obj2.vignet_dgeoy))
           dgeo_copy(dgeofield, NULL, outobj2.vignet_dgeoy,
 		prefs.vignet_dgeoysize[0],prefs.vignet_dgeoysize[1], ix,iy);
       }

src/back.c

@@ -52,8 +52,10 @@ 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;
+   OFF_T2	currentElement, wcurrentElement, currentElement2, wcurrentElement2;
    size_t	bufsize, bufsize2,
 		size,meshsize;
    int		i,j,k,m,n, step, nlines,
@@ -65,7 +67,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 +89,17 @@ void	makeback(picstruct *field, picstruct *wfield, int wscale_flag)
 
   wfcurpos = wfcurpos2 = 0;		/* to avoid gcc -Wall warnings */
   QFTELL(field->file, fcurpos, field->filename);
+  currentElement = (tab->currentElement == 0) ? 1 : tab->currentElement; // CFITSIO
+
   if (wfield)
+    {
     QFTELL(wfield->file, wfcurpos, wfield->filename);
+    wcurrentElement = (wtab->currentElement == 0) ? 1 : wtab->currentElement; // 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 +107,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 */
@@ -172,16 +183,19 @@ void	makeback(picstruct *field, picstruct *wfield, int wscale_flag)
       {
 /*---- Image size too big, we have to skip a few data !*/
       QFTELL(field->file, fcurpos2, field->filename);
-      if (wfield)
+      currentElement2 = (tab->currentElement == 0) ? 1 : tab->currentElement; // CFITSIO
+      if (wfield){
         QFTELL(wfield->file, wfcurpos2, wfield->filename);
+        wcurrentElement2 = (wtab->currentElement == 0) ? 1 : wtab->currentElement; // CFITSIO
+      }
       if (j == ny-1 && (n=field->height%field->backh))
         {
         meshsize = n*(size_t)w;
         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 +206,46 @@ 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 = currentElement2; // CFITSIO
+
       bm = backmesh;
       for (m=nx; m--; bm++)
         if (bm->mean <= -BIG)
@@ -230,6 +255,7 @@ void	makeback(picstruct *field, picstruct *wfield, int wscale_flag)
       if (wfield)
         {
         QFSEEK(wfield->file, wfcurpos2, SEEK_SET, wfield->filename);
+        wtab->currentElement = wcurrentElement2; // CFITSIO
         wbm = wbackmesh;
         for (m=nx; m--; wbm++)
           if (wbm->mean <= -BIG)
@@ -284,8 +310,11 @@ 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)
+  tab->currentElement = currentElement; // CFITSIO
+  if (wfield) {
     QFSEEK(wfield->file, wfcurpos, SEEK_SET, wfield->filename);
+    wfield->tab->currentElement =  wcurrentElement; // CFITSIO
+  }
 
 /* Median-filter and check suitability of the background map */
   NFPRINTF(OUTPUT, "Filtering background map(s)");

src/catout.c

@@ -1059,7 +1059,7 @@ void	reendcat()
   {
    keystruct	*key;
    tabstruct	*tab;
-   OFF_T	pos;
+   OFF_T2	pos;
    char		*head;
 
   switch(prefs.cat_type)
@@ -1127,4 +1127,3 @@ void	zerocat(void)
   }
 
 
-

src/field.c

@@ -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))

src/fits/fitsbody.c

@@ -85,7 +85,7 @@ 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->naxisn[0] * tab->naxisn[1];
   size = npix*sizeof(PIXTYPE);
   if (size < body_ramleft)
     {
@@ -93,6 +93,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 +121,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)
@@ -294,6 +296,97 @@ void	free_body(tabstruct *tab)
   return;
   }
 
+#ifdef HAVE_CFITSIO
+void writeCfitsio(tabstruct *tab, size_t spoonful, void *cbufdata0) {
+
+    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;
+}
+
+
+void readTileCompressed(tabstruct *tab,  size_t	spoonful, double* bufdata0) {
+
+    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
+    int datatype;
+    switch(tab->bitpix){
+        case BYTE_IMG:
+            datatype = TBYTE;
+            break;
+        case SHORT_IMG:
+            datatype = TSHORT;
+            break;
+        case LONG_IMG:
+            datatype = TLONG;
+            break;
+        case FLOAT_IMG:
+            datatype = TFLOAT;
+            break;
+        case DOUBLE_IMG:
+            datatype = TDOUBLE;
+            break;
+    }
+
+    int anynul;
+    double bscale = 1.0, bzero = 0.0, nulval = 0.;
+
+    // turn off any scaling so that we copy raw pixel values
+    status = 0; fits_set_bscale(tab->infptr,  bscale, bzero, &status);
+
+    // now read the image
+    status = 0; fits_read_img(tab->infptr, datatype,  tab->currentElement, spoonful, &nulval, bufdata0, &anynul, &status);
+
+    // 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;
+}
+#endif
+
+
+void readSpoonful(tabstruct *tab, size_t spoonful, double* bufdata0, char* bufdata) {
+#ifdef HAVE_CFITSIO
+    if(tab->isTileCompressed)
+        readTileCompressed(tab, spoonful, bufdata0);
+    else
+#endif
+        QFREAD(bufdata, spoonful*tab->bytepix, tab->cat->file, tab->cat->filename);
+}
+
+void writeSpoonful(tabstruct *tab, size_t spoonful, void *cbufdata0) {
+#ifdef HAVE_CFITSIO
+    if (0 && tab->infptr != NULL)
+        writeCfitsio(tab, spoonful, cbufdata0);
+    else
+#endif
+        QFWRITE(cbufdata0, spoonful * tab->bytepix, tab->cat->file, tab->cat->filename);
+
+}
+
 
 /******* read_body ************************************************************
 PROTO	read_body(tabstruct *tab, PIXTYPE *ptr, long size)
@@ -323,15 +416,21 @@ void	read_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
   int			curval, dval, blankflag, bswapflag, 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;
   bswapflag = *((char *)&ashort);	// Byte-swapping flag
 
@@ -346,7 +445,8 @@ void	read_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
         if (spoonful>size)
           spoonful = size;
         bufdata = (char *)bufdata0;
-        QFREAD(bufdata, spoonful*tab->bytepix, cat->file, cat->filename);
+        readSpoonful(tab, spoonful, bufdata0, bufdata);
+
         switch(tab->bitpix)
           {
           case BP_BYTE:
@@ -383,6 +483,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)
@@ -418,6 +519,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)
@@ -454,6 +556,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)
@@ -489,6 +592,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
@@ -499,6 +603,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))
@@ -660,7 +765,8 @@ void	read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size)
         if (spoonful>size)
           spoonful = size;
         bufdata = (char *)bufdata0;
-        QFREAD(bufdata, spoonful*tab->bytepix, cat->file, cat->filename);
+        readSpoonful(tab, spoonful, bufdata0, bufdata);
+
         switch(tab->bitpix)
           {
           case BP_BYTE:
@@ -670,6 +776,7 @@ void	read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size)
             break;
 
           case BP_SHORT:
+              if (!tab->isTileCompressed)
             if (bswapflag)
               swapbytes(bufdata, 2, spoonful);
 #pragma ivdep
@@ -678,6 +785,7 @@ void	read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size)
             break;
 
           case BP_LONG:
+              if (!tab->isTileCompressed)
             if (bswapflag)
               swapbytes(bufdata, 4, spoonful);
 #pragma ivdep
@@ -687,6 +795,7 @@ void	read_ibody(tabstruct *tab, FLAGTYPE *ptr, size_t size)
 
 #ifdef HAVE_LONG_LONG_INT
           case BP_LONGLONG:
+              if (!tab->isTileCompressed)
             if (bswapflag)
               swapbytes(bufdata, 8, spoonful);
 #pragma ivdep
@@ -931,6 +1040,10 @@ void	write_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
 #pragma ivdep
             for (i=spoonful; i--;)
               *(bufdata++) = (*(ptr++)-bz)/bs;
+
+
+            // TODO not yet writing CFitsIO from SExtractor. CFITSIO - only perform byte-swap if we are NOT writing a tile-compressed format using cfitsio
+            // if (tab->infptr == NULL) // TODO
             if (bswapflag)
               swapbytes(cbufdata0, 4, spoonful);
             }
@@ -952,7 +1065,8 @@ void	write_body(tabstruct *tab, PIXTYPE *ptr, size_t size)
                                 "read_body()");
             break;
           }
-        QFWRITE(cbufdata0, spoonful*tab->bytepix, cat->file, cat->filename);
+
+        writeSpoonful(tab, spoonful, cbufdata0);
         }
       break;
 

src/fits/fitscat.c

@@ -167,6 +167,64 @@ int	close_cat(catstruct *cat)
   return RETURN_OK;
   }
 
+int open_cfitsio(catstruct *cat)
+{
+#ifdef HAVE_CFITSIO
+    fitsfile *infptr;
+    int status;
+    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);
+    }
+    cat->tab->infptr = infptr;
+    return status;
+#else
+    return 1;
+#endif
+}
+
+int	close_cfitsio(catstruct *cat)
+{
+#ifdef HAVE_CFITSIO
+    int status = 0;
+	if (cat->tab->infptr) {
+
+		; fits_close_file(cat->tab->infptr, &status);
+		if (status != 0) {
+			fits_report_error(stderr, status);
+			printf("ERROR could not close FITS file with cfitsio: %s\n", cat->filename);
+		}
+		else {
+			cat->tab->infptr == NULL;
+		}
+
+	}
+	return status;
+#else
+    return 1;
+#endif
+}
+
+void hdu_seek(tabstruct *tab) {
+#ifdef HAVE_CFITSIO
+    if(!tab->infptr) {
+        tab->infptr = tab->prevtab->infptr;
+    }
+    int status = 0;
+    int hdutype;
+    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, tab->cat->filename);
+    if (tab->tabsize)
+        // IMPORTANT: moving to start of next header using fseek and cfitsio position rather than table size, as done previously
+        fseek(tab->cat->file, tab->infptr->Fptr->headstart[tab->hdunum], SEEK_SET);
+#else
+    if (tab->tabsize)
+        QFSEEK(tab->cat->file, PADTOTAL(tab->tabsize), SEEK_CUR, tab->cat->filename);
+#endif
+}
 
 /****** free_cat ***************************************************************
 PROTO	void free_cat(catstruct **cat, int ncat)
@@ -324,23 +382,24 @@ int	map_cat(catstruct *cat)
   prevtab = NULL;
   QCALLOC(tab, tabstruct, 1);
   tab->cat = cat;
+  cat->tab = tab;
   QFTELL(cat->file, tab->headpos, cat->filename);
+  open_cfitsio(cat);
+
   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);
+    tab->hdunum = ntab + 1;
     if (prevtab)
       {
       tab->prevtab = prevtab;
       prevtab->nexttab = tab;
       }
-    else
-      cat->tab = tab;
     prevtab = tab;
+    hdu_seek(tab);
     QCALLOC(tab, tabstruct, 1);
     tab->cat = cat;
     QFTELL(cat->file, tab->headpos, cat->filename);

src/fits/fitscat.h

@@ -36,6 +36,10 @@
 #include <sys/types.h>
 #endif
 
+#ifdef HAVE_CFITSIO
+#include FITSIO_H
+#endif
+
 #define	MAXCHARS	256	/* max. number of characters */
 #define WARNING_NMAX	1000	/* max. number of recorded warnings */
 
@@ -107,10 +111,11 @@ typedef long long			SLONGLONG;
 typedef union {int l[2];}		SLONGLONG;
 #endif
 
-#if defined(_FILE_OFFSET_BITS) && !defined(OFF_T)
-#define OFF_T	off_t
+// CFITSIO changed OFF_T to OFF_T2 due to clash with cfitsio lib
+#if defined(_FILE_OFFSET_BITS) && !defined(OFF_T2)
+#define OFF_T2	off_t
 #else
-#define OFF_T	long
+#define OFF_T2	long
 #endif
 
 /*------------------------------- constants ---------------------------------*/
@@ -180,8 +185,8 @@ typedef struct structtab
   char		*headbuf;		/* buffer containing the header */
   int		headnblock;		/* number of FITS blocks */
   char		*bodybuf;		/* buffer containing the body */
-  OFF_T		bodypos;		/* position of the body in the file */
-  OFF_T		headpos;		/* position of the head in the file */
+  OFF_T2		bodypos;		/* position of the body in the file */
+  OFF_T2		headpos;		/* position of the head in the file */
   struct structcat *cat;		/* (original) parent catalog */
   struct structtab *prevtab, *nexttab;	/* previous and next tab in chain */
   int		seg;			/* segment position */
@@ -191,6 +196,14 @@ typedef struct structtab
   int		swapflag;		/* mapped to a swap file ? */
   char		swapname[MAXCHARS];	/* name of the swapfile */
   unsigned int	bodysum;		/* Checksum of the FITS body */
+
+#ifdef HAVE_CFITSIO
+  fitsfile *infptr;                     /* a cfitsio pointer to the file */
+#endif
+  int hdunum;                           /* FITS HDU number for this 'table' */
+  int isTileCompressed;                 /* is this a tile compressed image?  */
+  long currentElement;                  /* tracks the current image pixel */
+
   }		tabstruct;
 
 
@@ -271,6 +284,7 @@ extern int	about_cat(catstruct *cat, FILE *stream),
 		add_tab(tabstruct *tab, catstruct *cat, int pos),
 		blank_keys(tabstruct *tab),
 		close_cat(catstruct *cat),
+		close_cfitsio(catstruct *cat),
 		copy_key(tabstruct *tabin, char *keyname, tabstruct *tabout,
 			int pos),
 		copy_tab(catstruct *catin, char *tabname, int seg,

src/fits/fitshead.c

@@ -114,15 +114,37 @@ void	readbasic_head(tabstruct *tab)
 
   filename = (tab->cat? tab->cat->filename : strcpy(name, "internal header"));
 
-  if (fitsread(tab->headbuf, "BITPIX  ", &tab->bitpix, H_INT, T_LONG)
+// CFITSIO
+  char NAXIS_KEYWORD[14];
+  char BITPIX_KEYWORD[14];
+
+  if (fitsread(tab->headbuf, "ZIMAGE  ", str, H_STRING, T_STRING) == RETURN_OK) {
+
+          tab->isTileCompressed = 1;
+#ifndef HAVE_CFITSIO
+          error(EXIT_FAILURE, "*Error*: No CFITSIO support but detected tile compression in ", filename);
+#endif
+          strcpy(NAXIS_KEYWORD, "ZNAXIS  ");
+          strcpy(BITPIX_KEYWORD, "ZBITPIX ");
+  }
+  else {
+          tab->isTileCompressed = 0;
+
+          strcpy(NAXIS_KEYWORD, "NAXIS   ");
+          strcpy(BITPIX_KEYWORD, "BITPIX  ");
+  }
+
+  if (fitsread(tab->headbuf, BITPIX_KEYWORD, &tab->bitpix, H_INT, T_LONG)
 	==RETURN_ERROR)
     error(EXIT_FAILURE, "*Error*: Corrupted FITS header in ", filename);
 
   tab->bytepix = tab->bitpix>0?(tab->bitpix/8):(-tab->bitpix/8);
 
-  if (fitsread(tab->headbuf, "NAXIS   ", &tab->naxis, H_INT, T_LONG)
-	==RETURN_ERROR)
+  if (fitsread(tab->headbuf, NAXIS_KEYWORD, &tab->naxis, H_INT, T_LONG)
+                  ==RETURN_ERROR) {
     error(EXIT_FAILURE, "*Error*: Corrupted FITS header in ", filename);
+  }
+
 
   tabsize = 0;
   if (tab->naxis>0)
@@ -133,7 +155,11 @@ void	readbasic_head(tabstruct *tab)
     tabsize = 1;
     for (i=0; i<tab->naxis && i<999; i++)
       {
-      sprintf(key,"NAXIS%-3d", i+1);
+       // CFITSIO
+       if (tab->isTileCompressed)
+               sprintf(key,"ZNAXIS%-2d", (i+1));
+       else
+               sprintf(key,"NAXIS%-3d", (i+1));
       if (fitsread(tab->headbuf, key, &tab->naxisn[i], H_INT, T_LONG)
 		==RETURN_ERROR)
         error(EXIT_FAILURE, "*Error*: incoherent FITS header in ", filename);
@@ -144,6 +170,10 @@ void	readbasic_head(tabstruct *tab)
 /*random groups parameters (optional)*/
   tab->pcount = 0;
   fitsread(tab->headbuf, "PCOUNT  ", &tab->pcount, H_INT, T_LONG);
+
+  // CFITSIO TODO HACK
+  if (tab->isTileCompressed) tab->pcount = 0;
+
   tab->gcount = 1;
   fitsread(tab->headbuf, "GCOUNT  ", &tab->gcount, H_INT, T_LONG);
 
@@ -284,7 +314,9 @@ int	readbintabparam_head(tabstruct *tab)
         key->htype = H_STRING;
         break;
       default:
-        error(EXIT_FAILURE, "*Error*: Unknown TFORM in ", cat->filename);
+        // CFITSIO TODO dodgy
+    	key->ttype = T_FLOAT;
+        //error(EXIT_FAILURE, "*Error*: Unknown TFORM in ", cat->filename);
       }
 
 /*--handle the special case of multimensional arrays*/

src/fits/fitswrite.c

@@ -365,7 +365,7 @@ void	end_writeobj(catstruct *cat, tabstruct *tab, char *buf)
 
   {
    keystruct	*key;
-   OFF_T	pos;
+   OFF_T2	pos;
    int		k;
 
 /* Make the table parameters reflect its content*/

src/makeit.c

@@ -282,6 +282,7 @@ void	makeit()
       if ((imatab->naxis < 2)
 	|| !strncmp(imatab->xtension, "BINTABLE", 8)
 	|| !strncmp(imatab->xtension, "ASCTABLE", 8))
+    	  if (!imatab->isTileCompressed)
         continue;
       next++;
       }
@@ -321,6 +322,7 @@ void	makeit()
     if (!forcextflag && ((imatab->naxis < 2)
 	|| !strncmp(imatab->xtension, "BINTABLE", 8)
 	|| !strncmp(imatab->xtension, "ASCTABLE", 8)))
+if (!imatab->isTileCompressed)
       continue;
 
     nok++;
@@ -720,6 +722,10 @@ static int	selectext(char *filename)
     if ((bracr=strrchr(bracl+1, ']')))
       *bracr = '\0';
     next = strtol(bracl+1, NULL, 0);
+
+    // VERY BAD HACK to check if this is tile-compressed, if so, add +1 to extension number requested
+    if (strstr(filename, ".fits.fz") != NULL) next = next + 1;
+
     return next;
     }
 

src/readimage.c

@@ -93,7 +93,12 @@ void	*loadstrip(picstruct *field, picstruct *wfield)
       else if (flags & INTERP_FIELD)
         copydata(field, 0, nbpix);
       else
+      {
+        tab->currentElement = 1;
         read_body(tab, data, nbpix);
+
+        //printf("read 1: data[0] = %f %f %f %f %f %f\n", data[0], data[1], data[222], data[777], data[field->width-1], data[nbpix-1]);
+      } 
       if (flags & (WEIGHT_FIELD|RMS_FIELD|BACKRMS_FIELD|VAR_FIELD))
         weight_to_var(field, data, nbpix);
       if ((flags & MEASURE_FIELD) && (check=prefs.check[CHECK_IDENTICAL]))
@@ -198,8 +203,12 @@ void	*loadstrip(picstruct *field, picstruct *wfield)
         backrmsline(field, field->ymax, data);
       else if (flags & INTERP_FIELD)
         copydata(field, field->stripylim*w, w);
-      else
+      else {
+
         read_body(tab, data, w);
+
+        //printf("read 2: data[0] = %f %f %f %f %f %f\n", data[0], data[1], data[222], data[777], data[field->width-1], data[w-1]);
+      }
       if (flags & (WEIGHT_FIELD|RMS_FIELD|BACKRMS_FIELD|VAR_FIELD))
         weight_to_var(field, data, w);
 

src/types.h

@@ -621,7 +621,7 @@ typedef struct pic
   int		interp_xtimeout;	/* interpolation timeout value in x */
   int		interp_ytimeout;	/* interpolation timeout value in y */
   struct pic	*reffield;	       	/* pointer to a reference field */
-  OFF_T		mefpos;			/* Position in a MEF file */
+  OFF_T2		mefpos;			/* Position in a MEF file */
   }	picstruct;