Astronomical image processing guide (How to list of values of a FITS image?)
FITS Fortran Guide
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:
- open of FITS
- get of its size (width, height)
- allocate memory for a matrix
- read data
- 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'