|
| 1 | +import astropy.stats |
| 2 | +import numpy as np |
| 3 | +from scipy import ndimage |
| 4 | +from attr import attributes, attr |
| 5 | + |
| 6 | + |
| 7 | +@attributes |
| 8 | +class IslandParams(object): |
| 9 | + """ |
| 10 | + Data structure for representing source 'islands' |
| 11 | +
|
| 12 | + Initialized with parent image, label index, and peak-pixel value. |
| 13 | +
|
| 14 | + Attributes: |
| 15 | + parent (SourceFindImage): The image the island was detected from. |
| 16 | + label_idx (int): Index of region in label-map of source image. |
| 17 | + peak_val (float): Peak pixel value |
| 18 | + peak_x_idx (int): Peak pixel x-index |
| 19 | + peak_y_idx (int): Peak pixel y-index |
| 20 | + xbar (float): Barycentric centre in x-pixel index |
| 21 | + ybar(float): Barycentric centre in y-pixel index |
| 22 | + """ |
| 23 | + parent = attr() |
| 24 | + label_idx = attr() |
| 25 | + peak_val = attr() |
| 26 | + peak_x_idx = attr(default=None) |
| 27 | + peak_y_idx = attr(default=None) |
| 28 | + xbar = attr(default=None) |
| 29 | + ybar = attr(default=None) |
| 30 | + |
| 31 | + def calculate_params(self): |
| 32 | + """ |
| 33 | + Analyses an 'island' to extract further parameters |
| 34 | + """ |
| 35 | + self.data = np.ma.MaskedArray( |
| 36 | + self.parent.data, |
| 37 | + mask=_label_mask(self.parent.label_map, self.label_idx), |
| 38 | + |
| 39 | + ) |
| 40 | + self.peak_y_idx, self.peak_x_idx = _max_pixel_index(self.data) |
| 41 | + sum = np.ma.sum(self.data) |
| 42 | + self.xbar = np.ma.sum(self.parent.xgrid * self.data) / sum |
| 43 | + self.ybar = np.ma.sum(self.parent.ygrid * self.data) / sum |
| 44 | + |
| 45 | + |
| 46 | +def _label_mask(labels_map, label_num): |
| 47 | + return ~(labels_map == label_num) |
| 48 | + |
| 49 | + |
| 50 | +def _max_pixel_index(masked_image): |
| 51 | + """ |
| 52 | + Returns max pixel index in np array ordering, i.e. (y_max, x_max) |
| 53 | + """ |
| 54 | + return np.unravel_index(np.ma.argmax(masked_image), |
| 55 | + masked_image.shape) |
| 56 | + |
| 57 | + |
| 58 | +class SourceFindImage(object): |
| 59 | + """ |
| 60 | + Data structure for collecting intermediate results from source-extraction. |
| 61 | +
|
| 62 | + This can be useful for verifying / debugging the sourcefinder results, |
| 63 | + and intermediate results can also be reused to save recalculation. |
| 64 | +
|
| 65 | + Args: |
| 66 | + data (array_like): numpy.ndarray or numpy.ma.MaskedArray containing |
| 67 | + image data. |
| 68 | + detection_n_sigma (float): Detection threshold as multiple of RMS |
| 69 | + analysis_n_sigma (float): Analysis threshold as multiple of RMS |
| 70 | + rms_est (float): RMS estimate (may be `None`, in which case RMS is |
| 71 | + estimated from the image data via sigma-clipping). |
| 72 | + """ |
| 73 | + |
| 74 | + def __init__(self, data, detection_n_sigma, analysis_n_sigma, |
| 75 | + rms_est=None): |
| 76 | + self.data = data |
| 77 | + self.detection_n_sigma = detection_n_sigma |
| 78 | + self.analysis_n_sigma = analysis_n_sigma |
| 79 | + |
| 80 | + self.rms_est = rms_est |
| 81 | + if self.rms_est is None: |
| 82 | + self.rms_est = _estimate_rms(self.data) |
| 83 | + self.bg_level = np.ma.median(self.data) |
| 84 | + self.analysis_thresh = self.bg_level + self.analysis_n_sigma * self.rms_est |
| 85 | + self.detection_thresh = self.bg_level + self.detection_n_sigma * self.rms_est |
| 86 | + |
| 87 | + self.ygrid, self.xgrid = np.indices(self.data.shape) |
| 88 | + |
| 89 | + # Label connected regions |
| 90 | + self.label_map, label_maxvals = _label_detection_islands( |
| 91 | + self.data, self.analysis_thresh, self.detection_thresh |
| 92 | + ) |
| 93 | + |
| 94 | + self.islands = [] |
| 95 | + for l_idx, l_maxval in label_maxvals.items(): |
| 96 | + self.islands.append( |
| 97 | + IslandParams(parent=self, label_idx=l_idx, peak_val=l_maxval)) |
| 98 | + for island in self.islands: |
| 99 | + island.calculate_params() |
| 100 | + |
| 101 | + |
| 102 | +def _estimate_rms(image): |
| 103 | + clipped_data = astropy.stats.sigma_clip(image) |
| 104 | + clipped_stddev = np.ma.std(clipped_data) |
| 105 | + return clipped_stddev |
| 106 | + |
| 107 | + |
| 108 | +def _label_detection_islands(data, analysis_thresh, detection_thresh): |
| 109 | + """ |
| 110 | + Find regions which are above analysis_thresh and peak above detection_thresh. |
| 111 | +
|
| 112 | + Args: |
| 113 | + data (array_like): Raw data (may be masked) |
| 114 | + analysis_thresh (float): Analysis ('single connected region') threshold |
| 115 | + detection_thresh (float): Detections threshold |
| 116 | +
|
| 117 | + Returns: |
| 118 | + tuple (array_like, dict): Tuple of `(label_map, valid_labels)`. |
| 119 | + `label_map` is an ndarray containing the pixel labels. |
| 120 | + `valid_labels` is a dict mapping valid label numbers to the maximum |
| 121 | + pixel-value in that label-region. |
| 122 | +
|
| 123 | + """ |
| 124 | + analysis_map = data > analysis_thresh |
| 125 | + label_map, n_labels = ndimage.label(analysis_map) |
| 126 | + |
| 127 | + label_maxvals = ndimage.maximum(data, label_map, |
| 128 | + index=range(1, n_labels + 1)) |
| 129 | + |
| 130 | + valid_label_maxvals = {} |
| 131 | + # Delabel islands that don't meet detection threshold: |
| 132 | + for zero_idx, maxval in enumerate(label_maxvals): |
| 133 | + label = zero_idx + 1 |
| 134 | + if maxval < detection_thresh: |
| 135 | + label_map[label_map == label] = 0. |
| 136 | + else: |
| 137 | + valid_label_maxvals[label] = maxval |
| 138 | + return label_map, valid_label_maxvals |
0 commit comments