1#!/usr/bin/env python
2"""Module to difference two FITS images."""
3
4import sys
5import os
6import numpy as np
7import pyfits
8
9
10def save_fits_image(filename, data, header, img1=None, img2=None):
11    """Save a FITS image."""
12    # Reshape to add the frequency axis
13    data = np.reshape(data, (1, 1, data.shape[0], data.shape[1]))
14    new_hdr = pyfits.header.Header()
15    for i, item in enumerate(header.items()):
16        if item[0] != 'HISTORY':
17            new_hdr.append(item)
18    if img1 and img2:
19        new_hdr.append(('HISTORY', '' * 60))
20        new_hdr.append(('HISTORY', '-' * 60))
21        new_hdr.append(('HISTORY', 'Diff created from image1 - image2:'))
22        new_hdr.append(('HISTORY', '- image1 : %s' % img1))
23        new_hdr.append(('HISTORY', '- image2 : %s' % img2))
24        new_hdr.append(('HISTORY', '-' * 60))
25        new_hdr.append(('HISTORY', '  - Max       : % .3e' % np.max(data)))
26        new_hdr.append(('HISTORY', '  - Min       : % .3e' % np.min(data)))
27        new_hdr.append(('HISTORY', '  - Mean      : % .3e' % np.mean(data)))
28        new_hdr.append(('HISTORY', '  - STD       : % .3e' % np.std(data)))
29        new_hdr.append(('HISTORY', '  - RMS       : % .3e' %
30                        np.sqrt(np.mean(data**2))))
31        new_hdr.append(('HISTORY', '-' * 60))
32        new_hdr.append(('HISTORY', '' * 60))
33
34    if (os.path.exists(filename)):
35        print '+ WARNING, output FITS file already exsits, overwriting.'
36        os.remove(filename)
37    pyfits.writeto(filename, data, new_hdr)
38
39
40def load_fits_image(filename):
41    """Load a FITS image."""
42
43    hdulist = pyfits.open(filename)
44    data = hdulist[0].data
45    hdr = hdulist[0].header
46    return np.squeeze(data), hdr
47
48
49def fits_diff(outname, file1, file2):
50    if not os.path.exists(file1):
51        print 'ERROR: missing input file', file1
52        return
53    if not os.path.exists(file2):
54        print 'ERROR: missing input file', file2
55        return
56    if os.path.exists(outname):
57        return
58
59    f1, h1 = load_fits_image(file1)
60    f2, h2 = load_fits_image(file2)
61    diff = f1 - f2
62    print '-' * 80
63    print '+ Image size  : %i x %i' % (f1.shape[0], f1.shape[1])
64    print '+ File 1      : %s' % file1
65    print '+ File 2      : %s' % file2
66    print '+ Diff        : file1 - file2'
67    print '+ Output name : %s' % (outname)
68    print '+ Diff stats:'
69    print '  - Max       : % .3e' % np.max(diff)
70    print '  - Min       : % .3e' % np.min(diff)
71    print '  - Mean      : % .3e' % np.mean(diff)
72    print '  - STD       : % .3e' % np.std(diff)
73    print '  - RMS       : % .3e' % np.sqrt(np.mean(diff**2))
74    print '-' * 80
75    save_fits_image(outname, diff, h1, file1, file2)
76
77
78if __name__ == '__main__':
79
80    if len(sys.argv) - 1 != 3:
81        print 'Usage: fits_diff.py <diff name> <file1> <file2>'
82        print ''
83        print 'Performs: diff = file1 - file2'
84        sys.exit(1)
85
86    outname = sys.argv[1]
87    file1 = sys.argv[2]
88    file2 = sys.argv[3]
89
90    f1, h1 = load_fits_image(file1)
91    f2, h2 = load_fits_image(file2)
92
93    diff = f1 - f2
94
95    print '-' * 80
96    print '+ Image size  : %i x %i' % (f1.shape[0], f1.shape[1])
97    print '+ File 1      : %s' % file1
98    print '+ File 2      : %s' % file2
99    print '+ Diff        : file1 - file2'
100    print '+ Output name : %s' % (outname)
101    print '+ Diff stats:'
102    print '  - Max       : % .3e' % np.max(diff)
103    print '  - Min       : % .3e' % np.min(diff)
104    print '  - Mean      : % .3e' % np.mean(diff)
105    print '  - STD       : % .3e' % np.std(diff)
106    print '  - RMS       : % .3e' % np.sqrt(np.mean(diff**2))
107    print '-' * 80
108
109    save_fits_image(outname, diff, h1, file1, file2)
110