Coverage for metadetect / interpolate.py: 5%
102 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-23 22:27 -0700
1"""
2interpolation utils - orig. from beckermr/pizza-cutter
3"""
4import numpy as np
5from scipy.interpolate import CloughTocher2DInterpolator
6import logging
8from numba import njit
10logger = logging.getLogger(__name__)
13@njit
14def _get_nearby_good_pixels(bad_msk, nbad, buff, iso_buff):
15 """
16 get the set of good pixels surrounding bad pixels.
18 Parameters
19 ----------
20 bad_msk : bool array
21 2d array of mask bits. True means it is a bad
22 pixel
23 nbad : int
24 The number of bad pixels.
25 buff : int
26 The size of the good pixel buffer around each bad pixel.
27 iso_buff : int
28 The size of the good pixel test buffer region around each bad pixel. If
29 a given bad pixel doesn't have any good pixels in this region, then it is
30 marked as isolated.
32 Returns
33 -------
34 bad_ind : array-like
35 The 1d indices of the pixels to interp in row*ncol + col.
36 bad_iso : array-like
37 An array of 1 if the bad pixel doesn't have any buffer pixels which are ok, 0
38 otherwise.
39 good_ind : array-like
40 The 1d indices of the good pixels to use in the interp in row*ncol + col.
41 """
43 nrows, ncols = bad_msk.shape
45 ngood = nbad*(2*buff+1)**2
46 bad_ind = np.zeros(ngood, dtype=np.int64)
47 bad_iso = np.zeros(ngood, dtype=np.int64)
48 good_ind = np.zeros(ngood, dtype=np.int64)
49 no_good = 1
51 ibad = 0
52 igood = 0
53 for row in range(nrows):
54 for col in range(ncols):
55 val = bad_msk[row, col]
56 if val:
57 bad_ind[ibad] = row * ncols + col
59 row_start = row - iso_buff
60 row_end = row + iso_buff
61 col_start = col - iso_buff
62 col_end = col + iso_buff
64 if row_start < 0:
65 row_start = 0
66 if row_end > (nrows-1):
67 row_end = nrows-1
68 if col_start < 0:
69 col_start = 0
70 if col_end > (ncols-1):
71 col_end = ncols-1
73 no_good = 1
74 for rc in range(row_start, row_end+1):
75 for cc in range(col_start, col_end+1):
76 tval = bad_msk[rc, cc]
77 if not tval:
78 no_good = 0
80 if no_good == 1:
81 bad_iso[ibad] = 1
83 if buff != iso_buff:
84 row_start = row - buff
85 row_end = row + buff
86 col_start = col - buff
87 col_end = col + buff
89 if row_start < 0:
90 row_start = 0
91 if row_end > (nrows-1):
92 row_end = nrows-1
93 if col_start < 0:
94 col_start = 0
95 if col_end > (ncols-1):
96 col_end = ncols-1
98 for rc in range(row_start, row_end+1):
99 for cc in range(col_start, col_end+1):
100 tval = bad_msk[rc, cc]
101 if not tval:
103 if igood == ngood:
104 raise RuntimeError('good_pix too small')
106 ind = rc*ncols + cc
107 good_ind[igood] = ind
108 igood += 1
110 ibad += 1
112 bad_ind = bad_ind[:ibad]
113 bad_iso = bad_iso[:ibad]
114 good_ind = good_ind[:igood]
116 return bad_ind, bad_iso, good_ind
119def interpolate_image_at_mask(
120 *, image, bad_msk, maxfrac=0.90, buff=4,
121 fill_isolated_with_noise=False, weight=None, rng=None, iso_buff=1,
122):
123 """
124 interpolate the bad pixels in an image
126 Parameters
127 ----------
128 image : array
129 the pixel data
130 bad_msk : array
131 boolean array, True means it is a bad pixel
132 maxfrac : float, optional
133 If the fraction of bad pixels is greater than this,
134 None is returned. Default is 0.90.
135 buff : int, optional
136 The buffer of good pixels around each bad pixel to keep for the interpolant.
137 weight : float, optional
138 The weight to use for generating noise when filling interiors of interpolated
139 regions with noise.
140 fill_isolated_with_noise : bool, optional
141 Fill isolated bad pixels with noise and then interp.
142 rng : np.random.RandomState, optional
143 An RNG to use if we are filling isolated bad pixels with noise.
144 iso_buff : int
145 The size of the good pixel test buffer region around each bad pixel. If
146 a given bad pixel doesn't have any good pixels in this region, then it is
147 marked as isolated.
149 Returns
150 -------
151 interp_image : array-like
152 The interpolated image.
153 """
154 nrows, ncols = image.shape
155 npix = bad_msk.size
157 nbad = bad_msk.sum()
158 bm_frac = nbad/npix
159 if bm_frac <= maxfrac and nbad < npix:
160 interp_image = image.copy()
162 bad_ind, bad_iso, _good_ind = _get_nearby_good_pixels(
163 bad_msk, nbad, buff, iso_buff,
164 )
165 good_ind = np.unique(_good_ind)
166 good_yx = np.unravel_index(good_ind, bad_msk.shape)
167 bad_yx = np.unravel_index(bad_ind, bad_msk.shape)
169 if fill_isolated_with_noise:
170 if rng is None:
171 raise RuntimeError(
172 "You must pass an RNG to fill an image with noise "
173 "when interpolating!"
174 )
176 if weight is None:
177 raise RuntimeError(
178 "You must pass a weight to fill an image with noise "
179 "when interpolating!"
180 )
182 msk = bad_iso == 1
183 if np.any(msk):
184 # mark them as ok pixels
185 bad_msk = bad_msk.copy()
186 bad_msk[bad_yx[0][msk], bad_yx[1][msk]] = False
188 # keep the ones we have to fill
189 noise_fill_yx = (bad_yx[0][msk], bad_yx[1][msk])
191 # recompute the good pixels so that they inlcude the ones we
192 # will noise fill
193 nbad = bad_msk.sum()
194 bad_ind, _, _good_ind = _get_nearby_good_pixels(
195 bad_msk, nbad, buff, iso_buff,
196 )
197 good_ind = np.unique(_good_ind)
198 bad_yx = np.unravel_index(bad_ind, bad_msk.shape)
199 good_yx = np.unravel_index(good_ind, bad_msk.shape)
201 shape = noise_fill_yx[0].shape
202 interp_image[noise_fill_yx[0], noise_fill_yx[1]] = rng.normal(
203 size=shape, scale=1.0/np.sqrt(weight)
204 )
206 good_pix = np.array(good_yx).T
207 bad_pix = np.array(bad_yx).T
208 good_im = interp_image[good_yx[0], good_yx[1]]
209 img_interp = CloughTocher2DInterpolator(
210 good_pix,
211 good_im,
212 fill_value=0.0,
213 )
214 interp_image[bad_msk] = img_interp(bad_pix)
216 return interp_image
218 else:
219 return None