From 57a8d8f12eebd8778d644366af668e881d967a41 Mon Sep 17 00:00:00 2001 From: John Cupitt Date: Tue, 26 Oct 2010 15:12:49 +0000 Subject: [PATCH] add fits read --- ChangeLog | 2 + TODO | 7 + configure.in | 19 +- libvips/format/Makefile.am | 1 + libvips/format/fits.c | 396 +++++++++++++++++++++++++++++++ libvips/format/format.c | 4 + libvips/morphology/im_cntlines.c | 56 ++--- libvips/morphology/im_dilate.c | 22 +- 8 files changed, 477 insertions(+), 30 deletions(-) create mode 100644 libvips/format/fits.c diff --git a/ChangeLog b/ChangeLog index dd42ad18..b4324789 100644 --- a/ChangeLog +++ b/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 diff --git a/TODO b/TODO index fa908e26..4fd3dd3d 100644 --- a/TODO +++ b/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 diff --git a/configure.in b/configure.in index 6d32b9e7..6eabda21 100644 --- a/configure.in +++ b/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 diff --git a/libvips/format/Makefile.am b/libvips/format/Makefile.am index 779b24b4..337dcf61 100644 --- a/libvips/format/Makefile.am +++ b/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@ diff --git a/libvips/format/fits.c b/libvips/format/fits.c new file mode 100644 index 00000000..aee29049 --- /dev/null +++ b/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 +#endif /*HAVE_CONFIG_H*/ +#include + +#ifdef HAVE_CFITSIO + +#include +#include +#include +#include + +#include +#include + +#include + +#ifdef WITH_DMALLOC +#include +#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*/ diff --git a/libvips/format/format.c b/libvips/format/format.c index 8d6a5231..bf401c7d 100644 --- a/libvips/format/format.c +++ b/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(); diff --git a/libvips/morphology/im_cntlines.c b/libvips/morphology/im_cntlines.c index 7b55110a..bd8e4359 100644 --- a/libvips/morphology/im_cntlines.c +++ b/libvips/morphology/im_cntlines.c @@ -1,16 +1,4 @@ -/* @(#) 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 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 +/* count lines * * Copyright: 1990, N. Dessipris. * @@ -20,6 +8,8 @@ * * 19/9/95 JC * - tidied up + * 23/10/10 + * - gtk-doc */ /* @@ -61,6 +51,23 @@ #include #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 ) ) + if( im_incheck( im ) || + im_check_uncoded( "im_cntlines", 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 ); diff --git a/libvips/morphology/im_dilate.c b/libvips/morphology/im_dilate.c index 1f0f1462..3df3c666 100644 --- a/libvips/morphology/im_dilate.c +++ b/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 )