Question about shape measurement in DM stack


I was trying to get the input semi-major axis, a, and semi-minor axis, b, parameters given the measured base_Sdssshape_xx,yy,xy (and their ext_HSMShape equivalents) for simulated sources, and I was wondering if these moments are already PSF corrected, i.e., they correspond to the intrinsic un-convolved source or, on the contrary, they have the PSF included.

The way I am trying to compute them is the following:

selection = cat_table['deblend_nChild']==0
Mxx = cat_table['base_SdssShape_xx'][selection]
Myy = cat_table['base_SdssShape_yy'][selection]
Mxy = cat_table['base_SdssShape_xy'][selection]
Muu_p_Mvv = Mxx + Myy
Muu_m_Mvv = numpy.sqrt((Mxx - Myy)**2 + 4*Mxy**2)
Muu = 0.5*(Muu_p_Mvv + Muu_m_Mvv)
Mvv = 0.5*(Muu_p_Mvv - Muu_m_Mvv)
theta = 0.5*np.arctan2(2*Mxy, Mxx - Myy)*180/np.pi
a = np.sqrt(Muu)
b = np.sqrt(Mvv)

Where cat_table is the input catalog table.

1 Like

The base_SdssShape algorithm does not correct for the PSF.

I strongly suggest using the lsst.afw.geom.ellipses module for converting between the various ellipticity conventions.

key = lsst.afw.table.QuadrupoleKey(catalog.schema["base_SdssShape"])
axes = [lsst.afw.geom.ellipses.Axes(src.get(key)) for src in catalog]
a = numpy.array([ax.getA() for ax in axes])
b = numpy.array([ax.getB() for ax in axes])
theta = numpy.array([ax.getTheta() for ax in axes])  # radians
1 Like

Thanks a lot @price!

If you want a shape measurement algorithm that does correct for the PSF, try meas_extensions_shapeHSM. The columns of interest would be ext_shapeHSM_HsmShapeRegauss_e1 and ext_shapeHSM_HsmShapeRegauss_e2.

e1, e2 = catalog["ext_shapeHSM_HsmShapeRegauss_e1"], catalog["ext_shapeHSM_HsmShapeRegauss_e2"]
ellipticities = (lsst.afw.geom.ellipses.Distortion(*ellip) for ellip in zip(e1, e2))
ellipses = (lsst.afw.geom.ellipses.SeparableDistortionDeterminantRadius(ellip, 1.0) for ellip in ellipticities if ellip.getE() <= 1)
axes = [lsst.afw.geom.ellipses.Axes(ellip) for ellip in ellipses]

After consultation with @jbosch, it seems that HSM doesn’t provide an unconvolved radius (“Deconvolved radii are really hard,” says Jim) so I’ve just plugged in a unit radius; scale it however you like.

1 Like

If you really need deconvolved sizes, here’s a way. When moments of the PSF are available (and @jbosch tells me they are), you can deconvolve analytically. I worked this out here:

This yields a deconvolved size with the same formula you used above (I’m not sure why deconvolved sizes should be harder than ellipticities from moments). However, moments are measured with some weight function to avoid noise and neighboring objects to dominate. In terms of accuracy, the larger the weight function used for measuring the moments, the better.

To clarify, I agree with @pmelchior on this - it’s just that people have historically spent more time on ellipticities, and the HSM algorithms already include (some attempt at) deconvolved ellipticities. Peter’s recommendation on how to deconvolve sizes is certainly at least as good as anything I can suggest.

Thanks a lot @pmelchior @jbosch. I will take a look at it.