diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,23 @@ +Copyright Zac Slade (c) 2018 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,24 @@ +[![Build Status](https://travis-ci.com/krakrjak/fits-parse.svg?branch=master)](https://travis-ci.com/krakrjak/fits-parse) +[![Documentation Status](https://readthedocs.org/projects/fits-parse/badge/?version=latest)](https://fits-parse.readthedocs.io/en/latest/?badge=latest) + +# fits-parse + +FITS - Flexible Image Transport System + +This project focuses on building a Haskell native parser for FITS files. This file format is used widely by Astronomers and those analyzing astronomical data. There is currently a C library called `fitsio` and a Haskell wrapper for the library. However, this library can be cumbersome to use and the Haskell bindings are very thin and close to the metal of the `fitsio` library. This project focuses on providing good Haskell bindings and datatypes by implementing the FITS spec directly in Haskell, without the need for a C library. The reference for this work comes from NASA in the US. Conformance is currently limited to the Version 4.0 specification with no extensions. See [the NASA site](https://fits.gsfc.nasa.gov/fits_standard.html) for more details and to download the spec. + +# Reading Documentation + +The [readthedocs](https://fits-parse.readthedocs.org) pages contain documentation for the *examples*. +The [API Haddock](https://krakrjak.github.io/fits-parse/index.html) docs for the *master* branch are +on github.io, release documentation will be on [Hackage](https://hackage.haskell.org/package/fits-parse). + +You can build the Haddocks for offline use with `stack build --haddock`. + +# RTD Building + +To build the [readthedocs](https://fits-parse.readthedocs.org) docs you need `sphinx`. To install it +locally for a single user use `pip` with the `--user` option to install `sphinx` and `sphinx-rtd-theme`. +Alternatively, on Ubuntu you can install the `python3-sphinx` and `python3-sphinx-rtd-theme` package. + + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = sphinx-build +SPHINXPROJ = fits-parse +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# fits-parse documentation build configuration file, created by +# sphinx-quickstart on Wed Apr 25 14:25:18 2018. +# +# This file is execfile()d with the current directory set to its +# containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +# import os +# import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- General configuration ------------------------------------------------ + +# If your documentation needs a minimal Sphinx version, state it here. +# +# needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# General information about the project. +project = 'fits-parse' +copyright = '2018, Zac Slade' +author = 'Zac Slade' + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +version = '0.0.1' +# The full version, including alpha/beta/rc tags. +release = '0.0.1' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This patterns also effect to html_static_path and html_extra_path +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = False + + +# -- Options for HTML output ---------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'alabaster' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Custom sidebar templates, must be a dictionary that maps document names +# to template names. +# +# This is required for the alabaster theme +# refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars +html_sidebars = { + '**': [ + 'relations.html', # needs 'show_related': True theme option to display + 'searchbox.html', + ] +} + + +# -- Options for HTMLHelp output ------------------------------------------ + +# Output file base name for HTML help builder. +htmlhelp_basename = 'fits-parsedoc' + + +# -- Options for LaTeX output --------------------------------------------- + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + # 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + # + 'pointsize': '12pt', + + # Additional stuff for the LaTeX preamble. + # + # 'preamble': '', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'fits-render.tex', 'fits-render Documentation', + 'Zac Slade', 'manual'), +] + + +# -- Options for manual page output --------------------------------------- + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'fits-render', 'fits-render Documentation', + [author], 1) +] + + +# -- Options for Texinfo output ------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'fits-render', 'fits-render Documentation', + author, 'fits-render', 'One line description of project.', + 'Miscellaneous'), +] diff --git a/docs/examples/omnibus.rst b/docs/examples/omnibus.rst new file mode 100644 --- /dev/null +++ b/docs/examples/omnibus.rst @@ -0,0 +1,46 @@ +Using the *omnibus* Example Program +=================================== + +The *omnibus* example program is a proving ground for the library by attempting +straightforward use case examples and tying them all into one program. + +The example *omnibus* program can be run as follows:: + + fits-parse$ stack exec -- omnibus + +It has a few options out might want to be aware of:: + + fits-parse$ stack exec -- omnibus --help + fits-parse - a FITS swiss army knife + + Usage: omnibus (FILE | --stdin) ([-o|--outfile ARG] | [--stdout]) + Parse, generate, and render FITS files + + Available options: + FILE Input FITS file to parse + --stdin Input FITS from stdin + -o,--outfile ARG Output filename to store + --stdout Output picture to stdout + -h,--help Show this help text + +Running *omnibus* on the example fits file *Spiral_2_30_0_300_10_0_NoGrad.fits* +gives an output like the following:: + + fits-parse$ stack exec omnibus example_fits_files/Spiral_2_30_0_300_10_0_NoGrad.fits + 29/Mar/2019:21:03:27 -0500: [DEBUG] input file size 1545444 bytes + 29/Mar/2019:21:03:27 -0500: [DEBUG] found 1 hdu record(s) + 29/Mar/2019:21:03:27 -0500: [DEBUG] Bit Format 32 bit IEEE single precision float + 29/Mar/2019:21:03:27 -0500: [DEBUG] data block size 1542564 bytes + 29/Mar/2019:21:03:27 -0500: [DEBUG] 2 Axes + 29/Mar/2019:21:03:27 -0500: [DEBUG] Axis: 1 count: 621 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Axis: 2 count: 621 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Total Pix Count: 385641 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Unwrapped Int Count: 0 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Unwrapped Double Count: 385641 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Vector Int Count: 0 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Vector Double Count: 385641 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Mean: 3.161230263379672e-2 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Range: 1.0 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Variance: 3.061296495598597e-2 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Std deviation: 0.17496583763186624 + 29/Mar/2019:21:03:27 -0500: [DEBUG] Std error of the mean: 2.817485308081582e-4 diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,25 @@ +.. image:: https://travis-ci.com/krakrjak/fits-parse.svg?branch=master + :target: https://travis-ci.com/krakrjak/fits-parse +.. image:: https://readthedocs.org/projects/fits-parse/badge/?version=latest + :target: https://fits-parse.readthedocs.io/en/latest/?badge=latest + :alt: Documentation Status + +Welcome to fits-parse's documentation! +====================================== + +If you use FITS files in your work or play and you use Hakell, or would like +to, then *fits-parse* just might be the library you've been looking for. This +project is still very young and only really works on the example FITS file it +ships with, but this will improve and may already be better than anticipated. + +In the **examples** folder in the source distribution you can find usage of +the library. The *omnibus* program glues together the basics and does logging +and is intended to be a kitchen sink sort of example program. + +`Haddock Docs `_ are the API docs for the project. + +.. toctree:: + :maxdepth: 2 + + examples/omnibus.rst + diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,36 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build +set SPHINXPROJ=fits-parse + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/examples/omnibus/Main.hs b/examples/omnibus/Main.hs new file mode 100644 --- /dev/null +++ b/examples/omnibus/Main.hs @@ -0,0 +1,147 @@ +module Main where + +-- base +import Control.Monad ( unless ) +import Data.Monoid( (<>) ) +import Control.Applicative( (<|>), (<**>) ) +import System.IO ( stdin ) + +import Options.Applicative( Parser, strOption, optional, short, long + , flag', metavar, help, execParser, info + , fullDesc, header, helper, progDesc + , argument, str ) + +import qualified Data.ByteString as BS +import qualified Data.ByteString.Lazy as LBS + +import qualified Data.Vector as V + +---- statistics +import Statistics.Sample ( range, mean, variance, stdDev, stdErrMean ) + +import System.Log.FastLogger( TimedFastLogger, ToLogStr, LogType( LogStderr ) + , defaultBufSize, newTimeCache, simpleTimeFormat + , toLogStr, newTimedFastLogger, withTimedFastLogger ) + + +-- local library imports +import Data.Fits ( HeaderDataUnit(..), HeaderData(..), Axis(..) + , parsePix, isBitPixInt, isBitPixFloat, pixsUnwrapI, pixsUnwrapD ) +import Data.Fits.MegaParser ( getAllHDUs ) + +-- | Paramaterized input type for files or standard input. +data Input = FileInput FilePath -- ^ The 'FileInput' constructor needs a path name + | StdInput -- ^ The `StdInput` constructor stands in for the obvious + +-- | Parameterized output type for files or standard output. +data Output = FileOutput FilePath -- ^ When you have a file use 'FileOutput' + | StdOutput -- ^ For a screen dump use 'StdOutput' + +fileInput :: Parser Input +fileInput = FileInput <$> argument str + ( metavar "FILE" + <> help "Input FITS file to parse" ) + +stdInput :: Parser Input +stdInput = flag' StdInput + ( long "stdin" + <> help "Input FITS from stdin" ) + +fileOutput :: Parser Output +fileOutput = FileOutput <$> strOption + ( long "outfile" + <> short 'o' + <> help "Output filename to store" ) + +stdOutput :: Parser Output +stdOutput = flag' StdOutput + ( long "stdout" + <> help "Output picture to stdout" ) + +inputParse :: Parser Input +inputParse = fileInput <|> stdInput + +outputParse :: Parser Output +outputParse = fileOutput <|> stdOutput + +data FitsConfig = FitsConfig + { infile :: Input + , outImage :: Maybe Output + } + +fitsConfigParser :: Parser FitsConfig +fitsConfigParser = FitsConfig <$> inputParse <*> optional outputParse + +main :: IO () +main = workOnFITS =<< execParser opts + where + opts = info (fitsConfigParser <**> helper) + ( fullDesc + <> progDesc "Parse, generate, and render FITS files" + <> header "fits-parse - a FITS swiss army knife" ) + +workOnFITS :: FitsConfig -> IO () +workOnFITS (FitsConfig i o) = do + timeCache <- newTimeCache simpleTimeFormat + withTimedFastLogger timeCache (LogStderr defaultBufSize) $ \logger -> do + fits <- bs i + myLog logger $ "[DEBUG] input file size " ++ show (BS.length fits) ++ " bytes\n" + hdus <- getAllHDUs fits + myLog logger ("[DEBUG] found " ++ show (length hdus) ++ " hdu record(s)\n") + mapM_ (processHDU logger) hdus + where + bs (FileInput f) = BS.readFile f + bs StdInput = BS.hGetContents stdin + +processHDU :: TimedFastLogger -> HeaderDataUnit -> IO () +processHDU logger hdu = do + myLog logger $ "[DEBUG] Bit Format " ++ show bpf ++ "\n" + myLog logger $ "[DEBUG] data block size " ++ show (BS.length pd) ++ " bytes\n" + myLog logger $ "[DEBUG] " ++ show (length ax) ++ " Axes\n" + unless (null ax) (mapM_ logAxes ax) + let pixCount = foldr (\x acc -> axisElementCount x * acc) 1 ax + pixs <- parsePix pixCount bpf (LBS.fromStrict pd) + let pxsI = if isBitPixInt bpf then pixsUnwrapI bpf pixs else [] + pxsD = if isBitPixFloat bpf then pixsUnwrapD bpf pixs else [] + pVI = V.fromList pxsI + pVD = V.fromList pxsD + myLog logger $ "[DEBUG] Total Pix Count: " ++ show (length pixs) ++ "\n" + myLog logger $ "[DEBUG] Unwrapped Int Count: " ++ show (length pxsI) ++ "\n" + myLog logger $ "[DEBUG] Unwrapped Double Count: " ++ show (length pxsD) ++ "\n" + myLog logger $ "[DEBUG] Vector Int Count: " ++ show (length pVI) ++ "\n" + myLog logger $ "[DEBUG] Vector Double Count: " ++ show (length pVD) ++ "\n" + if (length ax == 2) && isBitPixFloat bpf + then + bitMapProcess logger ax pVD + else + myLog logger "[DEBUG] skipping bitmap analysis.\n" + + where + hd = headerData hdu + pd = payloadData hdu + ax = axes hd + bpf = bitPixFormat hd + logAxes a = myLog logger $ "[DEBUG] Axis: " + ++ show (axisNumber a) + ++ " count: " + ++ show (axisElementCount a) ++ "\n" + +{- | If we happen to be working in floating point 2D, let's try the + following. +-} +bitMapProcess :: TimedFastLogger + -> [Axis] -- ^ Metadata about the column oriented axes + -> V.Vector Double -- ^ Data is stored in column-row major order + -> IO () +bitMapProcess logger [] _ = myLog logger "[ERROR] BitMap processing run with no axes.\n" +bitMapProcess logger [_] _ = myLog logger "[ERROR] BitMap processing run with only one axis.\n" +bitMapProcess logger (y:x) v = do + myLog logger $ "[DEBUG] Mean: " ++ show (mean v) ++ "\n" + myLog logger $ "[DEBUG] Range: " ++ show (range v) ++ "\n" + myLog logger $ "[DEBUG] Variance: " ++ show (variance v) ++ "\n" + myLog logger $ "[DEBUG] Std deviation: " ++ show (stdDev v) ++ "\n" + myLog logger $ "[DEBUG] Std error of the mean: " ++ show (stdErrMean v) ++ "\n" + + +myLog:: ToLogStr msg => TimedFastLogger -> msg -> IO () +myLog logger msg = logger $ \ft -> toLogStr ft <> toLogStr ": " <> toLogStr msg diff --git a/fits-parse.cabal b/fits-parse.cabal new file mode 100644 --- /dev/null +++ b/fits-parse.cabal @@ -0,0 +1,67 @@ +cabal-version: >=1.10 +name: fits-parse +version: 0.0.1 +license: BSD2 +license-file: LICENSE +copyright: Copyright (c) 2019 Zac Slade +maintainer: krakrjak@gmail.com +author: Zac Slade +homepage: https://github.com/krakrjak/fits-parse#readme +synopsis: Parse FITS files +description: + Parse and manipulate FITS data natively in Haskell +category: Science +build-type: Simple +extra-source-files: + docs/conf.py + docs/examples/omnibus.rst + docs/index.rst + docs/make.bat + docs/Makefile + fits_files/nonconformant/testkeys.fits + fits_files/nonconformant/testkeys2.fits + fits_files/Spiral_2_30_0_300_10_0_NoGrad.fits + README.md + +flag examples + description: + Do you want to build the examples? + default: False + manual: True + +library + exposed-modules: + Data.Fits + Data.Fits.MegaParser + hs-source-dirs: src + other-modules: + Paths_fits_parse + default-language: Haskell2010 + build-depends: + base >=4.7 && <5, + binary >=0.8.6.0 && <0.9, + bytestring >=0.10.8.2 && <0.11, + data-default >=0.7.1.1 && <0.8, + megaparsec >=7.0.4 && <7.1, + parser-combinators >=1.0.1 && <1.1, + text >=1.2.3.1 && <1.3, + text-latin1 >=0.3.1 && <0.4 + +executable omnibus + main-is: Main.hs + hs-source-dirs: examples/omnibus + other-modules: + Paths_fits_parse + default-language: Haskell2010 + build-depends: + JuicyPixels >=3.3.3 && <3.4, + base >=4.7 && <5, + bytestring >=0.10.8.2 && <0.11, + fast-logger >=2.4.13 && <2.5, + fits-parse -any, + optparse-applicative >=0.14.3.0 && <0.15, + statistics >=0.15.0.0 && <0.16, + vector >=0.12.0.2 && <0.13 + + if flag(examples) + buildable: True diff --git a/fits_files/Spiral_2_30_0_300_10_0_NoGrad.fits b/fits_files/Spiral_2_30_0_300_10_0_NoGrad.fits new file mode 100644 Binary files /dev/null and b/fits_files/Spiral_2_30_0_300_10_0_NoGrad.fits differ diff --git a/fits_files/nonconformant/testkeys.fits b/fits_files/nonconformant/testkeys.fits new file mode 100644 Binary files /dev/null and b/fits_files/nonconformant/testkeys.fits differ diff --git a/fits_files/nonconformant/testkeys2.fits b/fits_files/nonconformant/testkeys2.fits new file mode 100644 Binary files /dev/null and b/fits_files/nonconformant/testkeys2.fits differ diff --git a/src/Data/Fits.hs b/src/Data/Fits.hs new file mode 100644 --- /dev/null +++ b/src/Data/Fits.hs @@ -0,0 +1,364 @@ +{-| +Module : Data.Fits +Description : Types for FITS Data Units +Copyright : (c) Zac Slade, 2018 +License : BSD2 +Maintainer : krakrjak@gmail.com +Stability : experimental + +Definitions for the data types needed to parse an HDU in a FITS block. +-} + +{-# LANGUAGE PartialTypeSignatures, DataKinds, ExistentialQuantification + , ScopedTypeVariables, GADTs + , OverloadedStrings, TypeOperators, TypeFamilies #-} +module Data.Fits + ( -- * Data payload functions + parsePix + , pixsUnwrapI + , pixsUnwrapD + + -- * Main data types + , HeaderDataUnit(..) + , HeaderData(..) + , BitPixFormat(..) + , Pix(..) + + -- ** Header Data Types + , SimpleFormat(..) + , Axis(..) + , StringType(..) + , StringValue(..) + , NumberType(..) + , NumberModifier(..) + , NumberValue(..) + + -- * Utility + , isBitPixInt + , isBitPixFloat + , bitPixToWordSize + + -- ** Constants + , hduRecordLength + , hduMaxRecords + , hduBlockSize + + ) where + +---- text +import qualified Data.Text as T +---- bytestring +import qualified Data.ByteString as BS +import qualified Data.ByteString.Lazy as BL + +---- base +import Numeric.Natural ( Natural ) + +---- ghc +import GHC.TypeNats (KnownNat, Nat) + +---- text +import Data.Text ( Text ) + +---- bytestring +import Data.ByteString ( ByteString ) + +---- default +import Data.Default( Default, def ) + +---- binaryi +import Data.Binary +import Data.Binary.Get + +-- | A single record in the HDU is an eighty byte word. +{-@ type HDURecordLength = {v:Int | v = 80} @-} +{-@ hduRecordLength :: HDURecordLength @-} +hduRecordLength :: Int +hduRecordLength = 80 + +{-| The maximum amount of eighty byte records is thirty-six per the + standard. +-} +{-@ type HDUMaxRecords = {v:Int | v = 36} @-} +{-@ hduMaxRecords :: HDUMaxRecords @-} +hduMaxRecords :: Int +hduMaxRecords = 36 + +{-| The size of an HDU block is fixed at thirty-six eighty byte words. In + other words 2,880 bytes. These blocks are padded with zeros to this + boundary. +-} +{-@ type HDUBlockSize = {v:Int | v = 2880} @-} +{-@ hduBlockSize :: HDUBlockSize @-} +hduBlockSize :: Int +hduBlockSize = hduRecordLength * hduMaxRecords + +{-| There are many types of strings defined in the FITS documentation. + Refered to as "character string(s)" in the documentation, they can be + null, empty, undefined, or contain characters (printable ASCII only). +-} +data StringType = NullString -- ^ The NULL character string (e.g. AUTHOR=) + | EmptyString -- ^ An empty string (e.g. AUTHOR="") + | DataString -- ^ Plain ASCII data + +instance Show StringType where + show NullString = "Null String" + show EmptyString = "Empty quoted String" + show DataString = "String" + +-- | A 'StringValue' is a type paired with a possible value. +data StringValue = StringValue + { stringType :: StringType -- ^ Which 'StringType' is this value? + , stringValue :: Maybe Text -- ^ The payload of a character string + } + +{-| The default instance of 'StringValue' is an undefined string with no + text. +-} +instance Default StringValue where + def = StringValue NullString Nothing + +instance Show StringValue where + show (StringValue NullString _) = show NullString + show (StringValue EmptyString _) = show EmptyString + show (StringValue DataString Nothing) = "No good " ++ show DataString + show (StringValue DataString (Just s)) = show DataString ++ T.unpack s + +{-| The FITS standard allows for the encoding of unsigned integers, signed + integers, real numbers, and complex numbers. They are always ASCII + encoded. See 5.2 of the standard for more details. +-} +data NumberType = + IntegerType -- ^ HDU ASCII encoded integer number + | RealType -- ^ HDU ASCII encoded real number + | ComplexType -- ^ HDU ASCII encoded complex number + +-- | Utility data type to help with the ASCII representation of numbers +data NumberModifier = + Positive -- ^ HDU positive number value + | Negative -- ^ HDU negative number value + | Zero -- ^ HDU numeric value is zero, could be positive or negative + +{-| 'NumberValue' contains an encoded numeric record from a data field. + This data type still needs to be converted into more useful Haskell data + types. +-} +data NumberValue = NumberValue + { numberType :: NumberType + -- ^ Key to decoding the structure + , realModifier :: NumberModifier + -- ^ Encoding the sign of the real part + , realPart :: Text + -- ^ All 'NumberType' have a real part (sign stripped) + , imaginaryModifier :: Maybe NumberModifier + -- ^ Encoding the sign of the imaginary part + , imaginaryPart :: Maybe Text + -- ^ Only 'ComplexType' have an imaginary part (sign stripped) + , exponentModifier :: Maybe NumberModifier + -- ^ 'Positive', 'Negative', or 'Zero' + , exponent :: Maybe Int + -- ^ All 'NumberType' may have an exponent + } + +{-| The default instance for 'NumberValue' is a 32-bit integer of zero with + no imaginary or exponent parts. +-} +instance Default NumberValue where + def = NumberValue IntegerType Zero "0" Nothing Nothing Nothing Nothing + +{-| The standard defines two possible values for the SIMPLE keyword, T and + F. The T refers to a 'Conformant' format while F refers to + a 'NonConformant' format. At this time only the 'Conformant', T, format + is supported. +-} +data SimpleFormat = Conformant + -- ^ Value of SIMPLE=T in the header. /supported/ + | NonConformant + -- ^ Value of SIMPLE=F in the header. /unsupported/ + +-- | 'Axis' represents a single NAXIS record. +data Axis = Axis + { axisNumber :: Int -- ^ The axis number under consideration + , axisElementCount :: Int -- ^ The number of elements in this axis + } + +-- | The default instance for 'Axis' is NAXIS=0 with zero elements. +instance Default Axis where + def = Axis 0 0 + +{-| The 'BitPixFormat' is the nitty gritty of how the 'Axis' data is layed + out in the file. The standard recognizes six formats: unsigned 8 bit + integer, two's complement binary integers at 16, 32, and 64 bits along + with 32 and 64 bit IEEE floating point formats. +-} +data BitPixFormat = + EightBitInt -- ^ BITPIX = 8; unsigned binary integer of 8 bits + | SixteenBitInt -- ^ BITPIX = 16; two's complement binary integer of 16 bits + | ThirtyTwoBitInt -- ^ BITPIX = 32; two's complement binary integer of 32 bits + | SixtyFourBitInt -- ^ BITPIX = 64; two's complement binary integer of 64 bits + | ThirtyTwoBitFloat -- ^ BITPIX = -32; IEEE single precision floating point of 32 bits + | SixtyFourBitFloat -- ^ BITPIX = -64; IEEE double precision floating point of 64 bits + +instance Show BitPixFormat where + show EightBitInt = "8 bit unsigned integer" + show SixteenBitInt = "16 bit signed integer" + show ThirtyTwoBitInt = "32 bit signed integer" + show SixtyFourBitInt = "64 bit signed interger" + show ThirtyTwoBitFloat = "32 bit IEEE single precision float" + show SixtyFourBitFloat = "64 bit IEEE double precision float" + +{-| This utility function can be used to get the word count for data in an + HDU. +-} +bitPixToWordSize :: BitPixFormat -> Natural +bitPixToWordSize EightBitInt = 8 +bitPixToWordSize SixteenBitInt = 16 +bitPixToWordSize ThirtyTwoBitInt = 32 +bitPixToWordSize ThirtyTwoBitFloat = 32 +bitPixToWordSize SixtyFourBitInt = 64 +bitPixToWordSize SixtyFourBitFloat = 64 + +{-| This utility function can be used to get the size in bytes of the +- format. +-} +bitPixToByteSize :: BitPixFormat -> Natural +bitPixToByteSize EightBitInt = 1 +bitPixToByteSize SixteenBitInt = 2 +bitPixToByteSize ThirtyTwoBitInt = 4 +bitPixToByteSize ThirtyTwoBitFloat = 4 +bitPixToByteSize SixtyFourBitInt = 8 +bitPixToByteSize SixtyFourBitFloat = 8 + +{- | This utility functions quickly lets you know if you are dealing with + integer data. +-} +isBitPixInt :: BitPixFormat -> Bool +isBitPixInt EightBitInt = True +isBitPixInt SixteenBitInt = True +isBitPixInt ThirtyTwoBitInt = True +isBitPixInt SixtyFourBitInt = True +isBitPixInt _ = False + +{- | This utility functions quickly lets you know if you are dealing with + floating point data. +-} +isBitPixFloat :: BitPixFormat -> Bool +isBitPixFloat ThirtyTwoBitFloat = True +isBitPixFloat SixtyFourBitFloat = True +isBitPixFloat _ = False + +{- | Following `BitPixFormat` we have a tag for integer and floating point + values. We box them up to ease parsing. +-} +data Pix = PB Int | PI16 Int | PI32 Int | PI64 Int | PF Double | PD Double + +{- | Removes the `Pix` tag from an `Int` type within. -} +unPixI :: Pix -> Int +unPixI (PB b) = b +unPixI (PI16 i) = i +unPixI (PI32 i) = i +unPixI (PI64 i) = i +unPixI _ = error "Pix are not stored as integers, invalid unpacking" + +{- | Removes the `Pix` tag from a `Double` type within. -} +unPixD :: Pix -> Double +unPixD (PF d) = d +unPixD (PD d) = d +unPixD _ = error "Pix are not stored as floating point values, invalid unpacking" + +{- | Remove the Pix wrapper for integer `Pix` lists. -} +pixsUnwrapI :: BitPixFormat -> [Pix] -> [Int] +pixsUnwrapI EightBitInt pxs = map unPixI pxs +pixsUnwrapI SixteenBitInt pxs = map unPixI pxs +pixsUnwrapI ThirtyTwoBitInt pxs = map unPixI pxs +pixsUnwrapI SixtyFourBitInt pxs = map unPixI pxs +pixsUnwrapI _ _ = error "BitPixFormat is not an integer type" + +{- | Remove the `Pix` wrapper for floating point `Pix` lists. -} +pixsUnwrapD :: BitPixFormat -> [Pix] -> [Double] +pixsUnwrapD ThirtyTwoBitFloat pxs = map unPixD pxs +pixsUnwrapD SixtyFourBitFloat pxs = map unPixD pxs +pixsUnwrapD _ _ = error "BitPixFormat is not a floating point type" + +getPix :: BitPixFormat -> Get Pix +getPix EightBitInt = PB . fromIntegral <$> getInt8 +getPix SixteenBitInt = PI16 . fromIntegral <$> getInt16be +getPix ThirtyTwoBitInt = PI32 . fromIntegral <$> getInt32be +getPix SixtyFourBitInt = PI64 . fromIntegral <$> getInt64be +getPix ThirtyTwoBitFloat = PF . realToFrac <$> getFloatbe +getPix SixtyFourBitFloat = PD . realToFrac <$> getDoublebe + +getPixs :: Int -> BitPixFormat -> Get [Pix] +getPixs c bpf = do + empty <- isEmpty + if empty + then return [] + else do + p <- getPix bpf + ps <- getPixs (c - 1) bpf + return (p:ps) + +{- | This is the main low-level function which parses the data portion of an + HDU. You need and element count, a format and a bytestring. The resulting + list is produced in column-row major order as specified in the standard. +-} +parsePix :: Int -> BitPixFormat -> BL.ByteString -> IO [Pix] +parsePix c bpf bs = return $ runGet (getPixs c bpf) bs + +{- `pixDimsByCol` takes a list of Axis and gives a column-row major list of + axes dimensions. +-} +pixDimsByCol :: [Axis] -> [Int] +pixDimsByCol = map axisElementCount + +{- `pixDimsByRow` takes a list of Axis and gives a row-column major list of + axes dimensions. +-} +pixDimsByRow :: [Axis] -> [Int] +pixDimsByRow = reverse . pixDimsByCol + +{-| The header part of the HDU is vital carrying not only authorship + metadata, but also specifying how to make sense of the binary payload + that starts 2,880 bytes after the start of the 'HeaderData'. +-} + +data HeaderData = HeaderData + { simpleFormat :: SimpleFormat + -- ^ SIMPLE + , bitPixFormat :: BitPixFormat + -- ^ BITPIX + , axes :: [Axis] + -- ^ Axes metadata + , objectIdentifier :: StringValue + -- ^ OBJECT + , observationDate :: StringValue + -- ^ DATE + , originIdentifier :: StringValue + -- ^ OBJECT + , telescopeIdentifier :: StringValue + -- ^ TELESCOP + , instrumentIdentifier :: StringValue + -- ^ INSTRUME + , observerIdentifier :: StringValue + -- ^ OBSERVER + , authorIdentifier :: StringValue + -- ^ CREATOR + , referenceString :: StringValue + -- ^ REFERENC + } + +instance Default HeaderData where + def = HeaderData NonConformant EightBitInt [] + (def :: StringValue) (def :: StringValue) (def :: StringValue) + (def :: StringValue) (def :: StringValue) (def :: StringValue) + (def :: StringValue) (def :: StringValue) + +{-| The 'HeaderDataUnit' is the full HDU. Both the header information is + encoded alongside the 'Axis' payload. +-} +data HeaderDataUnit = HeaderDataUnit + { headerData :: HeaderData + -- ^ Just the header part of the HDU + , payloadData :: ByteString + -- ^ The actual data payload + } diff --git a/src/Data/Fits/MegaParser.hs b/src/Data/Fits/MegaParser.hs new file mode 100644 --- /dev/null +++ b/src/Data/Fits/MegaParser.hs @@ -0,0 +1,219 @@ +{-| +Module : Data.Fits.MegaParser +Description : MegaParsec based parser for an HDU. +Copyright : (c) Zac Slade, 2018 +License : BSD2 +Maintainer : krakrjak@gmail.com +Stability : experimental + +Parsing rules for an HDU in a FITS file. +-} + +{-# OPTIONS_HADDOCK hide #-} +{-# LANGUAGE OverloadedStrings #-} + +module Data.Fits.MegaParser where + +-- qualified imports +import qualified Data.ByteString as BS +import qualified Data.ByteString.Internal as BS ( c2w ) +import qualified Text.Megaparsec as M +import qualified Text.Megaparsec.Stream as M +import qualified Text.Megaparsec.Char as M +import qualified Text.Megaparsec.Char.Lexer as MCL +import qualified Data.Text as T +import qualified Data.Text.Encoding as TE + +-- explicit imports +import Control.Applicative ( (<$>) ) +import Control.Monad ( void ) +import Numeric.Natural ( Natural ) +import Data.ByteString ( ByteString ) +import Data.Text ( Text ) +import Data.Word ( Word8, Word16, Word32, Word64 ) +import Data.Void ( Void ) +import Data.Default ( def ) +import Text.Ascii ( isAscii ) +import Text.Megaparsec ( Parsec, ParseErrorBundle, (<|>), ()) + +-- full module imports +---- parser-combinators +import Control.Applicative.Permutations +---- base +import Data.Proxy + +-- local imports +import Data.Fits + +type Parser = Parsec Void Text +type ParseErr = ParseErrorBundle Text Void + +data DataUnitValues + = FITSUInt8 Word8 + | FITSInt16 Word16 + | FITSInt32 Word32 + | FITSInt64 Word64 + | FITSFloat32 Float + | FITSFloat64 Double + +headerBlockParse :: Parser HeaderData +headerBlockParse = do + -- TODO: Parse other sections of the neader for god's sake + -- HISTORY and COMMENT along with CONTINUE handling is missing + (simple, bitpix, axesDesc, bZero, bScale, ref, obs, instr, tele, object, pCreator, pDate, pEnd) <- + runPermutation $ + (,,,,,,,,,,,,) <$> toPermutation (parseSimple "simple") + <*> toPermutation (parseBitPix "bitpix") + <*> toPermutation ((parseAxisCount >>= parseNaxes) "axis parsing") + <*> toPermutationWithDefault 0 (parseBzero "bzero") + <*> toPermutationWithDefault 0 (parseBscale "bscale") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseReference "reference") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseObserver "observer") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseInstrument "instrument") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseTelescope "telescope") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseObject "object") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseCreator "creator") + <*> toPermutationWithDefault (StringValue NullString Nothing) (parseDate "date") + <*> toPermutation (parseEnd "end") + return defHeader { simpleFormat = simple + , bitPixFormat = bitpix + , axes = axesDesc + , referenceString = ref + , observerIdentifier = obs + , instrumentIdentifier = instr + , telescopeIdentifier = tele + , objectIdentifier = object + , observationDate = pDate + , authorIdentifier = pCreator } + where + defHeader = def :: HeaderData + +parseSimple :: Parser SimpleFormat +parseSimple = M.string' "simple" >> parseEquals + >> ((conformParse >> consumeDead >> return Conformant) + <|> (nonConformParse >> consumeDead >> return NonConformant)) + where + conformParse = M.char' 't' + nonConformParse = M.anySingle + +parseEquals :: Parser () +parseEquals = M.space >> M.char '=' >> M.space + +parseBitPix :: Parser BitPixFormat +parseBitPix = M.string' "bitpix" >> parseEquals + >> ((M.chunk "8" >> consumeDead >> return EightBitInt) + <|> (M.chunk "16" >> consumeDead >> return SixteenBitInt) + <|> (M.chunk "32" >> consumeDead >> return ThirtyTwoBitInt) + <|> (M.chunk "64" >> consumeDead >> return SixtyFourBitInt) + <|> (M.chunk "-32" >> consumeDead >> return ThirtyTwoBitFloat) + <|> (M.chunk "-64" >> consumeDead >> return SixtyFourBitFloat)) + +parseAxisCount :: Parser Natural +parseAxisCount = M.string' "naxis" >> parseEquals >> parseNatural + +parseNaxes :: Natural -> Parser [Axis] +parseNaxes n | n == 0 = return [] +parseNaxes n = do + axisNum <- M.string' "naxis" >> parseNatural + elemCount <- parseEquals >> parseNatural + ([buildAxis axisNum elemCount] ++) <$> parseNaxes (n - 1) + where + defAxis = def :: Axis + buildAxis an ec = defAxis { axisNumber = fromIntegral an + , axisElementCount = fromIntegral ec } + +parseBzero :: Parser Int +parseBzero = M.string' "bzero" >> parseEquals >> parseInteger + +parseBscale :: Parser Int +parseBscale = M.string' "bscale" >> parseEquals >> parseInteger + +parseReference :: Parser StringValue +parseReference = M.string' "referenc" >> parseEquals >> parseStringValue + +parseObserver :: Parser StringValue +parseObserver = M.string' "observer" >> parseEquals >> parseStringValue + +parseInstrument :: Parser StringValue +parseInstrument = M.string' "instrume" >> parseEquals >> parseStringValue + +parseTelescope :: Parser StringValue +parseTelescope = M.string' "telescop" >> parseEquals >> parseStringValue + +parseObject :: Parser StringValue +parseObject = M.string' "object" >> parseEquals >> parseStringValue + +parseCreator :: Parser StringValue +parseCreator = M.string' "creator" >> parseEquals >> parseStringValue + +parseDate :: Parser StringValue +parseDate = M.string' "date" >> parseEquals >> parseStringValue + +skipEmpty :: Parser () +skipEmpty = void (M.many $ M.satisfy ('\0' ==)) + +consumeDead :: Parser () +consumeDead = M.space >> skipEmpty + +parseEnd :: Parser () +parseEnd = M.string' "end" >> M.space <* M.eof + +parseNatural :: Parser Natural +parseNatural = do + v <- MCL.decimal + consumeDead + return $ fromIntegral v + +parseInteger :: Parser Int +parseInteger = do + v <- MCL.decimal + consumeDead + return v + +parseStringValue :: Parser StringValue +parseStringValue = do + -- The rules are weird, NULL means a NULL string, '' is an empty + -- string, a ' followed by a bunch of spaces and a close ' is + -- considered an empty string, and trailing whitespace is ignored + -- within the quotes, but not leading spaces. + ls <- M.between (M.char '\'') (M.char '\'') $ M.many $ M.anySingleBut '\'' + consumeDead + let v = M.tokensToChunk (Proxy :: Proxy Text) ls + if T.length v < 1 + then return (StringValue EmptyString Nothing) + else return (StringValue DataString (Just v)) + +countHeaderDataUnits :: ByteString -> IO Natural +countHeaderDataUnits bs = fromIntegral . length <$> getAllHDUs bs + +-- TODO: make the recursive case work. Currently limited to one HDU. +-- The current issue is that when the parser fails on an HDU parse, it +-- blows them all up instead of accepting the valid parsings. +getAllHDUs :: ByteString -> IO [HeaderDataUnit] +getAllHDUs bs = do + (hdu, rest) <- getOneHDU bs + return [hdu] +-- if BS.length rest < hduBlockSize then return [hdu] else return [hdu] + +getOneHDU :: ByteString -> IO (HeaderDataUnit, ByteString) +getOneHDU bs = + if isAscii header + then + case M.runParser headerBlockParse "FITS" (TE.decodeUtf8 header) of + Right mainHeader -> do + let (dataUnit, remainder) = BS.splitAt (fromIntegral $ dataSize mainHeader) rest + return (HeaderDataUnit mainHeader dataUnit, remainder) + Left e -> let err = M.errorBundlePretty e in error err + else error "Header data is not ASCII. Please Check your input file and try again" + where + (header, rest) = BS.splitAt hduBlockSize bs + +dataSize :: HeaderData -> Natural +dataSize h = paddedsize + where + axesCount = length $ axes h + wordCount = product $ map axisElementCount $ axes h + wordsize = fromIntegral . bitPixToWordSize $ bitPixFormat h + datasize = wordsize * wordCount + padding = if axesCount == 0 then 0 else fromIntegral hduBlockSize - datasize `mod` fromIntegral hduBlockSize + paddedsize = fromIntegral (datasize + padding)