Prompted in part by a recent ATLASGAL paper identifying pointing
offsets of about 3" in the BGPS, we revisit the BGPS pointing.
The ATLASGAL team compared the source locations in their catalog to
source locations in the Bolocam catalog by doing "nearest-match"
searches within a 40" radius (see their Figure 8, reproduced here)
Their comparison was over the range -10 < l < 21, so it only covered a
small fraction of the BGPS. It covered 13 fields with independent
pointing solutions, so it's possible that they have actually discovered
an offset only in some of our fields.
The catalog comparison, while interesting, is potentially quite flawed.
There's no guarantee that a source extraction algorithm will measure
source centers accurately when a "source" is just a local overdensity on
a complex background. Using source comparison will also lead to a bias
towards the most source-rich fields, e.g. l000 and l001, so an offset in
one of those fields would drastically affect the catalog offset.
There is a better way to compare pointing between two images that are
expected to be (nearly) identical. It is well-known that
cross-correlation is an effective technique for determining the offsets
between two identical images; I'll briefly summarize some of the
literature here.
Gratadour et al 2005 used a maximum likelihood estimator approach to
determine the "best-fit" offset between two images. This approach is
comparable to Guizar et al (2008), who implemented a fast solution
for (highly) sub-pixel image registration in matlab. In order for the
image registration to be fast, it must operate in fourier space, but to
get sub-pixel registration in fourier space, you need to either pad
(which is slow, and increases memory use drastically) or fit some
functional form around the peak of the cross-correlation image. The
alternative approach implemented by Guizar utilizes the Fourier scaling
theorem to create a zoomed-in image of the peak pixel, which allows
you to get much higher precision for a much lower computational cost.
My innovation is to use the minimum $\chi^2$ estimator to determine the
goodness of fit and therefore error bars on the best-fit offset.
Because the $\chi^2$ value for each offset is simply determined by
sums and multiplication ($\chi^2 = \sum
\frac{x_i-\mu_x}{\sigma_{x_i}^2}$), we can compute each term that
goes in to the $\chi^2$ value independently with fourier transforms,
then create goodness-of-fit contours around the $\chi^2$ minimum. The
statistical requirement for this approach to make sense is that the
errors on the data are gaussian distributed, which is an assumption we
inevitably make for astronomical images. I believe there is also a
requirement that the errors are independent, which may be more difficult
to satisfy, but in the Bolocam images it is satisfied, especially when
multiple independent observations are combined.
Strictly, this approach can only be used when the model data have the
same multiplicative scale as the fitted data. The peak will never be
wrong using this method, but the errors could be incorrect if the model
and data are multiplicatively offset. In principle, this can be
resolved in the future using a Mellin transform [see this site or
this for
a matlab approach and this for an academic paper
on it].
This is the approach I have implemented
at image-registration.rtfd.org. I used simulated test cases to
demonstrate that it is, indeed, effective and accurate. I used this
method to measure the offsets between the v2 data and the v1 data (which
should, in principle, be the same as the offsets between ATLASGAL and
v1) and the v2 vs Herschel Hi-Gal data (which should be zero).
There are actually a few methods implemented in image-registration, and
I compared those. There's a "dft" and a "$\chi^2$" approach, which are
the same (except $\chi^2$ includes realistic errors), a method where a
2D gaussian is fit to the peak of the cross-correlation image, and a
method where a 2nd-order Taylor expansion is performed around the peak
of the cross-correlation image. The latter two are biased. An
example comparison plot looks like this:
The grey dots are catalog centroid positions offsets measured between v1
and v2. The green cross represents the mean and standard deviation of
the grey points. The other data points, as labeled, show the offsets
between the l000 images in v1 and v2 as measured by the method shown.
They all have errorbars plotted, but the errorbars are generally smaller
than the points. The dark spot seen behind the purple point shows the
$\chi^2$ contours out to 8-$\sigma$: the error in the offset is tiny,
sub-arcsecond. In this case, the offsets nearly agree:
l000 catalog dx: -0.31 +/- 0.68 dy: 1.48 +/- 0.64
l000 $\chi^2$ dx: 1.74 +/- 0.03 dy: 1.41 +/- 0.03
This field agreed nicely between v1 and v2.
The comparison to Hi-Gal is perhaps more important; HiGal's pointing is
calibrated based on multi-wavelength observations, some of which include
actual stars. It's a space-based mission, so its pointing is more
stable. And finally, being a space mission, there's a large dedicated
team instead of a single, part-time individual working on the data.
Our offsets from Hi-Gal are pretty small in general, though not
trivially small.
And it turns out, the region that overlaps with ATLASGAL had more
serious pointing errors than the rest of the survey:
(note: both of the above plots are missing L=359 because I forgot it.
Fixing that now...)
The clearest problem field is l=15, with a longitude offset of -6"
between v2 and HiGal.... that's not the question, though. Somehow I've
lost the code that did the v1-HiGal offsets; I'll have to re-write that
tomorrow and let it run...
Update 12/13: I've spent the last couple days clearing up some issues
with the offsets. The error bars should be MUCH smaller than in the
above plots. The means are pretty similar, though.
Short story: the offsets between v1 and Hi-Gal are greater in the
ATLASGAL overlap regions than elsewhere, and in the right general
direction, but not quite as serious as they claimed. In v2, the
ATLASGAL overlap fields and the rest of the survey have the same mean
offsets, and those offsets are small (-0.5" in l, -1" in b).
The problem now is the table. If everything made sense,
(v1-v2)+(v2-higal)+(higal-v1) = 0. But that clearly isn't the case,
which implies an error in the method, which sucks since I'm claiming
this method is superior to alternatives. It's possible that I'm
actually underestimating the errors against Hi-Gal - that can be fixed
relatively easily - but the magnitude of the error won't affect the
centroid measurements. So I probably need to investigate one case very
carefully. l050 is a big problem case, with vector sums >1 pixel in
both directions. That will be my next line of investigation.
The approach will be:
-crop identical fields within l050 from v1, v2, herschel
-perform pointing comparison between them
-check that vector sum < sum of errors
I think - and hope - the trouble is just that I'm using inconsistent
sub-fields to compare Herschel with the two different Bolocam versions,
which is possible because of the way I selected these sub-fields. I'll
do more careful cropping, and probably re-do this analysis
degree-by-degree (with $512^2$ fields, in the hope that it speeds up the
FTs).
Update 12/14:
I've now cropped identical sections in each of the survey, 1 square
degree (512 pixels) each - which is great for speed. As a sidenote, a
little line profiling revealed that the make_cross_plots
code was the slow point in the process, and it is dominated by savefig
calls, not ffts.
I've run a careful examination of self consistency on the l=0 field,
with positive results: the offsets agree to well within the errorbars
(though there is some residual error at the 0.5" level).
However, a similar inspection of l=50 resulted in a major failure:
In this case, the problem is caused by W51 being exactly on the field
edge, leading to huge cross-correlation power at dx=0, but spread over a
large y range. My first thought is to try to downweight the edges,
which can be achieved by "zero-padding" the noise image, but with high
values instead of zero... or alternatively, by setting the edge region
to zero smoothly.
OK, first thought: Bad idea. Increasing the noise along the edges
drastically increases the small-shift autocorrelation for the noise,
which in turn ends up ruling out the small shifts as a fit possibility.
I don't think this really makes sense mathematically, but each step
does. Why would increasing the noise along the edges make the $\chi^2$
fit worse?
This revealed a serious bug in the code that, luckily, only affected
non-uniform error maps. Basically, I had decomposed the $\chi^2$
equation wrong (which is as bad as it sounds).
That total mess has been resolved now. The image edges are downweighted
with a gaussian of 12 pixels, error=100 outside and weight=0 outside
(with weight^2 inside... best to just view the source if you really want
to know the details). The new versions of the above diagrams:
Less than spectacular for l=50, but acceptable given the errors, which
are indeed significantly larger, as you might expect given the lower
total signal in l=50. Now I need to re-run the fits on every field.
OK, cool, last thing accomplished today (...by 8pm): offset comparison
by square degree for all fields. Again, I don't reproduce the magnitude
of the ATLASGAL-measured offsets, but the ATLASGAL fields are, on
average, more offset in longitude (to the negative) than the overall
average.
Curiously, for both v1 and v2, there appears to be a -1.5 deg shift in
latitude from Hi-Gal.
The vector sums are mostly sub-arcsecond, with most exceptions at l>50.
l=59,64, and 65 are particularly bad - but l=50 isn't so bad. So I
should do the "deep" examination of one or two of those fields... who
knows what new errors I'll turn up?
Here's the new v1-ATLASGAL offset plot:
