-
Notifications
You must be signed in to change notification settings - Fork 68
/
utils.py
1883 lines (1555 loc) · 70.4 KB
/
utils.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#
# Poppy utility functions
#
# These provide various utilities to measure the PSF's properties in certain ways, display it on screen etc.
#
import json
import logging
import os.path
import pickle
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate
import scipy.ndimage
import warnings
from astropy import config
import astropy.units as u
import astropy.io.fits as fits
import poppy
from importlib import reload
from . import accel_math
from .accel_math import xp, _scipy
try:
import pyfftw
except ImportError:
pyfftw = None
_log = logging.getLogger('poppy')
_loaded_fftw_wisdom = False
class FFTWWisdomWarning(RuntimeWarning):
pass
__all__ = ['display_psf', 'display_psf_difference', 'display_ee', 'measure_ee', 'measure_radius_at_ee',
'display_profiles', 'radial_profile',
'measure_radial', 'measure_fwhm', 'measure_sharpness', 'measure_centroid',
'measure_anisotropy', 'spectrum_from_spectral_type',
'specFromSpectralType', 'removePadding']
###########################################################################
#
# Display functions
#
def imshow(image, *args, **kwargs):
"""If needed, fetch array fromm GPU before imshow"""
return plt.imshow( accel_math.ensure_not_on_gpu(image),
*args, **kwargs)
def imshow_with_mouseover(image0, ax=None, *args, **kwargs):
"""Wrapper for matplotlib imshow that displays the value under the
cursor position
Wrapper for pyplot.imshow that sets up a custom mouseover display
formatter so that mouse motions over the image are labeled in the
status bar with pixel numerical value as well as X and Y coords.
"""
if ax is None:
ax = plt.gca()
image = accel_math.ensure_not_on_gpu(image0)
ax.imshow(image, *args, **kwargs)
aximage = ax.images[0].properties()['array']
# need to account for half pixel offset of array coordinates for mouseover relative to pixel center,
# so that the whole pixel from e.g. ( 1.5, 1.5) to (2.5, 2.5) is labeled with the coordinates of pixel (2,2)
# We use the extent and implementation to map back from the data coord to pixel coord
# There is probably an easier way to do this...
imext = ax.images[0].get_extent() # returns [-X, X, -Y, Y]
imsize = ax.images[0].get_size() # returns [sY, sX]g
def report_pixel(x, y):
# map data coords back to pixel coords
# and be sure to clip appropriately to avoid array bounds errors
img_y = np.floor((y - imext[2]) / (imext[3] - imext[2]) * imsize[0])
img_y = int(img_y.clip(0, imsize[0] - 1))
img_x = np.floor((x - imext[0]) / (imext[1] - imext[0]) * imsize[1])
img_x = int(img_x.clip(0, imsize[1] - 1))
return "(%6.3f, %6.3f) %-12.6g" % (x, y, aximage[img_y, img_x])
ax.format_coord = report_pixel
return ax
def display_psf(hdulist_or_filename, ext=0, vmin=1e-7, vmax=1e-1,
scale='log', cmap=None, title=None, imagecrop=None,
adjust_for_oversampling=False, normalize='None',
crosshairs=False, markcentroid=False, colorbar=True,
colorbar_orientation='vertical', pixelscale='PIXELSCL',
ax=None, return_ax=False, interpolation=None, cube_slice=None,
angular_coordinate_unit=u.arcsec):
"""Display nicely a PSF from a given hdulist or filename
This is extensively configurable. In addition to making an attractive display, for
interactive usage this function provides a live display of the pixel value at a
given (x,y) as you mouse around the image.
Parameters
----------
hdulist_or_filename : fits.hdulist or string
FITS file containing image to display.
ext : int
FITS extension. default = 0
vmin, vmax : float
min and max for image display scaling
scale : str
'linear' or 'log', default is log
cmap : matplotlib.cm.Colormap instance or None
Colormap to use. If not given, taken from user's
`poppy.conf.cmap_sequential` (Default: 'gist_heat').
title : string, optional
Set the plot title explicitly.
imagecrop : float
size of region to display (default is whole image)
adjust_for_oversampling : bool
rescale to conserve surface brightness for oversampled PSFs?
(Making this True conserves surface brightness but not
total flux.) Default is False, to conserve total flux.
normalize : string
set to 'peak' to normalize peak intensity =1, or to 'total' to
normalize total flux=1. Default is no normalization.
crosshairs : bool
Draw a crosshairs at the image center (0, 0)? Default: False.
markcentroid : bool
Draw a crosshairs at the image centroid location?
Centroiding is computed with the JWST-standard moving box
algorithm. Default: False.
colorbar : bool
Draw a colorbar on the image?
colorbar_orientation : 'vertical' (default) or 'horizontal'
How should the colorbar be oriented? (Note: Updating a plot and
changing the colorbar orientation is not supported. When replotting
in the same axes, use the same colorbar orientation.)
pixelscale : str or float
if str, interpreted as the FITS keyword name for the pixel scale in arcsec/pixels.
if float, used as the pixelscale directly.
ax : matplotlib.Axes instance
Axes to display into.
return_ax : bool
Return the axes to the caller for later use? (Default: False)
When True, this function returns a matplotlib.Axes instance, or a
tuple of (ax, cb) where the second is the colorbar Axes.
interpolation : string
Interpolation technique for PSF image. Default is None,
meaning it is taken from matplotlib's `image.interpolation`
rcParam.
cube_slice : int or None
if input PSF is a datacube from calc_datacube, which slice
of the cube should be displayed?
angular_coordinate_unit : astropy Unit
Coordinate unit to use for axes display. Default is arcseconds.
"""
if isinstance(hdulist_or_filename, str):
hdulist = fits.open(hdulist_or_filename)
elif isinstance(hdulist_or_filename, fits.HDUList):
hdulist = hdulist_or_filename
else:
raise ValueError("input must be a filename or FITS HDUList object")
# Get a handle on the input image
if hdulist[ext].data.ndim == 2:
im0 = hdulist[ext].data
psf_array_shape = hdulist[ext].data.shape
elif hdulist[ext].data.ndim == 3:
if cube_slice is None:
raise ValueError("To display a PSF datacube, you must set cube_slice=<#>.")
else:
im0 = hdulist[ext].data[cube_slice]
psf_array_shape = hdulist[ext].data.shape[1:]
else:
raise RuntimeError("Unsupported image dimensionality.")
# Normalization
if adjust_for_oversampling:
try:
scalefactor = hdulist[ext].header['OVERSAMP'] ** 2
except KeyError:
_log.error("Could not determine oversampling scale factor; "
"therefore NOT rescaling fluxes.")
scalefactor = 1
im = im0 * scalefactor
else:
# don't change normalization of actual input array, work with a copy!
im = im0.copy()
if normalize.lower() == 'peak':
_log.debug("Displaying image normalized to peak = 1")
im /= im.max()
elif normalize.lower() == 'total':
_log.debug("Displaying image normalized to PSF total = 1")
im /= im.sum()
if scale == 'linear':
norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
else:
norm = matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax)
if isinstance(pixelscale, str):
pixelscale = hdulist[ext].header[pixelscale]
else:
pixelscale = float(pixelscale)
if angular_coordinate_unit != u.arcsec:
coordinate_rescale = (1 * u.arcsec).to_value(angular_coordinate_unit)
pixelscale *= coordinate_rescale
else:
coordinate_rescale = 1
halffov_x = pixelscale * psf_array_shape[1] / 2.0
halffov_y = pixelscale * psf_array_shape[0] / 2.0
unit_label = str(angular_coordinate_unit)
extent = [-halffov_x, halffov_x, -halffov_y, halffov_y]
if cmap is None:
cmap = getattr(matplotlib.cm, poppy.conf.cmap_sequential)
# update and get (or create) image axes
ax = imshow_with_mouseover(
im,
extent=extent,
cmap=cmap,
norm=norm,
ax=ax,
interpolation=interpolation,
origin='lower'
)
ax.set_xlabel(unit_label)
if markcentroid:
_log.info("measuring centroid to mark on plot...")
ceny, cenx = measure_centroid(hdulist, ext=ext, units='arcsec', relativeto='center', boxsize=20, threshold=0.1)
ceny *= coordinate_rescale # if display coordinate unit isn't arcseconds, rescale the centroid accordingly
cenx *= coordinate_rescale
ax.plot(cenx, ceny, 'k+', markersize=15, markeredgewidth=1)
_log.info("centroid: (%f, %f) " % (cenx, ceny))
if imagecrop is not None:
halffov_x = min((imagecrop / 2.0, halffov_x))
halffov_y = min((imagecrop / 2.0, halffov_y))
ax.set_xbound(-halffov_x, halffov_x)
ax.set_ybound(-halffov_y, halffov_y)
if crosshairs:
ax.axhline(0, ls=':', color='k')
ax.axvline(0, ls=':', color='k')
if title is None:
try:
fspec = "%s, %s" % (hdulist[ext].header['INSTRUME'], hdulist[ext].header['FILTER'])
except KeyError:
fspec = str(hdulist_or_filename)
title = "PSF sim for " + fspec
ax.set_title(title)
if colorbar:
if ax.images[0].colorbar is not None:
# Reuse existing colorbar axes (Issue #21)
colorbar_axes = ax.images[0].colorbar.ax
cb = plt.colorbar(
ax.images[0],
ax=ax,
cax=colorbar_axes,
orientation=colorbar_orientation
)
else:
cb = plt.colorbar(
ax.images[0],
ax=ax,
orientation=colorbar_orientation
)
if scale.lower() == 'log':
ticks = np.logspace(np.log10(vmin), np.log10(vmax), int(np.round(np.log10(vmax / vmin) + 1)))
if colorbar_orientation == 'horizontal' and vmax == 1e-1 and vmin == 1e-8:
ticks = [1e-8, 1e-6, 1e-4, 1e-2, 1e-1] # looks better
cb.set_ticks(ticks)
cb.set_ticklabels(ticks)
if normalize.lower() == 'peak':
cb.set_label('Intensity relative to peak pixel')
else:
cb.set_label('Fractional intensity per pixel')
if return_ax:
if colorbar:
return ax, cb
else:
return ax
def display_psf_difference(hdulist_or_filename1=None, hdulist_or_filename2=None,
ext1=0, ext2=0, vmin=None, vmax=1e-4, title=None,
imagecrop=None, adjust_for_oversampling=False,
crosshairs=False, cmap=None, colorbar=True,
colorbar_orientation='vertical',
ax=None, return_ax=False,
normalize=False, normalize_to_second=False):
"""Display nicely the difference of two PSFs from given files
The two files may be FITS files on disk or FITS HDUList objects in memory. The two must have the same
shape and size.
Parameters
----------
hdulist_or_filename1, hdulist_or_filename2 : fits.HDUlist or string
FITS files containing images to difference
ext1, ext2 : int
FITS extension. default = 0
vmin, vmax : float
Image intensity scaling min and max.
title : string, optional
Title for plot.
imagecrop : float
Size of region to display (default is whole image).
adjust_for_oversampling : bool
Rescale to conserve surface brightness for oversampled PSFs?
(Making this True conserves surface brightness but not total flux.)
Default is False, to conserve total flux.
crosshairs : bool
Plot crosshairs over array center?
cmap : matplotlib.cm.Colormap instance or None
Colormap to use. If not given, use standard gray colormap.
colorbar : bool
Draw a colorbar on the image?
colorbar_orientation : 'vertical' (default) or 'horizontal'
How should the colorbar be oriented? (Note: Updating a plot and
changing the colorbar orientation is not supported. When replotting
in the same axes, use the same colorbar orientation.)
print\\_ : bool
Print RMS difference value for the images? (Default: False)
ax : matplotlib.Axes instance
Axes to display into.
return_ax : bool
Return the axes to the caller for later use? (Default: False)
When True, this function returns a matplotlib.Axes instance, or a
tuple of (ax, cb) where the second is the colorbar Axes.
normalize : bool
Display (difference image)/(mean image) instead of just
the difference image. Mutually exclusive to `normalize_to_second`.
(Default: False)
normalize_to_second : bool
Display (difference image)/(second image) instead of just
the difference image. Mutually exclusive to `normalize`.
(Default: False)
"""
if isinstance(hdulist_or_filename1, str):
hdulist1 = fits.open(hdulist_or_filename1)
elif isinstance(hdulist_or_filename1, fits.HDUList):
hdulist1 = hdulist_or_filename1
else:
raise ValueError("input must be a filename or HDUlist")
if isinstance(hdulist_or_filename2, str):
hdulist2 = fits.open(hdulist_or_filename2)
elif isinstance(hdulist_or_filename2, fits.HDUList):
hdulist2 = hdulist_or_filename2
else:
raise ValueError("input must be a filename or HDUlist")
if adjust_for_oversampling:
scalefactor = hdulist1[ext1].header['OVERSAMP'] ** 2
im1 = hdulist1[ext1].data * scalefactor
scalefactor = hdulist2[ext2].header['OVERSAMP'] ** 2
im2 = hdulist1[ext2].data * scalefactor
else:
im1 = hdulist1[ext1].data
im2 = hdulist2[ext2].data
diff_im = im1 - im2
if normalize and not normalize_to_second:
avg_im = (im1 + im2) / 2
diff_im /= avg_im
cbtitle = 'Image difference / average (per pixel)' # Relative intensity difference per pixel'
elif normalize_to_second and not normalize:
diff_im /= im2
cbtitle = 'Image difference / original (per pixel)' # Relative intensity difference per pixel'
else:
cbtitle = 'Intensity difference per pixel'
if vmin is None:
vmin = -vmax
norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
if cmap is None:
cmap = matplotlib.cm.gray
halffov_x = hdulist1[ext1].header['PIXELSCL'] * hdulist1[ext1].data.shape[1] / 2
halffov_y = hdulist1[ext1].header['PIXELSCL'] * hdulist1[ext1].data.shape[0] / 2
extent = [-halffov_x, halffov_x, -halffov_y, halffov_y]
ax = imshow_with_mouseover(diff_im, extent=extent, cmap=cmap, norm=norm, ax=ax,
origin='lower')
if imagecrop is not None:
halffov_x = min((imagecrop / 2, halffov_x))
halffov_y = min((imagecrop / 2, halffov_y))
ax.set_xbound(-halffov_x, halffov_x)
ax.set_ybound(-halffov_y, halffov_y)
if crosshairs:
ax.axhline(0, ls=":", color='k')
ax.axvline(0, ls=":", color='k')
if title is None:
title = "Difference of " + str(hdulist_or_filename1) + "-" + str(hdulist_or_filename2)
ax.set_title(title)
if colorbar:
cb = plt.colorbar(ax.images[0], ax=ax, orientation=colorbar_orientation)
# ticks = np.logspace(np.log10(vmin), np.log10(vmax), int(np.round(np.log10(vmax/vmin)+1)))
# if vmin == 1e-8 and vmax==1e-1:
# ticks = [1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1]
# ticks = [vmin, -0.5*vmax, 0, 0.5*vmax, vmax]
# cb.set_ticks(ticks)
# cb.set_ticklabels(ticks)
cb.set_label(cbtitle)
if return_ax:
if colorbar:
return ax, cb
else:
return ax
def display_ee(hdulist_or_filename=None, ext=0, overplot=False, ax=None, mark_levels=True, levels=None, **kwargs):
""" Display Encircled Energy curve for a PSF
The azimuthally averaged encircled energy is plotted as a function of radius.
Parameters
----------
hdulist_or_filename : fits.HDUlist or string
FITS file containing image to display encircled energy for.
ext : bool
FITS extension to use. Default is 0
overplot : bool
whether to overplot or clear and produce an new plot. Default false
ax : matplotlib Axes instance
axis to plot into. If not provided, current axis will be used.
mark_levels : bool
If set, mark and label on the plots the radii for 50%, 80%, 95% encircled energy.
Default is True
levels : list
if not None and mark_levels is true then this list specifies alternative levels
"""
if isinstance(hdulist_or_filename, str):
hdu_list = fits.open(hdulist_or_filename)
elif isinstance(hdulist_or_filename, fits.HDUList):
hdu_list = hdulist_or_filename
else:
raise ValueError("input must be a filename or HDUlist")
radius, profile, ee = radial_profile(hdu_list, ee=True, ext=ext, **kwargs)
if not overplot:
if ax is None:
plt.clf()
ax = plt.subplot(111)
ax.plot(radius, ee) # , nonposy='clip')
if not overplot:
ax.set_xlabel("Radius [arcsec]")
ax.set_ylabel("Encircled Energy")
if mark_levels:
if levels is None:
markers=[0.5, 0.8, 0.95]
else:
markers=levels
for level in markers:
ee_lev = radius[np.where(ee > level)[0][0]]
yoffset = 0 if level < 0.9 else -0.05
plt.text(ee_lev + 0.1, level + yoffset, 'EE=%2d%% at r=%.3f"' % (level * 100, ee_lev))
def display_profiles(hdulist_or_filename=None, ext=0, overplot=False, title=None, **kwargs):
""" Produce two plots of PSF radial profile and encircled energy
See also the display_ee function.
Parameters
----------
HDUlist_or_filename1,2 : fits.HDUlist or string
FITS files containing image to difference
ext : bool
FITS extension to use. Default is 0
overplot : bool
whether to overplot or clear and produce an new plot. Default false
title : string, optional
Title for plot
"""
if isinstance(hdulist_or_filename, str):
hdu_list = fits.open(hdulist_or_filename, ext=ext)
elif isinstance(hdulist_or_filename, fits.HDUList):
hdu_list = hdulist_or_filename
else:
raise ValueError("input must be a filename or HDUlist")
radius, profile, ee = radial_profile(hdu_list, ee=True, ext=ext, **kwargs)
if title is None:
try:
title = "%s, %s" % (hdu_list[ext].header['INSTRUME'], hdu_list[ext].header['FILTER'])
except KeyError:
title = str(hdulist_or_filename)
if not overplot:
plt.clf()
plt.title(title)
plt.xlabel("Radius [arcsec]")
plt.ylabel("PSF radial profile")
plt.subplot(2, 1, 1)
plt.semilogy(radius, profile)
fwhm = 2 * radius[np.where(profile < profile[0] * 0.5)[0][0]]
plt.text(fwhm, profile[0] * 0.5, 'FWHM = %.3f"' % fwhm)
plt.subplot(2, 1, 2)
# plt.semilogy(radius, ee, nonposy='clip')
plt.plot(radius, ee, color='r') # , nonposy='clip')
if not overplot:
plt.xlabel("Radius [arcsec]")
plt.ylabel("Encircled Energy")
for level in [0.5, 0.8, 0.95]:
if (ee > level).any():
ee_lev = radius[np.where(ee > level)[0][0]]
yoffset = 0 if level < 0.9 else -0.05
plt.text(ee_lev + 0.1, level + yoffset, 'EE=%2d%% at r=%.3f"' % (level * 100, ee_lev))
def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stddev=False, mad=False,
binsize=None, maxradius=None, normalize='None', pa_range=None, slice=0,
custom_function=None):
""" Compute a radial profile of the image.
This computes a discrete radial profile evaluated on the provided binsize. For a version
interpolated onto a continuous curve, see measure_radial().
Code taken pretty much directly from pydatatut.pdf
Parameters
----------
hdulist_or_filename : string
FITS HDUList object or path to a FITS file.
NaN values in the FITS data array are treated as masked and ignored in computing bin statistics.
ext : int
Extension in FITS file
ee : bool
Also return encircled energy (EE) curve in addition to radial profile?
center : tuple of floats
Coordinates (x,y) of PSF center, in pixel units. Default is image center.
binsize : float
size of step for profile. Default is pixel size.
stddev : bool
Compute standard deviation in each radial bin, not average?
mad : bool
Compute median absolute deviation (MAD) in each radial bin.
Cannot be used at same time as stddev; pick one or the other.
normalize : string
set to 'peak' to normalize peak intensity =1, or to 'total' to normalize total flux=1.
Default is no normalization (i.e. retain whatever normalization was used in computing the PSF itself)
maxradius : float, optional
Maximum radius to compute radial profile to. If not set, will be computed for all radii within the image.
pa_range : list of floats, optional
Optional specification for [min, max] position angles to be included in the radial profile.
I.e. calculate that profile only for some wedge, not the full image. Specify the PA in degrees
counterclockwise from +Y axis=0. Note that you can specify ranges across zero using negative numbers,
such as pa_range=[-10,10]. The allowed PA range runs from -180 to 180 degrees.
slice: integer, optional
Slice into a datacube, for use on cubes computed by calc_datacube. Default 0 if a
cube is provided with no slice specified.
custom_function : function
To evaluate an arbitrary function in each radial bin, provide some callable function that
takes a list of pixels and returns one float. For instance custome_function=np.min to compute
the minimum in each radial bin.
Returns
--------
results : tuple
Tuple containing (radius, profile) or (radius, profile, EE) depending on what is requested.
The radius gives the center radius of each bin, while the EE is given inside the whole bin
so you should use (radius+binsize/2) for the radius of the EE curve if you want to be
as precise as possible. The profile will be either the average within each bin, or the
standard or median absolute deviation within each bin if one of those options is selected.
"""
if isinstance(hdulist_or_filename, str):
hdu_list = fits.open(hdulist_or_filename)
elif isinstance(hdulist_or_filename, fits.HDUList):
hdu_list = hdulist_or_filename
else:
raise ValueError("input must be a filename or HDUlist")
if hdu_list[ext].header['NAXIS'] == 3:
# data cube, so pick out just one slice
image = hdu_list[ext].data[slice].copy() # don't change normalization of actual input array, work with a copy!
else:
image = hdu_list[ext].data.copy() # don't change normalization of actual input array, work with a copy!
if normalize.lower() == 'peak':
_log.debug("Calculating profile with PSF normalized to peak = 1")
image /= image.max()
elif normalize.lower() == 'total':
_log.debug("Calculating profile with PSF normalized to total = 1")
image /= image.sum()
pixelscale = hdu_list[ext].header['PIXELSCL']
if binsize is None:
binsize = pixelscale
y, x = np.indices(image.shape, dtype=float)
if center is None:
# get exact center of image
# center = (image.shape[1]/2, image.shape[0]/2)
center = tuple((a - 1) / 2.0 for a in image.shape[::-1])
x -= center[0]
y -= center[1]
r = np.sqrt(x ** 2 + y ** 2) * pixelscale / binsize # radius in bin size steps
if pa_range is None:
# Use full image
ind = np.argsort(r.flat)
sr = r.flat[ind] # sorted r
sim = image.flat[ind] # sorted image
else:
# Apply the PA range restriction
pa = np.rad2deg(np.arctan2(-x, y)) # Note the (-x,y) convention is needed for astronomical PA convention
mask = (pa >= pa_range[0]) & (pa <= pa_range[1])
ind = np.argsort(r[mask].flat)
sr = r[mask].flat[ind]
sim = image[mask].flat[ind]
ri = sr.astype(int) # sorted r as int
deltar = ri[1:] - ri[:-1] # assume all radii represented (more work if not)
rind = np.where(deltar)[0]
nr = rind[1:] - rind[:-1] # number in radius bin
csim = np.nan_to_num(sim).cumsum(dtype=float) # cumulative sum to figure out sums for each bin
# np.nancumsum is implemented in >1.12
tbin = csim[rind[1:]] - csim[rind[:-1]] # sum for image values in radius bins
radialprofile = tbin / nr
# pre-pend the initial element that the above code misses.
radialprofile2 = np.empty(len(radialprofile) + 1)
if rind[0] != 0:
radialprofile2[0] = csim[rind[0]] / (
rind[0] + 1) # if there are multiple elements in the center bin, average them
else:
radialprofile2[0] = csim[0] # otherwise if there's just one then just take it.
radialprofile2[1:] = radialprofile
# Compute radius values corresponding to the measured points in the radial profile.
# including handling the case where the innermost pixel may be more
# than one pixel from the center. This can happen if pa_range is not None, since for
# small ranges < 45 deg or so the innermost pixel that's valid in the mask may be
# more than one pixel from the center. It can also happen if we are computing a
# radial profile centered on an offset source outside of the FOV.
rr = np.arange(ri.min(), ri.min() + len(
radialprofile2)) * binsize + binsize * 0.5 # these should be centered in the bins, so add a half.
if maxradius is not None:
crop = rr < maxradius
rr = rr[crop]
radialprofile2 = radialprofile2[crop]
if stddev:
# Compute standard deviation in each radial bin
stddevs = np.zeros_like(radialprofile2)
r_pix = r * binsize
for i, radius in enumerate(rr):
if i == 0:
wg = np.where(r < radius + binsize / 2)
else:
wg = np.where((r_pix >= (radius - binsize / 2)) & (r_pix < (radius + binsize / 2)))
# wg = np.where( (r >= rr[i-1]) & (r <rr[i] )))
stddevs[i] = np.nanstd(image[wg])
return rr, stddevs
elif mad:
# Compute median absolute deviation in each radial bin
mads = np.zeros_like(radialprofile2)
r_pix = r * binsize
for i, radius in enumerate(rr):
if i == 0:
wg = np.where(r < radius + binsize / 2)
else:
wg = np.where((r_pix >= (radius - binsize / 2)) & (r_pix < (radius + binsize / 2)))
mads[i] = np.nanmedian(np.absolute(image[wg]-np.nanmedian(image[wg])))
return rr, mads
elif custom_function is not None:
# Compute some custom function in each radial bin
results = np.zeros_like(radialprofile2)
r_pix = r * binsize
for i, radius in enumerate(rr):
if i == 0:
wg = np.where(r < radius + binsize / 2)
else:
wg = np.where((r_pix >= (radius - binsize / 2)) & (r_pix < (radius + binsize / 2)))
if len(wg[0])==0: # Zero elements in this bin
results[i] = np.nan
else:
results[i] = custom_function(image[wg])
return rr, results
elif not ee:
# (Default behavior) Compute average in each radial bin
return rr, radialprofile2
else:
# also return the cumulative sum within each radial bin, i.e. the encircled energy
ee = csim[rind]
return rr, radialprofile2, ee
###########################################################################
#
# PSF evaluation functions
#
def measure_ee(hdulist_or_filename=None, ext=0, center=None, binsize=None, normalize='None'):
""" measure encircled energy vs radius and return as an interpolator
Returns a function object which when called returns the Encircled Energy inside a given radius,
for any arbitrary desired radius smaller than the image size.
Parameters
----------
hdulist_or_filename : string
Either a fits.HDUList object or a filename of a FITS file on disk
ext : int
Extension in that FITS file
center : tuple of floats
Coordinates (x,y) of PSF center. Default is image center.
binsize:
size of step for profile. Default is pixel size.
normalize : string
set to 'peak' to normalize peak intensity =1, or to 'total' to normalize total flux=1.
Default is no normalization (i.e. retain whatever normalization was used in computing the PSF itself)
Returns
--------
encircled_energy: function
A function which will return the encircled energy interpolated to any desired radius.
Examples
--------
>>> ee = measure_ee("someimage.fits") # doctest: +SKIP
>>> print("The EE at 0.5 arcsec is ", ee(0.5)) # doctest: +SKIP
"""
rr, radialprofile2, ee = radial_profile(hdulist_or_filename, ext, ee=True, center=center, binsize=binsize,
normalize=normalize)
# append the zero at the center
rr_ee = rr + (rr[1] - rr[0]) / 2.0 # add half a binsize to this, because the ee is measured inside the
# outer edge of each annulus.
rr0 = np.concatenate(([0], rr_ee))
ee0 = np.concatenate(([0], ee))
ee_fn = scipy.interpolate.interp1d(rr0, ee0, kind='cubic', bounds_error=False)
return ee_fn
def measure_radius_at_ee(hdulist_or_filename=None, ext=0, center=None, binsize=None, normalize='None'):
""" measure encircled energy vs radius and return as an interpolator
Returns a function object which when called returns the radius for a given Encircled Energy. This is the
inverse function of measure_ee
Parameters
----------
hdulist_or_filename : string
Either a fits.HDUList object or a filename of a FITS file on disk
ext : int
Extension in that FITS file
center : tuple of floats
Coordinates (x,y) of PSF center. Default is image center.
binsize:
size of step for profile. Default is pixel size.
normalize : string
set to 'peak' to normalize peak intensity =1, or to 'total' to normalize total flux=1.
Default is no normalization (i.e. retain whatever normalization was used in computing the PSF itself)
Returns
--------
radius: function
A function which will return the radius of a desired encircled energy.
Examples
--------
>>> ee = measure_radius_at_ee("someimage.fits") # doctest: +SKIP
>>> print("The EE is 50% at {} arcsec".format(ee(0.5))) # doctest: +SKIP
"""
rr, radialprofile2, ee = radial_profile(hdulist_or_filename, ext, ee=True, center=center, binsize=binsize,
normalize=normalize)
# append the zero at the center
rr_ee = rr + (rr[1] - rr[0]) / 2.0 # add half a binsize to this, because the EE is measured inside the
# outer edge of each annulus.
rr0 = np.concatenate(([0], rr_ee))
ee0 = np.concatenate(([0], ee))
radius_at_ee_fn = scipy.interpolate.interp1d(ee0, rr0, kind='cubic', bounds_error=False)
return radius_at_ee_fn
def measure_radial(hdulist_or_filename=None, ext=0, center=None, binsize=None):
""" measure azimuthally averaged radial profile of a PSF.
Returns a function object which when called returns the mean value at a given radius.
Parameters
----------
hdulist_or_filename : string
what it sounds like.
ext : int
Extension in FITS file
center : tuple of floats
Coordinates (x,y) of PSF center. Default is image center.
binsize:
size of step for profile. Default is pixel size.
Returns
--------
radial_profile: function
A function which will return the mean PSF value at any desired radius.
Examples
--------
>>> rp = measure_radial("someimage.fits") # doctest: +SKIP
>>> radius = np.linspace(0, 5.0, 100)
>>> plt.plot(radius, rp(radius), label="PSF") # doctest: +SKIP
"""
rr, radialprofile, ee = radial_profile(hdulist_or_filename, ext, ee=True, center=center, binsize=binsize)
radial_fn = scipy.interpolate.interp1d(rr, radialprofile, kind='cubic', bounds_error=False)
return radial_fn
def measure_fwhm(hdulist_or_filename, ext=0, center=None, plot=False, threshold=0.1):
""" Improved version of measuring FWHM, without any binning of image data.
Method: Pick out the image pixels which are above some threshold relative to the
peak intensity, then fit a Gaussian to those. Infer the FWHM based on the width of
the Gaussian.
Parameters
----------
hdulist_or_filename : string
what it sounds like.
ext : int
Extension in FITS file
center : tuple of floats
Coordinates (x,y) of PSF center, in pixel units. Default is image center.
threshold : float
Fraction relative to the peak pixel that is used to select the bright peak pixels
used in fitting the Gaussian. Default is 0.1, i.e. pixels brighter that 0.1 of
the maximum will be included. This is chosen semi-arbitrarily to include most of
the peak but exclude the first Airy ring for typical cases.
plot : bool
Display a diagnostic plot.
Returns
-------
fwhm : float
FWHM in arcseconds
"""
from astropy.modeling import models, fitting
if isinstance(hdulist_or_filename, str):
hdulist = fits.open(hdulist_or_filename)
elif isinstance(hdulist_or_filename, fits.HDUList):
hdulist = hdulist_or_filename
else:
raise ValueError("input must be a filename or HDUlist")
image = hdulist[ext].data.copy() # don't change normalization of actual input array; work with a copy
image = image / image.max() # Normalize the copy to peak=1
pixelscale = hdulist[ext].header['PIXELSCL']
_log.debug("Pixelscale is {} arcsec/pix.".format(pixelscale, ))
# Prepare array r with radius in arcseconds
y, x = np.indices(image.shape, dtype=float)
if center is None:
# get exact center of image
center = tuple((a - 1) / 2.0 for a in image.shape[::-1])
_log.debug("Using PSF center = {}".format(center))
x -= center[0]
y -= center[1]
r = np.sqrt(x ** 2 + y ** 2) * pixelscale # radius in arcseconds
# Select pixels above that threshold
wpeak = np.where(image > threshold) # note, image is normalized to peak=1 above
_log.debug("Using {} pixels above {} of peak".format(len(wpeak[0]), threshold))
rpeak = r[wpeak]
impeak = image[wpeak]
# Determine initial guess for Gaussian parameters
if 'DIFFLMT' in hdulist[ext].header:
std_guess = hdulist[ext].header['DIFFLMT'] / 2.354
else:
std_guess = measure_fwhm_radprof(hdulist, ext=ext, center=center, nowarn=True) / 2.354
_log.debug("Initial guess Gaussian sigma= {} arcsec".format(std_guess))
# Determine best fit Gaussian parameters
g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=std_guess)
g_init.mean.fixed = True
fit_g = fitting.LevMarLSQFitter()
g = fit_g(g_init, rpeak, impeak)
_log.debug("Fit results for Gaussian: {}, {}".format(g.amplitude, g.stddev))
# Convert from the fit result sigma parameter to FWHM.
# note, astropy fitting doesn't constrain the stddev to be positive for some reason.
# so take abs value here.
fwhm = 2 * np.sqrt(2 * np.log(2)) * np.abs(g.stddev)
if plot:
plt.loglog(rpeak, impeak, linestyle='none', marker='o', alpha=0.5)
rmin = rpeak[rpeak != 0].min()
plotr = np.linspace(rmin, rpeak.max(), 30)
plt.plot(plotr, g(plotr))
plt.xlabel("Radius [arcsec]")
plt.ylabel("Intensity relative to peak")
plt.axhline(0.5, ls=":")
plt.axvline(fwhm / 2, ls=':')
plt.text(0.1, 0.2, 'FWHM={:.4f} arcsec'.format(fwhm), transform=plt.gca().transAxes, )
plt.gca().set_ylim(threshold * .5, 2)
return fwhm
def measure_fwhm_radprof(HDUlist_or_filename=None, ext=0, center=None, level=0.5, nowarn=False):
""" Measure FWHM by interpolation of the radial profile.
This version is old/deprecated; see the new measure_fwhm instead.
However, this function is kept, for now, to provide a robust, simple backup
method which can be used to determine the initial guess for the model-fitting
approach in the newer measure_fwhm function.
This measures the full width at half maximum for the supplied PSF,
or optionally the full width at some other fraction of max.
Parameters
----------
HDUlist_or_filename, ext : string, int
Same as above
center : tuple
center to compute around. Default is image center.
level : float
Fraction of max to compute for; default is 0.5 for Half Max.
You can also measure widths at other levels e.g. FW at 10% max
by setting level=0.1
nowarn : bool
Set this to suppress the warning display that this function is deprecated.
But you probably shouldn't; only use this if you know what you're doing.
Returns
-------
fwhm : float
FWHM in arcseconds
"""
if not nowarn:
import warnings
warnings.warn("measure_fwhm_radprof uses a deprecated, older algorithm. "