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'