add fits read

ChangeLog

@@ -37,6 +37,8 @@
 - fix gtk-doc warnings
 - small mask load/save improvements
 - mask gtk-doc done
+- add cfitsio dependancy
+- add FITS reader
 
 12/5/10 started 7.22.2
 - the conditional image of ifthenelse can be any format, a (!=0) is added if

TODO

@@ -1,3 +1,10 @@
+- test fits reader more ... colour?
+
+- we get lots of valgrind errors with fits reading
+
+- why do we read_free so often during fits read?
+
+- are we triggering im_fits2vips() in each thread?? argh, try tracing
 
 
 - maybe im_draw_smudge() is too slow :-( also, we had a sanity failure with

configure.in

@@ -424,6 +424,20 @@ if test x"$with_matio" != "xno"; then
     ])
 fi
 
+# cfitsio
+AC_ARG_WITH([cfitsio], 
+  AS_HELP_STRING([--without-cfitsio], [build without cfitsio (default: test)]))
+
+if test x"$with_cfitsio" != "xno"; then
+  PKG_CHECK_MODULES(CFITSIO, cfitsio,
+    [AC_DEFINE(HAVE_CFITSIO,1,[define if you have cfitsio installed.])
+     with_cfitsio=yes
+     PACKAGES_USED="$PACKAGES_USED cfitsio"],
+    [AC_MSG_WARN([cfitsio not found; disabling cfitsio support])
+     with_cfitsio=no
+    ])
+fi
+
 # pangoft2
 AC_ARG_WITH([pangoft2], 
   AS_HELP_STRING([--without-pangoft2], 
@@ -556,14 +570,14 @@ fi
 # Gather all up for VIPS_CFLAGS, VIPS_INCLUDES, VIPS_LIBS and VIPS_CXX_LIBS
 # sort includes to get longer, more specific dirs first
 # helps, for example, selecting graphicsmagick over imagemagick
-VIPS_CFLAGS=`for i in $VIPS_CFLAGS $GTHREAD_CFLAGS $REQUIRED_CFLAGS $PANGOFT2_CFLAGS $FFTW3_CFLAGS $MAGICK_CFLAGS $PNG_CFLAGS $EXIF_CFLAGS $MATIO_CFLAGS $OPENEXR_CFLAGS  
+VIPS_CFLAGS=`for i in $VIPS_CFLAGS $GTHREAD_CFLAGS $REQUIRED_CFLAGS $PANGOFT2_CFLAGS $FFTW3_CFLAGS $MAGICK_CFLAGS $PNG_CFLAGS $EXIF_CFLAGS $MATIO_CFLAGS $CFITSIO_CFLAGS $OPENEXR_CFLAGS  
 do 
   echo $i 
 done | sort -ru`
 VIPS_CFLAGS=`echo $VIPS_CFLAGS`
 VIPS_CFLAGS="$VIPS_DEBUG_FLAGS $VIPS_CFLAGS"
 VIPS_INCLUDES="$PNG_INCLUDES $TIFF_INCLUDES $ZIP_INCLUDES $JPEG_INCLUDES $FFTW_INCLUDES $LCMS_INCLUDES" 
-VIPS_LIBS="$MAGICK_LIBS $PNG_LIBS $TIFF_LIBS $ZIP_LIBS $JPEG_LIBS $GTHREAD_LIBS $REQUIRED_LIBS $PANGOFT2_LIBS $FFTW3_LIBS $FFTW_LIBS $LCMS_LIBS $OPENEXR_LIBS $MATIO_LIBS $EXIF_LIBS -lm"
+VIPS_LIBS="$MAGICK_LIBS $PNG_LIBS $TIFF_LIBS $ZIP_LIBS $JPEG_LIBS $GTHREAD_LIBS $REQUIRED_LIBS $PANGOFT2_LIBS $FFTW3_LIBS $FFTW_LIBS $LCMS_LIBS $OPENEXR_LIBS $CFITSIO_LIBS $MATIO_LIBS $EXIF_LIBS -lm"
 # need -lstdc++ for (eg.) the C++ format loaders
 VIPS_CXX_LIBS="-lstdc++"
 
@@ -649,6 +663,7 @@ file import with libMagick: 		$with_magick
 ICC profile support with lcms: 		$with_lcms (version $with_lcms_ver)
 file import with OpenEXR: 		$with_OpenEXR
 file import with matio: 		$with_matio
+file import with cfitsio: 		$with_cfitsio
 text rendering with pangoft2: 		$with_pangoft2
 file import/export with libpng: 	$with_png
 file import/export with libtiff:	$with_tiff

libvips/format/Makefile.am

@@ -21,6 +21,7 @@ libformat_la_SOURCES = \
 	im_vips2tiff.c \
 	im_vips2raw.c \
 	matlab.c \
+	fits.c \
 	radiance.c 
 
 INCLUDES = -I${top_srcdir}/libvips/include @VIPS_CFLAGS@ @VIPS_INCLUDES@

libvips/format/fits.c

@@ -0,0 +1,396 @@
+/* Read FITS files with cfitsio
+ *
+ * 26/10/10
+ *	- from matlab.c
+ */
+
+/*
+
+    This file is part of VIPS.
+    
+    VIPS is free software; you can redistribute it and/or modify
+    it under the terms of the GNU Lesser General Public License as published by
+    the Free Software Foundation; either version 2 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 Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+ */
+
+/*
+
+    These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
+
+ */
+
+/*
+#define DEBUG
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif /*HAVE_CONFIG_H*/
+#include <vips/intl.h>
+
+#ifdef HAVE_CFITSIO
+
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <vips/vips.h>
+#include <vips/internal.h>
+
+#include <fitsio.h>
+
+#ifdef WITH_DMALLOC
+#include <dmalloc.h>
+#endif /*WITH_DMALLOC*/
+
+/* What we track during a cfitsio-file read.
+ */
+typedef struct {
+	char *filename;
+	IMAGE *out;
+
+	fitsfile *fptr;
+	int datatype;
+} Read;
+
+static void
+read_error( int status )
+{
+	char buf[80];
+
+	fits_get_errstatus( status, buf );
+	im_error( "fits", "%s", buf );
+}
+
+static void
+read_destroy( Read *read )
+{
+	IM_FREE( read->filename );
+	if( read->fptr ) {
+		int status;
+
+		if( fits_close_file( read->fptr, &status ) ) 
+			read_error( status );
+
+		read->fptr = NULL;
+	}
+
+	im_free( read );
+}
+
+static Read *
+read_new( const char *filename, IMAGE *out )
+{
+	Read *read;
+	int status;
+
+	if( !(read = IM_NEW( NULL, Read )) )
+		return( NULL );
+
+	read->filename = im_strdup( NULL, filename );
+	read->out = out;
+	read->fptr = NULL;
+
+	status = 0;
+
+	if( fits_open_file( &read->fptr, filename, READONLY, &status ) ) {
+		im_error( "fits", _( "unable to open \"%s\"" ), filename );
+		read_error( status );
+		read_destroy( read );
+		return( NULL );
+	}
+
+	return( read );
+}
+
+/* fits image types -> VIPS band formats. VIPS doesn't have 64-bit int, so no
+ * entry for LONGLONG_IMG (64).
+ */
+static int fits2vips_formats[][3] = {
+	{ BYTE_IMG, IM_BANDFMT_UCHAR, TBYTE },
+	{ SHORT_IMG,  IM_BANDFMT_USHORT, TUSHORT },
+	{ LONG_IMG,  IM_BANDFMT_UINT, TUINT },
+	{ FLOAT_IMG,  IM_BANDFMT_FLOAT, TFLOAT },
+	{ DOUBLE_IMG, IM_BANDFMT_DOUBLE, TDOUBLE }
+};
+
+static int
+fits2vips_get_header( Read *read )
+{
+	int status;
+	int bitpix;
+	int naxis;
+	long long int naxes[3];
+
+	int width, height, bands, format, type;
+	int keysexist;
+	int morekeys;
+	int i;
+
+	status = 0;
+
+	if( fits_get_img_paramll( read->fptr, 
+		3, &bitpix, &naxis, naxes, &status ) ) {
+		read_error( status );
+		return( -1 );
+	}
+	width = 1;
+	height = 1;
+	bands = 1;
+	switch( naxis ) {
+	case 3:
+		bands = naxes[2];
+
+	case 2:
+		height = naxes[1];
+
+	case 1:
+		width = naxes[0];
+		break;
+
+	default:
+		im_error( "fits", "bad number of axis %d", naxis );
+		return( -1 );
+	}
+
+	if( bands == 1 )
+		type = IM_TYPE_B_W;
+	else if( bands == 3 )
+		type = IM_TYPE_RGB;
+	else
+		type = IM_TYPE_MULTIBAND;
+
+	/* Get image format. We want the 'raw' format of the image, our caller
+	 * can convert using the meta info if they want.
+	 */
+	for( i = 0; i < IM_NUMBER( fits2vips_formats ); i++ )
+		if( fits2vips_formats[i][0] == bitpix )
+			break;
+	if( i == IM_NUMBER( fits2vips_formats ) ) {
+		im_error( "im_fits2vips", _( "unsupported bitpix %d\n" ),
+			bitpix );
+		return( -1 );
+	}
+	format = fits2vips_formats[i][1];
+	read->datatype = fits2vips_formats[i][2];
+
+	im_initdesc( read->out,
+		 width, height, bands,
+		 im_bits_of_fmt( format ), format,
+		 IM_CODING_NONE, type, 1.0, 1.0, 0, 0 );
+
+	/* Read all keys into meta.
+	 */
+	if( fits_get_hdrspace( read->fptr, &keysexist, &morekeys, &status ) ) {
+		read_error( status );
+		return( -1 );
+	}
+
+	for( i = 0; i < keysexist; i++ ) {
+		char key[81];
+		char value[81];
+		char comment[81];
+		char vipsname[100];
+
+		if( fits_read_keyn( read->fptr, i + 1, 
+			key, value, comment, &status ) ) {
+			read_error( status );
+			return( -1 );
+		}
+
+#ifdef DEBUG
+		printf( "fits: seen:\n" );
+		printf( " key == %s\n", key );
+		printf( " value == %s\n", value );
+		printf( " comment == %s\n", comment );
+#endif /*DEBUG*/
+
+		im_snprintf( vipsname, 100, "fits-%s", key );
+		if( im_meta_set_string( read->out, vipsname, value ) ) 
+			return( -1 );
+		im_snprintf( vipsname, 100, "fits-%s-comment", key );
+		if( im_meta_set_string( read->out, vipsname, comment ) ) 
+			return( -1 );
+	}
+
+	return( 0 );
+}
+
+static int
+fits2vips_header( const char *filename, IMAGE *out )
+{
+	Read *read;
+
+#ifdef DEBUG
+	printf( "fits2vips_header: reading \"%s\"\n", filename );
+#endif /*DEBUG*/
+
+	if( !(read = read_new( filename, out )) ) 
+		return( -1 );
+	if( fits2vips_get_header( read ) ) {
+		read_destroy( read );
+		return( -1 );
+	}
+	read_destroy( read );
+
+	return( 0 );
+}
+
+static int
+fits2vips_get_data( Read *read )
+{
+	IMAGE *im = read->out;
+	const int es = IM_IMAGE_SIZEOF_ELEMENT( im );
+
+	PEL *line_buffer;
+	PEL *band_buffer;
+	PEL *p, *q;
+	int x, y, b, z;
+	int status;
+
+	status = 0;
+
+	if( !(line_buffer = IM_ARRAY( im, IM_IMAGE_SIZEOF_LINE( im ), PEL )) ||
+		!(band_buffer = IM_ARRAY( im, es * im->Xsize, PEL )) ||
+		im_outcheck( im ) ||
+		im_setupout( im ) )
+		return( -1 );
+
+	for( y = 0; y < im->Ysize; y++ ) {
+		long int fpixel[3];
+
+		/* Start of scanline. We have to read top-to-bottom.
+		 */
+		fpixel[0] = 1;
+		fpixel[1] = im->Ysize - y;
+		fpixel[2] = 1;
+
+		for( b = 0; b < im->Bands; b++ ) {
+			fpixel[2] = b + 1;
+
+			/* Read one band of one scanline, then scatter-write
+			 * into the line buffer.
+			 */
+			if( fits_read_pix( read->fptr, 
+				read->datatype, fpixel, im->Xsize,
+				NULL, band_buffer, NULL, &status ) ) {
+				read_error( status );
+				return( -1 );
+			}
+
+			p = band_buffer;
+			q = line_buffer + b * es;
+			for( x = 0; x < im->Xsize; x++ ) {
+				for( z = 0; z < es; z++ )
+					q[z] = p[z];
+
+				p += es;
+				q += im->Bands * es;
+			}
+		}
+
+		if( im_writeline( y, im, line_buffer ) )
+			return( -1 );
+	}
+
+	return( 0 );
+}
+
+/**
+ * im_fits2vips:
+ * @filename: file to load
+ * @out: image to write to
+ *
+ * Read a FITS image file into a VIPS image. 
+ *
+ *
+ * See also: #VipsFormat.
+ *
+ * Returns: 0 on success, -1 on error.
+ */
+int
+im_fits2vips( const char *filename, IMAGE *out )
+{
+	Read *read;
+
+#ifdef DEBUG
+	printf( "fits2vips: reading \"%s\"\n", filename );
+#endif /*DEBUG*/
+
+	if( !(read = read_new( filename, out )) ) 
+		return( -1 );
+	if( fits2vips_get_header( read ) ||
+		fits2vips_get_data( read ) ) {
+		read_destroy( read );
+		return( -1 );
+	}
+
+	read_destroy( read );
+
+	return( 0 );
+}
+
+static int
+isfits( const char *filename )
+{
+	fitsfile *fptr;
+	int status;
+
+	status = 0;
+
+	if( fits_open_image( &fptr, filename, READONLY, &status ) ) {
+#ifdef DEBUG
+		printf( "isfits: error reading \"%s\"\n", filename );
+		read_error( status );
+#endif /*DEBUG*/
+
+		return( 0 );
+	}
+	fits_close_file( fptr, &status );
+
+	return( 1 );
+}
+
+static const char *fits_suffs[] = { ".fits", NULL };
+
+/* fits format adds no new members.
+ */
+typedef VipsFormat VipsFormatFits;
+typedef VipsFormatClass VipsFormatFitsClass;
+
+static void
+vips_format_fits_class_init( VipsFormatFitsClass *class )
+{
+	VipsObjectClass *object_class = (VipsObjectClass *) class;
+	VipsFormatClass *format_class = (VipsFormatClass *) class;
+
+	object_class->nickname = "fits";
+	object_class->description = _( "FITS" );
+
+	format_class->is_a = isfits;
+	format_class->header = fits2vips_header;
+	format_class->load = im_fits2vips;
+	format_class->save = NULL;
+	format_class->suffs = fits_suffs;
+}
+
+static void
+vips_format_fits_init( VipsFormatFits *object )
+{
+}
+
+G_DEFINE_TYPE( VipsFormatFits, vips_format_fits, VIPS_TYPE_FORMAT );
+
+#endif /*HAVE_CFITSIO*/

libvips/format/format.c

@@ -428,6 +428,10 @@ im__format_init( void )
 	extern GType vips_format_mat_get_type();
 	vips_format_mat_get_type();
 #endif /*HAVE_MATIO*/
+#ifdef HAVE_CFITSIO
+	extern GType vips_format_fits_get_type();
+	vips_format_fits_get_type();
+#endif /*HAVE_CFITSIO*/
 	vips_format_rad_get_type();
 #ifdef HAVE_MAGICK
 	extern GType vips_format_magick_get_type();

libvips/morphology/im_cntlines.c

@@ -1,16 +1,4 @@
-/* @(#)  Function which calculates the number of transitions
+/* count lines 
- * @(#) between black and white for the horizontal or the vertical
- * @(#) direction of an image.  black<128 , white>=128
- * @(#) The function calculates the number of transitions for all
- * @(#) Xsize or Ysize and returns the mean of the result
- * @(#) Input should be binary one channel
- * @(#)
- * @(#) int im_cntlines(im, nolines, flag)
- * @(#) IMAGE *im;
- * @(#) float *nolines;
- * @(#) int flag;	0 to cnt horizontal and 1 to count vertical lines
- * @(#)
- * @(#) Returns 0 on success and -1 on error
  *
  * Copyright: 1990, N. Dessipris.
  *
@@ -20,6 +8,8 @@
  *
  * 19/9/95 JC
  *	- tidied up
+ * 23/10/10
+ * 	- gtk-doc
  */
 
 /*
@@ -61,6 +51,23 @@
 #include <dmalloc.h>
 #endif /*WITH_DMALLOC*/
 
+/**
+ * im_cntlines:
+ * @im: input #IMAGE
+ * @nolines: output average number of lines
+ * @flag: 0 horizontal, 1 vertical
+ *
+ * Function which calculates the number of transitions
+ * between black and white for the horizontal or the vertical
+ * direction of an image.  black<128 , white>=128
+ * The function calculates the number of transitions for all
+ * Xsize or Ysize and returns the mean of the result
+ * Input should be one band, 8-bit.
+ *
+ * See also: im_erode().
+ *
+ * Returns: 0 on success, -1 on error.
+ */
 int 
 im_cntlines( IMAGE *im, double *nolines, int flag )
 {
@@ -68,16 +75,11 @@ im_cntlines( IMAGE *im, double *nolines, int flag )
 	PEL *line;
 	int cnt;
 
-	/* Is im valid? 
+	if( im_incheck( im ) ||
-	 */
+		im_check_uncoded( "im_cntlines", im ) ||
-	if( im_incheck( im ) )
+		im_check_mono( "im_cntlines", im ) ||
+		im_check_format( "im_cntlines", im, IM_BANDFMT_UCHAR ) )
 		return( -1 );
-	if( im->Coding != IM_CODING_NONE || im->BandFmt != IM_BANDFMT_UCHAR ||
-		im->Bands != 1 ) {
-		im_error( "im_cntlines", "%s", 
-			_( "1-band uchar uncoded only" ) );
-		return( -1 ); 
-	}
 	if( flag != 0 && flag != 1 ) {
 		im_error( "im_cntlines", "%s", 
 			_( "flag should be 0 (horizontal) or 1 (vertical)" ) );
@@ -92,9 +94,9 @@ im_cntlines( IMAGE *im, double *nolines, int flag )
 			PEL *p = line;
 			
 			for( x = 0; x < im->Xsize - 1; x++ ) {
-				if( p[0] < (PEL)128 && p[1] >= (PEL)128 )
+				if( p[0] < 128 && p[1] >= 128 )
 					cnt++;
-				else if( p[0] >= (PEL)128 && p[1] < (PEL)128 )
+				else if( p[0] >= 128 && p[1] < 128 )
 					cnt++;
 				
 				p++;
@@ -113,9 +115,9 @@ im_cntlines( IMAGE *im, double *nolines, int flag )
 			PEL *p2 = line + im->Xsize;
 			
 			for( x = 0; x < im->Xsize; x++ ) {
-				if( *p1 < (PEL)128 && *p2 >= (PEL)128 )
+				if( *p1 < 128 && *p2 >= 128 )
 					cnt++;
-				else if( *p1 >= (PEL)128 && *p2 < (PEL)128 )
+				else if( *p1 >= 128 && *p2 < 128 )
 					cnt++;
 
 				p1++;
@@ -125,7 +127,7 @@ im_cntlines( IMAGE *im, double *nolines, int flag )
 			line += im->Xsize;
 		}
 
-		*nolines = (float)cnt / (2.0 * im->Xsize);
+		*nolines = (float) cnt / (2.0 * im->Xsize);
 	}
 
 	return( 0 );

libvips/morphology/im_dilate.c

@@ -308,7 +308,27 @@ im_dilate_raw( IMAGE *in, IMAGE *out, INTMASK *m )
 	return( 0 );
 }
 
-/* The above, with a border to make out the same size as in.
+
+/**
+ * im_dilate:
+ * @in: input image
+ * @out: output image
+ * @m: dilate mask
+ *
+ * Dilate an image with a mask.
+ * The mask coefficients are either 255 (object) or 0 (black) or 128 (don't
+ * care).
+ * Input image are binary images with either 0 or 255 values, one channel 
+ * only. The program dilates a white object on a black background.
+ * The center of the mask is at location (m->xsize/2, m->ysize/2)
+ * integer division. The mask is expected to have an odd width and
+ * height.
+
+sets pixels in the output if
+ *
+ * See also: im_erode().
+ *
+ * Returns: 0 on success, -1 on error
  */
 int 
 im_dilate( IMAGE *in, IMAGE *out, INTMASK *m )