Coverage for metadetect / detect.py: 0%
253 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-19 22:39 -0800
1"""
3Code to take a set of images, for example coadds in different bands, run
4detection on them and produce a MEDS like interface to get postage stamps
6TODO factor out the detection code so we could, for example, drop in the LSST
7DM stack detection
9"""
10from __future__ import print_function
11import logging
12import numpy as np
13import esutil as eu
15from ngmix.medsreaders import NGMixMEDS, MultiBandNGMixMEDS
16from meds.util import get_image_info_struct
18from . import defaults
20logger = logging.getLogger(__name__)
22FWHM_FAC = 2*np.sqrt(2*np.log(2))
25class MEDSInterface(NGMixMEDS):
26 """
27 Wrap a full image with a MEDS-like interface
28 """
29 def __init__(self, obs, seg, cat):
30 self.obs = obs
31 self.seg = seg
32 self._image_types = (
33 'image', 'weight', 'seg', 'bmask', 'noise')
34 self._cat = cat
35 self._image_info = get_image_info_struct(1, 20)
37 def has_psf(self):
38 return True
40 def get_psf(self, iobj, icut):
41 """
42 get an image of the psf
43 """
44 return self.obs.psf.image.copy()
46 def get_cutout(self, iobj, icutout, type='image'):
47 """
48 Get a single cutout for the indicated entry
50 parameters
51 ----------
52 iobj:
53 Index of the object
54 icutout:
55 Index of the cutout for this object.
56 type: string, optional
57 Cutout type. Default is 'image'. Allowed
58 values are 'image','weight','seg','bmask'
60 returns
61 -------
62 The cutout image
63 """
65 self._check_indices(iobj, icutout=icutout)
67 if type == 'psf':
68 return self.get_psf(iobj, icutout)
70 im = self._get_type_image(type)
71 dims = im.shape
73 c = self._cat
74 orow = c['orig_start_row'][iobj, icutout]
75 ocol = c['orig_start_col'][iobj, icutout]
76 bsize = c['box_size'][iobj]
78 orow_box, row_box = self._get_clipped_boxes(dims[0], orow, bsize)
79 ocol_box, col_box = self._get_clipped_boxes(dims[1], ocol, bsize)
81 read_im = im[orow_box[0]:orow_box[1],
82 ocol_box[0]:ocol_box[1]]
84 subim = np.zeros((bsize, bsize), dtype=im.dtype)
85 subim += defaults.DEFAULT_IMAGE_VALUES[type]
87 subim[row_box[0]:row_box[1],
88 col_box[0]:col_box[1]] = read_im
90 return subim
92 def _get_type_image(self, type):
93 if type not in self._image_types:
94 raise ValueError("bad cutout type: '%s'" % type)
96 if type == 'image':
97 data = self.obs.image
98 elif type == 'weight':
99 data = self.obs.weight
100 elif type == 'bmask':
101 data = self.obs.bmask
102 elif type == 'noise':
103 data = self.obs.noise
104 elif type == 'seg':
105 data = self.seg
107 return data
109 def _get_clipped_boxes(self, dim, start, bsize):
110 """
111 get clipped boxes for slicing
113 If the box size goes outside the dimensions,
114 trim them back
116 parameters
117 ----------
118 dim: int
119 Dimension of this axis
120 start: int
121 Starting position in the image for this axis
122 bsize: int
123 Size of box
125 returns
126 -------
127 obox, box
129 obox: [start,end]
130 Start and end slice ranges in the original image
131 box: [start,end]
132 Start and end slice ranges in the output image
133 """
134 # slice range in the original image
135 obox = [start, start+bsize]
137 # slice range in the sub image into which we will copy
138 box = [0, bsize]
140 # rows
141 if obox[0] < 0:
142 obox[0] = 0
143 box[0] = 0 - start
145 im_max = dim
146 diff = im_max - obox[1]
147 if diff < 0:
148 obox[1] = im_max
149 box[1] = box[1] + diff
151 return obox, box
154class MEDSifier(object):
155 """
156 very simple MEDS maker for images. Assumes the images are perfectly
157 registered and are sky subtracted, with constant PSF and WCS.
159 The images are added together to make a detection image and sep, the
160 SExtractor wrapper, is run
162 parameters
163 ----------
164 mbobs: ngmix.MultiBandObsList
165 The data with the psf set.
166 sx_config: dict, optional
167 Dict holding sep extract parameters
168 meds_config: dict, optional
169 Dict holding MEDS parameters
170 nodet_flags: int
171 Integer representing bits to mask for detection. The results for all
172 bands are ored together if combining multiple bands into a detection
173 coadd. Default 0
174 """
175 def __init__(self, mbobs, sx_config, meds_config, nodet_flags=0):
176 self.mbobs = mbobs
177 self.nband = len(mbobs)
178 self.nodet_flags = nodet_flags
180 assert len(mbobs[0]) == 1, 'multi-epoch is not supported'
182 self._set_sx_config(sx_config)
183 self._set_meds_config(meds_config)
185 self._set_detim()
186 self._run_sep()
188 def get_multiband_meds(self):
189 """
190 get a MultiBandMEDS object holding all bands
191 """
193 mlist = []
194 for band in range(self.nband):
195 m = self.get_meds(band)
196 mlist.append(m)
198 return MultiBandNGMixMEDS(mlist)
200 def get_meds(self, band):
201 """
202 get fake MEDS interface to the specified band
203 """
204 obslist = self.mbobs[band]
205 obs = obslist[0]
206 return MEDSInterface(
207 obs=obs,
208 seg=self.seg,
209 cat=self.cat,
210 )
212 def _get_image_vars(self):
213 vars = []
214 for obslist in self.mbobs:
215 obs = obslist[0]
216 weight = obs.weight
217 w = np.where(weight > 0)
218 medw = np.median(weight[w])
219 vars.append(1/medw)
220 return np.array(vars)
222 def _set_detim(self):
224 detim = self.mbobs[0][0].image.copy()
225 detim *= 0
227 vars = self._get_image_vars()
228 weights = 1.0/vars
229 wsum = weights.sum()
230 detnoise = np.sqrt(1/wsum)
232 weights /= wsum
234 mask = np.zeros(detim.shape, dtype=bool)
236 for i, obslist in enumerate(self.mbobs):
237 obs = obslist[0]
238 detim += obs.image*weights[i]
239 if obs.has_bmask():
240 mask |= (obs.bmask & self.nodet_flags != 0)
242 self.detim = detim
243 self.detnoise = detnoise
244 self.detmask = mask
246 def _run_sep(self):
247 import sxdes
248 objs, seg = sxdes.run_sep(
249 image=self.detim,
250 noise=self.detnoise,
251 mask=self.detmask,
252 config=self.sx_config,
253 thresh=self.detect_thresh,
254 )
256 ncut = 2 # need this to make sure array
257 new_dt = [
258 ('iso_radius_arcsec', 'f4'),
259 ('id', 'i8'),
260 ('ncutout', 'i4'),
261 ('box_size', 'i4'),
262 ('file_id', 'i8', ncut),
263 ('orig_row', 'f4', ncut),
264 ('orig_col', 'f4', ncut),
265 ('orig_start_row', 'i8', ncut),
266 ('orig_start_col', 'i8', ncut),
267 ('orig_end_row', 'i8', ncut),
268 ('orig_end_col', 'i8', ncut),
269 ('cutout_row', 'f4', ncut),
270 ('cutout_col', 'f4', ncut),
271 ('psf_cutout_row', 'f4', ncut),
272 ('psf_cutout_col', 'f4', ncut),
273 ('dudrow', 'f8', ncut),
274 ('dudcol', 'f8', ncut),
275 ('dvdrow', 'f8', ncut),
276 ('dvdcol', 'f8', ncut),
277 ]
278 cat = eu.numpy_util.add_fields(objs, new_dt)
279 cat['id'] = np.arange(cat.size)
280 cat['ncutout'] = 1
282 jacob = self.mbobs[0][0].jacobian
283 cat['dudrow'][:, 0] = jacob.dudrow
284 cat['dudcol'][:, 0] = jacob.dudcol
285 cat['dvdrow'][:, 0] = jacob.dvdrow
286 cat['dvdcol'][:, 0] = jacob.dvdcol
288 cat['iso_radius_arcsec'] = (
289 cat['iso_radius'] * self.mbobs[0][0].jacobian.get_scale()
290 )
292 if cat.size > 0:
294 box_size = self._get_box_sizes(cat)
296 half_box_size = box_size//2
298 maxrow, maxcol = self.detim.shape
300 cat['box_size'] = box_size
302 cat['orig_row'][:, 0] = cat['y']
303 cat['orig_col'][:, 0] = cat['x']
305 orow = cat['orig_row'][:, 0].astype('i4')
306 ocol = cat['orig_col'][:, 0].astype('i4')
308 ostart_row = orow - half_box_size + 1
309 ostart_col = ocol - half_box_size + 1
310 oend_row = orow + half_box_size + 1 # plus one for slices
311 oend_col = ocol + half_box_size + 1
313 ostart_row.clip(min=0, out=ostart_row)
314 ostart_col.clip(min=0, out=ostart_col)
315 oend_row.clip(max=maxrow, out=oend_row)
316 oend_col.clip(max=maxcol, out=oend_col)
318 # could result in smaller than box_size above
319 cat['orig_start_row'][:, 0] = ostart_row
320 cat['orig_start_col'][:, 0] = ostart_col
321 cat['orig_end_row'][:, 0] = oend_row
322 cat['orig_end_col'][:, 0] = oend_col
323 cat['cutout_row'][:, 0] = \
324 cat['orig_row'][:, 0] - cat['orig_start_row'][:, 0]
325 cat['cutout_col'][:, 0] = \
326 cat['orig_col'][:, 0] - cat['orig_start_col'][:, 0]
328 psf_cen = (np.array(self.mbobs[0][0].psf.image.shape) - 1)/2
329 cat['psf_cutout_row'][:, 0] = psf_cen[0]
330 cat['psf_cutout_col'][:, 0] = psf_cen[1]
332 self.seg = seg
333 self.cat = cat
335 def _get_box_sizes(self, cat):
336 if cat.size == 0:
337 return []
339 mconf = self.meds_config
341 box_type = mconf['box_type']
342 if box_type == 'sigma_size':
343 sigma_size = self._get_sigma_size(cat)
344 row_size = cat['ymax'] - cat['ymin'] + 1
345 col_size = cat['xmax'] - cat['xmin'] + 1
347 # get max of all three
348 box_size = np.vstack((col_size, row_size, sigma_size)).max(axis=0)
350 elif box_type == 'iso_radius':
351 rad_min = mconf['rad_min'] # for box size calculations
352 rad_fac = mconf['rad_fac']
354 box_padding = mconf['box_padding']
355 rad = cat['iso_radius'].clip(min=rad_min)
357 box_rad = rad_fac*rad
359 box_size = (2*box_rad + box_padding).astype('i4')
360 else:
361 raise ValueError('bad box type: "%s"' % box_type)
363 box_size.clip(
364 min=mconf['min_box_size'],
365 max=mconf['max_box_size'],
366 out=box_size,
367 )
369 # now put in fft sizes
370 bins = [0]
371 bins.extend([sze for sze in defaults.ALLOWED_BOX_SIZES
372 if sze >= mconf['min_box_size']
373 and sze <= mconf['max_box_size']])
375 if bins[-1] != mconf['max_box_size']:
376 bins.append(mconf['max_box_size'])
378 bin_inds = np.digitize(box_size, bins, right=True)
379 bins = np.array(bins)
381 box_sizes = bins[bin_inds]
382 logger.debug('box sizes: %s' % box_sizes)
383 logger.debug('minmax: %s %s' % (box_sizes.min(), box_sizes.max()))
385 return box_sizes
387 def _get_sigma_size(self, cat):
388 """
389 "sigma" size, based on flux radius and ellipticity
390 """
391 mconf = self.meds_config
393 ellipticity = 1.0 - cat['b']/cat['a']
394 sigma = cat['flux_radius']*2.0/FWHM_FAC
395 drad = sigma*mconf['sigma_fac']
396 drad = drad*(1.0 + ellipticity)
397 drad = np.ceil(drad)
398 sigma_size = 2*drad.astype('i4') # sigma size is twice the radius
400 return sigma_size
402 def _set_sx_config(self, sx_config_in):
403 import sxdes
405 if sx_config_in is not None:
406 sx_config = {}
407 sx_config.update(sx_config_in)
408 sx_config['filter_kernel'] = np.array(sx_config['filter_kernel'])
410 # this isn't a keyword
411 self.detect_thresh = sx_config.pop('detect_thresh')
412 self.sx_config = sx_config
413 else:
414 self.sx_config = sxdes.SX_CONFIG
415 self.detect_thresh = sxdes.DETECT_THRESH
417 def _set_meds_config(self, meds_config):
418 self.meds_config = meds_config
421class CatalogMEDSifier(MEDSifier):
422 """
423 very simple MEDS maker for images. Assumes the images are perfectly
424 registered and are sky subtracted, with constant PSF and WCS.
426 Uses an input position and box size.
428 parameters
429 ----------
430 mbobs: ngmix.MultiBandObsList
431 The data with the psf set.
432 x : array-like
433 Input x/col in zero-indexed pixels.
434 y : array-like
435 Input y/row in zero-indexed pixels.
436 box_sizes : array-like
437 The box size for each stamp.
438 number : array, optional
439 If not None, the number in the seg map for the object.
440 seg : array, optional
441 If not None, the seg map for the objects.
442 Note the objects are assumed to be in numerical order
443 matching the seg map.
444 """
445 def __init__(self, mbobs, x, y, box_sizes, number=None, seg=None):
446 self.mbobs = mbobs
447 self.nband = len(mbobs)
448 assert len(mbobs[0]) == 1, 'multi-epoch is not supported'
449 self.x = x
450 self.y = y
451 self.box_sizes = box_sizes
453 self._set_cat_and_seg(number, seg)
455 def _set_cat_and_seg(self, number, seg):
456 ncut = 2 # need this to make sure array
457 new_dt = [
458 ('id', 'i8'),
459 ('number', 'i8'),
460 ('ncutout', 'i4'),
461 ('box_size', 'i4'),
462 ('file_id', 'i8', ncut),
463 ('orig_row', 'f4', ncut),
464 ('orig_col', 'f4', ncut),
465 ('orig_start_row', 'i8', ncut),
466 ('orig_start_col', 'i8', ncut),
467 ('orig_end_row', 'i8', ncut),
468 ('orig_end_col', 'i8', ncut),
469 ('cutout_row', 'f4', ncut),
470 ('cutout_col', 'f4', ncut),
471 ('psf_cutout_row', 'f4', ncut),
472 ('psf_cutout_col', 'f4', ncut),
473 ('dudrow', 'f8', ncut),
474 ('dudcol', 'f8', ncut),
475 ('dvdrow', 'f8', ncut),
476 ('dvdcol', 'f8', ncut),
477 ]
478 cat = np.zeros(self.x.shape[0], dtype=new_dt)
479 if number is None:
480 cat['id'] = np.arange(cat.size)
481 cat['number'] = cat['id'] + 1
482 else:
483 cat['id'] = number - 1
484 cat["number"] = number
485 cat['ncutout'] = 1
487 jacob = self.mbobs[0][0].jacobian
488 cat['dudrow'][:, 0] = jacob.dudrow
489 cat['dudcol'][:, 0] = jacob.dudcol
490 cat['dvdrow'][:, 0] = jacob.dvdrow
491 cat['dvdcol'][:, 0] = jacob.dvdcol
493 half_box_size = self.box_sizes//2
495 maxrow, maxcol = self.mbobs[0][0].image.shape
497 cat['box_size'] = self.box_sizes
499 cat['orig_row'][:, 0] = self.y
500 cat['orig_col'][:, 0] = self.x
502 orow = cat['orig_row'][:, 0].astype('i4')
503 ocol = cat['orig_col'][:, 0].astype('i4')
505 ostart_row = orow - half_box_size + 1
506 ostart_col = ocol - half_box_size + 1
507 oend_row = orow + half_box_size + 1 # plus one for slices
508 oend_col = ocol + half_box_size + 1
510 ostart_row.clip(min=0, out=ostart_row)
511 ostart_col.clip(min=0, out=ostart_col)
512 oend_row.clip(max=maxrow, out=oend_row)
513 oend_col.clip(max=maxcol, out=oend_col)
515 # could result in smaller than box_size above
516 cat['orig_start_row'][:, 0] = ostart_row
517 cat['orig_start_col'][:, 0] = ostart_col
518 cat['orig_end_row'][:, 0] = oend_row
519 cat['orig_end_col'][:, 0] = oend_col
520 cat['cutout_row'][:, 0] = \
521 cat['orig_row'][:, 0] - cat['orig_start_row'][:, 0]
522 cat['cutout_col'][:, 0] = \
523 cat['orig_col'][:, 0] - cat['orig_start_col'][:, 0]
525 psf_cen = (np.array(self.mbobs[0][0].psf.image.shape) - 1)/2
526 cat['psf_cutout_row'][:, 0] = psf_cen[0]
527 cat['psf_cutout_col'][:, 0] = psf_cen[1]
529 if seg is None:
530 self.seg = np.zeros_like(self.mbobs[0][0].image, dtype=np.int32)
531 else:
532 self.seg = seg
533 self.cat = cat