Possible method for handling trailed star images

My current astrometry software, and a lot of current astrometry software in general, models a "normal" star image as a Gaussian-plus-background. Thus, it assumes that, within a few pixels of the star (i.e., within the centroiding box), you can model the intensity of pixel (x, y) as

I(x, y) = B + A exp( -[(x-x₀)² + (y-y₀)²] / width²)

where A=amplitude, x₀, y₀ = star position, width=width of the star image, B=background. My software first estimates these five values using a center-of-mass method, then improves the values by using a least-squares fit. Pretty standard stuff. (Code for this is available in the Charon source code, as part of FINDSTAR.CPP. Specifically, look at the fit_psf_star( ) function. Except calling it a "PSF fit" is a slight exaggeration. Gory details as to how Charon does this are available here.)

OK, here's the less-standard part. Imagine that instead, over the duration of the image, we're tracking an asteroid so that the stars will trail dx pixels in x, dy pixels in y. (Note that I'm assuming here that we know the image orientation, so we can convert the raw math of "the object moved 3.4 arcseconds north and 8.3 arcsec east" into pixel-based coordinates such as "4.1 pixels in x and -1.4 pixels in y". Not a big assumption.)

In that case, we've got the same five parameters to fit, * but the intensity now is given as an integral of the above expression:

               .5
I(x,y) = B + A  exp( -[(x-x₀+t*dx)² + (y-y₀+t*dy)²] / width²) dt
              t=-.5

The integration runs from -.5 to +.5 so that we find the center, rather than the end of, one of the trails. As is shown below, it can be transformed a bit to result in finding out that it's the standard "area under a bell curve", a.k.a. "error function" or "erf". Do that, and you can evade doing a real numerical integration; there are several pieces of code that approximate erf very closely, using a variety of methods. You can click here for a C/C++ example implementation of the above function that does just that.

* David Tholen pointed out that the trail length is known only if we're tracking a target and letting stars trail (a common technique when dealing with faint minor planets). In the more usual case of measuring a trailed minor planet image, you will frequently not have any real clue what the trail length (i.e., dx and dy) are. In that situation, you'll have to do a least-squares fit of seven, not five, parameters.

Getting initial values for those five (or seven) parameters may be a little ugly, since the usual "simple centroiding" step won't work if the trails are very long. I plan to let users click on one end of the trail, then the other; that gives us initial values for x₀, y₀, dx, and dy. If the object's motion is known, we can compute dx and dy and "fix" them, least-squares-fitting only the remaining five parameters.

Jean-Claude Pelle suggested that this is related to a similar "trailing" problem: timing asteroid occultations by letting the target star move at a sidereal rate through a CCD camera. When it's occulted, the trail will stop (or at least get a lot thinner). When it emerges from occultation, the trail begins again. There's a discussion of this technique at

http://chez.mana.pf/~jcpelle/occultat_EN.htm

Transforming the above integral into the error function:

Set kx = x - x₀, ky = y - y₀ (constants during the integration); then

               .5
I(x,y) = B + A  exp( -[(kx+t*dx)² + (ky+t*dy)²] / width²) dt
              t=-.5

         .5
 = B + A  exp( -[ kx² + 2*kx*dx*t + (dx*t)² + ky² + 2*ky*dy*t + (dy*t)²] / width²) dt
        t=-.5

         .5
 = B + A  exp( -[ kx² + ky² + 2*(kx*dx+ky*dy)*t + (dx²+dy²)*t²] / width²) dt
        t=-.5

Set C=(kx²+ky²)/width², D=2*(kx*dx+ky*dy)/width², E=√(dx²+dy²)/width; note that all of these are constants during the integration:

         .5
   B + A  exp( -[C + D * t + E² * t²]) dt
        t=-.5

Substitute u = E t, du = E dt...

           .5 E
   B + (A/E)  exp( -[C + (D/E) * u + u²]) du
          u=-.5 E

Substitute v = u+(D/2E),

           .5 E + (D/2E)
   B + (A/E)  exp( -[C - (D/2E)² + v²]) dv
        v = -.5 E + (D/2E)

Since C - (D/2E)² is constant during the integration, we can pull it out of the integral to get...

                                  .5(E + D/E)
I(x, y) = B + (A/E) exp( -[C-(D/2E)²])   exp( -v²) dv
                                  v = -.5(-E + D/E)

where the integral is the error function erf().

C/C++ source code: The following code implements the above math. As you'll see, it includes test code for the error function; the rest is untested, but looks to be quite simple anyway.

#include <math.h>

/* Returns the integral of exp( -x^2) dx,  from x=0 to x=arg... the
"error function".  It does so using a lookup table of values of the
integral,  in .05 steps,  running from -.05 up to the point where the
values stopped changing (arg=4.35,  in the nine-place table I generated.)
(Not that it matters much,  but that asymptotic limit .866226925 equals
√pi / 2.)

The table was generated using an "older" err_func( ) that used Simpson's
method (see below).  Once I used that to generate a table,  it was
possible just to do a cubic interpolation within the table. A comparison
showed that the interpolated version is accurate to six places (worst case).
   */

#define N_LOOKUP 89
#define LOOKUP_STEP .05

static double err_func( const double arg)
{
static const double lookup_tbl[N_LOOKUP] = { -.049958365,
      .000000000, .049958365, .099667664, .148882553, .197365031,
      .244887887, .291237883, .336218591, .379652840, .421384703,
      .461281006, .499232335, .535153527, .568983677, .600685668,
      .630245271, .657669856, .682986783, .706241515, .727495536,
      .746824133, .764314099, .780061433, .794169078, .806744758,
      .817898943, .827742989, .836387469, .843940714, .850507572,
      .856188394, .861078228, .865266226, .868835241, .871861593,
      .874415002, .876558634, .878349281, .879837610, .881068494,
      .882081391, .882910749, .883586442, .884134195, .884576021,
      .884930627, .885213812, .885438833, .885616746, .885756711,
      .885866274, .885951610, .886017745, .886068745, .886107876,
      .886137752, .886160447, .886177601, .886190503, .886200159,
      .886207348, .886212675, .886216602, .886219483, .886221585,
      .886223112, .886224216, .886225009, .886225577, .886225981,
      .886226267, .886226469, .886226610, .886226709, .886226777,
      .886226825, .886226857, .886226879, .886226895, .886226905,
      .886226912, .886226916, .886226920, .886226922, .886226923,
      .886226924, .886226924, .886226925 };
   double table_loc, rval = 0., dx;
   int idx;

   if( arg < 0.)              /* let's make use of symmetry: */
      return( -err_func( -arg));
   table_loc = 1. + arg / LOOKUP_STEP;
   idx = (int)( table_loc - 1);
   if( idx + 3 >= N_LOOKUP)       /* past the range of the table */
      return( lookup_tbl[N_LOOKUP - 1]);
   dx = table_loc - (double)idx;
   rval = lookup_tbl[idx] * (dx - 1.) * (dx - 2.) * (dx - 3.) / -6.
         + lookup_tbl[idx + 1] * dx * (dx - 2.) * (dx - 3.) / 2.
         + lookup_tbl[idx + 2] * dx * (dx - 1.) * (dx - 3.) / -2.
         + lookup_tbl[idx + 3] * dx * (dx - 1.) * (dx - 2.) / 6.;
   return( rval);
}

#ifdef SIMPSONS_FORMULA_METHOD

/* I used the following version of err_func( ) to generate the lookup
table used in the above function.  It's not all that useful anymore,
and isn't really all that good or fast.  That's because all it had to
do was to work _once_,  so I could use it to generate the table.
It's included just for reference. */

static double err_func( const double arg)
{
   const double max_step_size = .001;
   double step_size, rval = 0., x = 0.;
   int i, n_steps;

   if( arg < 0.)              /* let's make use of symmetry: */
      return( -err_func( -arg));
   n_steps = (int)ceil( arg / max_step_size);
   if( !n_steps)              /* gotta have at least one step */
      n_steps++;
   if( n_steps % 2)           /* gotta have an even number of steps */
      n_steps++;
   step_size = arg / (double)n_steps;
   for( i = 0; i <= n_steps; i++, x += step_size)
      {
      double value = exp( -x * x);

      if( !i || i == n_steps)         /* first and last steps */
         rval += value;
      else if( i % 2)                     /* on odd steps... */
         rval += 4. * value;
      else                                /* even steps... */
         rval += 2. * value;
      }
   return( rval * step_size / 3.);
}
#endif

double trailed_intensity( const double a, const double backgnd,
         const double x0, const double y0,
         const double x, const double y,
         const double dx, const double dy, const double width)
{
   const double kx = x - x0;
   const double ky = y - y0;
   const double c = (kx * kx + ky * ky) / (width * width);
   const double d = 2 * (kx * dx + ky * dy) / (width * width);
   const double e = sqrt( dx * dx + dy * dy) / width;
   const double f = d / (2. * e);
   const double coeff = (a / e) * exp( -c + f * f);

   return( backgnd + coeff * (err_func( .5 * e + f) - err_func( -.5 * e + f)));
}

#ifdef TEST_FUNC
#include <stdio.h>
#include <stdlib.h>

   /* Used to generate a table of error function values,  for use in the
table lookup/cubic spline version of err_func;  then I used it again to
verify that the results of the lookup/spline approach were good. */

void main( int argc, char **argv)
{
   if( argc == 2)       /* just looking for one value */
      printf( "%.9lf\n", err_func( atof( argv[1])));
   else                 /* making a _table_ of values */
      {
      double step = atof( argv[1]), x = 0.;
      int i;

      for( i = 0; i < atoi( argv[2]); i++, x += step)
         printf( "%lf: %.9lf\n", x, err_func( x));
      }
}
#endif