# Photometric calibration and astrometric Jacobians

This post is an attempt to summarize, clarify, and expand this Slack conversation.

# Definitions

The (calibrated) images our pipelines operate on are typically either fluence images or surface brightness images:

• Pixels in fluence images have units that are (multiples) of e.g. janskys. While this is what astronomers usually mean when they say “flux”, it’s also what others mean when they say “flux density”, so I’m going to use fluence here to keep definitions precise and avoid that terminological mess.

• Pixels in surface brightness images have units that are (multiples) of e.g. janskys/arcsec^2.

These differ by a factor of the pixel area, which is related to the “astrometric Jacobian” - the determinant of the Jacobian of the WCS for an image (the mapping from that image to the sky) is the pixel area.

# What We Do Now

## ISR

Flat-fielding science images using dome flats - as our `IsrTask` does today - actually yields surface brightness images. They aren’t very good surface-brightness images, because non-uniform illumination, vignetting, chromatic effects, and various other details mean they they aren’t really very flat, and of course the multiplier to get to physical units (i.e. the zeropoint) is completely unknown at this stage. But they’re definitely more surface-brightness than fluence. That’s because most of the variation between pixels seen in a dome flat is actually due to pixel size variations, not sensitivity (i.e. QE) variations, so dividing by a typical dome flat (or average thereof) transforms what was (very roughly) a fluence raw image into (less roughly) a surface brightness image. While there is some small-scale structure in those pixel-size variations (“tree rings”, “ragged gates”), I’m mostly going to ignore that here, and focus on the large-scale structure, which is usually radial, and can be quite significant (~5% in HSC, I think).

## Single-Epoch Measurements

When we make photometric measurements on a surface brightness image (e.g. `calexp`), it’s mostly safe to think of the pixel area being constant over the region of interest for any particular source, and hence what we measure is the fluence of that source divided by the pixel area for that source’s region of interest.

That’s not quite true - one might point out, for instance, that if we’re doing aperture photometry that’s defined in pixel units, the angular size of that aperture is going to be different on different parts of the focal plane. But that’s a much smaller effect, and it’s easy to convince yourself that it goes to zero as the aperture gets larger. And as long as we measure and apply aperture corrections that map small-aperture and point-source photometry to large-aperture photometry, the errors in those will be small as well - at least for point sources; very extended objects are another story.

When we compute a per-CCD mapping from that measured photometry to the fluence measurements in a reference catalog, that mapping automatically includes the pixel area (Jacobian) factor as well as a fluence-fluence mapping (i.e. from scaled electrons to (e.g.) janskys). When we use one scaling factor per CCD to represent that mapping, we’re essentially assuming that we can get away with treating the entire CCD as having a constant pixel area (in addition to a constant mapping from DN to janskys, of course). That’s not a terrible assumption for most CCDs, but it would absolutely be a terrible assumption if we were to attempt to do it over the full focal plane of the large-format cameras we work with. The variation in the pixel area from the center to the outside of the focal plane is usually at least as large as the variation in the mapping from DN to janskys over the focal plane, and sometimes much larger.

## Multi-Epoch Calibration: meas_mosaic

After fitting its astrometric model, meas_mosaic fits a photometric model that explicitly includes this pixel-area term as a fixed model component that’s fully constrained by the astrometric model. Its only free parameters at this stage are those that constrain the fluence-to-fluence mapping, and hence its outputs are explicitly separated into a uncalibrated-surface-brightness-to-uncalibrated-fluence mapping (the Jacobian) and an uncalibrated-fluence-to-calibrated-fluence mapping.

## Multi-Epoch Calibration: jointcal

Jointcal’s photometric model just goes straight from uncalibrated surface-brightness to calibrated fluence, so it naturally outputs a monolithic (in these terms; it does have separate CCD and visit-level components) uncalibrated-surface-brightness-to-calibrated-fluence mapping. This is probably not as well-motivated as the meas_mosaic model, because it doesn’t make use of all the available information, but because the uncalibrated-fluence-to-calibrated-fluence mapping has power on roughly the same spatial scales as the Jacobian, it may do just as well in practice (and we have no evidence that it doesn’t).

## Multi-Epoch Calibration: FGCM

FGCM includes the Jacobian as an explicit known component. It’s actually much more important to do so in this context, because FGCM relies on being able to separate the full calibration mapping into distinct components with more physically-motivated parameterizations, so it can be fit without a reference catalog.

## Background Estimation

Our focal-plane level background estimation code runs on surface-brightness images, and that’s exactly what we want for background estimation: sky backgrounds are locally constant only on surface-brightness images, and hence in a fluence image the sky would look like it varies much more dramatically.

Before we actually warp, we scale the single-epoch images according to the photometric calibration we determined in a previous step. As described above, that photometric calibration is going to be a mapping from uncalibrated surface-brightness to calibrated fluence regardless of how it’s stored (i.e. meas_mosaic two-component `fcr` vs. `Calib`) or how we came by it - though it may be a poor mapping in the case of single-epoch-only calibrations precisely because that single-value-per-CCD `Calib` doesn’t account for the Jacobian well. As a result, our post-scaling, pre-warp single-epoch images should be fluence images.

Because our low-level image resampling code preserves fluence, our warps and coadds will be fluence images, too. That’s exactly what we want if our goal is to use those warps and/or coadds for detection, photometry, and other source measurements. I’m not really sure what we want in order to do background matching; I think we need fluence warps to successfully subtract away the sources, but we may want to then convert back to surface-brightness for the actual background modeling (something for @yusra to think about!).

# What We Ought to Do, Near-Term

We’ve been slowly rolling out a new, spatially-variable `PhotoCalib` object to replace `Calib` and unify the outputs of `meas_mosaic`, `jointcal`, and `FGCM`. One big point of contention has been whether `PhotoCalib` should represent (though until recently we didn’t realize we could frame it like this):

1. A mapping from any kind of image - surface-brightness or fluence - to calibrated fluence. When attached to a surface-brightness image, a `PhotoCalib` would include the Jacobian component of that mapping, but when attached to a fluence image, it would not.
2. A mapping from only from fluence images to calibrated fluence; scaling surface-brightness images or calibrating photometry measured on them would then require using a `SkyWcs` to apply the Jacobian as a separate step.

Some related points that are not in question:

• It is useful to be able to split an uncalibrated-surface-brightness to calibrated-fluence mapping into a Jacobian component and a fluence-fluence mapping.

• The `SkyWcs` class needs APIs to make it easy to extract and manipulate the Jacobian (including getting images of it efficiently), not least so we can divide a monolithic uncalibrated-surface-brightness to calibrated-fluence mapping by the Jacobian to obtain a fluence-to-fluence mapping.

I think the pretty clear answer here is (1): most consumers of a `PhotoCalib` just want to apply a calibration, and don’t want to have to acquire a `SkyWcs`, extract a Jacobian, and apply it after determining whether their photometry or image is a surface-brightness or fluence image. This is also (conveniently) the approach that requires the least work: the `PhotoCalib`s produced by `meas_mosaic` and `jointcal` already include the full uncalibrated-surface-brightness to calibrated-fluence mapping. And it makes `PhotoCalib` more of a direct replacement for the current `Calib`, which is also used to represent the full calibration for both surface-brightness images (e.g. `calexp`) and fluence images (e.g. coadds).

In order to support the (rarer, but still very important) use case of being able to examine just the fluence-to-fluence mapping, I think the right approach is to make sure the `BoundedField` held by `PhotoCalib` can easily be divided by the Jacobian `BoundedField` `SkyWcs` will soon be able to return. As an optimization, we could add a new method to `PhotoCalib` that takes a `SkyWcs` and returns a fluence-to-fluence `BoundedField` that is equivalent to dividing its full-mapping `BoundedField` by that `SkyWcs`’ Jacobian.

We should probably also make it possible to learn whether an image is a surface-brightness image or a fluence image, and I think it makes sense to make that a simple boolean or enum flag on `PhotoCalib`.

(@parejkoj, @swinbank: this is broadly consistent with what we discussed last week, though there are some new small details here, such as that boolean/enum; @parejkoj, can I ask you to make sure this is all ticketed?)

# What We Ought to Do, Longer-Term

I’ve argued above that we’re probably only making very small errors by considering the Jacobian to be constant over the scale of an object when doing photometry on surface-brightness images. But it’s probably easier to do it correctly than rigorously convince ourselves of that, because we can just multiply and divide our single-epoch images by the Jacobian in order to follow three simple rules:

• Always photometer (and detect, centroid, etc) on fluence images.
• Always warp and coadd fluence images.
• Always estimate and subtract backgrounds on surface-brightness images.

Multiplying our images by the Jacobian image should be a pretty fast operation, so I don’t think we should be afraid of doing it when we switch between those operations - the real compute expense is evaluating it in the first place, and I imagine we can transform that compute cost into a less-relevant memory cost by doing it once and caching it throughout single-epoch processing.

This also ties in nicely with how we plan to utilize chromatic flats (i.e. flatten to different SEDs) in the future:

• Always photometer (etc) after flattening to a fiducial SED. (This is not necessary, as we can definitely fix photometry with some other SED in catalog-space later, but it would make the bookkeeping a bit easier).
• Always warp and coadd after flattening to a fiducial SED (this makes the bookkeeping much easier).
• Always estimate and subtract backgrounds after flattening to the SED of the sky for that particular epoch.

As a result, I think we’ll actually just have one (spatially-varying) mapping to multiply and divide our single-epoch images by: the product of the Jacobian and the ratio of the sky-SED flat and the fiducial-SED flat. That single mapping transforms an image suitable for photometry, warping, and coaddition (which is fluence and fiducial-SED flattened) into an image suitable for background estimation (surface brightness, sky-SED flattened). That mapping will probably be a thing we want to attach to `Exposure` somehow in the future (perhaps as a `BoundedField`).

Following those rules would also remove the need to think about the Jacobian in `jointcal`, `FGCM`, or single-CCD photometric calibration, because their input photometry will have been done on a fluence image. And at that point we could retire `PhotoCalib`'s ability to provide calibrations for surface-brightness images, because the only surface-brightness images our pipeline deals with will be transitory and never photometered.

3 Likes

This is an excellent summary, and a lot to digest. I do have a few additional questions/comments that I’m curious where they fit in.

First of all, there’s no mention of star-flats/illumination correction (which in general are not part of the DM stack, but should they be?). Being able separate out this information can be very valuable, but it may be that this is a distinction without a difference for the long-term discussion. Basically, I first was concerned in the near-term plan in that lumps the two corrections together … not that this is incorrect, but that I don’t think we have the technology in the `BoundedField`s currently to describe these small scale variations. But maybe this is something that the `SkyWcs` `BoundedField` will be able to take care of; and certainly, the longer-term plan describes what seems to be the right thing to do. However, while I agree that we’re only making small errors by considering the Jacobian to be constant over the scale of an object, I think it’s the object-to-object errors that are not insignificant that the long-term plan will actually fix (by doing things correctly). So I want to make sure that the longer term plan isn’t lost.

Related, where does the mapping as determined from the collimated beam projector fit into this framework?

Second, I have always been confused by the terminology “flattening to a fiducial SED”. What is the mathematical operation that this describes?

Finally, FGCM (or rather, `fgcmcal` in the stack) I would say “can include” the Jacobian as an explicit known component (and would actually prefer this, but it is configurable).

1 Like

I think being able to distinguish the illumination correction from per-visit grey corrections is definitely desirable, and I know @parejkoj agrees…but neither meas_mosaic or jointcal does this now, so it’s not something I was thinking about when writing this post (I also think it’s at least mostly orthogonal to the Jacobian question).

We don’t in the current implementations of `BoundedField`, but it’s an abstract base class, and I think the interface could support implementations with small-scale variations.

But for very sharp variations - anything that breaks the assumption that the distortion is locally affine on the scale of the PSF - I think you have to do something else entirely, not just for photometry but for astrometry and shapes as well. Doing that really rigorously is really hard (you can’t really define a single commutator for convolution and distortion, so you probably need a sequence of convolutions and distortions interleaved), but @RHL and I have talked about approximately handling ragged gates by doing some simple (e.g. linear) interpolation up front.

So I want to make sure that the longer term plan isn’t lost.

+1

Related, where does the mapping as determined from the collimated beam projector fit into this framework?

I’m probably not the right person to answer this question in detail, but I think it provides both a way to get the illumination correction and a way to get the chromatic dependency of the flats. So once it’s operational we may not need to use stars to get what we call star flats today.

Second, I have always been confused by the terminology “flattening to a fiducial SED”. What is the mathematical operation that this describes?

If you have a flat that is a function of wavelength (as we will with the CBP and auxtel), this operation integrates that over wavelength for that fiducial SED to yield a flat that is not a function of wavelength, and then applies it to a science image. That is obviously the wrong flat for all sources in the images (none of which actually have the fiducial SED), and we’ll need to provide a way to correct that at the catalog level later. But assuming a fiducial SED for the images lets us coadd them and know the coadd was flattened with that fiducial SED rather having to compute the effective SED the coadd was effectively flattened with.

There’s one other term coming from the Jacobian. When we quote aperture magnitudes on CCDs their sizes are given in pixels, and the size of the aperture therefore varies. For the “standard” fluxes like `apInstFlux` this is taken out in the aperture corrections, but for non-standardised aperture fluxes (the curve of growth) we do not attempt to make allowance for this.