From f32c3f837e61ea3b2bcfade210a03fa47e95a5f3 Mon Sep 17 00:00:00 2001 From: Fons Rademakers <Fons.Rademakers@cern.ch> Date: Tue, 27 Jul 2010 12:56:50 +0000 Subject: [PATCH] From Claudi Astres: New module graf2d/fitsio to read (atronomical) images and data from FITS files using the cfitsio library (see http://heasarc.gsfc.nasa.gov/fitsio/). To get cfitsio on Ubuntu do "apt-get install cfitsio-dev" or on Mac install it from Fink or get the source from the NASA site. The module comes with some tutorials in the $ROOTSYS/tutorials/fitsio directory. git-svn-id: http://root.cern.ch/svn/root/trunk@34610 27541ba8-7e3a-0410-8455-c3a389f83636 --- Makefile | 9 +- config/Makefile.depend | 13 +- config/Makefile.in | 21 +- configure | 50 ++- graf2d/fitsio/Module.mk | 70 ++++ graf2d/fitsio/inc/LinkDef.h | 19 + graf2d/fitsio/inc/TFITS.h | 98 +++++ graf2d/fitsio/src/TFITS.cxx | 723 ++++++++++++++++++++++++++++++++++++ 8 files changed, 984 insertions(+), 19 deletions(-) create mode 100644 graf2d/fitsio/Module.mk create mode 100644 graf2d/fitsio/inc/LinkDef.h create mode 100644 graf2d/fitsio/inc/TFITS.h create mode 100644 graf2d/fitsio/src/TFITS.cxx diff --git a/Makefile b/Makefile index 5e4b6ff4a8b..576c010d1dc 100644 --- a/Makefile +++ b/Makefile @@ -117,7 +117,7 @@ endif ifeq ($(BUILDGLEW),yes) MODULES += graf3d/glew endif -MODULES += graf3d/gl graf3d/eve graf3d/gviz3d +MODULES += graf3d/gl graf3d/eve graf3d/gviz3d endif ifeq ($(BUILDMYSQL),yes) MODULES += sql/mysql @@ -173,6 +173,9 @@ endif ifeq ($(BUILDFFTW3),yes) MODULES += math/fftw endif +ifeq ($(BUILDFITSIO),yes) +MODULES += graf2d/fitsio +endif ifeq ($(BUILDGVIZ),yes) MODULES += graf2d/gviz endif @@ -211,7 +214,7 @@ ifeq ($(BUILDCINTEX),yes) MODULES += cint/cintex endif ifeq ($(BUILDCLING),yes) -# to be added to the unconditional MODULES list below once cling is in trunk +# to be added to the unconditional MODULES list above once cling is in trunk MODULES += cint/cling endif ifeq ($(BUILDROOFIT),yes) @@ -286,7 +289,7 @@ MODULES += core/unix core/winnt core/editline graf2d/x11 graf2d/x11ttf \ graf2d/qt gui/qtroot gui/qtgsi net/xrootd net/netx net/alien \ proof/proofd proof/proofx proof/clarens proof/peac proof/pq2 \ sql/oracle io/xmlparser math/mathmore cint/reflex cint/cintex \ - tmva io/hdfs \ + tmva io/hdfs graf2d/fitsio \ roofit/roofitcore roofit/roofit roofit/roostats \ math/minuit2 net/monalisa math/fftw sql/odbc math/unuran \ geom/gdml graf3d/eve net/glite misc/memstat \ diff --git a/config/Makefile.depend b/config/Makefile.depend index c8f8aaadb80..b4d31f332c9 100644 --- a/config/Makefile.depend +++ b/config/Makefile.depend @@ -57,7 +57,7 @@ TMVALIBDEPM = $(IOLIB) $(HISTLIB) $(MATRIXLIB) $(TREELIB) \ $(MINUITLIB) $(MATHCORELIB) $(XMLLIB) SPLOTLIBDEPM = $(MATRIXLIB) $(HISTLIB) $(TREELIB) $(TREEPLAYERLIB) \ $(G3DLIB) $(GRAFLIB) $(MATHCORELIB) -QUADPLIBDEPM = $(MATRIXLIB) +QUADPLIBDEPM = $(MATRIXLIB) GLLIBDEPM = $(G3DLIB) $(GUILIB) $(GRAFLIB) $(HISTLIB) $(GEDLIB) \ $(MATHCORELIB) $(TREELIB) $(TREEPLAYERLIB) HBOOKLIBDEPM = $(HISTLIB) $(MATRIXLIB) $(TREELIB) $(GRAFLIB) \ @@ -76,7 +76,7 @@ ALIENLIBDEPM = $(XMLLIB) $(NETXLIB) $(TREELIB) $(PROOFPLAYERLIB) \ ROOFITCORELIBDEPM = $(HISTLIB) $(GRAFLIB) $(MATRIXLIB) $(TREELIB) \ $(MINUITLIB) $(IOLIB) $(MATHCORELIB) $(FOAMLIB) ROOFITLIBDEPM = $(ROOFITCORELIB) $(TREELIB) $(IOLIB) $(MATRIXLIB) \ - $(MATHCORELIB) + $(MATHCORELIB) ROOSTATSLIBDEPM = $(ROOFITLIB) $(ROOFITCORELIB) $(TREELIB) $(IOLIB) \ $(MATRIXLIB) $(HISTLIB) $(MATHCORELIB) $(GRAFLIB) \ $(GPADLIB) @@ -115,7 +115,7 @@ GVIZ3DLIBDEPM = $(GUILIB) $(GEDLIB) $(GPADLIB) $(G3DLIB) $(GRAFLIB) \ MEMSTATLIBDEPM = $(TREELIB) $(GPADLIB) $(GRAFLIB) MEMSTATGUILIBDEPM = $(MEMSTATLIB) $(TREELIB) $(GUILIB) RECLIBDEPM = $(IOLIB) $(TREELIB) $(GUILIB) $(THREADLIB) - +FITSIOLIBDEPM = $(HISTLIB) $(GRAFLIB) $(MATRIXLIB) ifeq ($(EXPLICITLINK),yes) @@ -208,6 +208,7 @@ GVIZ3DLIBDEP = $(GVIZ3DLIBDEPM) MEMSTATLIBDEP = $(MEMSTATLIBDEPM) MEMSTATGUILIBDEP = $(MEMSTATGUILIBDEPM) RECLIBDEP = $(RECLIBDEPM) +FITSIOLIBDEP = $(FITSIOLIBDEPM) ifeq ($(PLATFORM),win32) @@ -222,7 +223,7 @@ G3DLIBEXTRA = lib/libGraf.lib lib/libHist.lib lib/libGpad.lib \ WIN32GDKLIBEXTRA = lib/libGraf.lib ASIMAGELIBEXTRA = lib/libGraf.lib lib/libMathCore.lib ASIMAGEGUILIBEXTRA = lib/libGraf.lib lib/libHist.lib lib/libGui.lib \ - lib/libASImage.lib lib/libRIO.lib + lib/libASImage.lib lib/libRIO.lib ASIMAGEGSLIBEXTRA = lib/libGraf.lib lib/libASImage.lib GEDLIBEXTRA = lib/libHist.lib lib/libGpad.lib lib/libGraf.lib \ lib/libTree.lib lib/libTreePlayer.lib lib/libGui.lib @@ -284,7 +285,7 @@ TMVALIBEXTRA = lib/libRIO.lib lib/libHist.lib lib/libMatrix.lib \ SPLOTLIBEXTRA = lib/libMatrix.lib lib/libHist.lib lib/libTree.lib \ lib/libTreePlayer.lib lib/libGraf3d.lib \ lib/libGraf.lib lib/libMathCore.lib -QUADPLIBEXTRA = lib/libMatrix.lib +QUADPLIBEXTRA = lib/libMatrix.lib GLLIBEXTRA = lib/libGraf3d.lib lib/libGui.lib lib/libGraf.lib \ lib/libHist.lib lib/libGed.lib lib/libMathCore.lib \ lib/libTree.lib lib/libTreePlayer.lib lib/libRIO.lib @@ -354,6 +355,7 @@ MEMSTATLIBEXTRA = lib/libTree.lib lib/libGpad.lib lib/libGraf.lib MEMSTATGUILIBEXTRA = lib/libMemStat.lib lib/libTree.lib lib/libGui.lib RECLIBEXTRA = lib/libGui.lib lib/libRIO.lib lib/libTree.lib \ lib/libThread.lib +FITSIOLIBEXTRA = lib/libHist.lib lib/libGraf.lib lib/libMatrix.lib else @@ -462,6 +464,7 @@ GVIZ3DLIBEXTRA = -Llib -lGui -lGed -lGpad -lGraf3d -lGraf -lGeom -lRGL MEMSTATLIBEXTRA = -Llib -lTree -lGpad -lGraf MEMSTATGUILIBEXTRA = -Llib -lMemStat -lTree -lGui RECLIBEXTRA = -Llib -lGui -lRIO -lTree -lThread +FITSIOLIBEXTRA = -Llib -lHist -lGraf -lMatrix endif diff --git a/config/Makefile.in b/config/Makefile.in index de62aabab28..ebd29359d6c 100644 --- a/config/Makefile.in +++ b/config/Makefile.in @@ -96,7 +96,7 @@ BUILDQTGSI := @buildqtgsi@ QTLIBDIR := @qtlibdir@ QTLIB := @qtlib@ QTINCDIR := $(filter-out /usr/include, @qtincdir@) -QTVERS := @qtvers@ +QTVERS := @qtvers@ QTMOCEXE := @qtmocexe@ BUILDRFIO := @buildrfio@ @@ -154,13 +154,13 @@ CHIRPLIBDIR := @chirplibdir@ CHIRPCLILIB := @chirplib@ CHIRPINCDIR := $(filter-out /usr/include, @chirpincdir@) -BUILDHDFS := @buildhdfs@ -HDFSLIBDIR := @hdfslibdir@ -HDFSCLILIB := @hdfslib@ -HDFSINCDIR := $(filter-out /usr/include, @hdfsincdir@) -JNIINCDIR := $(filter-out /usr/include, @jniincdir@) $(filter-out /usr/include, @jniincdir@/linux) -JVMCLILIB := @jvmlib@ -JVMLIBDIR := @jvmlibdir@ +BUILDHDFS := @buildhdfs@ +HDFSLIBDIR := @hdfslibdir@ +HDFSCLILIB := @hdfslib@ +HDFSINCDIR := $(filter-out /usr/include, @hdfsincdir@) +JNIINCDIR := $(filter-out /usr/include, @jniincdir@) $(filter-out /usr/include, @jniincdir@/linux) +JVMCLILIB := @jvmlib@ +JVMLIBDIR := @jvmlibdir@ BUILDALIEN := @buildalien@ ALIENLIBDIR := @alienlibdir@ @@ -194,6 +194,11 @@ FFTW3LIBDIR := @fftw3libdir@ FFTW3LIB := @fftw3lib@ FFTW3INCDIR := $(filter-out /usr/include, @fftw3incdir@) +BUILDFITSIO := @buildfitsio@ +CFITSIOLIBDIR := @cfitsiolibdir@ +CFITSIOLIB := @cfitsiolib@ +CFITSIOINCDIR := $(filter-out /usr/include, @cfitsioincdir@) + BUILDGVIZ := @buildgviz@ GRAPHVIZLIBDIR := @gvizlibdir@ GRAPHVIZLIB := @gvizlib@ diff --git a/configure b/configure index 73a2d0d4d8e..f70c3990958 100755 --- a/configure +++ b/configure @@ -55,6 +55,7 @@ options=" \ enable_exceptions \ enable_explicitlink \ enable_fftw3 \ + enable_fitsio \ enable_gviz \ enable_gdml \ enable_genvector \ @@ -177,6 +178,7 @@ HDFS \ PYTHIA6 \ PYTHIA8 \ FFTW3 \ +CFITSIO \ GVIZ \ PYTHONDIR \ DCACHE \ @@ -1180,12 +1182,13 @@ enable/disable options, prefix with either --enable- or --disable- exceptions Turn on compiler exception handling capability explicitlink Explicitly link with all dependent libraries fftw3 Fast Fourier Transform support, requires libfftw3 + fitsio Read images and data from FITS files, requires cfitsio gviz Graphs visualization support, requires graphviz gdml GDML writer and reader gfal GFAL support, requires libgfal globus Globus authentication support, requires Globus toolkit glite gLite support, requires libglite-api-wrapper v.3 from GSI (https://subversion.gsi.de/trac/dgrid/wiki) - gsl-shared enable linking against shared libraries for GSL (default no) + gsl-shared Enable linking against shared libraries for GSL (default no) hdfs HDFS support; requires libhdfs from HDFS >= 0.19.1 krb5 Kerberos5 support, requires Kerberos libs ldap LDAP support, requires (Open)LDAP libs @@ -1256,6 +1259,8 @@ with options, prefix with --with-, enables corresponding support ftgl-libdir FTGL support, location of libftgl fftw3-incdir FFTW3 support, location of fftw3.h fftw3-libdir FFTW3 support, location of libfftw3 (libfftw3-3 for windows) + cfitsio-incdir FITS support, location of fitsio2.h + cfitsio-libdir FITS support, location of libcfitsio gviz-incdir Graphviz support, location of gvc.h gviz-libdir Graphviz support, location of libgvplugin_core gfal-incdir GFAL support, location of gfal_api.h @@ -1515,8 +1520,10 @@ if test $# -gt 0 ; then --with-ftgl-libdir=*) ftgllibdir=$optarg ; enable_builtin_ftgl=no;; --with-fftw3-incdir=*) fftw3incdir=$optarg ; enable_fftw3="yes" ;; --with-fftw3-libdir=*) fftw3libdir=$optarg ; enable_fftw3="yes" ;; - --with-gviz-incdir=*) gvizincdir=$optarg ; enable_gviz="yes" ;; - --with-gviz-libdir=*) gvizlibdir=$optarg ; enable_gviz="yes" ;; + --with-cfitsio-incdir=*) cfitsioincdir=$optarg ; enable_fitsio="yes" ;; + --with-cfitsio-libdir=*) cfitsiolibdir=$optarg ; enable_fitsio="yes" ;; + --with-gviz-incdir=*) gvizincdir=$optarg ; enable_gviz="yes" ;; + --with-gviz-libdir=*) gvizlibdir=$optarg ; enable_gviz="yes" ;; --with-gfal-incdir=*) gfalincdir=$optarg ; enable_gfal="yes" ;; --with-gfal-libdir=*) gfallibdir=$optarg ; enable_gfal="yes" ;; --with-glew-incdir=*) glewincdir=$optarg ; enable_builtin_glew=no;; @@ -2944,6 +2951,38 @@ fi check_explicit "$enable_fftw3" "$enable_fftw3_explicit" \ "Explicitly required FFT in the West dependencies not fulfilled" +###################################################################### +# +### echo %%% FITSIO Support - Third party libraries +# +# (See http://heasarc.gsfc.nasa.gov/fitsio/) +# +# If the user has set the flags "--disable-fitsio", we don't check for +# FITSIO at all. +# +if test ! "x$enable_fitsio" = "xno"; then + # Check for cfitsio include and library + check_header "fitsio2.h" "$cfitsioincdir" $CFITSIO $CFITSIO/include \ + /usr/local/include /usr/include /opt/cfitsio/include \ + $finkdir/include + cfitsioinc=$found_hdr + cfitsioincdir=$found_dir + + # At this time, libcfitsio.a should always be prefered over .so, + # to avoid forcing users to install cfitsio. + check_library "libcfitsio" "no" "$cfitsiolibdir" \ + $CFITSIO $CFITSIO/lib $CFITSIO/.libs /usr/local/lib /usr/lib \ + /opt/cfitsio/lib $finkdir/lib + cfitsiolib=$found_lib + cfitsiolibdir=$found_dir + + if test "x$cfitsioincdir" = "x" || test "x$cfitsiolib" = "x"; then + enable_fitsio="no" + fi +fi +check_explicit "$enable_fitsio" "$enable_fitsio_explicit" \ + "Explicitly required cfitsio dependencies not fulfilled" + ###################################################################### # ### echo %%% graphviz Support - Third party libraries @@ -5260,6 +5299,7 @@ if test "x$show_pkglist" = "xyes" ; then test "x$enable_hdfs" = "xyes" && pl="$pl root-plugin-io-hdfs" test "x$enable_rfio" = "xyes" && pl="$pl root-plugin-io-rfio" test "x$enable_fftw3" = "xyes" && pl="$pl root-plugin-math-fftw3" + test "x$enable_fitsio" = "xyes" && pl="$pl root-plugin-graf2d-fitsio" test "x$enable_gviz" = "xyes" && pl="$pl root-plugin-graf2d-gviz" test "x$enable_minuit2" = "xyes" && pl="$pl root-plugin-math-minuit2" test "x$enable_pythia6" = "xyes" && pl="$pl root-plugin-montecarlo-pythia6" @@ -5823,6 +5863,10 @@ sed -e "s|@srcdir@|$srcdir|" \ -e "s|@curseshdr@|$curseshdr|" \ -e "s|@curseslib@|$curseslib|" \ -e "s|@curseslibdir@|$curseslibdir|" \ + -e "s|@cfitsioincdir@|$cfitsioincdir|" \ + -e "s|@cfitsiolib@|$cfitsiolib|" \ + -e "s|@cfitsiolibdir@|$cfitsiolibdir|" \ + -e "s|@buildfitsio@|$enable_fitsio|" \ -e "s|@buildgl@|$enable_opengl|" \ -e "s|@buildldap@|$enable_ldap|" \ -e "s|@buildmysql@|$enable_mysql|" \ diff --git a/graf2d/fitsio/Module.mk b/graf2d/fitsio/Module.mk new file mode 100644 index 00000000000..e8ff2fdcba8 --- /dev/null +++ b/graf2d/fitsio/Module.mk @@ -0,0 +1,70 @@ +# Module.mk for cfitsio module +# Copyright (c) 2010 Rene Brun and Fons Rademakers +# +# Author: Claudi Martinez, 24/07/2010 + +MODNAME := fitsio +MODDIR := graf2d/$(MODNAME) +MODDIRS := $(MODDIR)/src +MODDIRI := $(MODDIR)/inc + +FITSIODIR := $(MODDIR) +FITSIODIRS := $(FITSIODIR)/src +FITSIODIRI := $(FITSIODIR)/inc + +##### libFITSIO ##### +FITSIOL := $(MODDIRI)/LinkDef.h +FITSIODS := $(MODDIRS)/G__FITSIO.cxx +FITSIODO := $(FITSIODS:.cxx=.o) +FITSIODH := $(FITSIODS:.cxx=.h) + +FITSIOH := $(filter-out $(MODDIRI)/LinkDef%,$(wildcard $(MODDIRI)/*.h)) +FITSIOS := $(filter-out $(MODDIRS)/G__%,$(wildcard $(MODDIRS)/*.cxx)) +FITSIOO := $(FITSIOS:.cxx=.o) + +FITSIODEP := $(FITSIOO:.o=.d) $(FITSIODO:.o=.d) + +FITSIOLIB := $(LPATH)/libFITSIO.$(SOEXT) +FITSIOMAP := $(FITSIOLIB:.$(SOEXT)=.rootmap) + +# used in the main Makefile +ALLHDRS += $(patsubst $(MODDIRI)/%.h,include/%.h,$(FITSIOH)) +ALLLIBS += $(FITSIOLIB) +ALLMAPS += $(FITSIOMAP) + +# include all dependency files +INCLUDEFILES += $(FITSIODEP) + +##### local rules ##### +.PHONY: all-$(MODNAME) clean-$(MODNAME) distclean-$(MODNAME) + +include/%.h: $(FITSIODIRI)/%.h + cp $< $@ + +$(FITSIOLIB): $(FITSIOO) $(FITSIODO) $(ORDER_) $(MAINLIBS) $(FITSIOLIBDEP) + @$(MAKELIB) $(PLATFORM) $(LD) "$(LDFLAGS)" \ + "$(SOFLAGS)" libFITSIO.$(SOEXT) $@ "$(FITSIOO) $(FITSIODO)" \ + "$(FITSIOLIBEXTRA) $(CFITSIOLIBDIR) $(CFITSIOLIB)" + +$(FITSIODS): $(FITSIOH) $(FITSIOL) $(ROOTCINTTMPDEP) + @echo "Generating dictionary $@..." + $(ROOTCINTTMP) -f $@ -c $(FITSIOH) $(FITSIOL) + +$(FITSIOMAP): $(RLIBMAP) $(MAKEFILEDEP) $(FITSIOL) + $(RLIBMAP) -o $(FITSIOMAP) -l $(FITSIOLIB) \ + -d $(FITSIOLIBDEPM) -c $(FITSIOL) + +all-$(MODNAME): $(FITSIOLIB) $(FITSIOMAP) + +clean-$(MODNAME): + @rm -f $(FITSIOO) $(FITSIODO) + +clean:: clean-$(MODNAME) + +distclean-$(MODNAME): clean-$(MODNAME) + @rm -f $(FITSIODEP) $(FITSIODS) $(FITSIODH) $(FITSIOLIB) $(FITSIOMAP) + +distclean:: distclean-$(MODNAME) + +##### extra rules ###### +$(FITSIOO) $(FITSDO): CXXFLAGS += $(CFITSIOINCDIR:%=-I%) diff --git a/graf2d/fitsio/inc/LinkDef.h b/graf2d/fitsio/inc/LinkDef.h new file mode 100644 index 00000000000..ae2ac4a2ffd --- /dev/null +++ b/graf2d/fitsio/inc/LinkDef.h @@ -0,0 +1,19 @@ +/* @(#)root/asimage:$Id$ */ + +/************************************************************************* + * Copyright (C) 1995-2010, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms seeg $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class TFITSHDU+; + +#endif diff --git a/graf2d/fitsio/inc/TFITS.h b/graf2d/fitsio/inc/TFITS.h new file mode 100644 index 00000000000..bdf89f85805 --- /dev/null +++ b/graf2d/fitsio/inc/TFITS.h @@ -0,0 +1,98 @@ +// @(#)root/graf2d:$Id$ +// Author: Claudi Martinez, July 19th 2010 + +/************************************************************************* + * Copyright (C) 1995-2010, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef ROOT_TFITS +#define ROOT_TFITS + +////////////////////////////////////////////////////////////////////////// +// // +// TFITS // +// // +// Interface to FITS astronomical files. // +// // +////////////////////////////////////////////////////////////////////////// + +#ifndef ROOT_TNamed +#include "TNamed.h" +#endif +#ifndef ROOT_TMatrixDfwd +#include "TMatrixDfwd.h" +#endif +#ifndef ROOT_TVectorDfwd +#include "TVectorDfwd.h" +#endif + +class TArrayI; +class TArrayD; +class TH1; +class TASImage; +class TImagePalette; + + +class TFITSHDU : public TNamed { + +private: + void _release_resources(); + void _initialize_me(); + +public: + enum EHDUTypes { // HDU types + kImageHDU, + kTableHDU + }; + + struct HDURecord { // FITS HDU record + TString fKeyword; + TString fValue; + TString fComment; + }; + +protected: + TString fCleanFilePath; // Path to HDU's file (without filter) + struct HDURecord *fRecords; // HDU metadata records + Int_t fNRecords; // Number of records + enum EHDUTypes fType; // HDU type + TString fExtensionName; // Extension Name + Int_t fNumber; // HDU number (1=PRIMARY) + TArrayI *fSizes; // Image sizes in each dimension (when fType == kImageHDU) + TArrayD *fPixels; // Image pixels (when fType == kImageHDU) + + Bool_t LoadHDU(TString& filepath_filter); + static void CleanFilePath(const char *filepath_with_filter, TString &dst); + void PrintHDUMetadata(const Option_t *opt="") const; + void PrintFileMetadata(const Option_t *opt="") const; + +public: + TFITSHDU(const char *filepath_with_filter); + TFITSHDU(const char *filepath, Int_t extension_number); + TFITSHDU(const char *filepath, const char *extension_name); + ~TFITSHDU(); + + Int_t GetRecordNumber() const { return fNRecords; } + struct HDURecord *GetRecord(const char *keyword); + TString& GetKeywordValue(const char *keyword); + void Print(const Option_t *opt="") const; + + //Image readers + TH1 *ReadAsHistogram(); + TASImage *ReadAsImage(Int_t layer = 0, TImagePalette *pal = 0); + TMatrixD *ReadAsMatrix(Int_t layer = 0); + TVectorD *GetArrayRow(UInt_t row); + TVectorD *GetArrayColumn(UInt_t col); + + //Table readers + //TODO + + ClassDef(TFITSHDU,1) // Class interfacing FITS HDUs +}; + + +#endif diff --git a/graf2d/fitsio/src/TFITS.cxx b/graf2d/fitsio/src/TFITS.cxx new file mode 100644 index 00000000000..febd534030e --- /dev/null +++ b/graf2d/fitsio/src/TFITS.cxx @@ -0,0 +1,723 @@ +// @(#)root/graf2d:$Id$ +// Author: Claudi Martinez, July 19th 2010 + +/************************************************************************* + * Copyright (C) 1995-2010, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +// IMPLEMENTATION NOTE: +// CFITSIO library uses standard C types ('int', 'long', ...). Since these +// types may depend on the compiler machine (in 32-bit CPU 'int' and 'long' are 32-bit), +// we use the standard C types too. Using types as Long_t (which is defined to be 64-bit), +// may lead to type size mismatch. + +//______________________________________________________________________________ +/* Begin_Html +<center><h2>FITS file interface class</h2></center> +TFITS is a class that allows extracting images and data from FITS files and contains +several methods to manage them. +End_Html */ + +#include "TFITS.h" +#include "TROOT.h" +#include "TASImage.h" +#include "TArrayI.h" +#include "TArrayD.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TH3D.h" +#include "TVectorD.h" +#include "TMatrixD.h" + +#include "fitsio2.h" +#include <stdlib.h> + + +ClassImp(TFITSHDU) + +/**************************************************************************/ +// TFITSHDU +/**************************************************************************/ + +//______________________________________________________________________________ +void TFITSHDU::CleanFilePath(const char *filepath_with_filter, TString &dst) +{ + // Clean path from possible filter and put the result in 'dst'. + + dst = filepath_with_filter; + + Ssiz_t ndx = dst.Index("[", 1, 0, TString::kExact); + if (ndx != kNPOS) { + dst.Resize(ndx); + } +} + + +//______________________________________________________________________________ +TFITSHDU::TFITSHDU(const char *filepath_with_filter) +{ + // TFITSHDU constructor from file path with HDU selection filter. + + _initialize_me(); + TString finalpath = filepath_with_filter; + CleanFilePath(filepath_with_filter, fCleanFilePath); + + if (kFALSE == LoadHDU(finalpath)) { + _release_resources(); + throw -1; + } +} + +//______________________________________________________________________________ +TFITSHDU::TFITSHDU(const char *filepath, Int_t extension_number) +{ + // TFITSHDU constructor from filepath and extension number. + + _initialize_me(); + CleanFilePath(filepath, fCleanFilePath); + + //Add "by extension number" filter + TString finalpath; + finalpath.Form("%s[%d]", fCleanFilePath.Data(), extension_number); + + if (kFALSE == LoadHDU(finalpath)) { + _release_resources(); + throw -1; + } +} + +//______________________________________________________________________________ +TFITSHDU::TFITSHDU(const char *filepath, const char *extension_name) +{ + // TFITSHDU constructor from filepath and extension name. + + _initialize_me(); + CleanFilePath(filepath, fCleanFilePath); + + //Add "by extension number" filter + TString finalpath; + finalpath.Form("%s[%s]", fCleanFilePath.Data(), extension_name); + + + if (kFALSE == LoadHDU(finalpath)) { + _release_resources(); + throw -1; + } +} + +//______________________________________________________________________________ +TFITSHDU::~TFITSHDU() +{ + // TFITSHDU destructor. + + _release_resources(); +} + +//______________________________________________________________________________ +void TFITSHDU::_release_resources() +{ + // Release internal resources. + + if (fRecords) delete [] fRecords; + + if (fType == kImageHDU) { + if (fSizes) delete fSizes; + if (fPixels) delete fPixels; + } else { + //TODO + } +} + +//______________________________________________________________________________ +void TFITSHDU::_initialize_me() +{ + // Do some initializations. + + fRecords = 0; + fPixels = 0; + fSizes = 0; +} + +//______________________________________________________________________________ +Bool_t TFITSHDU::LoadHDU(TString& filepath_filter) +{ + // Load HDU from fits file satisfying the specified filter. + // Returns kTRUE if success. Otherwise kFALSE. + // If filter == "" then the primary array is selected + + fitsfile *fp=0; + int status = 0; + char errdescr[FLEN_STATUS+1]; + + // Open file with filter + fits_open_file(&fp, filepath_filter.Data(), READONLY, &status); + if (status) goto ERR; + + // Read HDU number + int hdunum; + fits_get_hdu_num(fp, &hdunum); + fNumber = Int_t(hdunum); + + // Read HDU type + int hdutype; + fits_get_hdu_type(fp, &hdutype, &status); + if (status) goto ERR; + fType = (hdutype == IMAGE_HDU) ? kImageHDU : kTableHDU; + + //Read HDU header records + int nkeys, morekeys; + char keyname[FLEN_KEYWORD+1]; + char keyvalue[FLEN_VALUE+1]; + char comment[FLEN_COMMENT+1]; + + fits_get_hdrspace(fp, &nkeys, &morekeys, &status); + if (status) goto ERR; + + fRecords = new struct HDURecord[nkeys]; + + for (int i = 1; i <= nkeys; i++) { + fits_read_keyn(fp, i, keyname, keyvalue, comment, &status); + if (status) goto ERR; + fRecords[i-1].fKeyword = keyname; + fRecords[i-1].fValue = keyvalue; + fRecords[i-1].fComment = comment; + } + + fNRecords = Int_t(nkeys); + + //Set extension name + fExtensionName = "PRIMARY"; //Default + for (int i = 0; i < nkeys; i++) { + if (fRecords[i].fKeyword == "EXTNAME") { + fExtensionName = fRecords[i].fValue; + break; + } + } + + //Read HDU's data + if (fType == kImageHDU) { + //Image + int param_ndims=0; + long *param_dimsizes; + + //Read image number of dimensions + fits_get_img_dim(fp, ¶m_ndims, &status); + if (status) goto ERR; + if (param_ndims > 0) { + //Read image sizes in each dimension + param_dimsizes = new long[param_ndims]; + fits_get_img_size(fp, param_ndims, param_dimsizes, &status); + if (status) goto ERR; + + fSizes = new TArrayI(param_ndims); + fSizes = new TArrayI(param_ndims); + for (int i = 0; i < param_ndims; i++) { //Use for loop to copy values instead of passing array to constructor, since 'Int_t' size may differ from 'long' size + fSizes->SetAt(param_dimsizes[i], i); + } + + delete [] param_dimsizes; + + //Read pixels + int anynul; + long *firstpixel = new long[param_ndims]; + double nulval = 0; + long npixels = 1; + + for (int i = 0; i < param_ndims; i++) { + npixels *= (long) fSizes->GetAt(i); //Compute total number of pixels + firstpixel[i] = 1; //Set first pixel to read from. + } + + double *pixels = new double[npixels]; + + fits_read_pix(fp, TDOUBLE, firstpixel, npixels, + (void *) &nulval, (void *) pixels, &anynul, &status); + + if (status) { + delete [] firstpixel; + delete [] pixels; + goto ERR; + } + + fPixels = new TArrayD(npixels, pixels); + + delete [] firstpixel; + delete [] pixels; + + } else { + //Null array + fSizes = new TArrayI(); + fPixels = new TArrayD(); + } + } else { + //Table + //TODO + + } + + // Close file + fits_close_file(fp, &status); + return kTRUE; + +ERR: + fits_get_errstatus(status, errdescr); + Warning("LoadHDU", "error opening FITS file. Details: %s", errdescr); + status = 0; + if (fp) fits_close_file(fp, &status); + return kFALSE; +} + +//______________________________________________________________________________ +struct TFITSHDU::HDURecord* TFITSHDU::GetRecord(const char *keyword) +{ + // Get record by keyword. + + for (int i = 0; i < fNRecords; i++) { + if (fRecords[i].fKeyword == keyword) { + return &fRecords[i]; + } + } + return 0; +} + +//______________________________________________________________________________ +TString& TFITSHDU::GetKeywordValue(const char *keyword) +{ + // Get the value of a given keyword. Return "" if not found. + + HDURecord *rec = GetRecord(keyword); + if (rec) { + return rec->fValue; + } else { + return *(new TString("")); + } +} + +//______________________________________________________________________________ +void TFITSHDU::PrintHDUMetadata(const Option_t *) const +{ + // Print records. + + for (int i = 0; i < fNRecords; i++) { + if (fRecords[i].fComment.Length() > 0) { + printf("%-10s = %s / %s\n", fRecords[i].fKeyword.Data(), fRecords[i].fValue.Data(), fRecords[i].fComment.Data()); + } else { + printf("%-10s = %s\n", fRecords[i].fKeyword.Data(), fRecords[i].fValue.Data()); + } + } +} + +//______________________________________________________________________________ +void TFITSHDU::PrintFileMetadata(const Option_t *opt) const +{ + // Print HDU's parent file's metadata. + + fitsfile *fp=0; + int status = 0; + char errdescr[FLEN_STATUS+1]; + int hducount, extnum; + int hdutype = IMAGE_HDU; + const char *exttype; + char extname[FLEN_CARD]="PRIMARY"; //room enough + int verbose = (opt[0] ? 1 : 0); + + // Open file with no filters: current HDU will be the primary one. + fits_open_file(&fp, fCleanFilePath.Data(), READONLY, &status); + if (status) goto ERR; + + // Read HDU count + fits_get_num_hdus(fp, &hducount, &status); + if (status) goto ERR; + printf("Total: %d HDUs\n", hducount); + + extnum = 0; + while(hducount) { + // Read HDU type + fits_get_hdu_type(fp, &hdutype, &status); + if (status) goto ERR; + + if (hdutype == IMAGE_HDU) { + exttype="IMAGE"; + } else if (hdutype == ASCII_TBL) { + exttype="ASCII TABLE"; + } else { + exttype="BINARY TABLE"; + } + + //Read HDU header records + int nkeys, morekeys; + char keyname[FLEN_KEYWORD+1]; + char keyvalue[FLEN_VALUE+1]; + char comment[FLEN_COMMENT+1]; + + fits_get_hdrspace(fp, &nkeys, &morekeys, &status); + if (status) goto ERR; + + struct HDURecord *records = new struct HDURecord[nkeys]; + + for (int i = 1; i <= nkeys; i++) { + fits_read_keyn(fp, i, keyname, keyvalue, comment, &status); + if (status) { + delete [] records; + goto ERR; + } + + records[i-1].fKeyword = keyname; + records[i-1].fValue = keyvalue; + records[i-1].fComment = comment; + + if (strcmp(keyname, "EXTNAME") == 0) { + //Extension name + strcpy(extname, keyvalue); + } + } + + //HDU info + printf(" [%d] %s (%s)\n", extnum, exttype, extname); + + //HDU records + if (verbose) { + for (int i = 0; i < nkeys; i++) { + if (comment[0]) { + printf(" %-10s = %s / %s\n", records[i].fKeyword.Data(), records[i].fValue.Data(), records[i].fComment.Data()); + } else { + printf(" %-10s = %s\n", records[i].fKeyword.Data(), records[i].fValue.Data()); + } + } + } + printf("\n"); + + delete [] records; + + hducount--; + extnum++; + if (hducount){ + fits_movrel_hdu(fp, 1, &hdutype, &status); + if (status) goto ERR; + } + } + + // Close file + fits_close_file(fp, &status); + return; + +ERR: + fits_get_errstatus(status, errdescr); + Warning("PrintFileMetadata", "error opening FITS file. Details: %s", errdescr); + status = 0; + if (fp) fits_close_file(fp, &status); +} + +//______________________________________________________________________________ +void TFITSHDU::Print(const Option_t *opt) const +{ + // Print metadata. + // Currently supported options: + // "" : print HDU record data + // "F" : print FITS file's extension names, numbers and types + // "F+": print FITS file's extension names and types and their record data + + if ((opt[0] == 'F') || (opt[0] == 'f')) { + PrintFileMetadata((opt[1] == '+') ? "+" : ""); + } else { + PrintHDUMetadata(""); + } +} + + +//______________________________________________________________________________ +TASImage *TFITSHDU::ReadAsImage(Int_t layer, TImagePalette *pal) +{ + // Read image HDU as a displayable image. Return 0 if conversion cannot be done. + // If the HDU seems to be a multilayer image, 'layer' parameter can be used + // to retrieve the specified layer (starting from 0) + + if (fType != kImageHDU) { + Warning("ReadAsImage", "this is not an image HDU."); + return 0; + } + + if ((fSizes->GetSize() != 2) && (fSizes->GetSize() != 3)) { + Warning("ReadAsImage", "could not convert image HDU to image because it has %d dimensions.", fSizes->GetSize()); + return 0; + } + + Int_t width, height; + UInt_t pixels_per_layer; + + width = Int_t(fSizes->GetAt(0)); + height = Int_t(fSizes->GetAt(1)); + + pixels_per_layer = UInt_t(width) * UInt_t(height); + + if (((fSizes->GetSize() == 2) && (layer > 0)) || ((fSizes->GetSize() == 3) && (layer > fSizes->GetAt(2)))) { + Warning("ReadAsImage", "layer out of bounds."); + return 0; + } + + // Get the maximum and minimum pixel values in the layer to auto-stretch pixels + Double_t maxval = 0, minval = 0; + register UInt_t i; + Double_t pixvalue; + Int_t offset = layer * pixels_per_layer; + + for (i = 0; i < pixels_per_layer; i++) { + pixvalue = fPixels->GetAt(offset + i); + + if (pixvalue > maxval) { + maxval = pixvalue; + } + + if ((i == 0) || (pixvalue < minval)) { + minval = pixvalue; + } + } + + //Build the image stretching pixels into a range from 0.0 to 255.0 + TASImage *im = new TASImage(width, height); + TArrayD *layer_pixels = new TArrayD(pixels_per_layer); + + Double_t factor = 255.0 / (maxval-minval); + for (i = 0; i < pixels_per_layer; i++) { + pixvalue = fPixels->GetAt(offset + i); + layer_pixels->SetAt(factor * (pixvalue-minval), i) ; + } + + if (pal == 0) { + // Default palette: grayscale palette + pal = new TImagePalette(256); + for (i = 0; i < 256; i++) { + pal->fPoints[i] = ((Double_t)i)/255.0; + pal->fColorAlpha[i] = 0xffff; + pal->fColorBlue[i] = pal->fColorGreen[i] = pal->fColorRed[i] = i << 8; + } + pal->fPoints[0] = 0; + pal->fPoints[255] = 1.0; + + im->SetImage(*layer_pixels, UInt_t(width), pal); + + delete pal; + + } else { + im->SetImage(*layer_pixels, UInt_t(width), pal); + } + + delete layer_pixels; + + return im; +} + +//______________________________________________________________________________ +TMatrixD* TFITSHDU::ReadAsMatrix(Int_t layer) +{ + // Read image HDU as a matrix. Return 0 if conversion cannot be done + // If the HDU seems to be a multilayer image, 'layer' parameter can be used + // to retrieve the specified layer (starting from 0) in matrix form + + if (fType != kImageHDU) { + Warning("ReadAsMatrix", "this is not an image HDU."); + return 0; + } + + + if ((fSizes->GetSize() != 2) && (fSizes->GetSize() != 3)) { + Warning("ReadAsMatrix", "could not convert image HDU to image because it has %d dimensions.", fSizes->GetSize()); + return 0; + } + + if (((fSizes->GetSize() == 2) && (layer > 0)) || ((fSizes->GetSize() == 3) && (layer > fSizes->GetAt(2)))) { + Warning("ReadAsMatrix", "layer out of bounds."); + return 0; + } + + Int_t width, height; + UInt_t pixels_per_layer; + Int_t offset; + + width = Int_t(fSizes->GetAt(0)); + height = Int_t(fSizes->GetAt(1)); + + pixels_per_layer = UInt_t(width) * UInt_t(height); + offset = layer * pixels_per_layer; + + + TMatrixD *mat = new TMatrixD(height, width); + double *layer_pixels = new double[pixels_per_layer]; + + for (UInt_t i = 0; i < pixels_per_layer; i++) { + layer_pixels[i] = fPixels->GetAt(offset + i); + } + + mat->Use(height, width, layer_pixels); + + return mat; +} + + +//______________________________________________________________________________ +TH1 *TFITSHDU::ReadAsHistogram() +{ + // Read image HDU as a histogram. Return 0 if conversion cannot be done. + // The returned object can be TH1D, TH2D or TH3D depending on data dimensionality. + // Please, check condition (returnedValue->IsA() == TH*D::Class()) to + // determine the object class. + // NOTE: do not confuse with image histogram! This function interprets + // the array as a histogram. It does not compute the histogram of pixel + // values of an image! Here "pixels" are interpreted as number of entries. + + if (fType != kImageHDU) { + Warning("ReadAsHistogram", "this is not an image HDU."); + return 0; + } + + TH1 *result=0; + + if ((fSizes->GetSize() != 1) && (fSizes->GetSize() != 2) && (fSizes->GetSize() != 3)) { + Warning("ReadAsHistogram", "could not convert image HDU to histogram because it has %d dimensions.", fSizes->GetSize()); + return 0; + } + + if (fSizes->GetSize() == 1) { + //1D + UInt_t Nx = UInt_t(fSizes->GetAt(0)); + UInt_t x; + + TH1D *h = new TH1D("", "", Int_t(Nx), 0, Nx-1); + + for (x = 0; x < Nx; x++) { + Long_t nentries = Long_t(fPixels->GetAt(x)); + if (nentries < 0) nentries = 0; //Crop negative values + h->Fill(x, nentries); + } + + result = h; + + } else if (fSizes->GetSize() == 2) { + //2D + UInt_t Nx = UInt_t(fSizes->GetAt(0)); + UInt_t Ny = UInt_t(fSizes->GetAt(1)); + UInt_t x,y; + + TH2D *h = new TH2D("", "", Int_t(Nx), 0, Nx-1, Int_t(Ny), 0, Ny-1); + + for (y = 0; y < Ny; y++) { + UInt_t offset = y * Nx; + for (x = 0; x < Nx; x++) { + Long_t nentries = Long_t(fPixels->GetAt(offset + x)); + if (nentries < 0) nentries = 0; //Crop negative values + h->Fill(x,y, nentries); + } + } + + result = h; + + } else if (fSizes->GetSize() == 3) { + //3D + UInt_t Nx = UInt_t(fSizes->GetAt(0)); + UInt_t Ny = UInt_t(fSizes->GetAt(1)); + UInt_t Nz = UInt_t(fSizes->GetAt(2)); + UInt_t x,y,z; + + TH3D *h = new TH3D("", "", Int_t(Nx), 0, Nx-1, Int_t(Ny), 0, Ny-1, Int_t(Nz), 0, Nz-1); + + + for (z = 0; z < Nz; z++) { + UInt_t offset1 = z * Nx * Ny; + for (y = 0; y < Ny; y++) { + UInt_t offset2 = y * Nx; + for (x = 0; x < Nx; x++) { + Long_t nentries = Long_t(fPixels->GetAt(offset1 + offset2 + x)); + if (nentries < 0) nentries = 0; //Crop negative values + h->Fill(x, y, z, nentries); + } + } + } + + result = h; + } + + return result; +} + +//______________________________________________________________________________ +TVectorD* TFITSHDU::GetArrayRow(UInt_t row) +{ + // Get a row from the image HDU when it's a 2D array. + + if (fType != kImageHDU) { + Warning("GetArrayRow", "this is not an image HDU."); + return 0; + } + + if (fSizes->GetSize() != 2) { + Warning("GetArrayRow", "could not get row from HDU because it has %d dimensions.", fSizes->GetSize()); + return 0; + } + + UInt_t i, offset, W,H; + + W = UInt_t(fSizes->GetAt(0)); + H = UInt_t(fSizes->GetAt(1)); + + + if (row >= H) { + Warning("GetArrayRow", "index out of bounds."); + return 0; + } + + offset = W * row; + double *v = new double[W]; + + for (i = 0; i < W; i++) { + v[i] = fPixels->GetAt(offset+i); + } + + TVectorD *vec = new TVectorD(W, v); + + delete [] v; + + return vec; +} + +//______________________________________________________________________________ +TVectorD* TFITSHDU::GetArrayColumn(UInt_t col) +{ + // Get a column from the image HDU when it's a 2D array. + + if (fType != kImageHDU) { + Warning("GetArrayColumn", "this is not an image HDU."); + return 0; + } + + if (fSizes->GetSize() != 2) { + Warning("GetArrayColumn", "could not get row from HDU because it has %d dimensions.", fSizes->GetSize()); + return 0; + } + + UInt_t i, W,H; + + W = UInt_t(fSizes->GetAt(0)); + H = UInt_t(fSizes->GetAt(1)); + + + if (col >= W) { + Warning("GetArrayColumn", "index out of bounds."); + return 0; + } + + double *v = new double[H]; + + for (i = 0; i < H; i++) { + v[i] = fPixels->GetAt(W*i+col); + } + + TVectorD *vec = new TVectorD(H, v); + + delete [] v; + + return vec; +} -- GitLab