Python notes 7/23/07 -- Sharpness Tools


Last weekend I wrote 4 Python programs which concern raw FITS image sharpness:

All programs, and examples, are available in the Sharpness project. This is an XCode project, however XCode is only being used as an editor and file manager: there are no targets or executables in the project. The programs are run under Python 2.5.1 (either with 'python program.py' or from the '#!/usr/bin/env python' line in the file), and you can use BBEdit (etc...) as the program editor.

sharpness.py

Usage: [python] sharpness.py directory [> output-file], or
       [python] sharpness.py fits-file [min] [max]

To run on a folder of files, use:

> sharpness.py contk1
contk1/data1050.a.fits    2004-07-13           16:17:23.972260           0.04539629     
contk1/data1047.a.fits    2004-07-13           16:17:16.547470           0.04383086     
contk1/data1048.a.fits    2004-07-13           16:17:19.018101           0.04270498     
contk1/data1049.a.fits    2004-07-13           16:17:21.525703           0.04204120     
contk1/data1045.a.fits    2004-07-13           16:17:11.548302           0.03994618     
contk1/data1042.a.fits    2004-07-13           16:17:03.969390           0.03979242     
contk1/data1043.a.fits    2004-07-13           16:17:06.517683           0.03971004     
contk1/data1044.a.fits    2004-07-13           16:17:08.984247           0.03949092     
contk1/data1046.a.fits    2004-07-13           16:17:13.994263           0.03925742     
contk1/data1041.a.fits    2004-07-13           16:17:01.523765           0.03888214     

The program reads all .fits files in the folder, estimates their sharpness, and prints the sorted results. You can also redirect this output to a text file. The program takes less than 0.5 seconds to process each image, on a 2x2GHz MacBook.

To estimate the sharpness of a single image, use:

> sharpness.py contk1/data1041.a.fits 
Min radius      Max radius      Max wavelength  Min wavelength  Count   Score
----------      ----------      --------------  --------------  -----   -----
0.00000000      8.00000000      512.00000000    64.00000000     197     0.73594810          
8.00000000      16.00000000     64.00000000     32.00000000     600     0.13945106          
16.00000000     24.00000000     32.00000000     21.33333333     996     0.05120139          
24.00000000     32.00000000     21.33333333     16.00000000     1416    0.02514732          
32.00000000     40.00000000     16.00000000     12.80000000     1816    0.01377158          
40.00000000     48.00000000     12.80000000     10.66666667     2188    0.00819796          
48.00000000     56.00000000     10.66666667     9.14285714      2632    0.00515955          
56.00000000     64.00000000     9.14285714      8.00000000      3008    0.00346567          
64.00000000     72.00000000     8.00000000      7.11111111      3388    0.00242458          
72.00000000     80.00000000     7.11111111      6.40000000      3840    0.00167302          
80.00000000     88.00000000     6.40000000      5.81818182      4232    0.00115471          
88.00000000     96.00000000     5.81818182      5.33333333      4604    0.00090070          
96.00000000     104.00000000    5.33333333      4.92307692      5032    0.00072927          
104.00000000    112.00000000    4.92307692      4.57142857      5432    0.00057161          
112.00000000    120.00000000    4.57142857      4.26666667      5844    0.00044305          
120.00000000    128.00000000    4.26666667      4.00000000      6208    0.00039044          
128.00000000    136.00000000    4.00000000      3.76470588      6656    0.00033251          
136.00000000    144.00000000    3.76470588      3.55555556      7012    0.00029439          
144.00000000    152.00000000    3.55555556      3.36842105      7432    0.00025345          
152.00000000    160.00000000    3.36842105      3.20000000      7848    0.00021591          
160.00000000    168.00000000    3.20000000      3.04761905      8252    0.00019948          
168.00000000    176.00000000    3.04761905      2.90909091      8640    0.00019078          
176.00000000    184.00000000    2.90909091      2.78260870      9080    0.00017908          
184.00000000    192.00000000    2.78260870      2.66666667      9428    0.00018725          
192.00000000    200.00000000    2.66666667      2.56000000      9848    0.00018150          
200.00000000    208.00000000    2.56000000      2.46153846      10248   0.00019118          
208.00000000    216.00000000    2.46153846      2.37037037      10628   0.00019263          
216.00000000    224.00000000    2.37037037      2.28571429      11080   0.00022758          
224.00000000    232.00000000    2.28571429      2.20689655      11464   0.00029652          
232.00000000    240.00000000    2.20689655      2.13333333      11868   0.00041801          
240.00000000    248.00000000    2.13333333      2.06451613      12288   0.00082729          
248.00000000    256.00000000    2.06451613      2.00000000      12654   0.00518241          
512 - 64: 0.73594810          
64 - 32:  0.13945106          
32 - 16:  0.07634871          
16 - 8:   0.03059476          
8 - 4:    0.00828738          
4 - 2:    0.00936998          
sharpness = 0.03888214          
contk1/data1041.a.fits    2004-07-13           16:17:01.523765           0.03888214     

The program also prints a table of frequency/wavelength values, and the percent energy in the image at each wavelength. The sharpness is the energy in the 16 - 4 pixel band.

To get the sharpness of a single image, plus plots, use:

> sharpness.py contk1/data1041.a.fits 0 3000

where the last two numbers on the line are the min and max scale values of the image. Note that the plots shown are of the log spectrum and average log spectrum, which are not the actual data used to estimate the sharpness. See the sharpness.py program text for more details. Also note that you must close all the windows to exit the program.

plotsharp.py

Usage: [python] plotsharp.py sharpness-file

To run on a .txt file of sharpness results, use:

> plotsharp.py sharpness.txt
2400 entries
Min time:  contk/data1041.a.fits     2004-07-13           16:17:01.523765           0.03888214     
Max time:  contk/data3540.a.fits     2004-07-13           18:30:32.704249           0.01595912     
Min sharp: contk/data3021.a.fits     2004-07-13           17:56:32.344494           0.00218952     
Max sharp: contk/data1161.a.fits     2004-07-13           16:22:02.456419           0.04652096     

Note that sharpness is estimated on the raw FITS image, full-frame. Therefore, any negative values, noise (readout, Poisson, or bad pixels), and the presence of the slit will affect this result. Also, bad images (usually meaning planet partially off-screen or multiple planet images) have not been removed from this set. Ideally, sharpness should be estimated after the image has been aligned and preprocessed, and should only be performed on a subset of the clouds which has been passed through a 2d Blackman filter. Once alignment and preprocessing have been performed, the sharpness can be recomputed in this manner.

Nevertheless, the results presented here correspond well with visual evaluation. Here are 6 images from best to worst, equally spaced in the sharpness.txt file output for 2400 7/13/07 images:


plotfits.py

Usage: [python] plotfits.py fits-file [min] [max] [header-flag (anything)]

To run on a .fits file, use:

> plotfits.py contk/data1041.a.fits 0 3000 1
SIMPLE  =                    T / DATA IS IN FITS FORMAT                         
BITPIX  =                   16 / 16 BITS TWOS COMPLEMENT INTEGERS               
NAXIS   =                    2 / NUMBER OF AXIS                                 
NAXIS1  =                  512 / PIXELS ON 1st MOST VARYING AXIS                
NAXIS2  =                  512 / PIXELS ON 2nd MOST VARYING AXIS                
DATAMIN =                -4960 / MIN DATA VALUE IN FILE                         
DATAMAX =                20779 / MAX DATA VALUE IN FILE                         
DATAMEAN=                 0.00 / MEAN DATA VALUE IN FILE                        
DIVISOR =                    4 / Normalization value                            
ORIGIN  = 'Institute for Astronomy'                                             
TELESCOP= 'NASA IRTF'                                                           
INSTRUME= 'SpeX, IRTF Guider/Camera'                                            
OBSERVER= 'Eliot Young and Mark Bullock'                                        
OBJECT  = 'Venus'                                                               
COMMENT = 'Morning Venus'                                                       
IRAFNAME= 'data1041.a.fits'                                                     
DSPTMFLE=               'none' / DSPTimingInfo File                             
BEAM    =                  'A' / Object(A) or sky(B)                            
TIME_OBS=    '16:17:01.523765' / UT TIME OF ACQISTION ('hh:mm:ss.ss')           
DATE_OBS=         '2004-07-13' / UT DATE OF ACQUISITION ('yyyy/mm/dd')          
TIME_GPS=       '17:01.594634' / UT TIMESTAMP from GPS ('mm:ss.usec')           
CAMMODE =                    0 / CameraMode is Basic                            
ITIME   =               0.2500 / INTEGRATION TIME IN SECONDS                    
CO_ADDS =                    4 / NUMBER OF INTEGRATIONS                         
CYCLES  =                  400 / Number of cycles                               
NUMARRAY=                    1 / Number of data sub-arrays                      
ARRAY0  =          0,0,512,512 / x,y,wid,hgt of data sub-arrays                 
BEAMPAT =                    0 / Beam Pattern                                   
BBMODE  =                    2 / BB Mode is BBMODE_D16                          
CBMODE  =                    1 / CB Mode is ARC_D                               
NDR     =                    1 / Number of Non-Destructive Reads                
SLOWCNT =                   12 / Number of Slow Counts                          
GRSTCNT =                  800 / Global Reset Count. 1 cnt = 25 nsec            
GPSTMFLG=                    1 / GPS_Time_Flag                                  
BGRESET =          1,2000,1600 / Backgroud Reset's flag, msec, cnt              
RMEX206 =                    0 / RemoveExtra206MuxData. MUX206=0                
TABLE_MS=              235.508 / msec to clock out array                        
CALMIR  =                  Out / menu=0,pos=252000                              
DIT     =                 open / menu=2,pos=160000                              
OSF     =                 Open / menu=0,pos=21333                               
ROTATOR =                90.00 / Rotator_Mechanical_Angle; pos=180000           
POSANGLE=                 0.00 / Sky_Position_Angle                             
SLIT    =               0.3x15 / menu=2,pos=89000                               
GRAT    =             LowRes15 / menu=5,pos=179350                              
GFLT    =                contK / menu=10,pos=448000                             
AFOC    =                50000 / pos=50000,volts= 0.24                          
PLATE_SC=                 0.12 / PlateScale arcsec/pixel                        
FLTZP   =                  0.0 / GFlt Filter Zero Point                         
EXT_COF =                0.000 / GFlt Extinction coefficient                    
TC330A  =  '30.00,74.13,30.00,Med30.00%' / Spectragraph ChA,ChB,SetPt,Heater    
TC330B  =  '30.00,11.89,30.00,Med10.00%' / Guider ChA,ChB,SetPt,Heater          
TC208   = ' 49.43  15.67  74.30  73.81  74.56  77.93  94.25 277.70'             
QTH_LAMP=                  Off / Quart Lamp                                     
INC_LAMP=                  Off / Incandesce Lamp                                
IR_SRC  =                  Off / IR Source                                      
ARG_SRC =                  Off / Argon Lamp                                     
SHUTTER =                Close / Cal Box shutter                                
TPD     = '  04:48:10.56  17:39:54.3 -03:25:29.47 1.511 2000.0 -O'              
RA      =          04:48:10.56 / Right Ascension                                
DEC     =           17:39:54.3 / Declination                                    
HA      =         -03:25:29.47 / Hour Angle                                     
AIRMASS =                1.511 / Air Mass                                       
EPOCH   =               2000.0 / Epoch                                          

To plot an image without showing the header, use:

> plotfits.py contk/data1041.a.fits 0 3000

To plot an image using auto-scale, omit the min and max values, or pass 0 for both:

> plotfits.py contk/data1041.a.fits

Note that you can display several images by using multiple plotfits.py commands, each terminated with the '&' character:

20: ~/Projects/Python/Sharpness > plotfits.py contk/data1041.a.fits &
[1] 283
21: ~/Projects/Python/Sharpness > plotfits.py contk/data1042.a.fits &
[2] 284
22: ~/Projects/Python/Sharpness > plotfits.py contk/data1043.a.fits &
[3] 285

mergesharp.py

Usage: [python] mergesharp.py sharpness-file1 sharpness-file2 [> output-file]

You can combine 2 sharpness .txt files with:

27: ~/Projects/Python/Sharpness > cat sharpness1.txt
contk/data1050.a.fits     2004-07-13                16:17:23.972260           0.04539629     
contk/data1047.a.fits     2004-07-13                16:17:16.547470           0.04383086     
contk/data1048.a.fits     2004-07-13                16:17:19.018101           0.04270498     
contk/data1049.a.fits     2004-07-13                16:17:21.525703           0.04204120     
contk/data1045.a.fits     2004-07-13                16:17:11.548302           0.03994618     
contk/data1042.a.fits     2004-07-13                16:17:03.969390           0.03979242     
contk/data1043.a.fits     2004-07-13                16:17:06.517683           0.03971004     
contk/data1044.a.fits     2004-07-13                16:17:08.984247           0.03949092     
contk/data1046.a.fits     2004-07-13                16:17:13.994263           0.03925742     
contk/data1041.a.fits     2004-07-13                16:17:01.523765           0.03888214     
28: ~/Projects/Python/Sharpness > cat sharpness2.txt
contk/data1051.a.fits     2004-07-13                16:17:26.529575           0.04340547     
contk/data1053.a.fits     2004-07-13                16:17:31.528383           0.04169112     
contk/data1052.a.fits     2004-07-13                16:17:28.970800           0.04164434     
contk/data1058.a.fits     2004-07-13                16:17:43.974738           0.04061172     
contk/data1060.a.fits     2004-07-13                16:17:49.016979           0.03716876     
contk/data1059.a.fits     2004-07-13                16:17:46.548851           0.03468087     
contk/data1057.a.fits     2004-07-13                16:17:41.512497           0.03426781     
contk/data1055.a.fits     2004-07-13                16:17:36.507754           0.03303268     
contk/data1056.a.fits     2004-07-13                16:17:38.970658           0.03005856     
contk/data1054.a.fits     2004-07-13                16:17:33.970351           0.02946241     
29: ~/Projects/Python/Sharpness > mergesharp.py sharpness1.txt sharpness2.txt
contk/data1050.a.fits     2004-07-13           16:17:23.972260           0.04539629     
contk/data1047.a.fits     2004-07-13           16:17:16.547470           0.04383086     
contk/data1051.a.fits     2004-07-13           16:17:26.529575           0.04340547     
contk/data1048.a.fits     2004-07-13           16:17:19.018101           0.04270498     
contk/data1049.a.fits     2004-07-13           16:17:21.525703           0.04204120     
contk/data1053.a.fits     2004-07-13           16:17:31.528383           0.04169112     
contk/data1052.a.fits     2004-07-13           16:17:28.970800           0.04164434     
contk/data1058.a.fits     2004-07-13           16:17:43.974738           0.04061172     
contk/data1045.a.fits     2004-07-13           16:17:11.548302           0.03994618     
contk/data1042.a.fits     2004-07-13           16:17:03.969390           0.03979242     
contk/data1043.a.fits     2004-07-13           16:17:06.517683           0.03971004     
contk/data1044.a.fits     2004-07-13           16:17:08.984247           0.03949092     
contk/data1046.a.fits     2004-07-13           16:17:13.994263           0.03925742     
contk/data1041.a.fits     2004-07-13           16:17:01.523765           0.03888214     
contk/data1060.a.fits     2004-07-13           16:17:49.016979           0.03716876     
contk/data1059.a.fits     2004-07-13           16:17:46.548851           0.03468087     
contk/data1057.a.fits     2004-07-13           16:17:41.512497           0.03426781     
contk/data1055.a.fits     2004-07-13           16:17:36.507754           0.03303268     
contk/data1056.a.fits     2004-07-13           16:17:38.970658           0.03005856     
contk/data1054.a.fits     2004-07-13           16:17:33.970351           0.02946241     

This program reads both files and resorts the combined results. You can also redirect this output to another .txt file. Note that although the sharpness.py program is fairly fast to run, you can get almost double the throughput on a dual-core machine by running two copies of sharpness.py, one in each of two terminal windows. You can then use mergesharp.py to combine the results.

Still To Do?


İSky Coyote 2007