## Python notes 7/20/07

This week I have tried to rewrite some of my C sharpness estimation code from a couple years ago in Python. Click here for a description and examples of how that code worked (section 2: 'selecting images'). One difference between then and now is that then the sharpness was estimated from a cloud-only subdisk extracted from aligned and preprocessed images and then run through a 2d Blackman filter before performing the FFT. Now, the estimate will need to be performed on raw images, full-frame. Hopefully it will work nearly as well.

I'm writing a Python program called 'sharpness.py' that reads either a single FITS image or a folder of images. If a single image, there is also an option to use PyQt to display the image and several intermediate results. Here is an example of reading an image (1041 on 7/13/04), taking the log(abs(fft)^2), and plotting the result wrt radial distance from the center of the transform:

Here is the radius of the transform and a plot of log(abs(fft)^2) x radius:

Here is the same thing plotting every 8th sample, and then a distribution made from 32 bins of the normalized average coefficients for spacial frequencies from 0 - 256 cycles/image (Nyquist limit):

Here is a printout of the distribution:

```Min radius      Max radius      Max wavelength  Min wavelength  Count   Score
----------      ----------      --------------  --------------  -----   -----
0.00000000      8.00000000      512.00000000    64.00000000     197     0.04307483
8.00000000      16.00000000     64.00000000     32.00000000     604     0.04029211
16.00000000     24.00000000     32.00000000     21.33333333     1000    0.03865690
24.00000000     32.00000000     21.33333333     16.00000000     1420    0.03746054
32.00000000     40.00000000     16.00000000     12.80000000     1820    0.03648501
40.00000000     48.00000000     12.80000000     10.66666667     2200    0.03562503
48.00000000     56.00000000     10.66666667     9.14285714      2636    0.03479370
56.00000000     64.00000000     9.14285714      8.00000000      3012    0.03402480
64.00000000     72.00000000     8.00000000      7.11111111      3392    0.03327356
72.00000000     80.00000000     7.11111111      6.40000000      3844    0.03253255
80.00000000     88.00000000     6.40000000      5.81818182      4244    0.03177532
88.00000000     96.00000000     5.81818182      5.33333333      4608    0.03104300
96.00000000     104.00000000    5.33333333      4.92307692      5036    0.03043207
104.00000000    112.00000000    4.92307692      4.57142857      5444    0.02991078
112.00000000    120.00000000    4.57142857      4.26666667      5848    0.02935766
120.00000000    128.00000000    4.26666667      4.00000000      6220    0.02908197
128.00000000    136.00000000    4.00000000      3.76470588      6660    0.02881528
136.00000000    144.00000000    3.76470588      3.55555556      7024    0.02866269
144.00000000    152.00000000    3.55555556      3.36842105      7436    0.02849598
152.00000000    160.00000000    3.36842105      3.20000000      7852    0.02832412
160.00000000    168.00000000    3.20000000      3.04761905      8264    0.02829784
168.00000000    176.00000000    3.04761905      2.90909091      8644    0.02823004
176.00000000    184.00000000    2.90909091      2.78260870      9084    0.02807906
184.00000000    192.00000000    2.78260870      2.66666667      9432    0.02800991
192.00000000    200.00000000    2.66666667      2.56000000      9852    0.02806763
200.00000000    208.00000000    2.56000000      2.46153846      10268   0.02812096
208.00000000    216.00000000    2.46153846      2.37037037      10640   0.02807495
216.00000000    224.00000000    2.37037037      2.28571429      11084   0.02808236
224.00000000    232.00000000    2.28571429      2.20689655      11468   0.02813657
232.00000000    240.00000000    2.20689655      2.13333333      11880   0.02808763
240.00000000    248.00000000    2.13333333      2.06451613      12300   0.02825625
248.00000000    256.00000000    2.06451613      2.00000000      12658   0.02843890
```

Here is the same thing for a blurrier image (3540 on 7/13/04):

```Min radius      Max radius      Max wavelength  Min wavelength  Count   Score
----------      ----------      --------------  --------------  -----   -----
0.00000000      8.00000000      512.00000000    64.00000000     197     0.04398423
8.00000000      16.00000000     64.00000000     32.00000000     604     0.04117782
16.00000000     24.00000000     32.00000000     21.33333333     1000    0.03929861
24.00000000     32.00000000     21.33333333     16.00000000     1420    0.03772854
32.00000000     40.00000000     16.00000000     12.80000000     1820    0.03641676
40.00000000     48.00000000     12.80000000     10.66666667     2200    0.03509932
48.00000000     56.00000000     10.66666667     9.14285714      2636    0.03386101
56.00000000     64.00000000     9.14285714      8.00000000      3012    0.03270312
64.00000000     72.00000000     8.00000000      7.11111111      3392    0.03153483
72.00000000     80.00000000     7.11111111      6.40000000      3844    0.03065055
80.00000000     88.00000000     6.40000000      5.81818182      4244    0.03022037
88.00000000     96.00000000     5.81818182      5.33333333      4608    0.02984994
96.00000000     104.00000000    5.33333333      4.92307692      5036    0.02962684
104.00000000    112.00000000    4.92307692      4.57142857      5444    0.02948804
112.00000000    120.00000000    4.57142857      4.26666667      5848    0.02931253
120.00000000    128.00000000    4.26666667      4.00000000      6220    0.02923771
128.00000000    136.00000000    4.00000000      3.76470588      6660    0.02905015
136.00000000    144.00000000    3.76470588      3.55555556      7024    0.02897174
144.00000000    152.00000000    3.55555556      3.36842105      7436    0.02886230
152.00000000    160.00000000    3.36842105      3.20000000      7852    0.02878383
160.00000000    168.00000000    3.20000000      3.04761905      8264    0.02869115
168.00000000    176.00000000    3.04761905      2.90909091      8644    0.02872059
176.00000000    184.00000000    2.90909091      2.78260870      9084    0.02867005
184.00000000    192.00000000    2.78260870      2.66666667      9432    0.02859581
192.00000000    200.00000000    2.66666667      2.56000000      9852    0.02855485
200.00000000    208.00000000    2.56000000      2.46153846      10268   0.02855297
208.00000000    216.00000000    2.46153846      2.37037037      10640   0.02850872
216.00000000    224.00000000    2.37037037      2.28571429      11084   0.02853153
224.00000000    232.00000000    2.28571429      2.20689655      11468   0.02866671
232.00000000    240.00000000    2.20689655      2.13333333      11880   0.02877451
240.00000000    248.00000000    2.13333333      2.06451613      12300   0.02892382
248.00000000    256.00000000    2.06451613      2.00000000      12658   0.02895105
```

Side-by-side:

Not too different. However, if instead of log(abs(fft)^2) one calculates just abs(fft)^2 and then computes the distribution of the normalized total (not average) coefficients at each radius, there is a difference in the percentage of energy in each frequency band of the image:

```Min radius      Max radius      Max wavelength  Min wavelength  Count   Score
----------      ----------      --------------  --------------  -----   -----
0.00000000      8.00000000      512.00000000    64.00000000     197     0.73503861
8.00000000      16.00000000     64.00000000     32.00000000     604     0.13970486
16.00000000     24.00000000     32.00000000     21.33333333     1000    0.05148516
24.00000000     32.00000000     21.33333333     16.00000000     1420    0.02524016
32.00000000     40.00000000     16.00000000     12.80000000     1820    0.01383987
40.00000000     48.00000000     12.80000000     10.66666667     2200    0.00826908
48.00000000     56.00000000     10.66666667     9.14285714      2636    0.00517407
56.00000000     64.00000000     9.14285714      8.00000000      3012    0.00347442
64.00000000     72.00000000     8.00000000      7.11111111      3392    0.00243512
72.00000000     80.00000000     7.11111111      6.40000000      3844    0.00169278
80.00000000     88.00000000     6.40000000      5.81818182      4244    0.00116907
88.00000000     96.00000000     5.81818182      5.33333333      4608    0.00091015
96.00000000     104.00000000    5.33333333      4.92307692      5036    0.00074590
104.00000000    112.00000000    4.92307692      4.57142857      5444    0.00057999
112.00000000    120.00000000    4.57142857      4.26666667      5848    0.00044773
120.00000000    128.00000000    4.26666667      4.00000000      6220    0.00039566
128.00000000    136.00000000    4.00000000      3.76470588      6660    0.00034138
136.00000000    144.00000000    3.76470588      3.55555556      7024    0.00029665
144.00000000    152.00000000    3.55555556      3.36842105      7436    0.00025688
152.00000000    160.00000000    3.36842105      3.20000000      7852    0.00021864
160.00000000    168.00000000    3.20000000      3.04761905      8264    0.00020143
168.00000000    176.00000000    3.04761905      2.90909091      8644    0.00019260
176.00000000    184.00000000    2.90909091      2.78260870      9084    0.00017961
184.00000000    192.00000000    2.78260870      2.66666667      9432    0.00018821
192.00000000    200.00000000    2.66666667      2.56000000      9852    0.00018359
200.00000000    208.00000000    2.56000000      2.46153846      10268   0.00019149
208.00000000    216.00000000    2.46153846      2.37037037      10640   0.00019281
216.00000000    224.00000000    2.37037037      2.28571429      11084   0.00022811
224.00000000    232.00000000    2.28571429      2.20689655      11468   0.00029883
232.00000000    240.00000000    2.20689655      2.13333333      11880   0.00042021
240.00000000    248.00000000    2.13333333      2.06451613      12300   0.00082809
248.00000000    256.00000000    2.06451613      2.00000000      12658   0.00517882

0.00000000      8.00000000      512.00000000    64.00000000     197     0.81966587
8.00000000      16.00000000     64.00000000     32.00000000     604     0.10535610
16.00000000     24.00000000     32.00000000     21.33333333     1000    0.03608526
24.00000000     32.00000000     21.33333333     16.00000000     1420    0.01417297
32.00000000     40.00000000     16.00000000     12.80000000     1820    0.00664905
40.00000000     48.00000000     12.80000000     10.66666667     2200    0.00330262
48.00000000     56.00000000     10.66666667     9.14285714      2636    0.00181718
56.00000000     64.00000000     9.14285714      8.00000000      3012    0.00112734
64.00000000     72.00000000     8.00000000      7.11111111      3392    0.00084038
72.00000000     80.00000000     7.11111111      6.40000000      3844    0.00062928
80.00000000     88.00000000     6.40000000      5.81818182      4244    0.00053931
88.00000000     96.00000000     5.81818182      5.33333333      4608    0.00040022
96.00000000     104.00000000    5.33333333      4.92307692      5036    0.00029603
104.00000000    112.00000000    4.92307692      4.57142857      5444    0.00025086
112.00000000    120.00000000    4.57142857      4.26666667      5848    0.00023486
120.00000000    128.00000000    4.26666667      4.00000000      6220    0.00018879
128.00000000    136.00000000    4.00000000      3.76470588      6660    0.00017594
136.00000000    144.00000000    3.76470588      3.55555556      7024    0.00016599
144.00000000    152.00000000    3.55555556      3.36842105      7436    0.00016438
152.00000000    160.00000000    3.36842105      3.20000000      7852    0.00017701
160.00000000    168.00000000    3.20000000      3.04761905      8264    0.00018789
168.00000000    176.00000000    3.04761905      2.90909091      8644    0.00019037
176.00000000    184.00000000    2.90909091      2.78260870      9084    0.00019563
184.00000000    192.00000000    2.78260870      2.66666667      9432    0.00022959
192.00000000    200.00000000    2.66666667      2.56000000      9852    0.00025868
200.00000000    208.00000000    2.56000000      2.46153846      10268   0.00023637
208.00000000    216.00000000    2.46153846      2.37037037      10640   0.00022866
216.00000000    224.00000000    2.37037037      2.28571429      11084   0.00023924
224.00000000    232.00000000    2.28571429      2.20689655      11468   0.00028155
232.00000000    240.00000000    2.20689655      2.13333333      11880   0.00039331
240.00000000    248.00000000    2.13333333      2.06451613      12300   0.00075834
248.00000000    256.00000000    2.06451613      2.00000000      12658   0.00456096
```

For example, compare the percent energy in the bands ending at wavelengths of 64 pixels, 32, 16, 8, 4, and 2:

```64:  0.73503861  0.81966587
32:  0.13970486  0.10535610
16:  0.02524016  0.01417297
8:   0.00347442  0.00112734
4:   0.00039566  0.00018879
2:   0.00517882  0.00456096
```

In particular, the sharper image has less energy in the 512-64 pixel band, and more energy in the bands ending at 16, 8, 4, and 2 pixels. If one also sums these bands by 'octave', the results are more striking:

```512-64:  0.73503861  0.81966587
64-32:   0.13970486  0.10535610
32-16:   0.07672532  0.05025823
16-8:    0.03075745  0.01289619
8-4:     0.00837640  0.00337972
4-2:     0.00939736  0.00844389
```

These results suggest that the 16-8 pixel and 8-4 pixel bands might be good indications of relative image sharpness. Obviously, the two images selected above are very different visually. But does this numerical difference also hold for 2 images closer together in apparent sharpness? Here are the reuslts from 2 images closer in apparent sharpness (1230 and 1236 on 7/13/04):

```512-64:  0.74809618  0.75003314
64-32:   0.12485443  0.12263589
32-16:   0.07275120  0.07601406
16-8:    0.03456287  0.03279663
8-4:     0.00993453  0.00890008
4-2:     0.00980079  0.00962019
```

Again, the 16-8 pixel and 8-4 pixel bands are consistent with the visual evaluation. I still need to test this program a bit more, and then set it up to run on a whole folder of files. I expect the calculation to take a few seconds per image, and the program syntax to be:

```> [python] sharpness.py directory [> output-file]
```

I might have this ready for you to test tomorrow or Sunday.

İSky Coyote 2007