IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12205


Ignore:
Timestamp:
Mar 2, 2007, 4:34:05 PM (19 years ago)
Author:
eugene
Message:

adding polynomial fits for chip and mosaic

Location:
branches/dvo-mods-2007-02/Ohana/src/relastro/src
Files:
3 added
6 edited

Legend:

Unmodified
Added
Removed

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitSimple.c

    r12068 r12205  
    88  for (i = 0; i < Nmatch; i++) {
    99    if (raw[i].mask) continue;
    10     fit_add (raw[i].X, raw[i].Y, ref[i].R, ref[i].D);
     10    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
    1111  }
    1212  fit_eval ();
    … …  
    2828   P,Q -> X,Y (polynomial transformation)
    2929*/
     30
     31/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
     32
     33   

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/ImageOps.c

    r12068 r12205  
    178178}
    179179
    180 StarData *getImageRaw (int im, int *Nstars) {
     180StarData *getImageRaw (int im, int *Nstars, int isMosaic) {
    181181
    182182  StarData *raw;
    183183 
    184184  ALLOCATE (raw, StarData, Nlist[im]);
     185
     186  if (isMosaic) {
     187    mosaic = getMosaicForImage (im);
     188  }
    185189
    186190  for (i = 0; i < Nlist[im]; i++) {
    … …  
    196200
    197201    /* note that for a Simple image, L,M = P,Q */
    198     XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].X, raw[i].Y, &image[im].coords);
     202    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
     203    if (isMosaic) {
     204      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, &mosaic);
     205      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &mosaic);
     206    } else {
     207      raw[i].P = raw[i].L;
     208      raw[i].Q = raw[i].M;
     209      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
     210    }
    199211  }
    200212 
    … …  
    203215}
    204216
    205 StarData *getImageRef (int im, int *Nstars) {
    206 
    207   StarData *raw;
    208  
    209   ALLOCATE (raw, StarData, Nlist[im]);
     217StarData *getImageRef (int im, int *Nstars, int isMosaic) {
     218
     219  StarData *ref;
     220 
     221  ALLOCATE (ref, StarData, Nlist[im]);
     222
     223  if (isMosaic) {
     224    mosaic = getMosaicForImage (im);
     225  }
    210226
    211227  for (i = 0; i < Nlist[im]; i++) {
    … …  
    217233    ref[i].R = catalog[c].average[n].R;
    218234    ref[i].D = catalog[c].average[n].D;
     235
     236    /* note that for a Simple image, L,M = P,Q */
     237    RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
     238    if (isMosaic) {
     239      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, &mosaic);
     240      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &mosaic);
     241    } else {
     242      ref[i].L = ref[i].P;
     243      ref[i].M = ref[i].Q;
     244      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
     245    }
    219246  }
    220247 
    221248  *Nstars = Nlist[im];
    222   return (raw);
    223 }
     249  return (ref);
     250}

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateChips.c

    r12068 r12205  
    1818
    1919    /* convert measure coordinates to raw entries */
    20     raw = getImageRaw (&image[i], &Nstars);
     20    raw = getImageRaw (&image[i], &Nstars, TRUE);
    2121
    2222    /* convert average coordinates to ref entries */
    23     ref = getImageRef (&image[i], &Nstars);
     23    ref = getImageRef (&image[i], &Nstars, TRUE);
    2424
    2525    /* XXX grab this code from mosastro */

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateMosaic.c

    r12068 r12205  
    1818
    1919    /* convert measure coordinates to raw entries */
    20     raw = getImageRaw (&image[i], &Nstars);
     20    raw = getImageRaw (&image[i], &Nstars, TRUE);
    2121
    2222    /* convert average coordinates to ref entries */
    23     ref = getImageRef (&image[i], &Nstars);
     23    ref = getImageRef (&image[i], &Nstars, TRUE);
    2424
    2525    /* XXX grab this code from mosastro */

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateSimple.c

    r12068 r12205  
    1919
    2020    /* convert measure coordinates to raw entries */
    21     raw = getImageRaw (&image[i], &Nstars);
     21    raw = getImageRaw (&image[i], &Nstars, FALSE);
    2222
    2323    /* convert average coordinates to ref entries */
    24     ref = getImageRef (&image[i], &Nstars);
     24    ref = getImageRef (&image[i], &Nstars, FALSE);
    2525
    2626    /* XXX grab this code from mosastro */

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/fitpoly.c

    r12068 r12205  
    11# include "relastro.h"
    22
    3 static int NTERM, NPOWER, NPARS, NORDER, Npts;
    4 static double **sum, **xsum, **ysum;
    5 static double **matrix, **vector;
    6 
    7 void fit_init (int order) {
     3/* these functions support simultaneous 2D fits to
     4   x2 = \sum a_i,j x1^i y1^j
     5   y2 = \sum b_i,j x1^i y1^j
     6
     7   the order of the fit (largest coefficient) is fixed to a single
     8   value for both x1,y1 terms and for x2,y2 fits
     9
     10   the code is currently confusing because we limit to i+j <= order.
     11   this could be cleaner if we used masks and allowed i <= order, j <= order
     12*/
     13
     14// XXX define a fit structure and drop the file static variables?
     15CoordFit *fit_init (int order) {
    816
    917  int i;
    1018
    11   Npts  = 0;
    12   NORDER = order;
    13   NPOWER = NORDER + 1;
    14   NTERM = 2*NORDER + 1;
    15   NPARS = (NORDER + 1)*(NORDER + 2) / 2;
    16 
    17   /* allocate arrays for fit solution */
    18   ALLOCATE (sum, double *, NTERM);
    19   ALLOCATE (xsum, double *, NTERM);
    20   ALLOCATE (ysum, double *, NTERM);
    21   for (i = 0; i < NTERM; i++) {
    22     ALLOCATE (sum[i], double, NTERM);
    23     bzero (sum[i], NTERM*sizeof(double));
    24     ALLOCATE (xsum[i], double, NTERM);
    25     bzero (xsum[i], NTERM*sizeof(double));
    26     ALLOCATE (ysum[i], double, NTERM);
    27     bzero (ysum[i], NTERM*sizeof(double));
    28   }
    29   ALLOCATE (matrix, double *, NPARS);
    30   ALLOCATE (vector, double *, NPARS);
    31   for (i = 0; i < NPARS; i++) {
    32     ALLOCATE (matrix[i], double, NPARS);
    33     ALLOCATE (vector[i], double, 2);
    34     bzero (vector[i], 2*sizeof(double));
    35     bzero (matrix[i], NPARS*sizeof(double));
    36   }
    37 
    38 }
    39 
    40 # define SCALE 1.0
    41 void fit_add (double x1, double y1, double x2, double y2) {
    42 
    43   int n, m;
     19  ALLOCATE (fit, CoordFit, 1);
     20
     21  fit[0].Npts   = 0;
     22  fit[0].Norder = order;
     23  fit[0].Nterms = order + 1;
     24  fit[0].Nsums  = 2*order + 1;
     25  fit[0].Nelems = SQ(order + 1);
     26
     27  /* summing arrays for fit solution */
     28
     29  // xsum[i][j] holds \sum (x2 wt x1^i y1^j)
     30  // ysum[i][j] holds \sum (y2 wt x1^i y1^j)
     31  ALLOCATE (fit[0].xsum, double *, fit[0].Nterms);
     32  ALLOCATE (fit[0].ysum, double *, fit[0].Nterms);
     33  for (i = 0; i < fit[0].Nterms; i++) {
     34    bzero (fit[0].xsum[i], fit[0].Nterms*sizeof(double));
     35    bzero (fit[0].ysum[i], fit[0].Nterms*sizeof(double));
     36  }     
     37
     38  // sum[i][j] holds \sum (wt x1^i y1^j)
     39  ALLOCATE (fit[0].sum, double *, fit[0].Nsums);
     40  for (i = 0; i < fit[0].Nsums; i++) {
     41    ALLOCATE (fit[0].sum[i], double, fit[0].Nsums);
     42    bzero (fit[0].sum[i], fit[0].Nsums*sizeof(double));
     43  }
     44
     45  // matrix, vector hold the final linear system
     46  ALLOCATE (fit[0].matrix, double *, fit[0].Nelems);
     47  ALLOCATE (fit[0].vector, double *, fit[0].Nelems);
     48  for (i = 0; i < fit[0].Nelems; i++) {
     49    ALLOCATE (fit[0].matrix[i], double, fit[0].Nelems);
     50    ALLOCATE (fit[0].vector[i], double, fit[0].Nelems);
     51    bzero (fit[0].matrix[i], fit[0].Nelems*sizeof(double));
     52    bzero (fit[0].vector[i], 2*sizeof(double));
     53  }
     54  return (fit);
     55}
     56
     57// XXX use implicit masks as below or explicit masks (with function to set?)
     58// XXX eg, add a global mask to this file and
     59void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) {
     60
     61  int ix, iy;
    4462  double xterm, yterm, term;
    4563
    4664  xterm = 1;
    47   for (n = 0; n < NTERM; n++) {
     65  for (ix = 0; ix < fit[0].Nsums; ix++) {
    4866    yterm = 1;
    49     for (m = 0; m < NTERM; m++) {
    50       term = xterm*yterm;
    51       if (n+m < NTERM) {
    52         sum[n][m] += term;
    53       }
    54       if (n+m < NPOWER) {
    55         xsum[n][m] += x2*term;
    56         ysum[n][m] += y2*term;
    57       }
    58       yterm *= y1/SCALE;
    59     }
    60     xterm *= x1/SCALE;
    61   }
    62   Npts ++;
    63 }
    64 
    65 /** I am renormalizing here by the max pivots to keep gaussj sane **
    66  ** would not be needed if the fit used scaled ind. variables **/
    67 void fit_eval () {
     67    for (iy = 0; iy < fit[0].Nsums; iy++) {
     68      term = xterm*yterm*wt;
     69      fit[0].sum[ix][iy] += term;
     70      if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) {
     71        fit[0].xsum[ix][iy] += x2*term;
     72        fit[0].ysum[ix][iy] += y2*term;
     73      }
     74      yterm *= y1;
     75    }
     76    xterm *= x1;
     77  }
     78  fit[0].Npts ++;
     79}
     80
     81/* convert the xsum,ysum,sum terms into vector,matrix and solve */
     82void fit_eval (CoordFit *fit) {
    6883
    6984  int i, j, n, m, M, N;
    7085
    71   if (Npts == 0) {
     86  if (fit[0].Npts == 0) {
    7287    fprintf (stderr, "warning: no valid pts\n");
    7388  }
    7489
    75   i = 0;
    76   for (m = 0; m < NPOWER; m++) {
    77     for (n = 0; n < NPOWER - m; n++, i++) {
    78       vector[i][0] = xsum[n][m];
    79       vector[i][1] = ysum[n][m];
    80     }   
    81   }
    82   j = 0;
    83   for (M = 0; M < NPOWER; M++) {
    84     for (N = 0; N < NPOWER - M; N++, j++) {
    85       i = 0;
    86       for (m = 0; m < NPOWER; m++) {
    87         for (n = 0; n < NPOWER - m; n++, i++) {
    88           matrix[i][j] = sum[n+N][m+M];
    89         }       
    90       }
    91     }
    92   }       
    93 # if (0)
    94   max = 0.0;
    95   for (i = 0; i < NPARS; i++) {
    96     for (j = 0; j < NPARS; j++) {
    97       max = MAX (max, fabs(matrix[i][j]));
    98     }
    99     max = MAX (max, fabs(vector[i][0]));
    100     max = MAX (max, fabs(vector[i][1]));
    101   }
    102   for (i = 0; i < NPARS; i++) {
    103     for (j = 0; j < NPARS; j++) {
    104       matrix[i][j] /= max;
    105     }
    106     vector[i][0] /= max;
    107     vector[i][1] /= max;
    108   }
    109 # endif
     90  /* remap the xsum,ysum terms into the vector */
     91  for (i = 0; i < fit[0].Nelems; i++) {
     92    ix = i % fit[0].Nterms;
     93    iy = i / fit[0].Nterms;
     94    vector[i][0] = xsum[ix][iy];
     95    vector[i][1] = ysum[ix][iy];
     96
     97    for (j = 0; j < fit[0].Nelems; j++) {
     98      jx = j % fit[0].Nterms;
     99      jy = j / fit[0].Nterms;
     100      matrix[i][j] = sum[ix+jx][iy+jy];
     101    }
     102  }
    110103
    111104  dgaussj (matrix, NPARS, vector, 2);
    112105
    113 # if (0)
    114   i = 0;
    115   for (m = 0; m < NPOWER; m++) {
    116     for (n = 0; n < NPOWER - m; n++, i++) {
    117       fprintf (stderr, "RA x^%dy^%d: %10.4g    DEC x^%dy^%d: %10.4g \n",
    118                n, m, vector[i][0], n, m, vector[i][1]);
    119     }   
    120   }
    121 # endif
     106  for (i = 0; i < fit[0].Nelems; i++) {
     107    ix = i % fit[0].Nterms;
     108    iy = i / fit[0].Nterms;
     109    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n",
     110             ix, iy, vector[i][0], ix, iy, vector[i][1]);
     111  }     
    122112}
    123113
    124114/* linear portion of fit : NORDER is 1 */
    125 void fit_apply_coords (Coords *coords) {
     115/* this should only apply to the polynomial, not the projection terms */
     116/* compare with psastro supporting code */
     117void fit_apply_coords (CoordFit *fit, Coords *coords) {
    126118
    127119  int i, j, Np, Nv, N;
    … …  
    130122  double R;
    131123
    132   /* update the higher order terms */
    133   if (NORDER > 1) {
    134     for (i = 0; i < NPOWER; i++) {
    135       for (j = 0; j < (NPOWER - i); j++) {
    136         if (i + j < 2) continue;
    137         Np = mkpolyterm (i, j);
    138         Nv = mkvector (i, j, NORDER);
    139         coords[0].polyterms[Np][0] = vector[Nv][0];
    140         coords[0].polyterms[Np][1] = vector[Nv][1];
    141       }
    142     }
    143   }
     124  /* I have L,M = fit(X,Y). set corresponding terms for coords */
     125
     126  // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ...
     127  // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
     128
     129  coords[0].crpix1,crpix2 = CoordsGetCenter (fit);
     130  modfit = CoordsSetCenter (fit, xo, yo);
     131
     132  // set pc1_1, pc1_2, pc2_1, pc2_2
     133  fit_to_coordterms (coords, modfit, 0, 1);
     134  fit_to_coordterms (coords, modfit, 1, 0);
     135
     136  fit_to_coordterms (coords, modfit, 0, 2);
     137  fit_to_coordterms (coords, modfit, 1, 1);
     138  fit_to_coordterms (coords, modfit, 2, 0);
     139
     140  fit_to_coordterms (coords, modfit, 0, 3);
     141  fit_to_coordterms (coords, modfit, 1, 2);
     142  fit_to_coordterms (coords, modfit, 2, 1);
     143  fit_to_coordterms (coords, modfit, 3, 0);
    144144
    145145  /* get the correct vector entries for the linear terms */
    146   N = mkvector (0, 0, NORDER);
    147146  coords[0].crval1 = vector[N][0]; 
    148147  coords[0].crval2 = vector[N][1];
Note: See TracChangeset for help on using the changeset viewer.