These notes are in reference to a chat I had with Naoki Yasuda of HSC, where we attempted to sort out the notational differences between jointcal and his documents describing meas_mosaic. I’ve attached said documents below, for reference.
A star s in exposure e is measured at position $x^{s,e},y^{s,e}$ (MeasuredStar $X_{\gamma,i}$) on CCD c, which corresponds to position $\alpha^s,\delta^s$ on the sky (FittedStar $F_k$), and which may be associated with a star in a reference catalog (RefStar $R_j$ with associated FittedStar $F_j$).
meas_mosiac
Stars with no reference catalog object are internal stars, while those with reference objects are external stars.
Sky positions $\alpha^s,\delta^s$ are mapped onto a tangent plane $\xi^{s,e},\eta^{s,e}$, taking into account the boresight $A^e,D^e$.
The CCD positions $X_c,Y_c$ and rotations $\theta_c$ are explicitly fit parameters.
The mapping is an nth order polynomial, with kth polynomial coefficients $a_k^e,b_k^e$ for the $\xi,\eta$ parts, respectively.
The fit is globally for the entire focal plane.
jointcal
A CcdImage $\gamma$ is the catalog, initial WCS, telescope pointing, etc. from 1 CCD c and exposure e.
$M_\gamma$ maps from Pixel to Tangent plane: $x,y\stackrel{M}{\rightarrow}\xi,\eta$.
$P_\gamma$ maps from sky to tangent plane (generally taken to be a gnomonic projection): $\alpha,\delta\stackrel{P}{\rightarrow}\xi,\eta$.
The weight $W_{\gamma,i}$ is the measurement weight of $M_\gamma(X_{\gamma,i})$.
RefStars $R_j$ and their associated FittedStars $F_j$ are mapped to the tangent plane via projection $P$ (explicitly gnomonic projection), with weight $W_j$.
CCD position/rotation are subsumed into the fitted mapping, and not dealt with explicitly.
Fits are per-CCD per-visit, not globally over the focal plane, but a full optical distortion could be included in the $P_\gamma$ term.
Parameters
Fitted
Mapping from Pixel to Tangent plane: if mapping $M_\gamma$ is a polynomial, it’s coefficients should be similar to $a_k^e,b_k^e$.
$X_c,Y_c$ and $\theta_c$ are only explicitly fit in meas_mosaic: they are part of $M_\gamma$ in jointcal.
Star positions on the sky: $F_k$ == $\alpha^s,\delta^s$ for internal stars.
Fixed
CCD positions: $X_{\gamma,i} == x^{s,e},y^{s,e}$
Reference star positions on the sky: $F_j = R_j == \alpha^s,\delta^s$ for external stars.
Projected coordinates: $\alpha,\delta\stackrel{P}{\rightarrow}\xi,\eta$ for both $P$ (reference terms) and $P_\gamma$ (measured terms).
Thoughts
Though explicitly fitting the CCD positions and rotations ($X_c,Y_c,\theta_c$) and the telescope boresight ($A^e,D^e$) is a nice feature of meas_mosaic, it would probably be best to include them in the “known” set of Transforms that are corrected for before we solve for the “perturbation”, since they shouldn’t change from one exposure to the next.
Meas_mosiac Taylor expands $a_k^e,b_k^e$ series and fits the full $\chi^2$ matrix, while jointcal Taylor expands and fits the second derivative matrix $d^2\chi^2/d\theta^2$, which is much sparser.
I didn’t quite follow the discussion of CCD positions. meas_mosaic fits for them explicitly, while jointcal doesn’t. What I didn’t pick up was whether jointcal has the same degrees of freedom – if a chip moves, can jointcal handle that?
I propose that we allow the CCDs to move (and rotate), but that we’d add some penalty function (e.g. k \sum (dx^2 + dy^2) where dx and dy are the differential offset between neighbouring chips). Is this what jointcal does? In this parameterisation, as we set k to large values the array becomes rigid. It’s a very flexible way to explore whether we need to free up these degrees of freedom.
If this is implemented via an EM-like iteration between updating the camera model and the astrometry that’s fine by me.
I think @boutigny or @PierreAstier would be best placed to answer these questions, but I’ll give it a shot: jointcal outputs essentially a FITS SIP polynomial (up to 9th order), and a SimplePolyMap applied to each CcdImage can probably take out a small amount of rotation. But there is no actual rotation fit in jointcal, so there’s nowhere to add a penalty function to. I could see (conceptually) how one might add that, but I don’t know how feasible it is in the actual code.
I still think a better choice would be to take out CCD rotation and shift before we do the astrometric shift, since I’d think we’d know those very well, and they shouldn’t change from exposure to exposure (should they?).
EM-like?
Do you mean “optical model” or do you mean “fitted polynomial”? The optical model is not fitted for (in the sense of some sequence of “real” optical distortions), just a polynomial, which I’d like to eventually be a perturbation on top of the “known” camera and optical distortions.
EM-like: EM (Expectation Maximisation) switches between estimating the latent and observable variables, and I wondered if you were doing something analogous here, solving for the CCD positions and other parameters.
I mean I want to ask for the as-observed optical model.
It sounds as if we’re fitting the polynomial correction to the optical model (saved as SIP – is this the external representation of a Chebyshev?) and a SimplePolyMap. Why aren’t we solving for an “affineMap” per chip instead? Would it work as well? Anyway, at the end of running jointcal I want to be able to ask for the as-delivered optical model.
I was assuming that we’ll run jointcal on a chunk of dithered data to determine the updates to the optical distortion model and the chip positions. If the structure of the code/outputs doesn’t permit this, what is our strategy?
The answer depends on what you call the “best-fit optical model”. My practical definition is a set of transformations that allow one to map pixel coordinates in any chip to a common space (typically a projection plane, tangent at a well defined point in the focal plane).
If you agree with this definition, then it is already there. My understanding is that a fit to a set of dithered exposures (and some external astrometric frame) is not going to deliver more than that, unless you have other inputs (e.g. camera blueprints). There could be a well defined way of splitting these pixels -> tangent plane transformations into transformations from CCD space to focal plane, followed by a single transformation from positions in the focal plane to angles on the sky. I have not tried to work it out. It is not clear to me if small acoplanarities of the chips can be accurately described by assuming this two-step (chips-> focal plane -> tangent plane) transformation pipe if one enforces the first step to be made only from shifts and rotations.
if the need is to be able to place objects detected in the same exposure but in different chips on a common coordinate system, then what jointcal (and scamp) fit is sufficient: I have done the exercise of fitting the model for Megacam and SuprimeCam and using the fitted model to astrometrically match a whole exposure at once.
The python scripts in jointcal do not fit the camera model by default.
The switch is simple, but we have not made it, because it is not needed
to just stack exposures.
One last remark: the AstromFit class manipulates an abstract model class, and so any non degenerate model can in principle be fitted using jointcal, as long as somebody is brave enough to code the derivatives of the residuals w.r.t the model parameters. There are two such models in the package at the moment.