Revision 1, 2000 Feb 7
Revision 2, 2000 Mar 1
Revsion 3, 2000 Apr 28
Photometric Precision
Measuring Positions
Error Estimation
Effects of coma on "Star"
Conclusions
What Does "Star" Actually Do?
Aopt = sum(b_i p_i)
Where b_i = Neff a_i
Neff = 1/sum(a_j^2)
and the optimum variance is Neff v where v is the variance for a
background pixel. The assumption that the source is weak involves
neglecting the additional pixel noise of the source.
This is the optimum solution. Any other solution, say
A2 = sum(c_i p_i)
has a higher variance var(A2) = v sum(c_j^2)
= v(Neff + sum((c_j - b_j)^2))
For a concrete example, I have measured the PSF for a particular MK III image (G0493946) and get Neff = 58. "Star" uses
c_i = 1 if a_i > some value
c_i = 0 otherwise, giving about 30 non-zero points in the 7x7 PSF.
I don't know exactly how many as "Star" doesn't
record the actual number - it ought to!
This is not an unbiased estimate and must be scaled by
1/sum(a_i c_i) = 2.76 - that's the measured Aopt/Astar - in order to give sum(a_i c_i) = 1
Yes! Star is already missing 60% of the flux for the MK III. Maybe the MK IV, where we know we have a problem because of the obvious coma, isn't actually much worse than the MK III - just different.
In this case, var(Aopt) = 58v
var(Astar) = 30*2.76^2v = 228v approximately.
The "Star" standard deviation would thus be expected to be about
double the optimum. This is the penalty one pays for the
simpler aperture method.
[It seems that I picked some images suffering partial obscuration by tree branches. The following results therefore need to be reworked with better data!]
The above theory claims that the precision for an optimum least squares fit should be about twice as good as for the aperture fit used by "Star". Let's see if it is!
I have processed 5 selected V-band image from CD5. These were chosen for their large degree of overlap so the same star is measured in more or less the same place on each image. I compared magnitudes only for those stars found on every image. The images had from 561 to 1132 stars found by Star (20 sigma limit) of which only 373 made it to the final analysis. 68 of those were Tycho stars.
Some stars were presumably below Star's limit on some images.
Some bright stars were saturated on one or more images.
Some non-saturated stars were lost in the process of removing
the saturated ones.
Some stars were close enough to other stars that the identification
was clearly in error on one or more images.
Some may even be variable!
I then subtracted a mean magnitude correction for each image: this is a more reliable process than the method of matching the Tycho magnitudes which I used elsewhere, since exactly the same stars are used for each image and there are lots of them. The correction was quite similar to that found previously using the Tycho magnitudes. Each star gives an estimate of the internal standard deviation of the magnitude measurement. Allowing for the degrees of freedom used up in the processing, the results are:
Standard deviations Star Least Squares Fit All 373 0.071 0.053 magnitudes Tycho stars (68) 0.037 0.034 Brightest 50% 0.051 0.042
The improvements are statistically significant. Not quite the hoped for factor of two but useful all the same.
Stone (Ronald C. Stone. A Comparison of Digital Centering Algorithms. 1989. AJ 97 pp1227-1237) promotes one-dimensional fitting methods such as the median fit used by Star. He quotes the conclusions of previous work as
"two-dimensional centering methods are much slower than the processing of marginal distributions and only slightly more accurate."It should be noted that Stone's review was in the context of real-time data reduction using hardware much less powerful than the machines most of us now use to store our e-mail. I think that as a result of his overriding requirement for speed, the accuracy advantages of 2-D fitting have been downplayed. The hardware constraint has now been removed and we need all the accuracy we can get. I would prefer to turn the conclusion round to read:
"two-dimensional centering methods are somewhat more accurate than one-dimensional methods using marginal distributions. When best accuracy is required, they are worth using in spite of the much heavier computational load."
Stone also points out that median methods (such as those used by Star) give increased position errors in cases where the PSF is narrow (such as the MK IV, at least near the centre of the image). He suggests a minimum FWHM of 4 pixels for the median method while his gaussian fit works down to a 2-pixel FWHM.
standard deviations in degrees.
File X1 Star X1 Fit X2 Star X2 Fit df H3R1438.882 0.000235 0.000118 0.000264 0.000156 90 H3R1439.880 0.000272 0.000118 0.000292 0.000162 99 H3R1440.875 0.000247 0.000122 0.000306 0.000151 97 H3R1442.871 0.000269 0.000128 0.000313 0.000166 88 H3R1446.853 0.000311 0.000142 0.000361 0.000180 86
or in mas (for those unable to multiply by 3,600,000 in their heads)
File X1 Star X1 Fit X2 Star X2 Fit H3R1438.882 846 425 951 572 H3R1439.880 989 425 1052 579 H3R1440.875 889 436 1101 543 H3R1442.871 968 461 1127 595 H3R1446.853 1120 508 1300 648
Yes - that's a significant improvement using the 2-D fit. Pretty close to the expected factor of two.
My "toy" model predicts an error in source amplitude of var(Aopt) = Neff v. This applies only in the limit of a weak source. For the brighter sources actually measured, one must allow for the photon noise of the source itself. Making this allowance, estimating the variance, v, from the residuals of the fit and using the fitted PSF, this gives an estimated standard deviation around 0.01 magnitudes for a typical Tycho source, varying rather slowly with magnitude. This is far less than the observed scatter and quite close to my original estimate. The much larger estimate of around 0.06 magnitudes that I posted recently was the result of a plain old-fashioned mistake that I introduced at the same time I put in the code to calculate photon noise properly. Allowing for the higher variance of the "Star" method, my estimates of scatter agree well with the estimates produced by "Star"; both sets of estimates are less than their respective observed scatters by a similar factor.
Similarly, I arrived (with no great confidence, but the result looks plausible!) at an estimate of the position error somewhat smaller than the measurements averaged over the whole image.
The agreement is not perfect. The estimated errors are a
funtion of radius (because of the variation of the PSF):
the observed scatter is not. That is, the theory predicts lower
scatter for sources near the middle of the image. This implies
that there are other sources of error - not too big a surprise!
These other errors must be fairly uniformly distributed across
the image so as to mask the better precision near the middle.
A likely candidate is small scale variations in sky transparency.
[Tree branches!]
I have been analysing CD5 data using the MK III version of "Star". This has involved chopping the MK IV image (2032x2032) into three strips each 768x2032, the X1 dimension of 768 matching the image dimension for the MK III. This gives an overlap between the strips. The overlap area gets analysed twice and the results are significantly different.
I assume that the differences are the result of the coma. "Star" computes a mean PSF and, because of coma, this differs for the different strips, resulting in different pixels being used in the two cases to compute the position and flux.
The effect results in both scatter and sytematic error in position:
One pixel equals 7500 mas. The systematic errors of 0.05 to 0.1 pixels (375 to 750 mas) are not small compared to the random errors investigated above and are well above Arne Henden's requirement of 100 mas. This is Bad News. The random scatter affects a minority of the points, giving an error of around half a pixel in either direction. Some, but not all, of the points showing this large effect are near the edges of the overlap region, where such things are to be expected.
There is also scatter in the measured magnitudes.
As expected, there is an effect for the PSF based magnitudes. There is also an effect for the 11x11 magnitudes. For both sets, the extreme points (2 for each set) arise from sources near the edges of the overlap region and are therefore not unexpected. The remaining scatter in the 11x11 magnitudes was unexpected since, according to my initial reading of the C code, the 11x11 box used depended only on the position of the brightest pixel. I initially thought that this was entirely independent of the PSF estimate. However, I now believe (subject to correction) that the selection of the brightest pixel is made in a filtered image and is therefore dependent on the estimated PSF.
The aperture fitting method used by "Star" has theoretically higher errors than those of an optimum least squares fit. This has been demonstrated to be the case for actual MK IV images. The errors found in practise are greater than the theory predicts: this is probably the result of other sources of error [including tree branches].
The aperture fitting method used by "Star" also introduces rather large systematic errors when used with comatic MK IV images.
Adoption of a new analysis method involving fitting the variations of the PSF across the image would result in a very large increase in computation requirements. I think we should at least consider this. Computing power gets cheaper every day and the potential improvement in precision is substantial.
My comments are based on the last version of "Star" that I downloaded. Some problems have been fixed since then.
In Stardoc.cpp:
LoadCatalog Get the current catalog.
DetectStars (Starlist.cpp)
MatchStars (Starlist.cpp)
MatchCatalog (Starlist.cpp)
In Starlist.cpp:
DetectStars:
IPPSF (Psf.cpp) Set PSF from bright stars
CorrNormalize (Psf.cpp) Make zero based PSF.
Correlate (1)
IPAperture (Aperture.cpp) Set up apertures.
FindNextPeak Find peaks in convolved image.
IPMeasuredStar (Star.cpp)
MatchStars Match catalog: set up coordinate transformation. (2)
MatchCatalog Match all the stars in the catalog. Optionally
adjust magnitudes.
(1) Convolution with the PSF would give a Wiener optimum amplitude (magnitude) estimate if all points had the same noise variance. Because of photon noise, the variance is not constant so the estimate is not quite optimum for the brighter stars. The use of the zero based PSF is discussed in Technical Note 12.
(2) The routines for matching catalogs are really quite amazing! They don't care which way up the image is and worry about the scaling factors only if the user tells them to. The coordinate transformations they set up currently only go as far as the linear terms (3 coefficients for each axis). For the MK IV, one really needs the quadratic terms (6 coefficients for each axis). This is relatively easy to change.
In Psf.cpp
IPPSF
FindNextPeak Find peaks in unconvolved image.
FWHM Uses linear interpolation starting from the peak to
determine the FWHM. (3)
If FWHM's are close enough to medians, normalize PSF (sum = 1)
and average with all the other successful candidates.(4)
CorrNormalize Subtracts a constant to bring mean or sum to zero. (5)
Scales. (6)
(3) This is different from the method used to determine the position of the source in that only a single row or column of pixels is used rather than the marginal sum and linear interpolation is used instead of cubic. Linear interpolation gives systematic errors in estimating small FWHM's but at least it doesn't run wild as cubic interpolation can!
(4)Note that the normalized response is simply averaged. It is not shifted to centre the source and the source can be off-centre by half a pixel in either coordinate. This is not a good idea if the FWHM is at all small - which it is.
The box size used for the PSF is calculated from the median
FWHM's:
Box size = 2*Round(median*2.5*0.42466) + 1
A centred gaussian would be down to 4% at the edge of the box
and would have a bit over 1% of its total flux outside the box.
Unfortunately, our PSF's are far from gaussian and have a great
deal more flux outside the box.
No note is made of how many sources are rejected because they have out-of-range FWHM's. Pity. It would be nice to know.
(5) The problem with this subtraction of a constant over the area of the box used to estimate the PSF is that boxes have edges. These edges give rise to quite serious artefacts near bright sources. The solution to this problem is to subtract something which goes smoothly to zero at the edges of the box. Unfortunately, this requires the use of a much larger box, resulting in much larger loss of data around the edges of the image. The larger number of points in the convolution also increase the computing time but this is not a serious problem.
(6) The scaling is such as to make the peak values in the convolved image match the total flux (sum of values above background) in the original image. The weighting factor is closely related to my NEff factor. Or at least it would be if the PSF box were large enough to catch most of the flux.
In Aperture.cpp
IPAperture (PSF option) Set up to use points for which the mean PSF of
the bright sources exceeds some value involving a
factor denomAdjust = 2.0 (7)
(other apertures) Set up to use points within an ellipse
filling the bounding rectangle.
(7) The PSF values are sorted into numerical order. Starting with the largest, points are accepted down to and including the first point for which the value of the PSF is less than
(the average of all the larger values)/denomAdjustDoing some sums on the backs of old envelopes, I reckon that for a centred gaussian with FWHM < 2.0 this algorithm accepts exactly two points - the brightest plus whichever of the four adjacent points is tested first. This is obviously asymmetrical in a random sort of way; the method was never intended for systems with such a small FWHM.
The method automatically eliminates the halo in the MK III and the comatic tail in the MK IV together with around 60% of the flux. Apart from the waste of potential data, this is not too serious for the MK III: the fraction of flux lost is quite constant across the image and therefore simple to calibrate out. The problem was much more serious for the original comatic MK IV lenses but should now be quite manageable.
In Star.cpp
IPMeasuredStar
BaseValue (Image.cpp) Median of pixels outside bounds. (8)
Sets up marginal and cumulative marginal distributions using
PSF aperture.
Interpolate Calculates x, y (9)
Calculates ra and dec
FWHM (Psf.cpp; covered above)
Sums flux over apertures. (10)
(8) According to my sums, the median is of 336 points except near the edges of the image. Routine BaseValue apparently uses an array dimensioned 100 to compute this. I presume that coping with this sort of thing is a feature of C. Elsewhere, the variance of this estimate is taken as 4*Sqr(baseNoise)/NumBasePixels. I think this is an overestimate; the 4 should be replaced (assuming gaussian noise) by something like Pi/2. In any case, in spite of the comment in the code that the contribution of the base estimate to the total variance of the flux estimate increases with the square of the number of source pixels, the lack of accuracy of the base estimate is not the major contributor to the error budget.
(9) The method used is akin to that proposed by Stone (Summary, p1237) combining the median method with image thresholding except that the data trimming takes place using a fixed window centred on the peak pixel. This is not the best way; Stone says (p1230)
The use of trimming requires that each centering solution be iterated until the trimming process is convergent.
The method suffers the disadvantages of median methods. Stone shows that only the Derivative Search method is worse affected by a small FWHM. He obtained an error of 0.008 arc sec for an FWHM of 2 pixels, deteriorating rapidly for smaller FWHM. Unfortunately, in his haste to convert from pixels to arc sec, Stone neglects to say how big his pixel was at this point. Elsewhere, he uses pixels of 1/8, 1/4 and 1/2 arc sec so, scaling appropriately for the MK IV, his 0.008 arc sec probably corresponds to a MK IV error of 0.48, 0.24 or 0.12 arc sec: in any case larger than Arne's specification of 100 mas.
The last version of Interpolate that I downloaded still contains the bug that returns a result in the range 0..1 if the correct result is in the range last-1..last. This bug never gets exercised by real data ... I think.
(10) Contrary to what some posters to the list seem to believe, all accepted points receive equal weight and all rejected points get a weight of zero. The flux determined using the PSF aperture is the sum of all the accepted pixels, not a weighted sum using the average PSF.
This method is not intended to be used with an FWHM as small as two pixels.
The flux determined using the PSF aperture is compared with that obtained from the convolved image, specifically with the peak value used to select the source in the first place. If these disagree by too much, the source is rejected. Again, it would be nice to know how many sources are rejected in this way.