Astronomical image processing guide (How to list of info keywords?)

A great advantage of the FITS format is support of any arbitrary information included into an image file. A place for the info is reserved in header of the FITS file. Every FITS header must contain a set of records (lines) containing certain mandatory keywords. The set may be followed by records with any other keywords.

The mandatory keywords are: SIMPLE, EXTEND and END. All others are optional. The records of header has the fixed structure:

1-8 - keyword
9 - = equal's sign
10-x - value (various lengths)
x -80 - comment (every text following slash '/')

All records are 80-bytes long and there is no separator between ones. The special keywords COMMENT and HISTORY has no defined a structure and introduce any text string (continuation lines are not allowed). Recent recommendations to contents of the records adds physical units to numerical values as a part of comment. Any set of keywords is not widely used. It means that various utilities may used different keyword for the same entity.

The most common keywords are included in following (artificial) example:

SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 382 / length of data axis 1
NAXIS2 = 255 / length of data axis 2
EXTEND = T / FITS dataset may contain extensions
EXPTIME = 60.000 / [s] Exposure time
DATE-OBS= '2008-06-06T01:02:41.390' / UTC of exposure start
FILTER = 'R ' / filter
OBJECT = 'Star' / Object name
COMMENT This file was written by XXX.
END

We have more ways to handle with FITS header keywords. Practically, we will need to list of all keywords (as a variant of more or less command on unix's shell) or list of value of some specified record. Both situations can be easy coded. We introduce of utility FITSless which will print all keywords (full header) without any switches and a specified value when any keyword will presented, its value will be printed.

The values are usually different kinds (from computer point of view). For example, the observer's name will coded as a string, an exposure time will a real number, the date of start will a specially formatted string. To print the information, we will use only string representation, but in real code it will be better use of appropriate data type. Unfortunately, there is no simple way how to do it in Fortran, C and cFITSIO under recent versions of ones.

program FITSless

implicit none

integer :: status
! status ... FITS status (0=no error)

integer :: j,blocksize

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to open

character(len=666) :: keyword
! keyword to print the record

character(len=666) :: record
! header's record

character(len=666) :: value, comment
! record's value & comment

integer :: nhead, hpos
! number of records in header and current position

! get first command line parameter
call get_command_argument(1, keyword)

status = 0
call ftopen(25,name,0,blocksize,status)
if( status /= 0 ) stop 'File not found.'

! list specified keyword or full header
if( keyword /= ' ' ) then

write(*,*) 'List of a record specified by your keyword:'
write(*,*) '-------------------------------------------'

call ftgkys(25,keyword,value, comment, status)
! checking status of the operation
if( status == 0 ) then
write(*,*) trim(keyword),' = ',trim(value),' /',trim(comment)
else
write(*,*) 'Keyword "',trim(keyword),'" not found.'
end if

else

write(*,*) 'Full header list:'
write(*,*) '-----------------'

call ftghps(25,nhead,hpos,status)
do j = 1, nhead
call ftgrec(25,j,record,status)
write(*,*) trim(record)
enddo

end if

call ftclos(25,status)

end program FITSless

The code is self-explanatory, I recommend try this following experiments:

• try give as parameter an arbitrary text string or upper/lower case keywords
• process any printed information by some another text tool (sed, grep,..)
• play with listing of a record on given position in the header
• modify code to read an input file name by command-line parameters

Astronomical image processing guide (How to save of an image to a FITS file?)

One of most common operation done on FITS files is creation of a new image. The creation process is straightforward. One is required to have an image (result of a mathematical algorithm) with known dimensions and that's all.

The FITS format supports up to seven-dimensional images without any strict limitation of its axis ranges. Practically, we've some limits due to physical properties of present computers. The most important is a size of the output image. The size is determined (approximately) as a product of image dimensions multiplied by number of bytes occupied by one pixel. For two-dimensional image of 100x100 pixels in both axis, represented by real numbers (4-byte, BITPIX=-32), we got size about 40 kilo-bytes. As we know from previous lecture, the images produced by an algorithm are usually saved as real data to preserve its numerical precision.

As example of a test image, I choose Bessel's function of zero kind J₀ which represents light's diffraction on cylindrical aperture. We can observe of square of Bessel's function in an ideal telescope as the image of a star (of course, we use only J₀ for simplicity, a star will looks differently). J₀ is non-standard Fortran function and may be not supported by all compilers (gfortran does it) so one is optional only and you can play with cosine under a strict-standard compiler.


An algorithm to save an image to FITS file is straightforward:

1. create image as an array
2. fill the array by an image
3. initiate a FITS
4. setup image parameters
5. save the data

program FITSsave

implicit none

integer :: status, bitpix, naxis, naxes(2)
! status ... FITS status (0=no error)
! naxis .. number of axes in the image (we set =2)
! naxes .. dimensions of the image (2-element array)

real, dimension(:,:), allocatable :: d
! data matrix

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to create

! aux
real :: x,y,r
integer :: i,j

! set dimensions of the new image
naxis = 2
naxes = (/ 100,100 /)

! set bipix of the image (try: 8,16,32 and -32)
bitpix = -32

! create a data storage (allocate memory) for the image
allocate(d(naxes(1),naxes(2)))

! fill image with values
do i = 1, naxes(1)
do j = 1, naxes(2)

! rectangular coordinates
! the left bottom pixel has 1,1
x = i
y = j

! distance from origin
r = sqrt(x**2 + y**2)

! set value
d(i,j) = cos(r/5.0)
!d(i,j) = besj0(r/5.0)
! uncomment for J0 (Bessel function of zero kind)
! J0 represents diffraction on cylindrical aperture
end do
end do

! save the data
status = 0
call ftinit(26,name,1,status)
call ftphps(26,bitpix,naxis,naxes,status)
call ftp2de(26,1,naxes(1),naxes(1),naxes(2),d,status)
call ftclos(26,status)

! free allocated memory
deallocate(d)

end program FITSsave

The code can be compiled and run by

host$ gfortran -Wall -o FITSsave FITSsave.f90 -L/usr/local/lib -lcfitsio
host$ ./FITSsave

the output file is named as image.fits and can be viewed by any FITS viewer, for example by ds9:

host$ ds9 image.fits


Notes.

Any handle with numerical operations in many computer languages may be little bit confusing. Fortran strictly distinguish between integer and real numbers. The notation 3/4 (both integers) products result 0 (reminder is forget), but 4/3 gives 1. Opposite with this, the notation 3.0/4.0 gives 0.75 (two significant places).

The function ftinit (the initial function for FITS file) can create only a new file. There is no way how to replace any existing file. This is simply a feature of cfitsio, not a bug. It means that you must remove the older file image.fits before run FITSsave again.

Astronomical image processing guide (How to list of values of a FITS image?)

Every FITS file is representation of an image usually created by an optical device. The image is quantized (sampled) to elementary cells called pixels. An information knows for any pixel are: a pixel coordinate (Cartesian x,y) and its captured optical flux (CCD's flux is a linear function of captured photons).

The pixels which represents an image, are rearranged and a captured flux is digitalised to a matrix. The matrix is saved to a FITS file by a defined algorithm. The pixels coordinates (integers) are arranged in that matrix by the way:

[M,1] [M,2] .. [M,N]
.. ..
[2,1] [2,2] .. [2,N]
[1,1] [1,2] .. [1,N]

The matrix represents an image of width of N pixels and M pixels of height. The origin and orientation is in usual mathematical fashion.

Every pixel is represented by a number. The kind of the number may be an integer and a real (real numbers are with fractional part). Raw images (pure product of a device) are usually represented by integer numbers in interval from 0 to 65535 (2^16) for CCD and from 0 to 4096 (2^12) for a digital camera. A processed images as result of mathematical operations are saved as a real numbers with floating point. (Arithmetical operations may reduce of its precision). The method (complicated on first sight) to store of data reflects many of astronomer's needs and save your disk space. The data representation is coded in parameter BITPIX in the fits header by the way (not all possibilities are included):

BITPIX bytes type range of values
8 1 integer 0 .. 255
16 2 integer 0 .. 65535
32 4 integer 0 .. 4294836225
-32 4 real -1e38 .. 1e38 (7 digits)

An operation on a image included in a FITS file is relative easy. Follow the instruction:

1. open of FITS
2. get of its size (width, height)
3. allocate memory for a matrix
4. read data
5. play with data

Fortran offers a very effective way for manipulation with matrixes. For a matrix D, we can select an i,j-element as D(i,j), a i-row D(i,:),i-column D(:,j) or submatrix D(1:10,50:60).

program FITSlist

implicit none

integer :: status, bitpix, naxis, naxes(2)
! status ... FITS status (0=no error)
! naxis .. number of axes in image (we require =2)
! naxes .. dimension of the image (2-element array)

integer :: i,j,blocksize,pcount,gcount,minvalue
logical :: extend, simple,anyf
! required by cFITSIO

real, dimension(:,:), allocatable :: d
! data matrix

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to open

status = 0
call ftopen(25,name,0,blocksize,status)
call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
allocate(d(naxes(1),naxes(2)))
call ftg2de(25,1,minvalue,naxes(1),naxes(1),naxes(2),d,anyf,status)
call ftclos(25,status)

! print value of a random pixel
write(*,*) '# d(1,1)=',d(1,1)

! print last tree values of first row
write(*,*) '# d(1,-10:)=',d(1,size(d,2)-3:)

! print of a submatrix with indexes
do i = 1,10
do j = 50,60
write(*,*) i,j,d(i,j)
end do
end do

end program FITSlist

The output of the code can be saved to a file by sequence of commands:

host$ gfortran -Wall -o FITSlist FITSlist.f90 -L/usr/local/lib -lcfitsio
host$ ./FITSlist > pixels

and easy plotted in a gnuplot with the command:

gnuplot> splot 'pixels'

Astronomical image processing guide (How to open a FITS image?)

A structure of a FITS file is a little bit complicated, but not too much. The FITS itself includes a header and a data part. The header of an image consists from meta-information about the image. The most important are the size (width and height) and data representation of the image.

The header is set of 80-byte (character) length records. Every record is represented by text line with structure:
KEYWORD = VALUE
The KEYWORD must be no longer than 8 characters. The '=' must be in 9 column. The width and height of an image is coded by the way:
NAXIS1 = 1628
NAXIS2 = 1236

How to get an image size?

To get this basic information by use of the cfitsio library, we can use of the piece of the code (source code):

! to compile: gfortran -Wall -o FITSsize o.f90 -L/usr/local/lib -lcfitsio

program FITSsize

implicit none

integer :: status, bitpix, naxis, naxes(2)
! status ... FITS status (0=no error)
! naxis .. number of axes in image (we require =2)
! naxes .. dimension of the image (2-element array)

integer :: blocksize,pcount,gcount
logical :: extend, simple
! required by cFITSIO

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to open
status = 0
call ftopen(25,name,0,blocksize,status)
call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
call ftclos(25,status)

if( status == 0 ) then
write(*,*) 'The image ',trim(name),' has the size:',naxes
else
write(*,*) 'The image "',trim(name),'" not found or not accessible.'
end if

end program FITSsize

The code may by compiled by command:

host$ gfortran -Wall -o FITSsize FITSsize.f90 -L/usr/local/lib -lcfitsio

where switch -Wall prints some warnings, -o specify name of the generated binary file (name of the routine), -L points path to cfitsio library (may be omitted, usually any system directory) and -lcfitsio links cFITSIO library (libcfitsio.a).

Type:

host$ ./FITSsize

to run. The utility will try to open file named as 'image.fits' (can be changed in declaration name = 'image.fits'). If this file is accessible it will print the size of the image. If the name can't be open, an error will appeared.

This code demonstrates basic idioms of FITS-specific ones:

• its look horribly
• there is a lot of declarations which meaning is too hard to remember
• the variable status must be set to zero before calling of any of cFITSIO routines
• the name of the image to open is a second argument to ftopen
• the routine ftghpr reads important parameters (coded by KEYWORDS as above) of the image

Astronomical image processing guide (Intro)

From time to time, I'm thinking about the most useful way to process of an astronomical image data. An use of wide-known software package utilities like Gaia, ds9, IRAF may be better or worse in comparison of a direct coding of a self-made routine? Both approaches has its own advantages and disadvantages so, I think, there is no an universal way to process its. For example, a lot of work can be done inside IRAF's (http://iraf.noao.edu/) environment, but sometimes may be faster and more suitable of coding of an own utility. Both approaches may be possible useful and important especially for a particular processing. Also, I thinks the develop of any own routine may be much, much better and simpler then use of prepared ones. While there is a lot of documents describing of various software packages, the processing in a computer language is described poorly. That's why I'm starting write this guide.

The basic operation to work with data is data file handling. A FITS format (http://heasarc.gsfc.nasa.gov/docs/heasarc/fits.html) is a wide-used format to storing of an astronomical data including both an images and a data tables. A general structure of FITS files may be really very complicated. Fortunately, W. Pence and etc. made the cFITSIO library to provide of a standard way to create, open and modify of any FITS files. The library offers interfaces for C and Fortran languages. There is also a lot of wrappers to others (Perl, Python, C++, ..). I'll use of Fortran (fortran 90/95/2003 dialect) in this guide to create a simple code but the change to C (for example) is straightforward.

A computer library (like cfitsio) is a file with a special structure which includes a set of (useful) functions. The file is usually created by compiling of a computer language (cfitsio is in C) to a machine-specific code.

How to install an environment to open FITS files?

1. Use of your package system and install gfortran (Fortran compiler).
2. Use of your package system and install cfitsio.

ad 1). The gfortran is included in most modern Linuxes. Its is possible to use also Intel's ifort (apparently faster). Solaris offers f90/f95. The GNU g95 provides binaries for other systems (BSD's, Mac OSX, ..).

ad 2) Many modern distributions (Debian, Ubuntu, Fedora..) offers cfitsio as a package. If your doesn't, download directly tarball. Move it to /tmp and execute:
host:/tmp$ tar zxf cfitsio-x.y.z.tar.gz
host:/tmp$ cd cfistio-x.y.z
host:/tmp$ ./configure # check output
host:/tmp$ make # compile
host:/tmp$ su # switch to root account
host:/tmp$ make install # install to /usr/local

To install it in an another place use --prefix parameter for configure (./configure --help). To uninstall it, type (as root) make uninstall. A successful installation is indicated by presence of drvrsmem.h,fitsio.h, fitsio2.h and longnam.h in /usr/local/include and libcfitsio.a in /usr/local/lib.