IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12068


Ignore:
Timestamp:
Feb 26, 2007, 5:00:31 PM (19 years ago)
Author:
eugene
Message:

working on simple image fitting

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

Legend:

Unmodified
Added
Removed

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

    r12050 r12068  
    22
    33/* do we want an init function which does the alloc and a clear function to free? */
    4 FitPM (double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
     4int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
    55
    6   PFfit *fit;
    76  double **A, **B;
    87
    … …  
    5756  dgaussj (A, 4, B, 1);
    5857
    59   ALLOCATE (fit, PMfit, 1);
    60 
    6158  fit[0].Ro = B[0];
    6259  fit[0].uR = B[1];
    … …  
    7875  free (B);
    7976
    80   return (fit);
     77  return (TRUE);
    8178}

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

    r12050 r12068  
    22
    33/* do we want an init function which does the alloc and a clear function to free? */
     4int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
    45
    5 FitPM (double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
    6 
    7   PFfit *fit;
    86  double **A, **B;
    97
    … …  
    8078  dgaussj (A, 5, B, 1);
    8179
    82   ALLOCATE (fit, PMfit, 1);
    83 
    8480  fit[0].Ro = B[0];
    8581  fit[0].uR = B[1];
    … …  
    10197  free (B);
    10298
    103   return (fit);
     99  /* get the chisq from the matrix values */
     100
     101  return (TRUE);
    104102}

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

    r12048 r12068  
    1212static Image        *image;
    1313static int          Nimage;
     14
     15Image *getimages (int *N) {
     16  *N = Nimage;
     17  return (image);
     18}
     19
     20Image *getimage (int N) {
     21  return (&image[N]);
     22}
    1423
    1524void initImages (Image *input, int N) {
    … …  
    8291}
    8392
     93/* modify this function to use the measure->imageID field */
    8494void matchImage (Catalog *catalog, int meas, int cat) {
    8595
    … …  
    96106    if (measure[0].t > stop[i]) continue;
    97107   
    98 # ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
    99     if (USE_GRID) {
    100 
    101       /* identify the ccd on the basis of the photcode name */
    102       pname = GetPhotcodeNamebyCode (image[i].photcode);
    103       filter = photcode[0].name;
    104       sprintf (base, "%s.%s.", MOSAICNAME, filter);
    105       if (strncmp (pname, base, strlen (base))) continue;
    106       p = pname + strlen(base);
    107       if (*p == 0) continue;
    108       ccdnum = atoi (p);
    109       /* some ambiguity here between seq number and id number */
    110 
    111       /* add this measurement to the grid cell for this chip */
    112       ave = measure[0].averef;
    113       ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
    114       dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
    115        
    116       /* X,Y always positive-definite in range 0,0 - dX, dY */
    117       RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    118       setGridMeasure (meas, cat, X, Y, ccdnum);
    119     }
    120 # endif
    121 
    122108    bin[cat][meas] = i;
    123109
    … …  
    133119    return;
    134120  }
    135   /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
    136 }
    137 
    138 float getMcal (int meas, int cat) {
    139 
    140   int i;
    141   float value;
    142 
    143   i = bin[cat][meas];
    144   if (i == -1) return (NO_MAG);
    145 
    146   if (image[i].code & IMAGE_BAD)  return (NO_MAG); 
    147   value = image[i].Mcal;
    148   return (value);
    149121}
    150122
    … …  
    156128  if (i == -1) return (NULL);
    157129  return (&image[i].coords);
    158 }
    159 
    160 /* determine Mcal values for all images */
    161 void setMcal (Catalog *catalog, int PoorImages) {
    162 
    163   int i, j, m, c, n, N, Nmax, mark, bad;
    164   float Msys, Mrel, Mmos, Mgrid;
    165   double *list, *dlist;
    166   StatType stats;
    167 
    168   if (FREEZE_IMAGES) return;
    169 
    170   if (PoorImages) {
    171     IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    172   }
    173 
    174   Nmax = 0;
    175   for (i = 0; i < Nimage; i++) {
    176     Nmax = MAX (Nmax, Nlist[i]);
    177   }
    178   ALLOCATE (list, double, Nmax);
    179   ALLOCATE (dlist, double, Nmax);
    180 
    181   for (i = 0; i < Nimage; i++) {
    182    
    183     /* on PoorImages run, skip good images */
    184     if (PoorImages) {
    185       bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
    186       if (!bad) continue;
    187     }     
    188 
    189     N = 0;
    190     for (j = 0; j < Nlist[i]; j++) {
    191      
    192       m = mlist[i][j];
    193       c = clist[i][j];
    194      
    195       if (catalog[c].measure[m].flags & MEAS_BAD) continue;
    196       if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
    197       if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
    198       if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
    199      
    200       n = catalog[c].measure[m].averef;
    201       Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
    202       list[N] = Msys - Mrel - Mmos - Mgrid;
    203       dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
    204       N++;
    205     }
    206     /* Nlist[i] is all measurements, N is good measurements */
    207 
    208     /* too few good measurements or too many bad measurements */
    209     if (!PoorImages) {
    210       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
    211       if (mark) {
    212         image[i].code |= ID_IMAGE_FEW;
    213       } else {
    214         image[i].code &= ~ID_IMAGE_FEW;
    215       }     
    216     }
    217    
    218     liststats (list, dlist, N, &stats);
    219     image[i].Mcal  = stats.mean;
    220     image[i].dMcal = stats.sigma;
    221     image[i].Xm    = 100.0*log10(stats.chisq);
    222   }
    223   free (list);
    224   free (dlist);
    225   if (PoorImages) {
    226     IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
    227     STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    228     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
    229   }
    230   return;
    231 }
    232 
    233 /* mark image if: abs(Mcal) too large, dMcal too large */
    234 void clean_images () {
    235 
    236   int i, N, mark, Nmark;
    237   double *mlist, *slist, *dlist;
    238   double MaxOffset, MaxScatter, MedOffset;
    239   StatType stats;
    240 
    241   if (FREEZE_IMAGES) return;
    242 
    243   if (VERBOSE) fprintf (stderr, "marking poor images\n");
    244 
    245   ALLOCATE (mlist, double, Nimage);
    246   ALLOCATE (slist, double, Nimage);
    247   ALLOCATE (dlist, double, Nimage);
    248 
    249   for (i = N = 0; i < Nimage; i++) {
    250     if (image[i].code & IMAGE_BAD) continue;
    251     mlist[N] = fabs (image[i].Mcal);
    252     slist[N] = image[i].dMcal;
    253     dlist[N] = 1;
    254     N++;
    255   }
    256   initstats ("MEAN");
    257   liststats (mlist, dlist, N, &stats);
    258   MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
    259   MedOffset = stats.median;
    260   liststats (slist, dlist, N, &stats);
    261   MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    262   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    263  
    264   Nmark = 0;
    265   for (i = 0; i < Nimage; i++) {
    266     mark = FALSE;
    267     image[i].code &= ~ID_IMAGE_POOR;
    268     mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset);
    269     if (mark) {
    270       Nmark ++;
    271       image[i].code |= ID_IMAGE_POOR;
    272     } else {
    273       image[i].code &= ~ID_IMAGE_POOR;
    274     }
    275   }
    276 
    277   fprintf (stderr, "%d images marked poor\n", Nmark);
    278   initstats (STATMODE);
    279   free (mlist);
    280   free (slist);
    281   free (dlist);
    282130}
    283131
    … …  
    330178}
    331179
    332 StatType statsImageN (Catalog *catalog) {
    333 
    334   int i, j, m, c, n, N;
    335   double *list, *dlist;
    336   StatType stats;
    337 
    338   bzero (&stats, sizeof (StatType));
    339   if (FREEZE_IMAGES) return (stats);
    340 
    341   ALLOCATE (list, double, Nimage);
    342   ALLOCATE (dlist, double, Nimage);
    343 
    344   n = 0;
    345   for (i = 0; i < Nimage; i++) {
    346     if (image[i].code & IMAGE_BAD)  continue;
    347 
    348     N = 0;
    349     for (j = 0; j < Nlist[i]; j++) {
    350 
    351       m = mlist[i][j];
    352       c = clist[i][j];
    353 
    354       if (getMcal  (m, c) == NO_MAG) continue;
    355       if (getMmos  (m, c) == NO_MAG) continue;
    356       if (getMgrid (m, c) == NO_MAG) continue;
    357       N++;
    358     }
    359     list[n] = N;
    360     dlist[n] = 1;
    361     n++;
    362   }
    363 
    364   liststats (list, dlist, n, &stats);
    365   free (list);
    366   free (dlist);
    367   return (stats);
    368 }
    369 
    370 StatType statsImageX (Catalog *catalog) {
    371 
    372   int i, n;
    373   double *list, *dlist;
    374   StatType stats;
    375 
    376   bzero (&stats, sizeof (StatType));
    377   if (FREEZE_IMAGES) return (stats);
    378 
    379   ALLOCATE (list, double, Nimage);
    380   ALLOCATE (dlist, double, Nimage);
    381 
    382   n = 0;
    383   for (i = 0; i < Nimage; i++) {
    384 
    385     if (image[i].code & IMAGE_BAD)  continue;
    386 
    387     list[n] = pow (10.0, 0.01*image[i].Xm);
    388     dlist[n] = 1;
    389     n++;
    390   }
    391 
    392   liststats (list, dlist, n, &stats);
    393   free (list);
    394   free (dlist);
    395   return (stats);
    396 }
    397 
    398 StatType statsImageM (Catalog *catalog) {
    399 
    400   int i, n;
    401   double *list, *dlist;
    402   StatType stats;
    403 
    404   bzero (&stats, sizeof (StatType));
    405   if (FREEZE_IMAGES) return (stats);
    406 
    407   ALLOCATE (list, double, Nimage);
    408   ALLOCATE (dlist, double, Nimage);
    409 
    410   n = 0;
    411   for (i = 0; i < Nimage; i++) {
    412 
    413     if (image[i].code & IMAGE_BAD)  continue;
    414 
    415     list[n] = image[i].Mcal;
    416     dlist[n] = 1;
    417     n++;
    418   }
    419 
    420   liststats (list, dlist, n, &stats);
    421   free (list);
    422   free (dlist);
    423   return (stats);
    424 }
    425 
    426 StatType statsImagedM (Catalog *catalog) {
    427 
    428   int i, n;
    429   double *list, *dlist;
    430   StatType stats;
    431 
    432   bzero (&stats, sizeof (StatType));
    433   if (FREEZE_IMAGES) return (stats);
    434 
    435   ALLOCATE (list, double, Nimage);
    436   ALLOCATE (dlist, double, Nimage);
    437 
    438   n = 0;
    439   for (i = 0; i < Nimage; i++) {
    440 
    441     if (image[i].code & IMAGE_BAD)  continue;
    442 
    443     list[n] = image[i].dMcal;
    444     dlist[n] = 1;
    445     n++;
    446   }
    447 
    448   liststats (list, dlist, n, &stats);
    449   free (list);
    450   free (dlist);
    451   return (stats);
    452 }
    453 
    454 Image *getimages (int *N) {
    455 
    456   *N = Nimage;
    457   return (image);
    458 }
    459 
    460 Image *getimage (int N) {
    461   return (&image[N]);
    462 }
    463 
     180StarData *getImageRaw (int im, int *Nstars) {
     181
     182  StarData *raw;
     183 
     184  ALLOCATE (raw, StarData, Nlist[im]);
     185
     186  for (i = 0; i < Nlist[im]; i++) {
     187    m = mlist[im][i];
     188    c = clist[im][i];
     189
     190    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     191    raw[i].X = catalog[c].measure[m].Xccd;
     192    raw[i].Y = catalog[c].measure[m].Yccd;
     193   
     194    raw[i].Mag  = catalog[c].measure[m].M;
     195    raw[i].dMag = catalog[c].measure[m].dM;
     196
     197    /* 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);
     199  }
     200 
     201  *Nstars = Nlist[im];
     202  return (raw);
     203}
     204
     205StarData *getImageRef (int im, int *Nstars) {
     206
     207  StarData *raw;
     208 
     209  ALLOCATE (raw, StarData, Nlist[im]);
     210
     211  for (i = 0; i < Nlist[im]; i++) {
     212    m = mlist[im][i];
     213    c = clist[im][i];
     214    n = catalog[c].measure[m].averef;
     215
     216    /* apply the current image transformation or use the current value of R+dR, D+dD? */
     217    ref[i].R = catalog[c].average[n].R;
     218    ref[i].D = catalog[c].average[n].D;
     219  }
     220 
     221  *Nstars = Nlist[im];
     222  return (raw);
     223}

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

    r12050 r12068  
    22
    33static int Nmax;
    4 static double *listX, *dlistX;
    5 static double *listY, *dlistY;
    6 static double *listR, *dlistR;
    7 static double *listD, *dlistD;
    8 static double *listT, *dlistT;
     4static double *X, *dX;
     5static double *Y, *dY;
     6static double *R, *dR;
     7static double *D, *dD;
     8static time_t *T;
     9static double *dT;
    910
    1011void initObjectData (Catalog *catalog, int Ncatalog) {
    … …  
    1920  }
    2021
    21   ALLOCATE (list, double, MAX (1, Nmax));
    22   ALLOCATE (dlist, double, MAX (1, Nmax));
     22  ALLOCATE (R, double, MAX (1, Nmax));
     23  ALLOCATE (D, double, MAX (1, Nmax));
     24  ALLOCATE (T, time_t, MAX (1, Nmax));
     25  ALLOCATE (X, double, MAX (1, Nmax));
     26  ALLOCATE (Y, double, MAX (1, Nmax));
     27  ALLOCATE (pR, double, MAX (1, Nmax));
     28  ALLOCATE (pD, double, MAX (1, Nmax));
     29
     30  ALLOCATE (dR, double, MAX (1, Nmax));
     31  ALLOCATE (dD, double, MAX (1, Nmax));
     32  ALLOCATE (dT, double, MAX (1, Nmax));
     33  ALLOCATE (dX, double, MAX (1, Nmax));
     34  ALLOCATE (dY, double, MAX (1, Nmax));
    2335} 
    2436
    … …  
    2941  StatType statsR, statsD;
    3042  Coords coords;
     43  PMFit fit;
    3144
    3245  coords.crval1 = 0;
    … …  
    5164        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    5265       
    53         listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    54         listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
    55         listT[N] = catalog[i].measure[m].t;
     66        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     67        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     68        T[N] = catalog[i].measure[m].t;
    5669
    5770        /* the astrometric errors are not being carried yet (but should be!) */
    5871        /* we use the photometric mag error as a weighting term */
    59         dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    60         dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    61         dlistT[N] = 1.0; /* XXX hard-wire this for now */
     72        dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     73        dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
     74        dT[N] = catalog[i].measure[m].dt;
    6275
    6376        N++;
    … …  
    6578      if (N < STAR_TOOFEW) { /* too few measurements */
    6679        catalog[i].average[j].code |= ID_STAR_FEW;
     80        continue;
    6781      } else {
    6882        catalog[i].average[j].code &= ~ID_STAR_FEW;
    … …  
    7084
    7185      /* we need to do the fit in a locally linear space; choose a ref coordinate */
    72       coords.crval1 = listR[0];
    73       coords.crval2 = listD[0];
     86      coords.crval1 = R[0];
     87      coords.crval2 = D[0];
    7488     
    7589      /* project all of the R,D coordinates to a plane centered on this coordinate */
    7690      /* calculate pR[i], pD[i] for each point */
    7791      for (k = 0; k < N; k++) {
    78           RD_to_XY (&listX[k], &listY[k], listR[k], listD[k], &coords);
    79           ParFactor (&pX[k], &pY[k], listR[k], listD[k]);
     92          RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
    8093      }   
    8194
    82       /* in here, we should be fitting the parallax and proper-motion components */
    83       /* XXX if I fit R,D,T as above, I will be getting the wrong answers near the pole */
    84       if (FitProperMotion) {
    85           fit = FitPM (listX, dlistX, listY, dlistY, listT, pR, pD, N);
    86       }
     95      /* fit the model components as needed */
     96      switch (FITTING_MODE) {
     97        case FIT_AVERAGE:
     98          liststats (R, dR, N, &statsR);
     99          liststats (D, dD, N, &statsD);
    87100
    88       if (FitProperMotion) {
    89           fit = FitPMandPar (listX, dlistX, listY, dlistY, listT, pR, pD, N);
    90       }
     101          fit.Ro = statsR.mean;
     102          fit.dR = statsR.sigma;
    91103
    92       if (AverageCoords) {
    93           liststats (listR, dlistR, N, &statsR);
    94           liststats (listD, dlistD, N, &statsD);
     104          fit.Do = statsD.mean;
     105          fit.dD = statsD.sigma;
    95106
    96           catalog[i].average[j].R = statsR.mean;
    97           catalog[i].average[j].dR = statsR.sigma;
     107          fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
     108          fit.Nfit = N;
     109          break;
    98110
    99           catalog[i].average[j].D = statsD.mean;
    100           catalog[i].average[j].dD = statsD.sigma;
     111        case FIT_PM_ONLY:
     112          FitPM (&fit, X, dX, Y, dY, T, N);
     113          break;
    101114
    102           chisq = 0.5*(statsR.chisq + statsD.chisq);
    103       }
     115        case FIT_PM_AND_PAR:
     116          for (k = 0; k < N; k++) {
     117            ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
     118          }
     119          FitPMandPar (&fit, X, dX, Y, dY, T, pR, pD, N);
     120          break;
     121        default:
     122          fprintf (stderr, "invalid fitting mode %d\n", FITTING_MODE);
     123          exit (2);
     124      }   
    104125
    105       catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
     126      catalog[i].average[j].R   = fit.Ro;
     127      catalog[i].average[j].D   = fit.Do;
     128      catalog[i].average[j].dR  = fit.dRo;
     129      catalog[i].average[j].dD  = fit.dDo;
     130
     131      catalog[i].average[j].uR  = fit.uR;
     132      catalog[i].average[j].uD  = fit.uD;
     133      catalog[i].average[j].duR = fit.duR;
     134      catalog[i].average[j].duD = fit.duD;
     135
     136      catalog[i].average[j].P   = fit.p;
     137      catalog[i].average[j].dP  = fit.dp;
     138
     139      catalog[i].average[j].Xp  = (fit.Npts > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
    106140    }
    107141  }

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

    r12048 r12068  
    3434  initMosaics (subset, Nsubset);
    3535 
     36  /* unlock, if we can (else, unlocked below) */
     37  if (!UPDATE) dvo_image_unlock (db);
     38
    3639  return (skylist);
    3740}

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

    r11742 r12068  
    2727  skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
    2828
    29   /* unlock, if we can (else, unlocked below) */
    30   if (!UPDATE) dvo_image_unlock (&db);
    31 
    3229  /* load catalog data from region files */
    3330  catalog = load_catalogs (skylist, &Ncatalog);
    … …  
    3532  /* match measurements with images, mosaics */
    3633  initImageBins  (catalog, Ncatalog);
    37   initMosaicBins (catalog, Ncatalog);
    38   initGridBins   (catalog, Ncatalog);
    39   initMrel (catalog, Ncatalog);
     34  // initMosaicBins (catalog, Ncatalog);
    4035
    4136  findImages (catalog, Ncatalog);
    42   findMosaics (catalog, Ncatalog);  /* also sets Grid values */
    43   SAVEPLOT = FALSE;
    44 
    45   setExclusions (catalog, Ncatalog);
     37  // findMosaics (catalog, Ncatalog);  /* also sets Grid values */
    4638
    4739  if (PLOTSTUFF) {
    … …  
    5446    UpdateObjects (catalog, Ncatalog);
    5547  }
    56   if (UpdateImages) {
    57     UpdateImages (catalog, Ncatalog);
     48  if (UpdateSimple) {
     49    UpdateSimple ();
     50  }
     51  if (UpdateChips) {
     52    UpdateChips ();
     53  }
     54  if (UpdateMosaics) {
     55    UpdateMosaics ();
    5856  }
    5957

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

    r12048 r12068  
    11# include "relastro.h"
    22
     3/* XXX I think all of these, except for X,Y selection, are done / can be done in bcatalog */
    34int setExclusions (Catalog *catalog, int Ncatalog) {
    45
Note: See TracChangeset for help on using the changeset viewer.