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.