26 import matplotlib.pyplot
as pyplot
38 from .
import algorithmsLib
42 fluxMin = pexConfig.Field(
43 doc =
"specify the minimum psfFlux for good Psf Candidates",
47 check =
lambda x: x >= 0.0,
49 fluxMax = pexConfig.Field(
50 doc =
"specify the maximum psfFlux for good Psf Candidates (ignored if == 0)",
53 check =
lambda x: x >= 0.0,
55 kernelSize = pexConfig.Field(
56 doc =
"size of the Psf kernel to create",
60 borderWidth = pexConfig.Field(
61 doc =
"number of pixels to ignore around the edge of PSF candidate postage stamps",
65 badFlags = pexConfig.ListField(
66 doc =
"List of flags which cause a source to be rejected as bad",
68 default = [
"base_PixelFlags_flag_edge",
69 "base_PixelFlags_flag_interpolatedCenter",
70 "base_PixelFlags_flag_saturatedCenter",
71 "base_PixelFlags_flag_crCenter",
72 "base_PixelFlags_flag_bad",
75 widthMin = pexConfig.Field(
76 doc =
"minimum width to include in histogram",
79 check =
lambda x: x >= 0.0,
81 widthMax = pexConfig.Field(
82 doc =
"maximum width to include in histogram",
85 check =
lambda x: x >= 0.0,
87 sourceFluxField = pexConfig.Field(
88 doc =
"Name of field in Source to use for flux measurement",
90 default =
"base_GaussianFlux_flux",
92 widthStdAllowed = pexConfig.Field(
93 doc =
"Standard deviation of width allowed to be interpreted as good stars",
96 check =
lambda x: x >= 0.0,
98 nSigmaClip = pexConfig.Field(
99 doc =
"Keep objects within this many sigma of cluster 0's median",
102 check =
lambda x: x >= 0.0,
106 pexConfig.Config.validate(self)
108 raise pexConfig.FieldValidationError(
"widthMin (%f) > widthMax (%f)"
112 """A class to handle key strokes with matplotlib displays"""
113 def __init__(self, axes, xs, ys, x, y, frames=[0]):
121 self.
cid = self.axes.figure.canvas.mpl_connect(
'key_press_event', self)
124 if ev.inaxes != self.
axes:
127 if ev.key
and ev.key
in (
"p"):
128 dist = numpy.hypot(self.
xs - ev.xdata, self.
ys - ev.ydata)
129 dist[numpy.where(numpy.isnan(dist))] = 1e30
131 which = numpy.where(dist == min(dist))
136 ds9.pan(x, y, frame=frame)
137 ds9.cmdBuffer.flush()
144 """Return a vector of centerIds based on their distance to the centers"""
145 assert len(centers) > 0
147 minDist = numpy.nan*numpy.ones_like(yvec)
148 clusterId = numpy.empty_like(yvec)
149 clusterId.dtype = int
151 for i, mean
in enumerate(centers):
152 dist = abs(yvec - mean)
154 update = dist == dist
156 update = dist < minDist
158 minDist[update] = dist[update]
159 clusterId[update] = i
163 def _kcenters(yvec, nCluster, useMedian=False, widthStdAllowed=0.15):
164 """A classic k-means algorithm, clustering yvec into nCluster clusters
166 Return the set of centres, and the cluster ID for each of the points
168 If useMedian is true, use the median of the cluster as its centre, rather than
171 Serge Monkewitz points out that there other (maybe smarter) ways of seeding the means:
172 "e.g. why not use the Forgy or random partition initialization methods"
173 however, the approach adopted here seems to work well for the particular sorts of things
174 we're clustering in this application
179 mean0 = sorted(yvec)[len(yvec)//10]
180 delta = mean0 * widthStdAllowed * 2.0
181 centers = mean0 + delta * numpy.arange(nCluster)
183 func = numpy.median
if useMedian
else numpy.mean
185 clusterId = numpy.zeros_like(yvec) - 1
186 clusterId.dtype = int
188 oclusterId = clusterId
191 if numpy.all(clusterId == oclusterId):
194 for i
in range(nCluster):
195 centers[i] = func(yvec[clusterId == i])
197 return centers, clusterId
199 def _improveCluster(yvec, centers, clusterId, nsigma=2.0, nIteration=10, clusterNum=0, widthStdAllowed=0.15):
200 """Improve our estimate of one of the clusters (clusterNum) by sigma-clipping around its median"""
202 nMember =
sum(clusterId == clusterNum)
205 for iter
in range(nIteration):
206 old_nMember = nMember
208 inCluster0 = clusterId == clusterNum
209 yv = yvec[inCluster0]
211 centers[clusterNum] = numpy.median(yv)
212 stdev = numpy.std(yv)
215 stdev_iqr = 0.741*(syv[int(0.75*nMember)] - syv[int(0.25*nMember)])
216 median = syv[int(0.5*nMember)]
218 sd = stdev
if stdev < stdev_iqr
else stdev_iqr
221 print "sigma(iqr) = %.3f, sigma = %.3f" % (stdev_iqr, numpy.std(yv))
222 newCluster0 = abs(yvec - centers[clusterNum]) < nsigma*sd
223 clusterId[numpy.logical_and(inCluster0, newCluster0)] = clusterNum
224 clusterId[numpy.logical_and(inCluster0, numpy.logical_not(newCluster0))] = -1
226 nMember =
sum(clusterId == clusterNum)
228 if nMember == old_nMember
or sd < widthStdAllowed * median:
233 def plot(mag, width, centers, clusterId, marker="o", markersize=2, markeredgewidth=0, ltype='-',
234 magType=
"model", clear=
True):
238 fig = pyplot.figure()
243 axes = fig.add_axes((0.1, 0.1, 0.85, 0.80));
245 xmin = sorted(mag)[int(0.05*len(mag))]
246 xmax = sorted(mag)[int(0.95*len(mag))]
248 axes.set_xlim(-17.5, -13)
249 axes.set_xlim(xmin - 0.1*(xmax - xmin), xmax + 0.1*(xmax - xmin))
252 colors = [
"r", "g", "b", "c", "m", "k",]
253 for k, mean
in enumerate(centers):
255 axes.plot(axes.get_xlim(), (mean, mean,),
"k%s" % ltype)
258 axes.plot(mag[l], width[l], marker, markersize=markersize, markeredgewidth=markeredgewidth,
259 color=colors[k%len(colors)])
261 l = (clusterId == -1)
262 axes.plot(mag[l], width[l], marker, markersize=markersize, markeredgewidth=markeredgewidth,
266 axes.set_xlabel(
"Instrumental %s mag" % magType)
267 axes.set_ylabel(
r"$\sqrt{(I_{xx} + I_{yy})/2}$")
275 A measurePsfTask star selector
277 ConfigClass = ObjectSizeStarSelectorConfig
281 Construct a star selector that looks for a cluster of small objects in a size-magnitude plot.
283 \param[in] config An instance of ObjectSizeStarSelectorConfig
297 """!Return a list of PSF candidates that represent likely stars
299 A list of PSF candidates may be used by a PSF fitter to construct a PSF.
301 \param[in] exposure the exposure containing the sources
302 \param[in] catalog a SourceCatalog containing sources that may be stars
303 \param[in] matches astrometric matches; ignored by this star selector
305 \return psfCandidateList a list of PSF candidates.
314 logger =
pexLogging.Log(pexLogging.getDefaultLog(),
"meas.algorithms.objectSizeStarSelector")
316 detector = exposure.getDetector()
317 pixToTanXYTransform =
None
318 if detector
is not None:
319 tanSys = detector.makeCameraSys(cameraGeom.TAN_PIXELS)
320 pixToTanXYTransform = detector.getTransformMap().get(tanSys)
326 xx = numpy.empty(len(catalog))
327 xy = numpy.empty_like(xx)
328 yy = numpy.empty_like(xx)
329 for i, source
in enumerate(catalog):
330 Ixx, Ixy, Iyy = source.getIxx(), source.getIxy(), source.getIyy()
331 if pixToTanXYTransform:
333 linTransform = pixToTanXYTransform.linearizeForwardTransform(p).getLinear()
334 m = geomEllip.Quadrupole(Ixx, Iyy, Ixy)
335 m.transform(linTransform)
336 Ixx, Iyy, Ixy = m.getIxx(), m.getIyy(), m.getIxy()
338 xx[i], xy[i], yy[i] = Ixx, Ixy, Iyy
340 width = numpy.sqrt(0.5*(xx + yy))
344 bad = reduce(
lambda x, y: numpy.logical_or(x, catalog.get(y)), badFlags,
False)
345 bad = numpy.logical_or(bad, flux < self.
_fluxMin)
346 bad = numpy.logical_or(bad, numpy.logical_not(numpy.isfinite(width)))
347 bad = numpy.logical_or(bad, numpy.logical_not(numpy.isfinite(flux)))
348 bad = numpy.logical_or(bad, width < self.
_widthMin)
349 bad = numpy.logical_or(bad, width > self.
_widthMax)
351 bad = numpy.logical_or(bad, flux > self.
_fluxMax)
352 good = numpy.logical_not(bad)
354 if not numpy.any(good):
355 raise RuntimeError(
"No objects passed our cuts for consideration as psf stars")
357 mag = -2.5*numpy.log10(flux[good])
364 import os, cPickle
as pickle
367 pickleFile = os.path.expanduser(os.path.join(
"~",
"widths-%d.pkl" % _ii))
368 if not os.path.exists(pickleFile):
372 with open(pickleFile,
"wb")
as fd:
373 pickle.dump(mag, fd, -1)
374 pickle.dump(width, fd, -1)
376 centers, clusterId =
_kcenters(width, nCluster=4, useMedian=
True,
379 if display
and plotMagSize
and pyplot:
380 fig =
plot(mag, width, centers, clusterId, magType=self._sourceFluxField.split(
".")[-1].title(),
381 marker=
"+", markersize=3, markeredgewidth=
None, ltype=
':', clear=
True)
388 if display
and plotMagSize
and pyplot:
389 plot(mag, width, centers, clusterId, marker=
"x", markersize=3, markeredgewidth=
None, clear=
False)
391 stellar = (clusterId == 0)
398 if display
and displayExposure:
399 ds9.mtv(exposure.getMaskedImage(), frame=frame, title=
"PSF candidates")
402 eventHandler =
EventHandler(fig.get_axes()[0], mag, width,
403 catalog.getX()[good], catalog.getY()[good], frames=[frame])
411 reply = raw_input(
"continue? [c h(elp) q(uit) p(db)] ").strip()
420 We cluster the points; red are the stellar candidates and the other colours are other clusters.
421 Points labelled + are rejects from the cluster (only for cluster 0).
423 At this prompt, you can continue with almost any key; 'p' enters pdb, and 'h' prints this text
425 If displayExposure is true, you can put the cursor on a point and hit 'p' to see it in ds9.
427 elif reply[0] ==
"p":
428 import pdb; pdb.set_trace()
429 elif reply[0] ==
'q':
434 if display
and displayExposure:
435 mi = exposure.getMaskedImage()
437 with ds9.Buffering():
438 for i, source
in enumerate(catalog):
444 ds9.dot(
"+", source.getX() - mi.getX0(),
445 source.getY() - mi.getY0(), frame=frame, ctype=ctype)
449 with ds9.Buffering():
450 psfCandidateList = []
451 for isStellar, source
in zip(stellar, [s
for g, s
in zip(good, catalog)
if g]):
456 psfCandidate = algorithmsLib.makePsfCandidate(source, exposure)
460 if psfCandidate.getWidth() == 0:
465 im = psfCandidate.getMaskedImage().getImage()
467 if not numpy.isfinite(vmax):
469 psfCandidateList.append(psfCandidate)
471 if display
and displayExposure:
472 ds9.dot(
"o", source.getX() - mi.getX0(), source.getY() - mi.getY0(),
473 size=4, frame=frame, ctype=ds9.CYAN)
474 except Exception
as err:
475 logger.logdebug(
"Failed to make a psfCandidate from source %d: %s" % (source.getId(), err))
477 return psfCandidateList
479 starSelectorRegistry.register(
"objectSize", ObjectSizeStarSelector)
def selectStars
Return a list of PSF candidates that represent likely stars.
A measurePsfTask star selector.
a place to record messages and descriptions of the state of processing.
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
def __init__
Construct a star selector that looks for a cluster of small objects in a size-magnitude plot...
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.